From 9bdf03cb7fcb28c1874f916cd115117c4ef42665 Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Wed, 20 Aug 2025 16:10:22 -0700 Subject: [PATCH] refactor dtopotools.create_dtopography to avoid memory issues As originally written, the deformation dtopo.dZ was computed for each subfault and stored, and then summed up over all subfaults. In a CSZ problem with 8000 subfaults and fine resolution dtopo grids, this ran out of memory. The refactored version computes and uses each subfault dtopo in a way that it can be discarded before computing the next subfault deformation. Also added a `slip_tol` parameter so that the deformation is not computed for any subfaults with smaller slip than this tolerance (which could be a large fraction of the subfaults in some situations, e.g. a small localized earthquake on a large subduction zone fault surface). Also the `verbose` parameter can now be an integer, in which case it prints the subfault number `k` during processing only if `numpy.mod(k,verbose) == 0`. --- src/python/geoclaw/dtopotools.py | 99 ++++++++++++++++++++++---------- 1 file changed, 70 insertions(+), 29 deletions(-) diff --git a/src/python/geoclaw/dtopotools.py b/src/python/geoclaw/dtopotools.py index 141bdfde9..368963bf2 100644 --- a/src/python/geoclaw/dtopotools.py +++ b/src/python/geoclaw/dtopotools.py @@ -760,7 +760,8 @@ def Mw(self): r"""Calculate the moment magnitude for a fault composed of subfaults.""" return Mw(self.Mo()) - def create_dtopography(self, x, y, times=[0., 1.], verbose=False): + def create_dtopography(self, x, y, times=[0., 1.], slip_tol=0.001, + verbose=False): r"""Compute change in topography and construct a dtopography object. Use subfaults' `okada` routine and add all @@ -768,9 +769,17 @@ def create_dtopography(self, x, y, times=[0., 1.], verbose=False): Raises a ValueError exception if the *rupture_type* is an unknown type. + Subfaults are ignored if abs(subfault.slip) < slip_tol + + If verbose == True then the subfault number is printed as each is + processed to show progress. If verbose is an integer, then it is + printed for the k'th subfault only if mod(k,verbose) == 0. + returns a :class`DTopography` object. """ + nignore = 0 # to keep track of how many have abs(slip) < slip_tol + dtopo = DTopography() dtopo.x = x dtopo.y = y @@ -783,21 +792,31 @@ def create_dtopography(self, x, y, times=[0., 1.], verbose=False): print("Making Okada dz for each of %s subfaults" \ % len(self.subfaults)) - for k,subfault in enumerate(self.subfaults): - if verbose: - sys.stdout.write("%s.." % k) - sys.stdout.flush() - subfault.okada(x,y) # sets subfault.dtopo with times=[0] - # and subfault.dtopo.dZ.shape[0] == 1 - if verbose: - sys.stdout.write("\nDone\n") + dz = numpy.zeros(X.shape) # to accumulate dZ for static rupture if self.rupture_type == 'static': + + for k,subfault in enumerate(self.subfaults): + if verbose: + if (type(verbose) is int) and (numpy.mod(k,verbose) != 0): + pass + else: + sys.stdout.write("%s.." % k) + sys.stdout.flush() + + if abs(subfault.slip) < slip_tol: + nignore += 1 + else: + dtopo_sub = subfault.okada(x,y,set_dtopo=False) + # returns dtopo with times=[0] and dtopo.dZ.shape[0] == 1 + + dz += dtopo_sub.dZ[0,:,:] + + if verbose: + sys.stdout.write("\nDone\n") + if len(times) > 2: raise ValueError("For static deformation, need len(times) <= 2") - dz = numpy.zeros(X.shape) - for subfault in self.subfaults: - dz += subfault.dtopo.dZ[0,:,:] if len(times) == 1: # only final deformation stored: @@ -811,33 +830,47 @@ def create_dtopography(self, x, y, times=[0., 1.], verbose=False): elif self.rupture_type in ['dynamic','kinematic']: - t_prev = -1.e99 - dzt = numpy.zeros(X.shape) - dZ = None - for t in times: - for k,subfault in enumerate(self.subfaults): + # array subfault_frac[k,jt] will store fraction of subfault k + # deformation to add to dZ[jt,:,:] at times[jt]: + subfault_frac = numpy.zeros((len(self.subfaults),len(times))) - rf = rise_fraction([t_prev,t], + # compute these values from rupture and rise times, without yet + # computing deformations: + for k,subfault in enumerate(self.subfaults): + subfault_frac[k,:] = rise_fraction(times, subfault.rupture_time, subfault.rise_time, subfault.rise_time_starting, subfault.rise_shape) - dfrac = rf[1] - rf[0] - if dfrac > 0.: - dzt = dzt + dfrac * subfault.dtopo.dZ[0,:,:] + # now loop over subfaults, applying Okada to get + # deformation dtopo_sub due to single subfault, then loop + # over times to apply this at each time: - dzt = numpy.array(dzt, ndmin=3) # convert to 3d array - if dZ is None: - dZ = dzt.copy() + dZ = numpy.zeros([len(times)] + list(X.shape)) + + for k,subfault in enumerate(self.subfaults): + if abs(subfault.slip) < slip_tol: + nignore += 1 else: - dZ = numpy.append(dZ, dzt, axis=0) - t_prev = t + dtopo_sub = subfault.okada(x,y,set_dtopo=False) + # returns dtopo with times=[0] and dtopo.dZ.shape[0] == 1 + + for jt,t in enumerate(times): + dfrac = subfault_frac[k,jt] # fraction to add at t + + if dfrac > 0.: + dZ[jt,:,:] += dfrac * dtopo_sub.dZ[0,:,:] + dtopo.dZ = dZ else: raise ValueError("Unrecognized rupture_type: %s" % self.rupture_type) + if verbose and nignore > 0: + print('Ignored %i subfaults with abs(slip) < slip_tol = %.3fm' \ + % (nignore,slip_tol)) + # Store for user self.dtopo = dtopo @@ -1682,12 +1715,18 @@ def _llz2utm(self,lon,lat,projection_zone=None): return x,y - def okada(self, x, y): + def okada(self, x, y, set_dtopo=True): r""" Apply Okada to this subfault and return a DTopography object. :Input: - x,y are 1d arrays + - set_dtopo (bool): If True, also set self.dtopo = dtopo, the + DTopography object created. + Defaults to True for backward compatibility, but create_dtopography + has been refactored so it is no longer necessary to save dtopo + for each subfault separately, which may use too much memory if + there are many subfaults with fine resolution dtopo files. :Output: - DTopography object with dZ array of shape (1,len(x),len(y)) with single static displacement and times = [0.]. @@ -1774,7 +1813,8 @@ def okada(self, x, y): dtopo.Y = Y dtopo.dZ = numpy.array(dz, ndmin=3) dtopo.times = [0.] - self.dtopo = dtopo + if set_dtopo: + self.dtopo = dtopo elif self.coordinate_specification == 'triangular': @@ -1866,7 +1906,8 @@ def okada(self, x, y): dtopo.dY = numpy.array(dY, ndmin=3) dtopo.dZ = numpy.array(dZ, ndmin=3) dtopo.times = [0.] - self.dtopo = dtopo + if set_dtopo: + self.dtopo = dtopo return dtopo