Open rjleveque opened 4 years ago
@mjberger and I also just discussed the possibility of specifying that the finest level(s) should be frozen beyond some time, after which re-gridding is never done for these level(s). This would be useful e.g. in GeoClaw simulations where the finest levels cover a fixed coastal region at very fine resolution, and frequent re-gridding is both expensive and unnecessary. This may take a bit of work to implement but should be possible.
See https://github.com/clawpack/geoclaw/issues/454 for a related suggestion in GeoClaw, eliminating the automatic creation of regions corresponding to each topo or dtopo file.
Another suggestion: I think in setrun.py, in addition to e.g.
amrdata.amr_levels_max = 3
we should also have something like:
amrdata.amr_levels_max_outside_flagregions = 2
as the default maximum refinement level in portions of the computational domain not covered by any flagregion.
This would be equivalent to specifying a flagregion with spatial extent covering the entire domain and with minlevel = 1, maxlevel = 2
(in the case shown above), but without having to specify a flagregion. This would be convenient and also perhaps clearer to the user.
Currently the default is that any level up to amr_levels_max
is allowed for any cell that is not in any flagregion (whether the cell is flagged or not then depends on other criteria such as flag2refine or Richardson).
Sometimes this is what one wants, but often not. In particular, if the finest levels are only intended to be used in a small part of the domain (e.g. one coastal site for a tsunami simulation) then there might be a flagregion that forces refinement to this level in a small region, but another flagregion would have to be introduced to forbid refinement to this level elsewhere, regardless of the size of wave.
For the former suggestion were you thinking of specifying a time at which a level would be frozen?
@mandli: Yes, it would be good to be able to specify a time at which each level is frozen. Perhaps we could allow this to also vary with level, with default time infinity, and with the requirement that frozen_time[L+1] <= frozen_time[L]
since if a level is frozen then all finer levels must also be frozen.
For the tsunami application I mentioned, we often have a finest level with minlevel==maxlevel
but that doesn't turn on until some time time t1
just before we expect the tsunami to arrive (based on coarse exploratory runs) and so it would be natural to freeze at time t1
.
However, one potentially tricky issue is that the regridding that introduces this finest grid will not in general happen exactly at time t1
, but only at the next regridding time after t1
. So it has to do this regridding to introduce the fine level in the first place before it is frozen, and hence it's not enough to check if t > frozen_time
when deciding whether to regrid this level, at least not if the user sets frozen_time = t1
. If the user is aware of this they could specify a frozen_time
somewhat larger than t1
, if they know how much larger it needs to be (to guarantee that at least one regridding happened between t1
and this time), or we could try to automate this by keeping track of whether it has already regridded after frozen_time
and only freezing it after the first regridding following frozen_time
.
Or since in general there would probably only be one such time for each run (or very few at least), we could force the code to adjust the time step if need to hit that time exactly (as we do for output times) and then force a regridding at exactly that time and freeze afterwards.
None of these is a perfect solution, but could be workable.
For discussion, fleshing out an idea we've discussed briefly in the past...
See also http://www.clawpack.org/refinement.html
Current flagging strategy in amrclaw:
(ignoring adjoint flagging to start)
setrun.py sets values:
If it is regrid time on level L, the following is done in flagger.f, looping over all patches on this level:
call flagregions2 to set some values to DONTFLAG if forbidden or DOFLAG if required to refine further, done as follows: If the cell center is in the union of some set of regions then the cell is marked
where the max is over these regions in both cases (yes, both are max). The cell remains UNSET if it does not lie in any region, or if
in which case refinement is neither required nor forbidden
Then buffer flagged points and cluster into rectanglar level L+1 patches.
If adjoint_flagging == True then the above is still done, but now:
Proposed modification (ignoring adjoint):
setrun.py sets values:
In each case, a scalar value means that value applies to all levels L, while an array such as
would mean use Richardson flagging on levels L = 1 and 2 but not on level 3 or above, and the tolerance could also vary by level.
The above algorithm would be unchanged except that on level L the values of flag_richardson, richardson_tol, etc. would be the ones specified for that level. This should be a very minor change.
Note that flag_richardson and flag2refine can each be True or False independent of the other. If both True, then if a cell is UNSET by regions we first check flag2refine and if it remains UNSET then it is checked by flag_richardson. So it will end up flagged if either criterion says it should be flagged (unless the regions already determined it as DOFLAG or DONTFLAG).
Adjoint flagging:
Currently if adjoint_flagging = True then either flag2refine and/or flag_richardson can be applied, with the tolerance applied to the adjoint inner product rather than directly to q or the error estimate.
This could be extended to arrays as above by letting adjoint_flagging be an array and apply the current approach on each level based on adjoint_flagging[L].
A more general approach would be to allow adjoint_flagging in addition to one or both of flag2refine (based on gradient of q) and flag_richardson (error estimates of q). Ideally we would allow two new possibilities, adjoint_flagging based on magnitude of inner product of q and the adjoint, and adjoint_flagging based on the inner product of the error estimate with the adjoint. This does not seem worth fully implementing at this time because doing both Richardson flagging on q and on the adjoint simultaneously would require more refactoring and probably would not be used.
But I do suggest we introduce:
for more clarity of what is being done in setrun.py, and with the requirement that flag_adjoint_richardson and flag_richardson cannot both be True on the same level.
In particular, allow both flag2refine and flag_adjoint to be True on the same level, and remove the code that computes the adjoint inner product from flag2refine since it does not belong there anyway -- the original idea was that flag2refine.f90 should be a simple routine for gradient flagging that the user could easily modify to implement some other desired flagging based directly on q.
This set of modifications would in particular allow setting, e.g.
in order to use adjoint flagging only on levels L=1 and 2, and flagging based on q in all finer levels. This is a common use case in GeoClaw where we want to use adjoint flagging to track waves over the ocean on levels up to L+1=3, but not on the finer grids near the coastal location of interest, where wave_tolerance flagging is both much cheaper to apply and also properly identifies onshore locations that are flooded (whereas the adjoint of the linearized ocean at rest is identically 0 onshore).