mod_constrained_transport.t

Go to the documentation of this file.

module mod_constrained_transport
  implicit none
  public
 
contains
  !> re-calculate the magnetic field from the vector potential in a completely
  !> divergency free way
  subroutine recalculateb
    use mod_global_parameters
    use mod_fix_conserve
    use mod_physics
    use mod_ghostcells_update
 
    integer :: igrid,iigrid
 
    call init_comm_fix_conserve(1,ndim,^nd)
 
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ! Make zero the magnetic fluxes
       ! Fake advance, storing electric fields at edges
       block=>ps(igrid)
       call fake_advance(igrid,1,^nd,ps(igrid))
 
    end do
    !$OMP END PARALLEL DO
 
    ! Do correction
    call recvflux(1,ndim)
    call sendflux(1,ndim)
    call fix_conserve(ps,1,ndim,1,^nd)
 
    call fix_edges(ps,1,^nd)
 
    ! Now we fill the centers for the staggered variables
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
       call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
       ! update cell center magnetic field
       call phys_face_to_center(ixm^ll,ps(igrid))
       call phys_to_conserved(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
    end do
    !$OMP END PARALLEL DO
 
    call getbc(global_time,0.d0,ps,iwstart,nwgc)
 
  end subroutine recalculateb
 
  !> fake advance a step to calculate magnetic field
  subroutine fake_advance(igrid,idim^LIM,s)
    use mod_global_parameters
    use mod_fix_conserve
 
    integer       :: igrid,idim^LIM
    type(state)   :: s
 
    double precision             :: dx^D
    double precision             :: fC(ixG^T,1:nwflux,1:ndim)
    double precision             :: fE(ixG^T,sdim:3)
 
    dx^d=rnode(rpdx^d_,igrid);
 
    call fake_update(ixg^ll,s,fc,fe,dx^d)
 
    call store_flux(igrid,fc,idim^lim,^nd)
    call store_edge(igrid,ixg^ll,fe,idim^lim) 
 
  end subroutine fake_advance
 
  !>fake update magnetic field from vector potential
  subroutine fake_update(ixI^L,s,fC,fE,dx^D)
    use mod_global_parameters
 
    integer       :: ixI^L
    type(state)   :: s
    double precision             :: fC(ixI^S,1:nwflux,1:ndim)
    double precision             :: fE(ixI^S,sdim:3)
    double precision             :: dx^D
 
    double precision                   :: A(s%ixGs^S,1:3)
    integer                            :: ixIs^L,ixO^L,idir
 
    associate(ws=>s%ws,x=>s%x)
 
    a=zero
    ws=zero
 
    ixis^l=s%ixGs^l;
    ixo^l=ixi^l^lsubnghostcells;
 
    fc=0.d0
    call b_from_vector_potentiala(ixis^l, ixi^l, ixo^l, ws, x, a)
 
    ! This is important only in 3D
    do idir=sdim,3
       fe(ixi^s,idir) =-a(ixi^s,idir)
    end do
 
    end associate
  end subroutine fake_update
 
  !> calculate magnetic field from vector potential A at cell edges
  subroutine b_from_vector_potentiala(ixIs^L, ixI^L, ixO^L, ws, x, A)
    use mod_global_parameters
    use mod_usr_methods, only: usr_init_vector_potential
 
    integer, intent(in)                :: ixIs^L, ixI^L, ixO^L
    double precision, intent(inout)    :: ws(ixIs^S,1:nws),A(ixIs^S,1:3)
    double precision, intent(in)       :: x(ixI^S,1:ndim)
 
    double precision                   :: xC(ixIs^S,1:ndim),xCC(ixIs^S,1:ndim)
    double precision                   :: circ(ixIs^S,1:ndim)
    integer                            :: ixC^L, hxC^L, idim1, idim2, idir
 
    a=zero
    ! extend one layer of cell center locations in xCC
    xc=0.d0
    xcc=0.d0
    xcc(ixi^s,1:ndim)=x(ixi^s,1:ndim)
    {
    xcc(ixismin^d^d%ixI^s,1:ndim)=x(iximin^d^d%ixI^s,1:ndim)
    xcc(ixismin^d^d%ixIs^s,^d)=x({iximin^dd,},^d)-block%dx({iximin^dd,},^d)
    \}
    {^ifthreed
    xcc(iximin1:iximax1,ixismin2,ixismin3,1)=x(iximin1:iximax1,iximin2,iximin3,1)
    xcc(ixismin1,iximin2:iximax2,ixismin3,2)=x(iximin1,iximin2:iximax2,iximin3,2)
    xcc(ixismin1,ixismin2,iximin3:iximax3,3)=x(iximin1,iximin2,iximin3:iximax3,3)
    }
 
    do idir=sdim,3
      ixcmax^d=ixomax^d;
      ixcmin^d=ixomin^d-1+kr(idir,^d);
      do idim1=1,ndim
        ! Get edge coordinates
        if (idim1/=idir) then
          xc(ixc^s,idim1)=xcc(ixc^s,idim1)+half*block%dx(ixc^s,idim1)
        else
          xc(ixc^s,idim1)=xcc(ixc^s,idim1)
        end if
      end do
      ! Initialise vector potential at the edge
      call usr_init_vector_potential(ixis^l, ixc^l, xc, a(ixis^s,idir), idir)
      a(ixc^s,idir)=a(ixc^s,idir)*block%dsC(ixc^s,idir)
    end do
 
    ! Take the curl of the vector potential 
    circ=zero
    ! Calculate circulation on each face
    do idim1=1,ndim ! Coordinate perpendicular to face 
      ixcmax^d=ixomax^d;
      ixcmin^d=ixomin^d-kr(idim1,^d);
      do idim2=1,ndim
        do idir=sdim,3 ! Direction of line integral
          if(lvc(idim1,idim2,idir)==0) cycle
          ! Assemble indices
          hxc^l=ixc^l-kr(idim2,^d);
          ! Add line integrals in direction idir
          circ(ixc^s,idim1)=circ(ixc^s,idim1)&
                           +lvc(idim1,idim2,idir)* &
                            (a(ixc^s,idir)&
                            -a(hxc^s,idir))
        end do
      end do
      ! Divide by the area of the face to get B
      where(block%surfaceC(ixc^s,idim1)==0)
        circ(ixc^s,idim1)=zero
      elsewhere
        circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
      end where
      ws(ixc^s,idim1) = circ(ixc^s,idim1)
    end do
 
  end subroutine b_from_vector_potentiala
 
  !> Reconstruct scalar q within ixO^L to 1/2 dx in direction idir
  !> Return both left and right reconstructed values
  !> Only implemented for hll flavour of CT
  subroutine reconstruct(ixI^L,ixC^L,idir,q,qL,qR)
    use mod_limiter
    use mod_global_parameters
 
    integer, intent(in)                :: ixI^L, ixC^L, idir
    double precision, intent(in)       :: q(ixI^S)
    double precision, intent(out)      :: qL(ixI^S), qR(ixI^S)
 
    double precision                   :: qC(ixI^S)
    double precision,dimension(ixI^S)  :: dqC,ldq,rdq
    integer                            :: ixO^L,jxC^L,gxC^L,hxC^L
 
    jxc^l=ixc^l+kr(idir,^d);
    gxcmin^d=ixcmin^d-kr(idir,^d);gxcmax^d=jxcmax^d;
    hxc^l=gxc^l+kr(idir,^d);
 
    qr(gxc^s) = q(hxc^s)
    ql(gxc^s) = q(gxc^s)
 
    select case (type_limiter(block%level))
    case (limiter_ppm)
       ! the ordinary grid-index:
       ixomin^d=ixcmin^d+kr(idir,^d);
       ixomax^d=ixcmax^d;
       call ppmlimitervar(ixi^l,ixo^l,idir,q,q,ql,qr)
    case (limiter_mp5)
       call mp5limitervar(ixi^l,ixc^l,idir,q,ql,qr)
    case (limiter_weno3)
       call weno3limiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,1)
    case (limiter_wenoyc3)
       call weno3limiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,2)
    case (limiter_weno5)
       call weno5limiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,1)
    case (limiter_weno5nm)
       call weno5nmlimiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,1)
    case (limiter_wenoz5)
       call weno5limiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,2)
    case (limiter_wenoz5nm)
       call weno5nmlimiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,2)
    case (limiter_wenozp5)
       call weno5limiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,3)
    case (limiter_wenozp5nm)
       call weno5nmlimiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,3)
    case (limiter_weno5cu6)
       call weno5cu6limiter(ixi^l,ixc^l,idir,q,ql,qr)
    case (limiter_teno5ad)
       call teno5adlimiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr)
    case (limiter_weno7)
       call weno7limiter(ixi^l,ixc^l,idir,q,ql,qr,1)
    case (limiter_mpweno7)
       call weno7limiter(ixi^l,ixc^l,idir,q,ql,qr,2)
    case (limiter_venk)
       call venklimiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr)
    case default
       dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
       call dwlimiter2(dqc,ixi^l,gxc^l,idir,type_limiter(block%level),ldq,rdq)
       ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
       qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
    end select
 
  end subroutine reconstruct
 
end module mod_constrained_transport