mod_b0.t

Go to the documentation of this file.

module mod_b0
  implicit none
  private
  public :: set_b0_grid
 
contains
 
  subroutine set_b0_grid(igrid)
    use mod_global_parameters
    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(b0fieldalloccoarse) 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
end module mod_b0