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