Center points updates

src/python/geoclaw/center_gauges.py

@@ -0,0 +1,105 @@
+
+from clawpack.geoclaw import center_points
+
+def center_all_gaugedata(rundata, dx, dy, modify_rundata=True,
+                         print_appends=True, verbose=False):
+    """
+    Function that can be called from setrun after setting rundata.gaugedata,
+    to center all of the gauges using the domain edges and the specified dx,dy,
+    the resolution at which you want them centered (in cells that will be
+    an integer number of dx,dy away from the domain edges).
+
+    Note that specifying any dx,dy less than the finest resolution in the
+    simulation will at least insure that the gauges do not lie on grid
+    cell edges on any AMR level (which can create noisy output if rounding
+    errors cause values to be recorded at some times from one cell and at
+    other times from its neighbor).
+    """
+    import copy
+
+    if len(rundata.gaugedata.gauges) == 0:
+        print('No gauges found to center')
+        return rundata
+
+    # assume the x,y values specified in gaugedata.gauges are the desired
+    # locations that we will adjust a bit if necessary:
+    # recall gaugedata.gauges is a list of lists for each gauge of the form:
+    #    [gaugeno, x, y, t1, t2]
+
+    x_desired = [gauge[1] for gauge in rundata.gaugedata.gauges]       
+    y_desired = [gauge[2] for gauge in rundata.gaugedata.gauges]       
+
+    # edges of cells should be integer multiples of dx,dy from domain edge:
+    x_edge = rundata.clawdata.lower[0]
+    y_edge = rundata.clawdata.lower[1]
+
+    x_centered, y_centered = center_points.adjust_xy(x_desired, y_desired,
+                                       x_edge, y_edge, dx, dy,
+                                       verbose=verbose)
+
+    # reset the x,y values for each gauge:
+    gauges = copy.deepcopy(rundata.gaugedata.gauges)
+
+    for k,gauge in enumerate(gauges):
+        gauge[1] = float(x_centered[k])
+        gauge[2] = float(y_centered[k])
+
+    if print_appends:
+        print('Modifications for setrun:')
+        for k in range(len(gauges)):
+            print('gauges.append([%i, %.8f, %.8f, %g, %g])' \
+                    % tuple(gauges[k]))
+
+    if modify_rundata:
+        rundata.gaugedata.gauges = gauges
+        print('*** Shifted rundata.gaugedata.gauges if necessary to center in cells')
+    return rundata
+
+if __name__=='__main__':
+
+    # Sample code...
+
+    # Center gauges based on rundata in setrun.py, assuming they 
+    # should be centered at the finest AMR level specified by
+    # `amr_level_max` and the refinement ratios.
+
+    # Optional argument when executing for path to setrun file to use
+
+    import sys
+    from clawpack.clawutil.util import fullpath_import
+
+    if len(sys.argv) > 1:
+        setrun_file = sys.argv[1]
+    else:
+        setrun_file = 'setrun.py'
+
+    print(f'Will center gauges at finest resolution based on setrun in {setrun_file}')
+    setrun = fullpath_import(setrun_file)
+
+
+    rundata = setrun.setrun()
+    amr_levels_max = rundata.amrdata.amr_levels_max
+    rrx = rundata.amrdata.refinement_ratios_x
+    rry = rundata.amrdata.refinement_ratios_y
+    lower = rundata.clawdata.lower
+    upper = rundata.clawdata.upper
+    dx_level1 = (upper[0] - lower[0]) / rundata.clawdata.num_cells[0]
+    dy_level1 = (upper[1] - lower[1]) / rundata.clawdata.num_cells[1]
+
+    dx_finest = dx_level1
+    dy_finest = dy_level1
+    for k in range(amr_levels_max-1):
+        dx_finest = dx_finest / rrx[k]
+        dy_finest = dy_finest / rry[k]
+
+    print(f'Centering based on domain with:')
+    print(f'          xlower = {lower[0]:.6f}, ylower = {lower[1]:.6f}')
+    print(f'          dx_level1 = {dx_level1:.6e}, dy_level1 = {dy_level1:.6e}')
+    print(f'          amr_levels_max = {amr_levels_max} with finest resolution')
+    print(f'          dx_finest = {dx_finest:.6e}, dy_finest = {dy_finest:.6e}')
+
+    # center gauges and print out the revised lines to insert in setrun.py
+    center_all_gaugedata(rundata, dx_finest, dy_finest, modify_rundata=False,
+                         print_appends=True, verbose=False)
+
+

src/python/geoclaw/center_points.py

@@ -2,18 +2,23 @@
 The function adjust_xy will center a point in a finite volume grid cell of
 a given resolution on a given domain.
 
-This can be used in particular to take adjust an approximate point
+This can be used in particular to adjust an approximate point
 (x_desired, y_desired) where you want a computational gauge, to obtain
 a point that lies exactly in the center of a grid cell at a given resolution.
 This eliminates the need to interpolate between cell values in GeoClaw output,
 which has issues for cells near the shoreline as described at
 https://www.clawpack.org/nearshore_interp.html
 
+Can also be used to center fgmax or fgout points. In this case, center one
+point from the uniform grid. Then, assuming the grid has the same resolution as
+the specified dx,dy, all of the points on the grid should be aligned with
+centers of the finite volume grid at that resolution.
+
 Functions:
 
  - adjust: utility function to adjust in one dimension (x or y)
  - adjust_xy: the function to call to center a point in both x and y.
- - test: a simple example
+ - test_adjust_xy: a simple example
 
 """
 
@@ -22,13 +27,12 @@
 def adjust(z_desired, z_edge, dz, verbose=False):
     """
     Given a desired location z (either x or y) for a gauge or other point
-    of interest, adjust the location is it is offset (integer + 1/2)*dz
-    from z_edge, an arbitrary edge location (e.g. an edge of the
-    computational domain or any point offset integer*dz from the domain edge).
+    of interest, adjust the location it is offset (integer + 1/2)*dz
+    from z_edge, an arbitrary edge location (e.g. an edge of the computational
+    domain, or any point offset by integer*dz from the domain edge).
     This will put the new location in the center of a finite volume cell
     provided this point is in a patch with resolution dz.
     """
-    # z represents either x or y
     i = np.round((z_desired-z_edge - 0.5*dz)/dz)
     z_centered = z_edge + (i+0.5)*dz
     if verbose:
@@ -37,6 +41,9 @@ def adjust(z_desired, z_edge, dz, verbose=False):
         print("adjusted from %15.9f" % z_desired)
         print("           to %15.9f" % z_centered)
         print("   shifted by %15.9f = %.3f*dz" % (zshift,zfrac))
+        print("   (z_centered - z_edge)/dz should be int+0.5: ", \
+                (z_centered - z_edge)/dz)
+        print(f'   z_desired = {z_desired}; z_centered = {z_centered:.12f}; z_edge = {z_edge:.12f};  dz = {dz:.12f}')
     return z_centered
 
 def adjust_xy(x_desired, y_desired, x_edge, y_edge, dx, dy, verbose=False):
@@ -63,9 +70,14 @@ def adjust_xy(x_desired, y_desired, x_edge, y_edge, dx, dy, verbose=False):
 
     - xc,yc: centered points that lie within dx/2, dy/2 of the desired
       location(s) and with (xc - x_edge)/dx and (yc - y_edge)/dy
-      equal to an integer + 0.5,
+      equal to an integer + 0.5.
+
+      Returns numpy arrays if inputs are lists or arrays, 
+      floats if inputs are scalar
     """
 
+    return_scalar = np.isscalar(x_desired)
+
     # convert single values or lists/tuples to numpy arrays
     x_desired = np.array(x_desired, ndmin=1)
     y_desired = np.array(y_desired, ndmin=1)
@@ -81,32 +93,37 @@ def adjust_xy(x_desired, y_desired, x_edge, y_edge, dx, dy, verbose=False):
         x_centered.append(x)
         y_centered.append(y)
 
-    if len(x_centered) == 1:
+    if return_scalar:
         return float(x_centered[0]), float(y_centered[0])
     else:
         return np.array(x_centered), np.array(y_centered)
 
-def test():
+def test_adjust_xy():
 
     x_desired = [-122.01, -122.02]
     y_desired = [47.001, 47.002]
 
     # grid resolution on which to center point:
-    dx = 1/(3*3600.)
+    dx = dy = 1/(3*3600.)
 
     # lower left edge of computational domain:
     x_edge = -123.
     y_edge = 45.
 
+    print('x_edge = ',x_edge)
+    print('y_edge = ',y_edge)
+    print('dx = ',dx)
+    print('dy = ',dy)
     print('Desired x = ',x_desired)
     print('Desired y = ',y_desired)
 
-    xc,yc = adjust_xy(x_desired,y_desired,x_edge,y_edge,dx,dx,verbose=True)
+    xc,yc = adjust_xy(x_desired,y_desired,x_edge,y_edge,dx,dy,verbose=True)
 
     print('Centered x = ',xc)
     print('Offsets in x in units of 1/3 arcsec: ', (xc-x_edge)*3*3600)
     print('Centered y = ',yc)
     print('Offsets in y in units of 1/3 arcsec: ', (yc-y_edge)*3*3600)
 
+
 if __name__ == '__main__':
-    test()
+    test_adjust_xy()