gca_fill_gridvars Subroutine

private subroutine gca_fill_gridvars()

Uses

  • proc~~gca_fill_gridvars~~UsesGraph proc~gca_fill_gridvars gca_fill_gridvars module~mod_geometry mod_geometry proc~gca_fill_gridvars->module~mod_geometry module~mod_global_parameters mod_global_parameters proc~gca_fill_gridvars->module~mod_global_parameters module~mod_usr_methods mod_usr_methods proc~gca_fill_gridvars->module~mod_usr_methods module~mod_comm_lib mod_comm_lib module~mod_geometry->module~mod_comm_lib module~mod_basic_types mod_basic_types module~mod_global_parameters->module~mod_basic_types module~mod_connectivity mod_connectivity module~mod_global_parameters->module~mod_connectivity module~mod_constants mod_constants module~mod_global_parameters->module~mod_constants module~mod_physicaldata mod_physicaldata module~mod_global_parameters->module~mod_physicaldata module~mod_variables mod_variables module~mod_global_parameters->module~mod_variables mpi mpi module~mod_global_parameters->mpi module~mod_variables->module~mod_basic_types

Arguments

None

Calls

proc~~gca_fill_gridvars~~CallsGraph proc~gca_fill_gridvars gca_fill_gridvars phys_to_primitive phys_to_primitive proc~gca_fill_gridvars->phys_to_primitive proc~fill_gridvars_default fill_gridvars_default proc~gca_fill_gridvars->proc~fill_gridvars_default proc~gradient gradient proc~gca_fill_gridvars->proc~gradient proc~mpistop mpistop proc~gca_fill_gridvars->proc~mpistop proc~fields_from_mhd fields_from_mhd proc~fill_gridvars_default->proc~fields_from_mhd proc~gradient->proc~mpistop dsin dsin proc~gradient->dsin mpi_abort mpi_abort proc~mpistop->mpi_abort proc~fields_from_mhd->phys_to_primitive proc~particle_get_current particle_get_current proc~fields_from_mhd->proc~particle_get_current proc~curlvector curlvector proc~particle_get_current->proc~curlvector proc~curlvector->proc~mpistop proc~curlvector->dsin dcos dcos proc~curlvector->dcos

Source Code

  subroutine gca_fill_gridvars
    use mod_global_parameters
    use mod_usr_methods, only: usr_particle_fields
    use mod_geometry

    integer                                   :: igrid, iigrid, idir, idim
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:ndir) :: beta
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:nw)   :: w,wold
    double precision                          :: current(ixGlo1:ixGhi1,&
       ixGlo2:ixGhi2,ixGlo3:ixGhi3,7-2*ndir:3)
    integer                                   :: idirmin
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:ndir) :: vE, bhat
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
       ixGlo3:ixGhi3)        :: kappa, kappa_B, absB, tmp

    call fill_gridvars_default()

    do iigrid=1,igridstail; igrid=igrids(iigrid);
      w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         1:nw) = ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,1:nw)
      call phys_to_primitive(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,ixGlo1,&
         ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,w,ps(igrid)%x)

      ! grad(kappa B)
      absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) = &
         sqrt(sum(gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(:))**2,dim=ndim+1))
      vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         1) = gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(2)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(3)) - gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(3)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(2))
      vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         2) = gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(3)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(1)) - gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(1)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(3))
      vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         3) = gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(1)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(2)) - gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(2)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(1))
      do idir=1,ndir
        where (absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) .gt. 0.d0)
          vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             idir) = vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             idir) / absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3)**2
        elsewhere
          vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,idir) = 0.d0
        end where
      end do
      if (any(sum(vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,:)**2,&
         dim=ndim+1) .ge. c_norm**2) .and. relativistic) then
        call mpistop("GCA FILL GRIDVARS: vE>c! ABORTING...")
      end if
      if (any(vE .ne. vE)) then
        call mpistop("GCA FILL GRIDVARS: NaNs IN vE! ABORTING...")
      end if

      if (relativistic) then
        kappa(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
           ixGlo3:ixGhi3) = 1.d0/sqrt(1.0d0 - sum(vE(ixGlo1:ixGhi1,&
           ixGlo2:ixGhi2,ixGlo3:ixGhi3,:)**2,dim=ndim+1)/c_norm**2)
      else
        kappa(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) = 1.d0
      end if
      kappa_B(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) = absB(ixGlo1:ixGhi1,&
         ixGlo2:ixGhi2,ixGlo3:ixGhi3) / kappa(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
         ixGlo3:ixGhi3)

      if (any(kappa_B .ne. kappa_B)) then
        call mpistop("GCA FILL GRIDVARS: NaNs IN kappa_B! ABORTING...")
      end if

      tmp=0.d0
      do idim=1,ndim
        call gradient(kappa_B,ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
           ixGlo1+1,ixGlo2+1,ixGlo3+1,ixGhi1-1,ixGhi2-1,ixGhi3-1,idim,tmp)
        gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
           grad_kappa_B(idim)) = tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
           ixGlo3:ixGhi3)
      end do

      ! bhat
      do idir=1,ndir
        where (absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) .gt. 0.d0)
          bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             idir) = gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,bp(idir)) / absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
        elsewhere
          bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,idir) = 0.d0
        end where
      end do
      if (any(bhat .ne. bhat)) then
        call mpistop("GCA FILL GRIDVARS: NaNs IN bhat! ABORTING...")
      end if

      do idir=1,ndir
        ! (b dot grad) b and the other directional derivatives
        do idim=1,ndim
          call gradient(bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,idir),&
             ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,ixGlo1+1,ixGlo2+1,&
             ixGlo3+1,ixGhi1-1,ixGhi2-1,ixGhi3-1,idim,tmp)
          gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             b_dot_grad_b(idir)) = gridvars(igrid)%w(ixGlo1:ixGhi1,&
             ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             b_dot_grad_b(idir)) + bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,idim) * tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
          gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             vE_dot_grad_b(idir)) = gridvars(igrid)%w(ixGlo1:ixGhi1,&
             ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             vE_dot_grad_b(idir)) + vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,idim) * tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
          call gradient(vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,idir),&
             ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,ixGlo1+1,ixGlo2+1,&
             ixGlo3+1,ixGhi1-1,ixGhi2-1,ixGhi3-1,idim,tmp)
          gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             b_dot_grad_vE(idir)) = gridvars(igrid)%w(ixGlo1:ixGhi1,&
             ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             b_dot_grad_vE(idir)) + bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,idim) * tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
          gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             vE_dot_grad_vE(idir)) = gridvars(igrid)%w(ixGlo1:ixGhi1,&
             ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             vE_dot_grad_vE(idir)) + vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,idim) * tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
        end do
      end do

    end do

  end subroutine gca_fill_gridvars