mod_refine.t

Go to the documentation of this file.

module mod_refine
 
  implicit none
  private
 
  public :: refine_grids
 
contains
 
 
  !> refine one block to its children blocks
  subroutine refine_grids(child_igrid,child_ipe,igrid,ipe,active)
    use mod_global_parameters
    use mod_initialize_amr, only: initial_condition
    use mod_amr_solution_node, only: alloc_node, dealloc_node
  
    integer, dimension(2^D&), intent(in) :: child_igrid, child_ipe
    integer, intent(in) :: igrid, ipe
    logical, intent(in) :: active
  
    integer :: ic^d
  
    ! allocate solution space for new children
    {do ic^db=1,2\}
       call alloc_node(child_igrid(ic^d))
    {end do\}
  
    if ((time_advance .and. active).or.convert.or.reset_grid) then
       ! prolong igrid to new children
       call prolong_grid(child_igrid,child_ipe,igrid,ipe)
    else
       ! Fill new created children with initial condition
       {do ic^db=1,2\}
          call initial_condition(child_igrid(ic^d))
       {end do\}
    end if
  
    ! remove solution space of igrid to save memory when converting data
    if(convert) call dealloc_node(igrid)
  end subroutine refine_grids
  
  !> prolong one block
  subroutine prolong_grid(child_igrid,child_ipe,igrid,ipe)
    use mod_physics, only: phys_to_primitive, phys_to_conserved, &
         phys_to_prolong, phys_from_prolong
    use mod_global_parameters
    use mod_amr_fct, only: old_neighbors
  
    integer, dimension(2^D&), intent(in) :: child_igrid, child_ipe
    integer, intent(in) :: igrid, ipe
  
    double precision :: dxco^d, xcomin^d, dxfi^d, xfimin^d
    integer :: ix^l, ichild, ixco^l, ic^d
  
    ix^l=ixm^ll^ladd1;
 
    ! Convert parent to prolongation space:
    !   phys_to_prolong: EoS-aware (rho, v, T) — avoids Jensen's inequality
    !   phys_to_primitive: standard (rho, v, p) — used with prolongprimitive
    ! the conversions may read block%B0
    block=>ps(igrid)
    if(associated(phys_to_prolong)) then
      call phys_to_prolong(ixg^ll,ix^l,ps(igrid)%w,ps(igrid)%x)
    else if(prolongprimitive) then
      call phys_to_primitive(ixg^ll,ix^l,ps(igrid)%w,ps(igrid)%x)
    end if
  
    xcomin^d=rnode(rpxmin^d_,igrid)\
    dxco^d=rnode(rpdx^d_,igrid)\
  
    if(stagger_grid) call old_neighbors(child_igrid,child_ipe,igrid,ipe)
  
    {do ic^db=1,2\}
      ichild=child_igrid(ic^d)
  
      ixcomin^d=ixmlo^d+(ic^d-1)*block_nx^d/2\
      ixcomax^d=ixmhi^d+(ic^d-2)*block_nx^d/2\
  
      xfimin^d=rnode(rpxmin^d_,ichild)\
      dxfi^d=rnode(rpdx^d_,ichild)\
      ^d&dxlevel(^d)=dxfi^d;
      call prolong_2nd(ps(igrid),ixco^l,ps(ichild), &
           dxco^d,xcomin^d,dxfi^d,xfimin^d,igrid,ichild)
      !call prolong_1st(ps(igrid)%w,ixCo^L,ps(ichild)%w,ps(ichild)%x)
    {end do\}
  
    ! Convert parent back from prolongation space
    block=>ps(igrid)
    if(associated(phys_to_prolong)) then
      call phys_from_prolong(ixg^ll,ix^l,ps(igrid)%w,ps(igrid)%x)
    else if (prolongprimitive) then
      call phys_to_conserved(ixg^ll,ix^l,ps(igrid)%w,ps(igrid)%x)
    end if
 
  end subroutine prolong_grid
  
  !> do 2nd order prolongation
  subroutine prolong_2nd(sCo,ixCo^L,sFi,dxCo^D,xComin^D,dxFi^D,xFimin^D,igridCo,igridFi)
    use mod_physics, only: phys_to_conserved, phys_to_primitive, &
         phys_handle_small_values, phys_wb_prolong, &
         phys_to_prolong, phys_from_prolong
    use mod_global_parameters
    use mod_amr_fct, only: already_fine, prolong_2nd_stg
  
    integer, intent(in) :: ixco^l, igridfi, igridco
    double precision, intent(in) :: dxco^d, xcomin^d, dxfi^d, xfimin^d
    type(state), intent(in)      :: sco
    type(state), intent(inout)   :: sfi
  
    double precision :: slopel, sloper, slopec, signc, signr
    double precision :: slope(nw,ndim)
    double precision :: eta^d
    integer :: ixco^d, jxco^d, hxco^d, ixfi^d, ix^d, idim, iw, ixcg^l, el
    logical :: fine_^l
  
    associate(wco=>sco%w, wfi=>sfi%w)
    ixcg^l=ixco^l;
    {do ixco^db = ixcg^lim^db
       ! lower left grid index in finer child block
       ixfi^db=2*(ixco^db-ixcomin^db)+ixmlo^db\}
  
       do idim=1,ndim
          hxco^d=ixco^d-kr(^d,idim)\
          jxco^d=ixco^d+kr(^d,idim)\
  
          do iw=1,nw
             slopel=wco(ixco^d,iw)-wco(hxco^d,iw)
             sloper=wco(jxco^d,iw)-wco(ixco^d,iw)
             slopec=half*(sloper+slopel)
  
             ! get limited slope
             signr=sign(one,sloper)
             signc=sign(one,slopec)
             !select case(prolong_limiter)
             !case(1)
             !  ! unlimited
             !  slope(iw,idim)=slopeC
             !case(2)
             !  ! minmod
             !  slope(iw,idim)=signR*max(zero,min(dabs(slopeR), &
             !                                    signR*slopeL))
             !case(3)
             !  ! woodward
             !  slope(iw,idim)=two*signR*max(zero,min(dabs(slopeR), &
             !                     signR*slopeL,signR*half*slopeC))
             !case(4)
             !  ! koren
             !  slope(iw,idim)=signR*max(zero,min(two*signR*slopeL, &
             !   (dabs(slopeR)+two*slopeL*signR)*third,two*dabs(slopeR)))
             !case default
               slope(iw,idim)=signc*max(zero,min(dabs(slopec), &
                                 signc*slopel,signc*sloper))
             !end select
          end do
       end do
       ! cell-centered coordinates of coarse grid point
       !^D&xCo^D=xCo({ixCo^DD},^D)
       {do ix^db=ixfi^db,ixfi^db+1 \}
          ! cell-centered coordinates of fine grid point
          !^D&xFi^D=xFi({ix^DD},^D)
          if(slab_uniform) then
            ! normalized distance between fine/coarse cell center
            ! in coarse cell: ranges from -0.5 to 0.5 in each direction
            ! (origin is coarse cell center)
            ! hence this is +1/4 or -1/4 on cartesian mesh
            !eta^D=(xFi^D-xCo^D)*invdxCo^D;
            eta^d=0.5d0*(dble(ix^d-ixfi^d)-0.5d0);
          else
            {! forefactor is -0.5d0 when ix=ixFi and +0.5d0 for ixFi+1
            eta^d=(dble(ix^d-ixfi^d)-0.5d0)*(one-sfi%dvolume(ix^dd) &
                  /sum(sfi%dvolume(ixfi^d:ixfi^d+1^d%ix^dd)))  \}
          end if
          wfi(ix^d,1:nw) = wco(ixco^d,1:nw) &
                                + {(slope(1:nw,^d)*eta^d)+}
       {end do\}
    {end do\}
    if(stagger_grid) then
      call already_fine(sfi,igridfi,fine_^l)
      call prolong_2nd_stg(sco,sfi,ixco^l,ixm^ll,dxco^d,xcomin^d,dxfi^d,xfimin^d,.false.,fine_^l)
    end if
  
    ! Convert child cells from prolongation space to conserved.
    block=>ps(igridfi)
    if(associated(phys_from_prolong) .and. associated(phys_wb_prolong)) then
      ! EoS-aware interpolation produced (rho, v, T) in the child cells.
      ! WB is also active: re-balance pressure via the HSE recurrence so
      ! the prolonged state preserves the discrete WB reference exactly.
      !   (rho,v,T) -[from_prolong]-> (rho,m,E) -[to_primitive]-> (rho,v,p)
      !   -[wb_prolong]-> (rho_eq, v, p_eq) -[to_conserved]-> (rho_eq, m, E_eq)
      call phys_from_prolong(ixg^ll,ixm^ll,wfi,sfi%x)
      call phys_to_primitive(ixg^ll,ixm^ll,wfi,sfi%x)
      call phys_wb_prolong(ixg^ll,ixm^ll,wfi,sfi%x)
      call phys_to_conserved(ixg^ll,ixm^ll,wfi,sfi%x)
    else if(associated(phys_from_prolong)) then
      ! EoS-aware only: (rho, v, T) -> (rho, rho*v, E) using EoS tables
      call phys_from_prolong(ixg^ll,ixm^ll,wfi,sfi%x)
    else if(associated(phys_wb_prolong)) then
      ! WB only: rebuild p from HSE recurrence using interpolated T = p/rho
      if(.not. prolongprimitive .and. .not. associated(phys_to_prolong)) then
        call phys_to_primitive(ixg^ll,ixm^ll,wfi,sfi%x)
      end if
      call phys_wb_prolong(ixg^ll,ixm^ll,wfi,sfi%x)
      call phys_to_conserved(ixg^ll,ixm^ll,wfi,sfi%x)
    else
      if(prolongprimitive) call phys_to_conserved(ixg^ll,ixm^ll,wfi,sfi%x)
    end if
 
    ! Small-value repair runs here, after the conversion above, not before it. wFi is in
    ! EoS prolongation space until then -- (rho, v, T), with T in the e_ slot and velocity in
    ! the momentum slots -- whereas phys_handle_small_values is told the layout by
    ! prolongprimitive, which is .false. here. Called earlier it therefore read T as total
    ! energy and formed 0.5*v^2/rho instead of 0.5*m^2/rho, and with
    ! small_values_method='replace' overwrote cells on the strength of that misreading.
    ! Symptom: only the energy slot of newly refined blocks was destroyed (rho/m/Ne within 4%
    ! of untouched blocks, e 5.7e6x too large), dt collapsing 1.9E-04 -> 2.7E-20 at the first
    ! regrid. Invisible whenever fix_small_values=.false., because that skips this call.
    if(fix_small_values) call phys_handle_small_values(prolongprimitive,wfi,sfi%x,ixg^ll,ixm^ll,'prolong_2nd')
    end associate
  
  end subroutine prolong_2nd
  
  !> do 1st order prolongation
  subroutine prolong_1st(wCo,ixCo^L,wFi,xFi)
    use mod_global_parameters
  
    integer, intent(in) :: ixco^l
    double precision, intent(in) :: wco(ixg^t,nw), xfi(ixg^t,1:ndim)
    double precision, intent(out) :: wfi(ixg^t,nw)
  
    integer :: ixco^d, ixfi^d, iw
    integer :: ixfi^l
  
    {do ixco^db = ixco^lim^db
       ixfi^db=2*(ixco^db-ixcomin^db)+ixmlo^db\}
       forall(iw=1:nw) wfi(ixfi^d:ixfi^d+1,iw)=wco(ixco^d,iw)
    {end do\}
  
  end subroutine prolong_1st
 
end module mod_refine