set_B0.t

Go to the documentation of this file.

subroutine set_b0_grid(igrid)
  use mod_global_parameters
  use mod_mhd_phys, only: mhd_semirelativistic
 
  integer, intent(in) :: igrid
 
  integer :: ixCoG^L
 
  ixcogmin^d=1;
  ixcogmax^d=(ixghi^d-2*nghostcells)/2+2*nghostcells;
 
  call set_b0_cell(ps(igrid)%B0(ixg^t,:,0),ps(igrid)%x,ixg^ll,ixg^ll)
  if(mhd_semirelativistic) call set_b0_cell(psc(igrid)%B0(ixcog^s,:,0),psc(igrid)%x,ixcog^l,ixcog^l)
  call set_j0_cell(igrid,ps(igrid)%J0,ixg^ll,ixm^ll^ladd1)
  call set_b0_face(igrid,ps(igrid)%x,ixg^ll,ixm^ll)
 
end subroutine set_b0_grid
 
subroutine set_b0_cell(wB0,x,ixI^L,ix^L)
  use mod_usr_methods, only: usr_set_b0
  use mod_global_parameters
  use mod_geometry
  
  integer, intent(in):: ixI^L,ix^L
  double precision, intent(inout) :: wB0(ixI^S,1:ndir)
  double precision, intent(in) :: x(ixI^S,1:ndim)
 
  wb0(ix^s,1:ndir)=zero
 
  ! approximate cell-averaged B0 as cell-centered B0
  select case (coordinate)
  case (spherical)
     {^nooned
     if (dabs(bdip)>smalldouble) then
        wb0(ix^s,1)=2.0d0*bdip*dcos(x(ix^s,2))/x(ix^s,1)**3
        wb0(ix^s,2)=bdip*dsin(x(ix^s,2))/x(ix^s,1)**3
     end if
  
     if (abs(bquad)>smalldouble) then
        wb0(ix^s,1)=wb0(ix^s,1) &
             +bquad*0.5d0*(1.0d0+3.0d0*dcos(2.0d0*x(ix^s,2)))/x(ix^s,1)**4
        wb0(ix^s,2)=wb0(ix^s,2)+bquad*dsin(2.0d0*x(ix^s,2))/x(ix^s,1)**4
     end if
     if (abs(boct)>smalldouble) then
        wb0(ix^s,1)=wb0(ix^s,1) &
                     +boct*(10.0d0*dcos(2.0d0*x(ix^s,2))-2.0d0) &
                          *dcos(x(ix^s,2))/x(ix^s,1)**5
        wb0(ix^s,2)=wb0(ix^s,2) &
                     +boct*1.5d0*(3.0d0+5.0d0*dcos(2.0d0*x(ix^s,2))) &
                          *dsin(x(ix^s,2))/x(ix^s,1)**5
     end if
     }
  end select
  
  if (associated(usr_set_b0)) call usr_set_b0(ixi^l,ix^l,x,wb0)
 
end subroutine set_b0_cell
 
subroutine set_j0_cell(igrid,wJ0,ixI^L,ix^L)
  use mod_usr_methods, only: usr_set_j0
  use mod_global_parameters
  use mod_geometry
 
  integer, intent(in):: igrid,ixI^L,ix^L
  double precision, intent(inout) :: wJ0(ixI^S,7-2*ndir:3)
  integer :: idirmin0, idirmin
 
  if(associated(usr_set_j0)) then
    call usr_set_j0(ixi^l,ix^l,ps(igrid)%x,wj0)
  else
    idirmin0 = 7-2*ndir
    call curlvector(ps(igrid)%B0(ixi^s,:,0),ixi^l,ix^l,wj0,idirmin,idirmin0,ndir)
  end if
 
end subroutine set_j0_cell
 
subroutine set_b0_face(igrid,x,ixI^L,ixO^L)
  use mod_global_parameters
 
  integer, intent(in) :: igrid, ixI^L, ixO^L
  double precision, intent(in) :: x(ixI^S,1:ndim)
 
  double precision :: delx(ixI^S,1:ndim)
  double precision :: xC(ixI^S,1:ndim),xshift^D
  integer :: idims, ixC^L, hxO^L, ix, idims2
 
  if(slab_uniform)then
    ^d&delx(ixi^s,^d)=rnode(rpdx^d_,igrid)\
  else
    ! for all non-cartesian and stretched cartesian coordinates
    delx(ixi^s,1:ndim)=ps(igrid)%dx(ixi^s,1:ndim)
  endif
 
 
  do idims=1,ndim
    hxo^l=ixo^l-kr(idims,^d);
    if(stagger_grid) then
      ! ct needs all transverse cells
      ixcmax^d=ixomax^d+nghostcells-nghostcells*kr(idims,^d); ixcmin^d=hxomin^d-nghostcells+nghostcells*kr(idims,^d);
    else
      ! ixC is centered index in the idims direction from ixOmin-1/2 to ixOmax+1/2
      ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
    end if
    ! always xshift=0 or 1/2
    xshift^d=half*(one-kr(^d,idims));
    do idims2=1,ndim
      select case(idims2)
      {case(^d)
        do ix = ixc^lim^d
          ! xshift=half: this is the cell center coordinate
          ! xshift=0: this is the cell edge i+1/2 coordinate
          xc(ix^d%ixC^s,^d)=x(ix^d%ixC^s,^d)+(half-xshift^d)*delx(ix^d%ixC^s,^d)
        end do\}
      end select
    end do
    call set_b0_cell(ps(igrid)%B0(ixi^s,:,idims),xc,ixi^l,ixc^l)
  end do
 
end subroutine set_b0_face