mod_errest.t

Go to the documentation of this file.

module mod_errest
  use mod_comm_lib, only: mpistop
  implicit none
  private
 
  public :: errest, forcedrefine_grid_io 
 
contains
 
  !> Do all local error estimation which determines (de)refinement
  subroutine errest
    use mod_forest, only: refine, buffer
    use mod_global_parameters
 
    double precision :: factor
    integer :: igrid, iigrid, ixcog^l
    logical, dimension(:,:), allocatable :: refine2
 
    if (igridstail==0) return
 
    select case (refine_criterion)
    case (1) 
       ! all refinement solely based on user routine usr_refine_grid
    case (2)
       ! Error estimation is based on Lohner's original scheme
    !$OMP PARALLEL DO PRIVATE(igrid)
       do iigrid=1,igridstail; igrid=igrids(iigrid);
          block=>ps(igrid)
          call lohner_orig_grid(igrid)
       end do
    !$OMP END PARALLEL DO
    case (3)
       ! Error estimation is based on Lohner's scheme
    !$OMP PARALLEL DO PRIVATE(igrid)
       do iigrid=1,igridstail; igrid=igrids(iigrid);
          block=>ps(igrid)
          call lohner_grid(igrid)
       end do
    !$OMP END PARALLEL DO
    case default
       call mpistop("Unknown error estimator")
    end select
 
    ! enforce additional refinement on e.g. coordinate and/or time info here
    if (nbufferx^d/=0|.or.) then
       allocate(refine2(max_blocks,npe))
       call mpi_allreduce(refine,refine2,max_blocks*npe,mpi_logical,mpi_lor, &
                          icomm,ierrmpi)
       refine=refine2
    end if
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       block=>ps(igrid)
       call forcedrefine_grid(igrid,ps(igrid)%w)
    end do
    !$OMP END PARALLEL DO
 
    if (nbufferx^d/=0|.or.) &
    buffer=.false.
 
  end subroutine errest
 
  subroutine lohner_grid(igrid)
    use mod_usr_methods, only: usr_var_for_errest, usr_refine_threshold
    use mod_forest, only: coarsen, refine
    use mod_physics, only: phys_energy
    use mod_global_parameters
 
    integer, intent(in) :: igrid
 
    double precision                   :: epsilon, threshold, wtol(1:nw), xtol(1:ndim)
    double precision, dimension(ixM^T) :: numerator, denominator, error
    double precision, dimension(ixG^T) :: tmp, tmp1, tmp2
    double precision                   :: w(ixg^t,1:nw)
    integer                            :: iflag, idims, idims2, level
    integer                            :: ix^l, hx^l, jx^l, h2x^l, j2x^l, ix^d
    logical, dimension(ixG^T)          :: refineflag, coarsenflag
 
    epsilon = 1.0d-6
    level   = node(plevel_,igrid)
    ix^l=ixm^ll^ladd1;
 
    error=zero
 
    w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
 
    if(b0field .and. allocated(iw_mag)) then
      if(phys_energy) &
        w(ixg^t,iw_e)=w(ixg^t,iw_e)+0.5d0*sum(ps(igrid)%B0(ixg^t,:,0)**2,dim=ndim+1) &
                        + sum(w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=ndim+1)
      w(ixg^t,iw_mag(:))=w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
    end if
 
    do iflag=1,nw+1
 
       if(w_refine_weight(iflag)==0.d0) cycle
       numerator=zero
 
       if (iflag > nw) then
          if (.not. associated(usr_var_for_errest)) then
             call mpistop("usr_var_for_errest not defined")
          else
             call usr_var_for_errest(ixg^ll,ixg^ll,iflag,ps(igrid)%w,ps(igrid)%x,tmp1)
          end if
       end if
 
       do idims=1,ndim
          hx^l=ix^l-kr(^d,idims);
          jx^l=ix^l+kr(^d,idims);
          if (iflag<=nw) then
            if (logflag(iflag)) then
              tmp(ix^s)=dlog10(w(jx^s,iflag))-dlog10(w(hx^s,iflag))
            else
              tmp(ix^s)=w(jx^s,iflag)-w(hx^s,iflag)
            end if
          else
            if (logflag(iflag)) then
              tmp(ix^s)=dlog10(tmp1(jx^s))-dlog10(tmp1(hx^s))
            else
              tmp(ix^s)=tmp1(jx^s)-tmp1(hx^s)
            end if
          end if
          do idims2=1,ndim
             h2x^l=ixm^ll-kr(^d,idims2);
             j2x^l=ixm^ll+kr(^d,idims2);
             numerator=numerator+(tmp(j2x^s)-tmp(h2x^s))**2
          end do
       end do
       denominator=zero
       do idims=1,ndim
          if (iflag<=nw) then
             if (logflag(iflag)) then
              ! With logflag the differenced quantity is u = log10(w), for which
              ! Lohner's eps*|u| filter is not meaningful: it is a relative floor,
              ! scale-invariant for a linear variable but not for a log, where a
              ! change of unit is additive and |u| becomes the distance from
              ! w_code = 1. The fractional character is already carried by
              ! d(log w), so the floor is constant and amr_wavefilter is in dex.
              tmp=one
             else
              tmp=dabs(w(ixg^t,iflag))
             end if
          else
             if (logflag(iflag)) then
              tmp=one                      ! dex floor, as in the w-slot branch above
             else
              tmp=dabs(tmp1(ixg^t))
             end if
          end if
          hx^l=ix^l-kr(^d,idims);
          jx^l=ix^l+kr(^d,idims);
          tmp2(ix^s)=tmp(jx^s)+tmp(hx^s)
          hx^l=ixm^ll-2*kr(^d,idims);
          jx^l=ixm^ll+2*kr(^d,idims);
          if (iflag<=nw) then
            if (logflag(iflag)) then
              tmp(ixm^t)=dabs(dlog10(w(jx^s,iflag))&
                         -dlog10(w(ixm^t,iflag))) &
                         +dabs(dlog10(w(ixm^t,iflag))&
                         -dlog10(w(hx^s,iflag)))
            else
               tmp(ixm^t)=dabs(w(jx^s,iflag)-w(ixm^t,iflag)) &
                          +dabs(w(ixm^t,iflag)-w(hx^s,iflag))
            end if
          else
            if (logflag(iflag)) then
              tmp(ixm^t)=dabs(dlog10(tmp1(jx^s))-dlog10(tmp1(ixm^t))) &
                        +dabs(dlog10(tmp1(ixm^t))-dlog10(tmp1(hx^s)))
            else
               tmp(ixm^t)=dabs(tmp1(jx^s)-tmp1(ixm^t)) &
                          +dabs(tmp1(ixm^t)-tmp1(hx^s))
            end if
          end if
          do idims2=1,ndim
             h2x^l=ixm^ll-kr(^d,idims2);
             j2x^l=ixm^ll+kr(^d,idims2);
             denominator=denominator &
                        +(tmp(ixm^t)+amr_wavefilter(level)*(tmp2(j2x^s)+tmp2(h2x^s)))**2
          end do
       end do
       error=error+w_refine_weight(iflag)*dsqrt(numerator/max(denominator,epsilon))
    end do
    
    refineflag=.false.
    coarsenflag=.false.
    threshold=refine_threshold(level)
    {do ix^db=ixmlo^db,ixmhi^db\}
    
       if (associated(usr_refine_threshold)) then
          wtol(1:nw)   = w(ix^d,1:nw)
          xtol(1:ndim) = ps(igrid)%x(ix^d,1:ndim)
          call usr_refine_threshold(wtol, xtol, threshold, global_time, level)
       end if
    
       if (error(ix^d) >= threshold) then
          refineflag(ix^d) = .true.
       else if (error(ix^d) <= derefine_ratio(level)*threshold) then
          coarsenflag(ix^d) = .true.
       end if
    {end do\}
 
    if (any(refineflag(ixm^t)).and.level<refine_max_level) refine(igrid,mype)=.true.
    if (all(coarsenflag(ixm^t)).and.level>1) coarsen(igrid,mype)=.true.
 
  end subroutine lohner_grid
 
  subroutine lohner_orig_grid(igrid)
    use mod_usr_methods, only: usr_var_for_errest, usr_refine_threshold
    use mod_forest, only: coarsen, refine
    use mod_global_parameters
 
    integer, intent(in) :: igrid
 
    double precision :: epsilon, threshold, wtol(1:nw), xtol(1:ndim)
    double precision, dimension(ixM^T) :: numerator, denominator, error
    double precision, dimension(ixG^T) :: dp, dm, dref, tmp1
    integer :: iflag, idims, level
    integer :: ix^l, hx^l, jx^l, ix^d
    logical, dimension(ixG^T) :: refineflag, coarsenflag
 
    epsilon=1.0d-6
    level=node(plevel_,igrid)
    ix^l=ixm^ll;
 
    error=zero
    do iflag=1,nw+1
       if(w_refine_weight(iflag)==0.d0) cycle
       numerator=zero
       denominator=zero
 
       if (iflag > nw) then
          if (.not. associated(usr_var_for_errest)) then
             call mpistop("usr_var_for_errest not defined")
          else
             call usr_var_for_errest(ixg^ll,ixg^ll,iflag,ps(igrid)%w,ps(igrid)%x,tmp1)
          end if
       end if
 
       do idims=1,ndim
          hx^l=ix^l-kr(^d,idims);
          jx^l=ix^l+kr(^d,idims);
          if (iflag<=nw) then
            if (logflag(iflag)) then
              dp(ix^s)=dlog10(ps(igrid)%w(jx^s,iflag))-dlog10(ps(igrid)%w(ix^s,iflag))
              dm(ix^s)=dlog10(ps(igrid)%w(ix^s,iflag))-dlog10(ps(igrid)%w(hx^s,iflag))
              dref(ixm^t)=dabs(dlog10(ps(igrid)%w(jx^s,iflag)))&
                         + 2.0d0 * dabs(dlog10(ps(igrid)%w(ixm^t,iflag))) &
                         + dabs(dlog10(ps(igrid)%w(hx^s,iflag)))
            else
              dp(ix^s)=ps(igrid)%w(jx^s,iflag)-ps(igrid)%w(ix^s,iflag)
              dm(ix^s)=ps(igrid)%w(ix^s,iflag)-ps(igrid)%w(hx^s,iflag)
              dref(ixm^t)=dabs(ps(igrid)%w(jx^s,iflag))+2.0d0*dabs(ps(igrid)%w(ixm^t,iflag)) &
                          +dabs(ps(igrid)%w(hx^s,iflag))
            end if
          else
            if (logflag(iflag)) then
              dp(ix^s)=dlog10(tmp1(jx^s))-dlog10(tmp1(ix^s))
              dm(ix^s)=dlog10(tmp1(ix^s))-dlog10(tmp1(hx^s))
              dref(ix^s)=dabs(dlog10(tmp1(jx^s)))&
                         + 2.0d0 * dabs(dlog10(tmp1(ix^s))) &
                         + dabs(dlog10(tmp1(hx^s)))
            else
              dp(ix^s)=tmp1(jx^s)-tmp1(ix^s)
              dm(ix^s)=tmp1(ix^s)-tmp1(hx^s)
              dref(ix^s)=dabs(tmp1(jx^s))+2.0d0*dabs(tmp1(ix^s)) &
                          +dabs(tmp1(hx^s))
            end if
          end if
 
          numerator(ixm^t)=numerator+(dp(ixm^t)-dm(ixm^t))**2
          denominator(ixm^t)=denominator &
               + (dabs(dp(ixm^t)) + dabs(dm(ixm^t)) + amr_wavefilter(level)*dref(ixm^t))**2
  
       end do
       error=error+w_refine_weight(iflag)*dsqrt(numerator/max(denominator,epsilon))
    end do
 
    refineflag=.false.
    coarsenflag=.false.
 
    threshold=refine_threshold(level)
    {do ix^db=ixmlo^db,ixmhi^db\}
 
       if (associated(usr_refine_threshold)) then
          wtol(1:nw)   = ps(igrid)%w(ix^d,1:nw)
          xtol(1:ndim) = ps(igrid)%x(ix^d,1:ndim)
          call usr_refine_threshold(wtol, xtol, threshold, global_time, level)
       end if
 
       if (error(ix^d) >= threshold) then
          refineflag(ix^d) = .true.
       else if (error(ix^d) <= derefine_ratio(level)*threshold) then
          coarsenflag(ix^d) = .true.
       end if
    {end do\}
 
    if (any(refineflag(ixm^t)).and.level<refine_max_level) refine(igrid,mype)=.true.
    if (all(coarsenflag(ixm^t)).and.level>1) coarsen(igrid,mype)=.true.
 
  end subroutine lohner_orig_grid
 
  subroutine forcedrefine_grid(igrid,w)
    use mod_usr_methods, only: usr_refine_grid
    use mod_forest, only: coarsen, refine, buffer
    use mod_global_parameters
 
    integer, intent(in) :: igrid
    double precision, intent(in) :: w(ixg^t,nw)
 
    double precision :: qt
    integer :: level
    integer :: my_refine, my_coarsen
    logical, dimension(ixG^T) :: refineflag
 
    level=node(plevel_,igrid)
 
    ! initialize to 0
    my_refine   = 0
    my_coarsen  = 0
 
    if (time_advance) then
       qt=global_time+dt
    else
       qt=global_time
    end if
 
    if (associated(usr_refine_grid)) then
       call usr_refine_grid(igrid,level,ixg^ll,ixm^ll,qt,w,ps(igrid)%x, &
            my_refine,my_coarsen)
    end if
 
    if (my_coarsen==1) then
       if (level>1) then
          refine(igrid,mype)=.false.
          coarsen(igrid,mype)=.true.
       else
          refine(igrid,mype)=.false.
          coarsen(igrid,mype)=.false.
       end if
    end if
 
    if (my_coarsen==-1)then
       coarsen(igrid,mype)=.false.
    end if
 
    if (my_refine==1) then
       if (level<refine_max_level) then
          refine(igrid,mype)=.true.
          coarsen(igrid,mype)=.false.
       else
          refine(igrid,mype)=.false.
          coarsen(igrid,mype)=.false.
       end if
    end if
 
    if (my_refine==-1) then
      refine(igrid,mype)=.false.
    end if
  
    if (nbufferx^d/=0|.or.) then
       if (refine(igrid,mype) .and. .not.buffer(igrid,mype)) then
          refineflag(ixm^t)=.true.
          call refinebuffer(igrid,refineflag)
       end if
    end if
  
  end subroutine forcedrefine_grid
  
  subroutine forcedrefine_grid_io(igrid,w)
    use mod_forest, only: coarsen, refine
    use mod_global_parameters
 
    integer, intent(in)          :: igrid
    double precision, intent(in) :: w(ixg^t,nw)
 
    logical, dimension(ixG^T) :: refineflag
    integer                   :: level, my_levmin, my_levmax
 
    level=node(plevel_,igrid)
 
    if (level_io > 0) then
       my_levmin = level_io
       my_levmax = level_io
    else
       my_levmin = max(1,level_io_min)
       my_levmax = min(refine_max_level,level_io_max)
    end if
 
    if (level>my_levmax) then
      refine(igrid,mype)=.false.
      coarsen(igrid,mype)=.true.
    elseif (level<my_levmin) then
      refine(igrid,mype)=.true.
      coarsen(igrid,mype)=.false.
    end if
 
    if (level==my_levmin .or. level==my_levmax) then
      refine(igrid,mype)=.false.
      coarsen(igrid,mype)=.false.
    end if
 
    if(refine(igrid,mype).and.level>=refine_max_level)refine(igrid,mype)=.false.
    if(coarsen(igrid,mype).and.level<=1)coarsen(igrid,mype)=.false.
 
  end subroutine forcedrefine_grid_io
 
  subroutine refinebuffer(igrid,refineflag)
    use mod_forest, only: refine, buffer
    use mod_global_parameters
 
    integer, intent(in) :: igrid
    logical, dimension(ixG^T), intent(in) :: refineflag
 
    integer :: ishiftbuf^d, i^d, ix^l, ineighbor, ipe_neighbor, level
 
    ishiftbuf^d=ixmhi^d-ixmlo^d-nbufferx^d+1;
    {do i^db=-1,1\}
       ixmin^d=max(ixmlo^d,ixmlo^d+i^d*ishiftbuf^d);
       ixmax^d=min(ixmhi^d,ixmhi^d+i^d*ishiftbuf^d);
       if (ixmax^d<ixmin^d|.or.) cycle
       if (any(refineflag(ix^s))) then
          select case (neighbor_type(i^d,igrid))
          case (neighbor_coarse)
             ineighbor=neighbor(1,i^d,igrid)
             ipe_neighbor=neighbor(2,i^d,igrid)
             if (.not.refine(ineighbor,ipe_neighbor)) then
                buffer(ineighbor,ipe_neighbor)=.true.
                refine(ineighbor,ipe_neighbor)=.true.
             end if
          case (neighbor_sibling)
             level=node(plevel_,igrid)
             if (level<refine_max_level) then
                ineighbor=neighbor(1,i^d,igrid)
                ipe_neighbor=neighbor(2,i^d,igrid)
                if (.not.refine(ineighbor,ipe_neighbor)) then
                   buffer(ineighbor,ipe_neighbor)=.true.
                   refine(ineighbor,ipe_neighbor)=.true.
                end if
             end if
          end select
       end if
    {end do\}
 
  end subroutine refinebuffer
 
end module mod_errest