Allow different AMR flagging on each level and flagregion changes

[rjleveque]

The changes I proposed in clawpack/amrclaw#262 would naturally be carried over to GeoClaw, allowing different flagging methods and/or tolerances at each AMR level.

In addition, I propose the following:

Currently in setrun.py one can specify flagregions (either as traditional rectangular "regions" or as Ruled Rectangles). But regions are also automatically created to go along with every topo or dtopo file, since the way a topofile is specified, for example includes parameters minlevel, maxlevel, t1, t2. These are used to define a rectangular region whose spatial extent is exactly that of the topofile. For each dtopo file a region is created with spatial extent and time interval determined by the dtopo data and only minlevel, maxlevel are specified along with the path to the dtopo file. This was originally intended to make it convenient in cases where the user wants to define a certain level of refinement around the earthquake source region, for example, over the duration of the earthquake.

However, this capability seems to be little used in practice and much more frequently is the cause of confusion for new users. And even experienced users often get messed up with refinement not happening where desired because of the parameters set for a topo or dtopo file were over-riding the regions that were explicitly set.

I propose we eliminate this option, so that specifying a topo or dtopo file only specifies this input data and does not implicitly generate a region. If the user wants a region with exactly the same spatial extent as the topo file, then they should explicitly set a flagregion with this spatial extent.

Another possibility for providing this convenience: the new way of specifying flagregions allows setting either:

flagregion.spatial_region_type = 1  # Rectangle
flagregion.spatial_region = [x1,x2,y1,y2]

for a rectangle, or e.g.

flagregion.spatial_region_type = 2  # Ruled Rectangle
flagregion.spatial_region_file = \
    os.path.abspath('RuledRectangle_Trapezoid.data')

for a Ruled rectangle.

We could also allow something like:

  flagregion.spatial_region_type = 3  # topofile extent
  flagregion.spatial_region_file = os.path.abspath('etopo1.asc')

in which case it would read the spatial extent [x1,x2,y1,y2] from the topofile header and set up a rectangular region. This would allow the convenience originally intended but in a much clearer manner than the current system.

[rjleveque]

I added another suggestion to clawpack/amrclaw#262 that would be useful in particular for GeoClaw, namely introducing amrdata.amr_levels_max_outside_flagregions.

That reminded me of another GeoClaw parameter that often causes confusion and/or unexpected behavior, the parameters deep_depth and max_level_deep. These are often confusing to people and don't work so well anyway and I propose eliminating them.

The original idea was to try to make sure that fine grids only needed near the coast should not be introduced in deeper water when it can be avoided, because the wave speed sqrt(g*h) is larger in deeper water and so small fine grids extending farther than needed into deep water potentially require much smaller time steps than if the fine grid is only very close to shore, due to the CFL restriction.

So the idea is that if the water depth is greater than deep_depth then a cell would not be flagged for refinement beyond level max_level_deep, even if it's in a region where it would be flagged otherwise. This is confusing because it's not clear how this works in conjunction with maxlevel of a flagregion, for example. Also it doesn't always work all that well -- even if a cell isn't flagged it may get refined anyway if it is part of the buffer or if it is in a rectangular patch after clustering the flagged and buffered cells. So there might still be cells refined with water depth much greater than deep_depth.

It is probably better to try to keep refinement away from deep water by specifying flagregions that hug the coast, which is now much easier to do with Ruled Rectangles, e.g. this region that follows the continental shelf, or by specifying a flagregion that covers a set of fgmax points.

[mandli]

Agreed.

[donnaaboise]

As somebody who recently stumbled onto the dtopo/region issue, I can attest to the confusion it initially caused. One source of confusion for me was that if a region was created for a dtopo or topo file, it never showed up in the regions.data file, and so it was unclear how the logic was applied to these hidden regions.

I have since taken @rjleveque 's advice, and created a region that explicitly does what adding a dtopo/topo file was doing implicitly. This still seems like a useful in some cases, and so I am wondering if it makes sense to have an option to explicitly create a region from a topo/dtopo file. Something like :

regions.append(create_region_from_topo(topofile, t1,t2,minlevel,maxlevel))

or

regions.append(create_region_from_dtopo(dtopofile, minlevel,maxlevel))

where in the dtopo case, the time interval is taken from the dtopo info.

These regions would then show up as usual regions in the regions. data file, and so it would then be clear to the user what was going on. This would also save the user the need to parse the topo/dtopo files for extent and time intervals.

[mandli]

My thought was mostly along the lines of what @donnaaboise is suggesting, namely processing through the topo/dtopo such that

1. for all old versions of the topo/dtopo specification that regions are added to the regions array and
2. for a new version of the topo/dtopo specification that these regions are not created.
   This would serve the purpose of putting all region specifications in one place (regions.data) and allow the specification of topo/dtopo to be decoupled from the older specification format but still support the old one. I would suggest that we handle this in the appropriate ClawData objects so that this is transparent to a user. One nice consequence of containing all regions inside of regions.data is that we could write a plotting function that could plot all the regions and show (or color) the level bounds over a topo file provided.

[rjleveque]

@donnaaboise - That is exactly what I was suggesting in proposing we allow the user to specify a flagregion with something like:

flagregion.spatial_region_type = 3  # topofile extent
flagregion.spatial_region_file = os.path.abspath('etopo1.asc')

which would set up an ordinary flagregion but taking the extent from the specified topo file. In the case of a dtopo file it could also take t1, t2 from the dtopo file by default or we could allow specifying a different t1, t2 (as would often be desired since even if the dtopo is instantaneous and finished after 1 second, say, we may want to continue resolving this region for a longer time period).

Note that in v5.7.0 there is a new flagregions.data file for all the flagregions specified. For backward compatibility there is also still a regions.data file, for regions specified the pre-5.7.0 way, but I would like to also deprecate that when we make these other changes so that all flag region info will be in flagregions.data.

[mandli]

This was to accommodate the adjoint flagging as well as the ruled rectangles?

[rjleveque]

@mandli: You mean the new flagregions.data? It was to accommodate the different set of input data needed for ruled rectangles and other options that can now be set related to flagregions.

Also I think flagregions is a more descriptive name than just regions since people talk about regions of the computational domain in other contexts and it's never been so clear to some users what we mean by regions.

The old regions.data still has the old format and only allows rectangles, but rectangles can also be specified in the flagregions.data format so I'd like to get rid of the old ones and have all flagregions specified in one file (including any that might be generated from topo or dtopo extents).

Regarding plotting all regions, this is already supported in GeoClaw as kml files, both for rectangles and ruled rectangles. Just add this to the bottom of setrun.py after rundata is computed in __main__:

 from clawpack.geoclaw import kmltools
 kmltools.make_input_data_kmls(rundata)

which in particular calls kmltools.regions2kml which is currently set up to read in both regions.data and flagregions.data.

One reason I added a flagregion.name attribute to a clawpack.amrclaw.data.FlagRegion object is so that it shows up with an appropriate name in the Google Earth index rather than just as region0 etc.

Note that kmltools.make_input_data_kmls also makes a kml for each topo file and dtopo file, and ones showing all the gauge locations, fgmax grids, etc. This is very useful (and should probably be better documented and advertised).

[mandli]

@rjleveque Thanks for the clarification. I was thinking we would also have something apart from KML to visualize these but I had forgotten about that capability.

[rjleveque]

Yes, it would be nice to have other plotting tools too, particularly for amrclaw simulations that have nothing to do with geosciences!