errest.t

Go to the documentation of this file.

!> Do all local error estimation which determines (de)refinement
subroutine errest
  use mod_forest, only: refine, buffer
  use mod_global_parameters
 
  integer :: igrid, iigrid, ixCoG^L
  double precision :: factor
  logical, dimension(:,:), allocatable :: refine2
 
  if (igridstail==0) return
 
  select case (refine_criterion)
  case (0) 
     ! all refinement solely based on user routine usr_refine_grid
  case (1) 
     ! simply compare w_n-1 with w_n and trigger refinement on relative
     ! differences
  !$OMP PARALLEL DO PRIVATE(igrid)
     do iigrid=1,igridstail; igrid=igrids(iigrid);
        block=>ps(igrid)
        call compare1_grid(igrid,pso(igrid)%w,ps(igrid)%w)
     end do
  !$OMP END PARALLEL DO
  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
 
  integer                            :: iflag, idims, idims2, level
  integer                            :: ix^L, hx^L, jx^L, h2x^L, j2x^L, ix^D
  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)
  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) 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.0d0
        end do
     end do
     denominator=zero
     do idims=1,ndim
        if (iflag<=nw) then
           if (logflag(iflag)) then
            tmp=dabs(dlog10(w(ixg^t,iflag)))
           else
            tmp=dabs(w(ixg^t,iflag))
           end if
        else
           if (logflag(iflag)) then
            tmp=dabs(dlog10(tmp1(ixg^t)))
           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
 
  integer :: iflag, idims, level
  integer :: ix^L, hx^L, jx^L, ix^D
  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
  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.0d0
 
        denominator(ixm^t)=denominator &
             + (dabs(dp(ixm^t)) + dabs(dm(ixm^t)) + amr_wavefilter(level)*dref(ixm^t))**2.0d0
 
     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 compare1_grid(igrid,wold,w)
  use mod_usr_methods, only: usr_refine_threshold
  use mod_forest, only: coarsen, refine
  use mod_global_parameters
 
  integer, intent(in) :: igrid
  double precision, intent(in) :: wold(ixG^T,1:nw), w(ixG^T,1:nw)
 
  integer :: ix^D, iflag, level
  double precision :: epsilon, threshold, wtol(1:nw), xtol(1:ndim)
  double precision :: average, error
  double precision :: averages(nw)
  logical, dimension(ixG^T) :: refineflag, coarsenflag
 
  ! identify the points to be flagged in two steps:
  !  step I: compare w_n-1 with w_n solution, store w_for_refine in auxiliary
  !  step II: transfer w_for_refine from auxiliary to refine and coarsen
 
  epsilon=1.0d-6
 
  refineflag(ixm^t) = .false.
  coarsenflag(ixm^t) = .false.
  level=node(plevel_,igrid)
  threshold=refine_threshold(level)
  {do ix^db=ixmlo^db,ixmhi^db \}
     average=zero
     error=zero
     do iflag=1,nw+1
        if(w_refine_weight(iflag)==0) cycle
        averages(iflag) = w(ix^d,iflag)-wold(ix^d,iflag)
        average=average+w_refine_weight(iflag)*abs(averages(iflag))
        if (abs(wold(ix^d,iflag))<smalldouble)then
           error=error+w_refine_weight(iflag)* &
              abs(averages(iflag))/(abs(wold(ix^d,iflag))+epsilon)
        else
           error=error+w_refine_weight(iflag)* &
              abs(averages(iflag))/(abs(wold(ix^d,iflag)))
        end if
     end do
 
     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 >= threshold) then
        refineflag(ix^d) = .true.
     else if (error <= derefine_ratio(level)*threshold) then
        coarsenflag(ix^d) = .true.
     end if
  {end do\}
 
  if (any(refineflag(ixm^t))) then
     if (level<refine_max_level) refine(igrid,mype)=.true.
  end if
  if (time_advance) then
     if (all(coarsenflag(ixm^t)).and.level>1) coarsen(igrid,mype)=.true.
  end if
 
end subroutine compare1_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)
 
  integer :: level
  integer :: my_refine, my_coarsen
  double precision :: qt
  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
     if (refine_criterion==1) then
        qt=global_time+dt
     else
        qt=global_time
     end if
  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
  endif
 
  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)
 
  integer                   :: level, my_levmin, my_levmax
  logical, dimension(ixG^T) :: refineflag
 
  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