From 8ccb6704640f47b210026874f9c711b4312a475c Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Sun, 21 Dec 2025 18:13:13 -0800 Subject: [PATCH 1/2] fix center_points.adjust_xy to work for single gauge, improve comments and printing --- src/python/geoclaw/center_points.py | 41 ++++++++++++++++++++--------- 1 file changed, 29 insertions(+), 12 deletions(-) diff --git a/src/python/geoclaw/center_points.py b/src/python/geoclaw/center_points.py index da4e7f1db..1ac895635 100644 --- a/src/python/geoclaw/center_points.py +++ b/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() From 5d19757707865c2c7dc0b672b8f8d73d1a1d3c3b Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Sat, 27 Dec 2025 11:17:14 -0800 Subject: [PATCH 2/2] utility function center_gauges.py to center gauges that have already been specified in setrun.py --- src/python/geoclaw/center_gauges.py | 105 ++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 src/python/geoclaw/center_gauges.py diff --git a/src/python/geoclaw/center_gauges.py b/src/python/geoclaw/center_gauges.py new file mode 100644 index 000000000..5dbf17b45 --- /dev/null +++ b/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) + +