MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_errest.t
Go to the documentation of this file.
1 module mod_errest
2  use mod_comm_lib, only: mpistop
3  implicit none
4  private
5 
6  public :: errest, forcedrefine_grid_io
7 
8 contains
9 
10  !> Do all local error estimation which determines (de)refinement
11  subroutine errest
12  use mod_forest, only: refine, buffer
14 
15  integer :: igrid, iigrid, ixcog^l
16  double precision :: factor
17  logical, dimension(:,:), allocatable :: refine2
18 
19  if (igridstail==0) return
20 
21  select case (refine_criterion)
22  case (1)
23  ! all refinement solely based on user routine usr_refine_grid
24  case (2)
25  ! Error estimation is based on Lohner's original scheme
26  !$OMP PARALLEL DO PRIVATE(igrid)
27  do iigrid=1,igridstail; igrid=igrids(iigrid);
28  block=>ps(igrid)
29  call lohner_orig_grid(igrid)
30  end do
31  !$OMP END PARALLEL DO
32  case (3)
33  ! Error estimation is based on Lohner's scheme
34  !$OMP PARALLEL DO PRIVATE(igrid)
35  do iigrid=1,igridstail; igrid=igrids(iigrid);
36  block=>ps(igrid)
37  call lohner_grid(igrid)
38  end do
39  !$OMP END PARALLEL DO
40  case default
41  call mpistop("Unknown error estimator")
42  end select
43 
44  ! enforce additional refinement on e.g. coordinate and/or time info here
45  if (nbufferx^d/=0|.or.) then
46  allocate(refine2(max_blocks,npe))
47  call mpi_allreduce(refine,refine2,max_blocks*npe,mpi_logical,mpi_lor, &
48  icomm,ierrmpi)
49  refine=refine2
50  end if
51  !$OMP PARALLEL DO PRIVATE(igrid)
52  do iigrid=1,igridstail; igrid=igrids(iigrid);
53  block=>ps(igrid)
54  call forcedrefine_grid(igrid,ps(igrid)%w)
55  end do
56  !$OMP END PARALLEL DO
57 
58  if (nbufferx^d/=0|.or.) &
59  buffer=.false.
60 
61  end subroutine errest
62 
63  subroutine lohner_grid(igrid)
65  use mod_forest, only: coarsen, refine
66  use mod_physics, only: phys_energy
68 
69  integer, intent(in) :: igrid
70 
71  integer :: iflag, idims, idims2, level
72  integer :: ix^L, hx^L, jx^L, h2x^L, j2x^L, ix^D
73  double precision :: epsilon, threshold, wtol(1:nw), xtol(1:ndim)
74  double precision, dimension(ixM^T) :: numerator, denominator, error
75  double precision, dimension(ixG^T) :: tmp, tmp1, tmp2
76  double precision :: w(ixG^T,1:nw)
77  logical, dimension(ixG^T) :: refineflag, coarsenflag
78 
79  epsilon = 1.0d-6
80  level = node(plevel_,igrid)
81  ix^l=ixm^ll^ladd1;
82 
83  error=zero
84 
85  w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
86 
87  if(b0field .and. allocated(iw_mag)) then
88  if(phys_energy) &
89  w(ixg^t,iw_e)=w(ixg^t,iw_e)+0.5d0*sum(ps(igrid)%B0(ixg^t,:,0)**2,dim=ndim+1) &
90  + sum(w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=ndim+1)
91  w(ixg^t,iw_mag(:))=w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
92  end if
93 
94  do iflag=1,nw+1
95 
96  if(w_refine_weight(iflag)==0.d0) cycle
97  numerator=zero
98 
99  if (iflag > nw) then
100  if (.not. associated(usr_var_for_errest)) then
101  call mpistop("usr_var_for_errest not defined")
102  else
103  call usr_var_for_errest(ixg^ll,ixg^ll,iflag,ps(igrid)%w,ps(igrid)%x,tmp1)
104  end if
105  end if
106 
107  do idims=1,ndim
108  hx^l=ix^l-kr(^d,idims);
109  jx^l=ix^l+kr(^d,idims);
110  if (iflag<=nw) then
111  if (logflag(iflag)) then
112  tmp(ix^s)=dlog10(w(jx^s,iflag))-dlog10(w(hx^s,iflag))
113  else
114  tmp(ix^s)=w(jx^s,iflag)-w(hx^s,iflag)
115  end if
116  else
117  if (logflag(iflag)) then
118  tmp(ix^s)=dlog10(tmp1(jx^s))-dlog10(tmp1(hx^s))
119  else
120  tmp(ix^s)=tmp1(jx^s)-tmp1(hx^s)
121  end if
122  end if
123  do idims2=1,ndim
124  h2x^l=ixm^ll-kr(^d,idims2);
125  j2x^l=ixm^ll+kr(^d,idims2);
126  numerator=numerator+(tmp(j2x^s)-tmp(h2x^s))**2
127  end do
128  end do
129  denominator=zero
130  do idims=1,ndim
131  if (iflag<=nw) then
132  if (logflag(iflag)) then
133  tmp=dabs(dlog10(w(ixg^t,iflag)))
134  else
135  tmp=dabs(w(ixg^t,iflag))
136  end if
137  else
138  if (logflag(iflag)) then
139  tmp=dabs(dlog10(tmp1(ixg^t)))
140  else
141  tmp=dabs(tmp1(ixg^t))
142  end if
143  end if
144  hx^l=ix^l-kr(^d,idims);
145  jx^l=ix^l+kr(^d,idims);
146  tmp2(ix^s)=tmp(jx^s)+tmp(hx^s)
147  hx^l=ixm^ll-2*kr(^d,idims);
148  jx^l=ixm^ll+2*kr(^d,idims);
149  if (iflag<=nw) then
150  if (logflag(iflag)) then
151  tmp(ixm^t)=dabs(dlog10(w(jx^s,iflag))&
152  -dlog10(w(ixm^t,iflag))) &
153  +dabs(dlog10(w(ixm^t,iflag))&
154  -dlog10(w(hx^s,iflag)))
155  else
156  tmp(ixm^t)=dabs(w(jx^s,iflag)-w(ixm^t,iflag)) &
157  +dabs(w(ixm^t,iflag)-w(hx^s,iflag))
158  end if
159  else
160  if (logflag(iflag)) then
161  tmp(ixm^t)=dabs(dlog10(tmp1(jx^s))-dlog10(tmp1(ixm^t))) &
162  +dabs(dlog10(tmp1(ixm^t))-dlog10(tmp1(hx^s)))
163  else
164  tmp(ixm^t)=dabs(tmp1(jx^s)-tmp1(ixm^t)) &
165  +dabs(tmp1(ixm^t)-tmp1(hx^s))
166  end if
167  end if
168  do idims2=1,ndim
169  h2x^l=ixm^ll-kr(^d,idims2);
170  j2x^l=ixm^ll+kr(^d,idims2);
171  denominator=denominator &
172  +(tmp(ixm^t)+amr_wavefilter(level)*(tmp2(j2x^s)+tmp2(h2x^s)))**2
173  end do
174  end do
175  error=error+w_refine_weight(iflag)*dsqrt(numerator/max(denominator,epsilon))
176  end do
177 
178  refineflag=.false.
179  coarsenflag=.false.
180  threshold=refine_threshold(level)
181  {do ix^db=ixmlo^db,ixmhi^db\}
182 
183  if (associated(usr_refine_threshold)) then
184  wtol(1:nw) = w(ix^d,1:nw)
185  xtol(1:ndim) = ps(igrid)%x(ix^d,1:ndim)
186  call usr_refine_threshold(wtol, xtol, threshold, global_time, level)
187  end if
188 
189  if (error(ix^d) >= threshold) then
190  refineflag(ix^d) = .true.
191  else if (error(ix^d) <= derefine_ratio(level)*threshold) then
192  coarsenflag(ix^d) = .true.
193  end if
194  {end do\}
195 
196  if (any(refineflag(ixm^t)).and.level<refine_max_level) refine(igrid,mype)=.true.
197  if (all(coarsenflag(ixm^t)).and.level>1) coarsen(igrid,mype)=.true.
198 
199  end subroutine lohner_grid
200 
201  subroutine lohner_orig_grid(igrid)
203  use mod_forest, only: coarsen, refine
205 
206  integer, intent(in) :: igrid
207 
208  integer :: iflag, idims, level
209  integer :: ix^L, hx^L, jx^L, ix^D
210  double precision :: epsilon, threshold, wtol(1:nw), xtol(1:ndim)
211  double precision, dimension(ixM^T) :: numerator, denominator, error
212  double precision, dimension(ixG^T) :: dp, dm, dref, tmp1
213  logical, dimension(ixG^T) :: refineflag, coarsenflag
214 
215  epsilon=1.0d-6
216  level=node(plevel_,igrid)
217  ix^l=ixm^ll;
218 
219  error=zero
220  do iflag=1,nw+1
221  if(w_refine_weight(iflag)==0.d0) cycle
222  numerator=zero
223  denominator=zero
224 
225  if (iflag > nw) then
226  if (.not. associated(usr_var_for_errest)) then
227  call mpistop("usr_var_for_errest not defined")
228  else
229  call usr_var_for_errest(ixg^ll,ixg^ll,iflag,ps(igrid)%w,ps(igrid)%x,tmp1)
230  end if
231  end if
232 
233  do idims=1,ndim
234  hx^l=ix^l-kr(^d,idims);
235  jx^l=ix^l+kr(^d,idims);
236  if (iflag<=nw) then
237  if (logflag(iflag)) then
238  dp(ix^s)=dlog10(ps(igrid)%w(jx^s,iflag))-dlog10(ps(igrid)%w(ix^s,iflag))
239  dm(ix^s)=dlog10(ps(igrid)%w(ix^s,iflag))-dlog10(ps(igrid)%w(hx^s,iflag))
240  dref(ixm^t)=dabs(dlog10(ps(igrid)%w(jx^s,iflag)))&
241  + 2.0d0 * dabs(dlog10(ps(igrid)%w(ixm^t,iflag))) &
242  + dabs(dlog10(ps(igrid)%w(hx^s,iflag)))
243  else
244  dp(ix^s)=ps(igrid)%w(jx^s,iflag)-ps(igrid)%w(ix^s,iflag)
245  dm(ix^s)=ps(igrid)%w(ix^s,iflag)-ps(igrid)%w(hx^s,iflag)
246  dref(ixm^t)=dabs(ps(igrid)%w(jx^s,iflag))+2.0d0*dabs(ps(igrid)%w(ixm^t,iflag)) &
247  +dabs(ps(igrid)%w(hx^s,iflag))
248  end if
249  else
250  if (logflag(iflag)) then
251  dp(ix^s)=dlog10(tmp1(jx^s))-dlog10(tmp1(ix^s))
252  dm(ix^s)=dlog10(tmp1(ix^s))-dlog10(tmp1(hx^s))
253  dref(ix^s)=dabs(dlog10(tmp1(jx^s)))&
254  + 2.0d0 * dabs(dlog10(tmp1(ix^s))) &
255  + dabs(dlog10(tmp1(hx^s)))
256  else
257  dp(ix^s)=tmp1(jx^s)-tmp1(ix^s)
258  dm(ix^s)=tmp1(ix^s)-tmp1(hx^s)
259  dref(ix^s)=dabs(tmp1(jx^s))+2.0d0*dabs(tmp1(ix^s)) &
260  +dabs(tmp1(hx^s))
261  end if
262  end if
263 
264  numerator(ixm^t)=numerator+(dp(ixm^t)-dm(ixm^t))**2
265  denominator(ixm^t)=denominator &
266  + (dabs(dp(ixm^t)) + dabs(dm(ixm^t)) + amr_wavefilter(level)*dref(ixm^t))**2
267 
268  end do
269  error=error+w_refine_weight(iflag)*dsqrt(numerator/max(denominator,epsilon))
270  end do
271 
272  refineflag=.false.
273  coarsenflag=.false.
274 
275  threshold=refine_threshold(level)
276  {do ix^db=ixmlo^db,ixmhi^db\}
277 
278  if (associated(usr_refine_threshold)) then
279  wtol(1:nw) = ps(igrid)%w(ix^d,1:nw)
280  xtol(1:ndim) = ps(igrid)%x(ix^d,1:ndim)
281  call usr_refine_threshold(wtol, xtol, threshold, global_time, level)
282  end if
283 
284  if (error(ix^d) >= threshold) then
285  refineflag(ix^d) = .true.
286  else if (error(ix^d) <= derefine_ratio(level)*threshold) then
287  coarsenflag(ix^d) = .true.
288  end if
289  {end do\}
290 
291  if (any(refineflag(ixm^t)).and.level<refine_max_level) refine(igrid,mype)=.true.
292  if (all(coarsenflag(ixm^t)).and.level>1) coarsen(igrid,mype)=.true.
293 
294  end subroutine lohner_orig_grid
295 
296  subroutine forcedrefine_grid(igrid,w)
298  use mod_forest, only: coarsen, refine, buffer
300 
301  integer, intent(in) :: igrid
302  double precision, intent(in) :: w(ixG^T,nw)
303 
304  integer :: level
305  integer :: my_refine, my_coarsen
306  double precision :: qt
307  logical, dimension(ixG^T) :: refineflag
308 
309  level=node(plevel_,igrid)
310 
311  ! initialize to 0
312  my_refine = 0
313  my_coarsen = 0
314 
315  if (time_advance) then
316  qt=global_time+dt
317  else
318  qt=global_time
319  end if
320 
321  if (associated(usr_refine_grid)) then
322  call usr_refine_grid(igrid,level,ixg^ll,ixm^ll,qt,w,ps(igrid)%x, &
323  my_refine,my_coarsen)
324  end if
325 
326  if (my_coarsen==1) then
327  if (level>1) then
328  refine(igrid,mype)=.false.
329  coarsen(igrid,mype)=.true.
330  else
331  refine(igrid,mype)=.false.
332  coarsen(igrid,mype)=.false.
333  end if
334  end if
335 
336  if (my_coarsen==-1)then
337  coarsen(igrid,mype)=.false.
338  end if
339 
340  if (my_refine==1) then
341  if (level<refine_max_level) then
342  refine(igrid,mype)=.true.
343  coarsen(igrid,mype)=.false.
344  else
345  refine(igrid,mype)=.false.
346  coarsen(igrid,mype)=.false.
347  end if
348  end if
349 
350  if (my_refine==-1) then
351  refine(igrid,mype)=.false.
352  end if
353 
354  if (nbufferx^d/=0|.or.) then
355  if (refine(igrid,mype) .and. .not.buffer(igrid,mype)) then
356  refineflag(ixm^t)=.true.
357  call refinebuffer(igrid,refineflag)
358  end if
359  end if
360 
361  end subroutine forcedrefine_grid
362 
363  subroutine forcedrefine_grid_io(igrid,w)
364  use mod_forest, only: coarsen, refine
366 
367  integer, intent(in) :: igrid
368  double precision, intent(in) :: w(ixg^t,nw)
369 
370  integer :: level, my_levmin, my_levmax
371  logical, dimension(ixG^T) :: refineflag
372 
373  level=node(plevel_,igrid)
374 
375  if (level_io > 0) then
376  my_levmin = level_io
377  my_levmax = level_io
378  else
379  my_levmin = max(1,level_io_min)
380  my_levmax = min(refine_max_level,level_io_max)
381  end if
382 
383  if (level>my_levmax) then
384  refine(igrid,mype)=.false.
385  coarsen(igrid,mype)=.true.
386  elseif (level<my_levmin) then
387  refine(igrid,mype)=.true.
388  coarsen(igrid,mype)=.false.
389  end if
390 
391  if (level==my_levmin .or. level==my_levmax) then
392  refine(igrid,mype)=.false.
393  coarsen(igrid,mype)=.false.
394  end if
395 
396  if(refine(igrid,mype).and.level>=refine_max_level)refine(igrid,mype)=.false.
397  if(coarsen(igrid,mype).and.level<=1)coarsen(igrid,mype)=.false.
398 
399  end subroutine forcedrefine_grid_io
400 
401  subroutine refinebuffer(igrid,refineflag)
402  use mod_forest, only: refine, buffer
404 
405  integer, intent(in) :: igrid
406  logical, dimension(ixG^T), intent(in) :: refineflag
407 
408  integer :: ishiftbuf^D, i^D, ix^L, ineighbor, ipe_neighbor, level
409 
410  ishiftbuf^d=ixmhi^d-ixmlo^d-nbufferx^d+1;
411  {do i^db=-1,1\}
412  ixmin^d=max(ixmlo^d,ixmlo^d+i^d*ishiftbuf^d);
413  ixmax^d=min(ixmhi^d,ixmhi^d+i^d*ishiftbuf^d);
414  if (ixmax^d<ixmin^d|.or.) cycle
415  if (any(refineflag(ix^s))) then
416  select case (neighbor_type(i^d,igrid))
417  case (neighbor_coarse)
418  ineighbor=neighbor(1,i^d,igrid)
419  ipe_neighbor=neighbor(2,i^d,igrid)
420  if (.not.refine(ineighbor,ipe_neighbor)) then
421  buffer(ineighbor,ipe_neighbor)=.true.
422  refine(ineighbor,ipe_neighbor)=.true.
423  end if
424  case (neighbor_sibling)
425  level=node(plevel_,igrid)
426  if (level<refine_max_level) then
427  ineighbor=neighbor(1,i^d,igrid)
428  ipe_neighbor=neighbor(2,i^d,igrid)
429  if (.not.refine(ineighbor,ipe_neighbor)) then
430  buffer(ineighbor,ipe_neighbor)=.true.
431  refine(ineighbor,ipe_neighbor)=.true.
432  end if
433  end if
434  end select
435  end if
436  {end do\}
437 
438  end subroutine refinebuffer
439 
440 end module mod_errest
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
subroutine refinebuffer(igrid, refineflag)
Definition: mod_errest.t:402
subroutine forcedrefine_grid(igrid, w)
Definition: mod_errest.t:297
subroutine lohner_grid(igrid)
Definition: mod_errest.t:64
subroutine lohner_orig_grid(igrid)
Definition: mod_errest.t:202
subroutine, public errest
Do all local error estimation which determines (de)refinement.
Definition: mod_errest.t:12
subroutine, public forcedrefine_grid_io(igrid, w)
Definition: mod_errest.t:364
Module with basic grid data structures.
Definition: mod_forest.t:2
logical, dimension(:,:), allocatable, save refine
Definition: mod_forest.t:70
logical, dimension(:,:), allocatable, save buffer
Definition: mod_forest.t:70
logical, dimension(:,:), allocatable, save coarsen
AMR flags and grids-in-use identifier per processor (igrid,ipe)
Definition: mod_forest.t:70
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
double precision global_time
The global simulation time.
integer, dimension(3, 3) kr
Kronecker delta tensor.
logical, dimension(:), allocatable logflag
double precision, dimension(:), allocatable amr_wavefilter
refinement: lohner estimate wavefilter setting
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer, parameter plevel_
integer, dimension(:), allocatable, parameter d
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
integer npe
The number of MPI tasks.
logical b0field
split magnetic field as background B0 field
integer nbufferx
Number of cells as buffer zone.
double precision, dimension(:), allocatable w_refine_weight
Weights of variables used to calculate error for mesh refinement.
double precision, dimension(:), allocatable refine_threshold
Error tolerance for refinement decision.
integer refine_max_level
Maximal number of AMR levels.
double precision, dimension(:), allocatable derefine_ratio
integer max_blocks
The maximum number of grid blocks in a processor.
integer, dimension(:,:), allocatable node
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
logical phys_energy
Solve energy equation or not.
Definition: mod_physics.t:27
Module with all the methods that users can customize in AMRVAC.
procedure(a_refine_threshold), pointer usr_refine_threshold
procedure(refine_grid), pointer usr_refine_grid
procedure(var_for_errest), pointer usr_var_for_errest