mod_finite_difference.t

Go to the documentation of this file.

!> Module with finite difference methods for fluxes
module mod_finite_difference
 
  implicit none
  private
 
  public :: fd
  public :: centdiff
 
contains
 
  subroutine fd(qdt,dtfactor,ixI^L,ixO^L,idims^LIM,qtC,sCT,qt,snew,fC,fE,dxs,x)
    use mod_physics
    use mod_source, only: addsource2
    use mod_finite_volume, only: reconstruct_lr
    use mod_global_parameters
    use mod_usr_methods
 
    double precision, intent(in)                                     :: qdt, dtfactor, qtc, qt, dxs(ndim)
    integer, intent(in)                                              :: ixi^l, ixo^l, idims^lim
    double precision, dimension(ixI^S,1:ndim), intent(in)            :: x
 
    type(state)                                                      :: sct, snew
    double precision, dimension(ixI^S,1:nwflux,1:ndim), intent(out)  :: fc
    double precision, dimension(ixI^S,sdim:3)                        :: fe
 
    double precision, dimension(ixI^S,1:nwflux)                      :: fct
    double precision, dimension(ixI^S,1:nw)                          :: fm, fp, fmr, fpl, wprim
    ! left and right constructed status in conservative form
    double precision, dimension(ixI^S,1:nw) :: wlc, wrc
    ! left and right constructed status in primitive form, needed for better performance
    double precision, dimension(ixI^S,1:nw) :: wlp, wrp
    double precision, dimension(ixI^S)      :: cmaxc, cminc
    double precision, dimension(ixI^S)      :: hspeed
    double precision, dimension(ixO^S)      :: inv_volume
    double precision, dimension(1:ndim)     :: dxinv
    logical :: transport, active
    integer :: idims, iw, ixc^l, ix^l, hxo^l, kxc^l, kxr^l
    type(ct_velocity) :: vcts
 
    associate(wct=>sct%w,wnew=>snew%w)
 
    fc=0.d0
    wprim=wct
    call phys_to_primitive(ixi^l,ixi^l,wprim,x)
 
    b0i = 0 ! fd uses centered values in phys_get_flux  
    do idims= idims^lim
 
       !b0i=idims
 
       ! Get fluxes for the whole grid (mesh+nghostcells)
       {^d& ixmin^d = ixomin^d - nghostcells * kr(idims,^d)\}
       {^d& ixmax^d = ixomax^d + nghostcells * kr(idims,^d)\}
 
       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);
         ixmax^d=ixmax^d+nghostcells-nghostcells*kr(idims,^d); ixmin^d=ixmin^d-nghostcells+nghostcells*kr(idims,^d);
         kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
         kxr^l=kxc^l+kr(idims,^d);
         ! wRp and wLp are defined at the same locations, and will correspond to
         ! the left and right reconstructed values at a cell face. Their indexing
         ! is similar to cell-centered values, but in direction idims they are
         ! shifted half a cell towards the 'lower' direction.
         wrp(kxc^s,1:nw)=wprim(kxr^s,1:nw)
         wlp(kxc^s,1:nw)=wprim(kxc^s,1:nw)
       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
 
       call phys_get_flux(wct,wprim,x,ixg^ll,ix^l,idims,fct)
 
       do iw=iwstart,nwflux
          ! Lax-Friedrich splitting:
          fp(ix^s,iw) = half * (fct(ix^s,iw) + tvdlfeps * cmax_global * wct(ix^s,iw))
          fm(ix^s,iw) = half * (fct(ix^s,iw) - tvdlfeps * cmax_global * wct(ix^s,iw))
       end do ! iw loop
 
       ! now do the reconstruction of fp and fm:
       call reconstructl(ixi^l,ixc^l,idims,fp,fpl)
       call reconstructr(ixi^l,ixc^l,idims,fm,fmr)
 
       fc(ixc^s,iwstart:nwflux,idims) = fpl(ixc^s,iwstart:nwflux) + fmr(ixc^s,iwstart:nwflux)
 
       if(stagger_grid) then
         ! apply limited reconstruction for left and right status at cell interfaces
         call reconstruct_lr(ixi^l,ixc^l,ixc^l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
         if(h_correction) then
           call phys_get_h_speed(wprim,x,ixi^l,ixo^l,idims,hspeed)
         end if
         call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^l,ixc^l,idims,hspeed,cmaxc,cminc)
         call phys_get_ct_velocity(vcts,wlp,wrp,ixi^l,ixc^l,idims,cmaxc,cminc)
       end if
 
    end do !idims loop
    !b0i=0
 
    if(stagger_grid) call phys_update_faces(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,snew,vcts)
 
    if(slab_uniform) then
      dxinv=-qdt/dxs
      do idims= idims^lim
        hxo^l=ixo^l-kr(idims,^d);
        if(local_timestep) then
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims) = -block%dt(ixi^s)*dtfactor/dxs(idims) * fc(ixi^s,iw,idims)
            wnew(ixo^s,iw)=wnew(ixo^s,iw)+(fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))
          end do ! iw loop
        else
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims) = dxinv(idims) * fc(ixi^s,iw,idims)
            wnew(ixo^s,iw)=wnew(ixo^s,iw)+(fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))
          end do ! iw loop
        endif
      end do ! Next idims
    else
      inv_volume=1.d0/block%dvolume(ixo^s)
      do idims= idims^lim
        hxo^l=ixo^l-kr(idims,^d);
        if(local_timestep) then
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims)=-block%dt(ixi^s)*dtfactor*fc(ixi^s,iw,idims)*block%surfaceC(ixi^s,idims)
            wnew(ixo^s,iw)=wnew(ixo^s,iw)+ &
                 (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))*inv_volume(ixo^s)
          end do ! iw loop
        else
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims)=-qdt*fc(ixi^s,iw,idims)*block%surfaceC(ixi^s,idims)
            wnew(ixo^s,iw)=wnew(ixo^s,iw)+ &
                 (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))*inv_volume(ixo^s)
          end do ! iw loop
        endif ! local_timestep
      end do ! Next idims
    end if
 
    if (.not.slab.and.idimsmin==1) call phys_add_source_geom(qdt,dtfactor,ixi^l,ixo^l,wct,wprim,wnew,x)
 
    if(stagger_grid) call phys_face_to_center(ixo^l,snew)
 
    ! check and optionally correct unphysical values
    if(fix_small_values) then
       call phys_handle_small_values(.false.,wnew,x,ixi^l,ixo^l,'fd')
       if(crash) then
         ! replace erroneous values with values at previous step
         wnew=wct
         if(stagger_grid) snew%ws=sct%ws
       end if
    endif
 
    call addsource2(qdt*dble(idimsmax-idimsmin+1)/dble(ndim), &
         dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim), &
         ixi^l,ixo^l,1,nw,qtc,wct,wprim,qt,wnew,x,.false.,active)
 
    end associate
 
  end subroutine fd
 
  subroutine reconstructl(ixI^L,iL^L,idims,w,wLC)
    use mod_global_parameters
    use mod_mp5
    use mod_limiter
    use mod_comm_lib, only: mpistop
 
    integer, intent(in)             :: ixI^L, iL^L, idims
    double precision, intent(in)    :: w(ixI^S,1:nw)
 
    double precision, intent(out)   :: wLC(ixI^S,1:nw) 
 
    double precision                :: ldw(ixI^S), dwC(ixI^S)
    integer                         :: jxR^L, ixC^L, jxC^L, kxC^L, iw
    double precision                :: a2max
 
    select case (type_limiter(block%level))
    case (limiter_mp5)
       call mp5limiterl(ixi^l,il^l,idims,w,wlc)
    case (limiter_weno5)
       call weno5limiterl(ixi^l,il^l,idims,w,wlc,1)
    case (limiter_weno5nm)
       call weno5nmlimiterl(ixi^l,il^l,idims,w,wlc,1)
    case (limiter_wenoz5)
       call weno5limiterl(ixi^l,il^l,idims,w,wlc,2)
    case (limiter_wenoz5nm)
       call weno5nmlimiterl(ixi^l,il^l,idims,w,wlc,2)
    case (limiter_wenozp5)
       call weno5limiterl(ixi^l,il^l,idims,w,wlc,3)
    case (limiter_wenozp5nm)
       call weno5nmlimiterl(ixi^l,il^l,idims,w,wlc,3)
    case default 
 
       kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
 
       wlc(kxc^s,iwstart:nwflux) = w(kxc^s,iwstart:nwflux)
 
       jxr^l=il^l+kr(idims,^d);
 
       ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
       jxc^l=ixc^l+kr(idims,^d);
 
       do iw=iwstart,nwflux
          dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
 
           if(need_global_a2max) then 
             a2max=a2max_global(idims)
           else
             select case(idims)
             case(1)
               a2max=schmid_rad1
             {^iftwod
             case(2)
               a2max=schmid_rad2}
             {^ifthreed
             case(2)
               a2max=schmid_rad2
             case(3)
               a2max=schmid_rad3}
             case default
               call mpistop("idims is wrong in mod_limiter")
             end select
           end if
 
          call dwlimiter2(dwc,ixi^l,ixc^l,idims,type_limiter(block%level),ldw=ldw,a2max=a2max)
 
          wlc(il^s,iw)=wlc(il^s,iw)+half*ldw(il^s)
       end do
    end select
 
  end subroutine reconstructl
 
  subroutine reconstructr(ixI^L,iL^L,idims,w,wRC)
    use mod_global_parameters
    use mod_mp5
    use mod_limiter
    use mod_comm_lib, only: mpistop
 
    integer, intent(in)             :: ixI^L, iL^L, idims
    double precision, intent(in)    :: w(ixI^S,1:nw)
 
    double precision, intent(out)   :: wRC(ixI^S,1:nw) 
 
    double precision                :: rdw(ixI^S), dwC(ixI^S)
    integer                         :: jxR^L, ixC^L, jxC^L, kxC^L, kxR^L, iw
    double precision                :: a2max
 
    select case (type_limiter(block%level))
    case (limiter_mp5)
       call mp5limiterr(ixi^l,il^l,idims,w,wrc)
    case (limiter_weno5)
       call weno5limiterr(ixi^l,il^l,idims,w,wrc,1)
    case (limiter_weno5nm)
       call weno5nmlimiterr(ixi^l,il^l,idims,w,wrc,1)
    case (limiter_wenoz5)
       call weno5limiterr(ixi^l,il^l,idims,w,wrc,2)
    case (limiter_wenoz5nm)
       call weno5nmlimiterr(ixi^l,il^l,idims,w,wrc,2)
    case (limiter_wenozp5)
       call weno5limiterr(ixi^l,il^l,idims,w,wrc,3)
    case (limiter_wenozp5nm)
       call weno5nmlimiterr(ixi^l,il^l,idims,w,wrc,3)
    case default 
 
       kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
       kxr^l=kxc^l+kr(idims,^d);
 
       wrc(kxc^s,iwstart:nwflux)=w(kxr^s,iwstart:nwflux)
 
       jxr^l=il^l+kr(idims,^d);
       ixcmax^d=jxrmax^d; ixcmin^d=ilmin^d-kr(idims,^d);
       jxc^l=ixc^l+kr(idims,^d);
 
       do iw=iwstart,nwflux
          dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
 
           if(need_global_a2max) then
             a2max=a2max_global(idims)
           else
             select case(idims)
             case(1)
               a2max=schmid_rad1
             {^iftwod
             case(2)
               a2max=schmid_rad2}
             {^ifthreed
             case(2)
               a2max=schmid_rad2
             case(3)
               a2max=schmid_rad3}
             case default
               call mpistop("idims is wrong in mod_limiter")
             end select
           end if
 
          call dwlimiter2(dwc,ixi^l,ixc^l,idims,type_limiter(block%level),rdw=rdw,a2max=a2max)
 
          wrc(il^s,iw)=wrc(il^s,iw)-half*rdw(jxr^s)
       end do
    end select
 
  end subroutine reconstructr
 
  subroutine centdiff(method,qdt,dtfactor,ixI^L,ixO^L,idims^LIM,qtC,sCT,qt,s,fC,fE,dxs,x)
 
    ! Advance the flow variables from global_time to global_time+qdt within ixO^L by
    ! fourth order centered differencing in space 
    ! for the dw/dt+dF_i(w)/dx_i=S type equation.
    ! wCT contains the time centered variables at time qtC for flux and source.
    ! w is the old value at qt on input and the new value at qt+qdt on output.
    use mod_physics
    use mod_finite_volume, only: reconstruct_lr
    use mod_global_parameters
    use mod_source, only: addsource2
    use mod_usr_methods
    use mod_variables
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: method
    integer, intent(in) :: ixi^l, ixo^l, idims^lim
    double precision, intent(in) :: qdt, dtfactor, qtc, qt, dxs(ndim)
    type(state)      :: sct, s
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision :: fc(ixi^s,1:nwflux,1:ndim)
    double precision :: fe(ixi^s,sdim:3)
 
    double precision :: v(ixi^s,ndim), f(ixi^s, nwflux)
 
    double precision, dimension(ixI^S,1:nw) :: wprim, wlc, wrc
    ! left and right constructed status in primitive form, needed for better performance
    double precision, dimension(ixI^S,1:nw) :: wlp, wrp
    double precision, dimension(ixI^S)      :: vlc, cmaxlc, cmaxrc
    double precision, dimension(ixI^S,1:number_species)      :: cmaxc
    double precision, dimension(ixI^S,1:number_species)      :: cminc
    double precision, dimension(ixI^S)      :: hspeed
    double precision, dimension(ixO^S)      :: inv_volume
 
    double precision :: dxinv(1:ndim)
    integer :: idims, iw, ix^l, hxo^l, ixc^l, jxc^l, hxc^l, kxc^l, kkxc^l, kkxr^l
    type(ct_velocity) :: vcts
    logical :: transport, new_cmax, patchw(ixi^s), active
 
    associate(wct=>sct%w,w=>s%w)
    ! two extra layers are needed in each direction for which fluxes are added.
    ix^l=ixo^l;
    do idims= idims^lim
       ix^l=ix^l^ladd2*kr(idims,^d);
    end do
 
    if (ixi^l^ltix^l|.or.|.or.) then
       call mpistop("Error in centdiff: Non-conforming input limits")
    end if
 
    fc=0.d0
    wprim=wct
    call phys_to_primitive(ixi^l,ixi^l,wprim,x)
 
    ! get fluxes
    do idims= idims^lim
       b0i=idims
 
       ix^l=ixo^l^ladd2*kr(idims,^d); 
       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);
         ixmax^d=ixmax^d+nghostcells-nghostcells*kr(idims,^d);
         ixmin^d=ixmin^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
       hxc^l=ixc^l-kr(idims,^d); 
       jxc^l=ixc^l+kr(idims,^d); 
       kxc^l=ixc^l+2*kr(idims,^d); 
 
       kkxcmin^d=iximin^d; kkxcmax^d=iximax^d-kr(idims,^d);
       kkxr^l=kkxc^l+kr(idims,^d);
       wrp(kkxc^s,1:nwflux)=wprim(kkxr^s,1:nwflux)
       wlp(kkxc^s,1:nwflux)=wprim(kkxc^s,1:nwflux)
 
       ! apply limited reconstruction for left and right status at cell interfaces
       call reconstruct_lr(ixi^l,ixc^l,ixc^l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
 
       if(stagger_grid) then
         if(h_correction) then
           call phys_get_h_speed(wprim,x,ixi^l,ixo^l,idims,hspeed)
         end if
         call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^l,ixc^l,idims,hspeed,cmaxc,cminc)
         call phys_get_ct_velocity(vcts,wlp,wrp,ixi^l,ixc^l,idims,cmaxc,cminc)
       end if
 
       ! Calculate velocities from upwinded values
       call phys_get_cmax(wlp,x,ixg^ll,ixc^l,idims,cmaxlc)
       call phys_get_cmax(wrp,x,ixg^ll,ixc^l,idims,cmaxrc)
       ! now take the maximum of left and right states
       vlc(ixc^s)=max(cmaxrc(ixc^s),cmaxlc(ixc^s))
 
       call phys_get_flux(wct,wprim,x,ixi^l,ix^l,idims,f)
 
       ! Center flux to interface
       if(method==fs_cd) then
          fc(ixc^s,iwstart:nwflux,idims)=half*(f(ixc^s,iwstart:nwflux)+f(jxc^s,iwstart:nwflux))
       else
          ! f_i+1/2= (-f_(i+2) +7 f_(i+1) + 7 f_i - f_(i-1))/12
          fc(ixc^s,iwstart:nwflux,idims)=(-f(kxc^s,iwstart:nwflux)+7.0d0*(f(jxc^s,iwstart:nwflux) + &
               f(ixc^s,iwstart:nwflux))-f(hxc^s,iwstart:nwflux))/12.0d0
 
          do iw=iwstart,nwflux
             ! add rempel dissipative flux, only second order version for now
             fc(ixc^s,iw,idims)=fc(ixc^s,iw,idims)-tvdlfeps*half*vlc(ixc^s) &
                  *(wrc(ixc^s,iw)-wlc(ixc^s,iw))
          end do
       end if
 
    end do       !next idims
    b0i=0
 
    if(stagger_grid) call phys_update_faces(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,s,vcts)
 
    if(slab_uniform) then
      dxinv=-qdt/dxs
      do idims= idims^lim
        hxo^l=ixo^l-kr(idims,^d);
        if(local_timestep) then
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims)=-block%dt(ixi^s)*dtfactor/dxs(idims)*fc(ixi^s,iw,idims)
            ! result: f_(i+1/2)-f_(i-1/2) = [-f_(i+2)+8(f_(i+1)-f_(i-1))+f_(i-2)]/12
            w(ixo^s,iw)=w(ixo^s,iw)+(fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))
          end do    !next iw
        else
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims)=dxinv(idims)*fc(ixi^s,iw,idims)
            ! result: f_(i+1/2)-f_(i-1/2) = [-f_(i+2)+8(f_(i+1)-f_(i-1))+f_(i-2)]/12
            w(ixo^s,iw)=w(ixo^s,iw)+(fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))
          end do    !next iw
        endif
      end do ! Next idims
    else
      inv_volume=1.d0/block%dvolume
      do idims= idims^lim
        hxo^l=ixo^l-kr(idims,^d);
        if(local_timestep) then
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims)=-block%dt(ixi^s)*dtfactor*block%surfaceC(ixi^s,idims)*fc(ixi^s,iw,idims)
            w(ixo^s,iw)=w(ixo^s,iw)+ &
                 (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))*inv_volume(ixo^s)
          end do    !next iw
        else
          do iw=iwstart,nwflux
            fc(ixi^s,iw,idims)=-qdt*block%surfaceC(ixi^s,idims)*fc(ixi^s,iw,idims)
            w(ixo^s,iw)=w(ixo^s,iw)+ &
                 (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))*inv_volume(ixo^s)
          end do    !next iw
        endif 
      end do ! Next idims
    end if
 
    if (.not.slab.and.idimsmin==1) call phys_add_source_geom(qdt,dtfactor,ixi^l,ixo^l,wct,wprim,w,x)
 
    if(stagger_grid) call phys_face_to_center(ixo^l,s)
 
    ! check and optionally correct unphysical values
    if(fix_small_values) then
       call phys_handle_small_values(.false.,w,x,ixi^l,ixo^l,'centdiff')
    endif
 
    call addsource2(qdt*dble(idimsmax-idimsmin+1)/dble(ndim), &
         dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim), &
         ixi^l,ixo^l,1,nw,qtc,wct,wprim,qt,w,x,.false.,active)
 
    end associate
  end subroutine centdiff
 
end module mod_finite_difference