Skip to content

solver.py

LevelSetSolver dataclass 

LevelSetSolver(
    Δx: float,
    Δt: float = 0,
    n_total: int | None = None,
    i_step: int | None = None,
    t_total: float = 0,
    Δt_reinitialize: float = 0.1,
    κ_boost: float = 1,
    n_slabfailure: int | None = None,
    n_pad_pixels: int | None = None,
    band_width: float | None = None,
    fm_order: int | None = None,
    n_pixels_xz0: tuple[int, int] | None = None,
    n_pixels_xz: tuple[int, int] | None = None,
    do_extend_line: bool = True,
    do_φ_everywhere: bool | None = None,
    domain: Domain | None = None,
    raw_domain: Domain | None = None,
    surface: (
        StraightLineSurface
        | ErrorFunctionSurface
        | SemiCircularArcSurface
        | QuarterCircularArcSurface
        | None
    ) = None,
    substrate: (
        OneLayerSubstrate | TwoLayerSubstrate | MultiLayerSubstrate | None
    ) = None,
    model: Isotropic | Step | ExponentialActivation | None = None,
    gma: GeometricMechanicsAnalysis | None = None,
)


              flowchart TD
              erosionfront.levelset.solver.LevelSetSolver[LevelSetSolver]
              erosionfront.levelset.base.LevelSetSolverBase[LevelSetSolverBase]
              erosionfront.levelset.speed.LevelSetSpeedMethods[LevelSetSpeedMethods]
              erosionfront.levelset.grid.LevelSetGridMethods[LevelSetGridMethods]
              erosionfront.levelset.elementary.LevelSetElementaryMethods[LevelSetElementaryMethods]
              erosionfront.levelset.data.LevelSetData[LevelSetData]
              erosionfront.levelset.gradient.GradientMixin[GradientMixin]
              erosionfront.levelset.shock.ShocksMixin[ShocksMixin]
              erosionfront.levelset.export.ExportMixin[ExportMixin]

                              erosionfront.levelset.base.LevelSetSolverBase --> erosionfront.levelset.solver.LevelSetSolver
                                erosionfront.levelset.speed.LevelSetSpeedMethods --> erosionfront.levelset.base.LevelSetSolverBase
                                erosionfront.levelset.grid.LevelSetGridMethods --> erosionfront.levelset.speed.LevelSetSpeedMethods
                                erosionfront.levelset.elementary.LevelSetElementaryMethods --> erosionfront.levelset.grid.LevelSetGridMethods
                                erosionfront.levelset.data.LevelSetData --> erosionfront.levelset.elementary.LevelSetElementaryMethods
                




                erosionfront.levelset.gradient.GradientMixin --> erosionfront.levelset.solver.LevelSetSolver
                
                erosionfront.levelset.shock.ShocksMixin --> erosionfront.levelset.solver.LevelSetSolver
                
                erosionfront.levelset.export.ExportMixin --> erosionfront.levelset.solver.LevelSetSolver
                


              click erosionfront.levelset.solver.LevelSetSolver href "" "erosionfront.levelset.solver.LevelSetSolver"
              click erosionfront.levelset.base.LevelSetSolverBase href "" "erosionfront.levelset.base.LevelSetSolverBase"
              click erosionfront.levelset.speed.LevelSetSpeedMethods href "" "erosionfront.levelset.speed.LevelSetSpeedMethods"
              click erosionfront.levelset.grid.LevelSetGridMethods href "" "erosionfront.levelset.grid.LevelSetGridMethods"
              click erosionfront.levelset.elementary.LevelSetElementaryMethods href "" "erosionfront.levelset.elementary.LevelSetElementaryMethods"
              click erosionfront.levelset.data.LevelSetData href "" "erosionfront.levelset.data.LevelSetData"
              click erosionfront.levelset.gradient.GradientMixin href "" "erosionfront.levelset.gradient.GradientMixin"
              click erosionfront.levelset.shock.ShocksMixin href "" "erosionfront.levelset.shock.ShocksMixin"
              click erosionfront.levelset.export.ExportMixin href "" "erosionfront.levelset.export.ExportMixin"

Solver implementing level-set solution of geomorphic Hamilton-Jacobi eqn.

Subclasses LevelSetSolverBase and uses GradientMixin, ShocksMixin, ExportMixin.

Parameters

See LevelSetSolverBase and its antecedents.

Attributes

φ: MaskedArray Signed-distance function φ grid, masked to narrow band across surface. Δφ_x: MaskedArray Finite-difference approx of ∂φ/∂x. Δφ_z: MaskedArray Finite-difference approx of ∂φ/∂z. Δ2φ_x2: MaskedArray Finite-difference approx of ∂2φ/∂x2. Δ2φ_z2: MaskedArray Finite-difference approx of ∂2φ/∂z2. β: MaskedArray Surface tilt angle grid. ξ_fn_β: Callable Surface-normal erosion speed model function ξ(β). ξ_model: MaskedArray Grid of model erosion speed ξ_model for narrow-band across the surface. ξ0: MaskedArray Grid of reference erosion speed ξ_0 for narrow-band across the surface. ξ: MaskedArray Grid of actual erosion speed ξ for narrow-band across the surface, computed as ξ = ξ_0 * ξ_model H_ls: MaskedArray Level-set numerical Hamiltonian (not to be confused with the geomorphic Hamiltonian) used in the HJE solution of erosion-front motion. Again computed for a narrow-band across the surface. dHlsdφx: MaskedArray Finite-difference approximation of the derivative of the level-set Hamiltonian wrt ∂φ/∂x. This term is used to compute a numerical viscosity such that we obtain a 'viscosity solution' of the HJE. Again computed for a narrow-band across the surface. dHlsdφz: MaskedArray Finite-difference approximation of the derivative of the level-set Hamiltonian wrt ∂φ/∂z. This term is used to compute a numerical viscosity such that we obtain a 'viscosity solution' of the HJE. Again computed for a narrow-band across the surface. κx: float Numerical viscosity in the x direction. Again computed for a narrow-band across the surface. κz: float Numerical viscosity in the x direction. Again computed for a narrow-band across the surface. dφdt: MaskedArray Rate of change of signed-distance function φ wrt time. Obtained from the (level-set) Hamilton-Jacobi equation. Again computed for a narrow-band across the surface. φ_next: MaskedArray | None Signed-distance function φ grid for the next time step. Again computed for a narrow-band across the surface. φ_everywhere: NDArray Signed-distance function φ computed across the whole grid. This field is not extrapolated by rather simply extended from the narrow band by filling the remaining grid cells with max/min values. φ_aux: NDArray | None An auxiliary grid used to compute which pixels are subject to 'slab failure'.

Used by:

Methods:

  • __post_init__ –

    Complete instantiation by computing domain, pixel size, offset.

  • above_or_below_ji –

    Test whether pixel i,j lies above or below surface.

  • above_or_below_xz –

    Test whether point x,z lies above or below surface.

  • compute_contour –

    Generate a continuous (for now) contour of grid function f.

  • compute_domain –

    (Re)build domain after resizing.

  • compute_line_pixel_offsets –

    Measure signed offsets between rasterized line pixels and line curve.

  • compute_pixel_dimensions –

    Calculate numbers of pixels spanning the grid.

  • compute_pixel_origin –

    Calculate pixel offset of domain x,z position (0,0).

  • evolve –

    Iteratively propagate the erosion front surface.

  • export_data –

    Tool to export time T slice data to HDF5 file.

  • export_metadata –

    Tool to simulation metadata to JSON file.

  • find_shock –

    Locate main break-in-slope by smoothing & changepoint detection.

  • find_shocks –

    Find breaks-in-slope in the set of T-sliced surfaces.

  • first_update –

    Set up time-stepping solution with requisite initializations.

  • fm_distance –

    Use fast marching to compute signed distance from a ref grid (zeroset).

  • fm_distance_everywhere –

    FM compute a narrowband signed distance and then extend everywhere.

  • fwd_back_differences –

    Compute forward, backward, mean & double finite differences.

  • get_extent –

    Get domain (x,z) extent.

  • get_front_points –

    Fetch points along erosion front surface aka zeroset of f.

  • get_line_points –

    Convert Shapely geometry Line into (x,z) points.

  • get_points_xz –

    Get surface (x,z) coordinates.

  • get_substrate_ξ0 –

    Generate relative η grid (effective speed ξ0).

  • get_x_limits –

    Get domain x min, max.

  • get_z_limits –

    Get domain z min, max.

  • get_φ_next –

    Return next signed-distance grid φ_next or just φ if not available.

  • ij_to_xz –

    Convert pixel indices [i,j] into pixel center coords (x,z).

  • initialize –

    Initialize level-set solver.

  • ji_to_xz –

    Convert reversed pixel indices [j,i] into pixel center coords (x,z).

  • make_line –

    Convert surface (x,z) coordinates into a Shapely Line geometry.

  • make_polygon –

    Create a 'bounding' polygon from surface x,z points.

  • measure_normal_distance –

    Compute distance between line and point using Shapely.

  • measure_φ_ji –

    Compute signed distance between line and i,j pixel using Shapely.

  • measure_φ_xz –

    Compute signed distance between line and xz point using Shapely.

  • pad_domain –

    Expand the domain vertically by pixel width amounts.

  • pad_grid –

    Pad grid with specified number of pixels below and above.

  • rasterize_T_slices –

    Rasterize time T slices onto grids.

  • rasterize_line –

    Rasterize surface line as x,z points into grid pixels using Rasterio.

  • record_T_slice –

    Record selected sim data for a single time step.

  • redo_φ_everywhere –

    Reassign signed distance values to pixels beyond the narrow band.

  • update –

    Update grids for one time step of duration Δt.

  • xz_to_ij –

    Convert pixel center coords (x,z) into pixel indices [i,j].

__post_init__ 

__post_init__()

Complete instantiation by computing domain, pixel size, offset.

Attributes

raw_domain: Domain Model domain instance. n_pixels_xz: tuple[int, int]: Grid pixel size. n_pixels_xz0: tuple[int, int]: Origin location in pixels.

Source code in erosionfront/levelset/elementary.py 
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
def __post_init__(self):
    """
    Complete instantiation by computing domain, pixel size, offset.

    Attributes
    ----------
    raw_domain: Domain
        Model domain instance.
    n_pixels_xz: tuple[int, int]:
        Grid pixel size.
    n_pixels_xz0: tuple[int, int]:
        Origin location in pixels.
    """
    super().__post_init__()
    # Compute bounds, extent, pixel dimensions of grid.
    # These methods MUST be provided by any subclass.
    if self.surface is None:
        raise ValueError("Surface not specified")
    self.raw_domain \
        = self.compute_domain(
            self.surface.x, self.surface.z, self.Δx,
        )
    self.n_pixels_xz \
        = self.compute_pixel_dimensions(self.raw_domain, self.Δx,)
    self.n_pixels_xz0 \
        = self.compute_pixel_origin(self.raw_domain, self.Δx,)

above_or_below_ji 

above_or_below_ji(ji: tuple[int, int] | NDArray, in_rock_polygon: Polygon) -> int

Test whether pixel i,j lies above or below surface.

Uses "in_rock_polygon" to test if pixel is within the bedrock or in the air. In other words, it tests which side of the surface zeroset the pixel lies.

The pixel center is used in this computation.

Parameters

ji: tuple[float,float] | NDArray Test pixel with reversed indexes. in_rock_polygon: Polygon Shapely polygon defining "below ground".

Atrributes

Uses ji_to_xz method.

Returns

-1 if "in the air" and +1 if "in the rock".

Source code in erosionfront/levelset/grid.py 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
def above_or_below_ji(
        self,
        ji: tuple[int,int] | NDArray,
        in_rock_polygon: Polygon,
    ) -> int:
    """
    Test whether pixel i,j lies above or below surface.

    Uses "in_rock_polygon" to test if pixel is within the bedrock or in 
    the air. In other words, it tests which side of the surface zeroset
    the pixel lies.

    The pixel center is used in this computation.

    Parameters
    ----------
    ji: tuple[float,float] | NDArray
        Test pixel with reversed indexes.
    in_rock_polygon: Polygon
        Shapely polygon defining "below ground".

    Atrributes
    ----------
    Uses ji_to_xz method.

    Returns
    -------
    -1 if "in the air" and +1 if "in the rock".
    """
    # integer pixel indices convert to pixel *center* x,z coords
    return (
        -1 if shapely.contains_xy(in_rock_polygon, *self.ji_to_xz(ji)) 
        else +1
    )

above_or_below_xz staticmethod 

above_or_below_xz(xz: tuple[float, float] | NDArray, in_rock_polygon: Polygon) -> int

Test whether point x,z lies above or below surface.

Uses "in_rock_polygon" to test if point is within the bedrock or in the air. In other words, it tests which side of the surface zeroset the point lies.

Parameters

xz: tuple[float,float] | NDArray Test point. in_rock_polygon: Polygon Shapely polygon defining "below ground".

Returns

-1 if "in the air" and +1 if "in the rock".

Source code in erosionfront/levelset/grid.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
@staticmethod
def above_or_below_xz(
        xz: tuple[float,float] | NDArray,
        in_rock_polygon: Polygon,
    ) -> int:
    """
    Test whether point x,z lies above or below surface.

    Uses "in_rock_polygon" to test if point is within the bedrock or in 
    the air. In other words, it tests which side of the surface zeroset
    the point lies.

    Parameters
    ----------
    xz: tuple[float,float] | NDArray
        Test point.
    in_rock_polygon: Polygon
        Shapely polygon defining "below ground".

    Returns
    -------
    -1 if "in the air" and +1 if "in the rock".
    """
    return (
        -1 if shapely.contains_xy(in_rock_polygon, *xz,) else +1
    )

compute_contour 

compute_contour(f: NDArray | MaskedArray, level: float = 0) -> NDArray

Generate a continuous (for now) contour of grid function f.

BUG: the contouring function actually returns a set of disjoint contour segments, often with more than one segment. Here it's incorrectly assumed these segments are contiguous and can simply be concatenated. This is not true! But to do better downstream use of this method needs to handle the segment set properly, which is going to be inconvenient.

Parameters

f: NDArray | MaskedArray Smooth grid (e.g., signed distance) on which to compute contour. level: float=0 Chosen contour level.

Attributes

Uses ij_to_xz method.

Returns

NDArray Single sequence of x,z contour coordinates as (2,N) array.

Source code in erosionfront/levelset/grid.py 
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
def compute_contour(
        self, 
        f: NDArray | MaskedArray, 
        level: float=0,
    ) -> NDArray:
    """
    Generate a continuous (for now) contour of grid function f.

    BUG: the contouring function actually returns a set of
    disjoint contour segments, often with more than one segment.
    Here it's incorrectly assumed these segments are contiguous 
    and can simply be concatenated. This is not true! But to do
    better downstream use of this method needs to handle the segment set
    properly, which is going to be inconvenient.

    Parameters
    ----------
    f: NDArray | MaskedArray
        Smooth grid (e.g., signed distance) on which to compute contour.
    level: float=0
        Chosen contour level.

    Attributes
    ----------
    Uses ij_to_xz method.

    Returns
    -------
    NDArray
        Single sequence of x,z contour coordinates as (2,N) array.
    """
    contours: list = skimage.measure.find_contours(
        f, 
        level=level, 
        mask=~f.mask if type(f) is MaskedArray else None,
    )
    # print(len(contours))
    #HACK: drop second contour (presumably generated by front sliding
    #      off grid) if it's too short & thus likely erroneous
    # if len(contours)>1:
    #     # print(contours[0].shape, contours[1].shape)
    #     if contours[1].shape[0]<10:
    #         contours = contours[0]
    return self.ij_to_xz(np.fliplr(np.vstack(contours))).T

compute_domain staticmethod 

compute_domain(x: tuple | NDArray, z: tuple | NDArray, resolution: float) -> Domain

(Re)build domain after resizing.

Parameters

x: tuple | NDArray Vectors of x coordinates for which a domain geometry needs to be measured. z: tuple | NDArray Vectors of z coordinates for which a domain geometry needs to be measured. resolution: float Grid pixel size.

Returns

Domain Domain geometry class instance.

Source code in erosionfront/levelset/elementary.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
@staticmethod
def compute_domain(
        x: tuple | NDArray,
        z: tuple | NDArray,
        resolution: float,
    ) -> Domain:
    """
    (Re)build domain after resizing.

    Parameters
    ----------
    x: tuple | NDArray
        Vectors of x coordinates for which a domain geometry needs
        to be measured.
    z: tuple | NDArray
        Vectors of z coordinates for which a domain geometry needs
        to be measured.
    resolution: float
        Grid pixel size.

    Returns
    -------
    Domain
        Domain geometry class instance.
    """
    n_digits: int = int(round(-np.log10(resolution)+0.5))+1
    min: Callable = lambda xz: round(float(
        (np.min(np.array(xz)/resolution) - 0.5) 
    ) * resolution, n_digits)
    max: Callable = lambda xz: round(float(
        (np.max(np.array(xz)/resolution) + 0.5) 
    ) * resolution, n_digits)
    return Domain(min(x), max(x), min(z), max(z),)

compute_line_pixel_offsets 

compute_line_pixel_offsets(
    line: LineString, ρ: NDArray, do_φ_everywhere: bool | None = False
) -> NDArray

Measure signed offsets between rasterized line pixels and line curve.

Spatial aliasing during rasterization means each line pixel is likely offset a bit from the line curve, so this function computes these offsets.

Parameters

line: LineString ρ: NDArray

bool

Compute distances for all grid pixels, not just those not masked.

Attributes

Uses measure_φ_ji and above_or_below_ji methods.

Returns

NDArray Grid with signed-distance normal offsets from the line curve at each rasterized line pixel.

Source code in erosionfront/levelset/grid.py 
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
def compute_line_pixel_offsets(
        self, 
        line: LineString,
        ρ: NDArray,
        do_φ_everywhere: bool | None=False,
    ) -> NDArray:
    """
    Measure signed offsets between rasterized line pixels and line curve.

    Spatial aliasing during rasterization means each line pixel is likely
    offset a bit from the line curve, so this function computes these
    offsets.

    Parameters
    ----------
    line: LineString
    ρ: NDArray

    do_φ_everywhere: bool
        Compute distances for all grid pixels, not just those not masked.

    Attributes
    ----------
    Uses measure_φ_ji and above_or_below_ji methods.

    Returns
    -------
    NDArray
        Grid with signed-distance normal offsets from the line curve
        at each rasterized line pixel.
    """
    offsets_grid: NDArray \
        = np.zeros_like(ρ, dtype=np.float64,)
    for ji_ in np.argwhere(ρ!=0):
        # integer pixel indices convert to pixel *center* x,z coords
        offsets_grid[*ji_] = self.measure_φ_ji(line,ji_,)
    d_min: float = np.min(offsets_grid)
    d_max: float = np.max(offsets_grid)
    if do_φ_everywhere:
        for ji_ in np.argwhere(ρ==0):
            offsets_grid[*ji_] = (
                d_min 
                if self.above_or_below_ji(ji_, self.in_rock_polygon,)==-1 
                else d_max
            )
    return offsets_grid

compute_pixel_dimensions staticmethod 

compute_pixel_dimensions(domain: Domain, resolution: float) -> tuple[int, int]

Calculate numbers of pixels spanning the grid.

Parameters

domain: Domain Grid domain geometry. resolution: float Grid pixel size.

Returns

tuple[int, int] Width and height of grid in pixel numbers.

Source code in erosionfront/levelset/elementary.py 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
@staticmethod
def compute_pixel_dimensions(
        domain: Domain,
        resolution: float,
    ) -> tuple[int, int]:
    """
    Calculate numbers of pixels spanning the grid.

    Parameters
    ----------
    domain: Domain
        Grid domain geometry.
    resolution: float
        Grid pixel size.

    Returns
    -------
    tuple[int, int]
        Width and height of grid in pixel numbers.
    """
    n_i_pixels: int = int(
        (domain.x_max-domain.x_min)/resolution + 0.5
    )
    n_j_pixels: int = int(
        (domain.z_max-domain.z_min)/resolution + 0.5
    )
    return (n_i_pixels, n_j_pixels,)

compute_pixel_origin staticmethod 

compute_pixel_origin(domain: Domain, resolution: float) -> tuple[int, int]

Calculate pixel offset of domain x,z position (0,0).

Parameters

domain: Domain Grid domain geometry. resolution: float Grid pixel size.

Returns

tuple[int, int] Pixel indexes of bottom-left origin of grid.

Source code in erosionfront/levelset/elementary.py 
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
@staticmethod
def compute_pixel_origin(
        domain: Domain,
        resolution: float,
    ) -> tuple[int, int]:
    """
    Calculate pixel offset of domain x,z position (0,0).

    Parameters
    ----------
    domain: Domain
        Grid domain geometry.
    resolution: float
        Grid pixel size.

    Returns
    -------
    tuple[int, int]
        Pixel indexes of bottom-left origin of grid.        
    """
    i_origin: int = int(domain.x_min/resolution+0.5)
    j_origin: int = int(domain.z_min/resolution+0.5)
    return (i_origin, j_origin,)

evolve 

evolve(
    Δt: float = 0,
    t_span: float | None = None,
    n_steps: int | None = None,
    Δt_reinitialize: float | None = None,
    κ_boost: float = 1,
    report_pc: int = 5,
    n_record_steps: int = 5,
    do_initial: bool = False,
    do_timing: bool = False,
) -> None

Iteratively propagate the erosion front surface.

Parameters

Δt: float=0 Simulation time step. t_span: float | None = None Duration of simulation. n_steps: int | None = None Number of time steps to take in the simulation. If neither t_span nor n_steps is specified, If neither are, a single step is assumed. Δt_reinitialize: float | None = None Period of reinitialization of φ signed-distance function grid. κ_boost: float = 1 Scale factor for boosting numerical viscosity. report_pc: int = 5 Reporting interval as percentage of runtime. n_record_steps: int = 5 Number of steps between T slicings. do_initial: bool=False Specify if this is the first update, and act accordingly. do_timing: bool=False Choose whether to perform timings of sim stages. Attributes

Most of the sim grids are modified.

Source code in erosionfront/levelset/base.py 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
def evolve(
        self, 
        Δt: float=0, 
        t_span: float | None = None,
        n_steps: int | None = None,
        Δt_reinitialize: float | None = None, 
        κ_boost: float = 1,
        report_pc: int = 5,
        n_record_steps: int = 5,
        do_initial: bool=False,
        do_timing: bool=False,
    ) -> None:
    """
    Iteratively propagate the erosion front surface.

    Parameters
    ----------
    Δt: float=0
        Simulation time step.
    t_span: float | None = None
        Duration of simulation.
    n_steps: int | None = None
        Number of time steps to take in the simulation.
        If neither `t_span` nor `n_steps` is specified, 
        If neither are, a single step is assumed.
    Δt_reinitialize: float | None = None
        Period of reinitialization of φ signed-distance function grid.
    κ_boost: float = 1
        Scale factor for boosting numerical viscosity.
    report_pc: int = 5
        Reporting interval as percentage of runtime.
    n_record_steps: int = 5
        Number of steps between T slicings.
    do_initial: bool=False
        Specify if this is the first update, and act accordingly.
    do_timing: bool=False
        Choose whether to perform timings of sim stages.
    Attributes
    ----------
    Most of the sim grids are modified.
    """
    self.Δt = Δt
    if self.Δt>0 and self.model is not None:
        self.n_slabfailure = int(self.model.Δt_slabfailure/self.Δt)
    self.κ_boost = κ_boost
    if Δt_reinitialize is None:
        self.Δt_reinitialize = Δt
    else:
        self.Δt_reinitialize = Δt_reinitialize
    n_steps_: int
    t_span_: float
    if n_steps is None:
        if t_span is None:
            n_steps_ = 1
            t_span_ = Δt
        else:
            n_steps_ = int(t_span/Δt+0.5)
            t_span_ = t_span
    else:
        n_steps_ = n_steps
        t_span_ = n_steps*Δt
    print(
        rf"Run time T={t_span_} for n_steps={n_steps_} "
        + rf"with Δt={Δt} and Δx={self.Δx}"
    )
    do_report: Callable \
        = lambda step: (
            False if int(n_steps_*report_pc/100)==0
            else (step+1)%int(n_steps_*report_pc/100)==0
        )
    # do_record: Callable \
    #     = lambda step: (
    #         False if int(n_steps_*record_pc/100)==0 
    #         else (step+1)%int(n_steps_*record_pc/100)==0
    #     )
    rnd: int = int(np.round(-np.log10(Δt if Δt>0 else 1))+1)
    def record_report(step) -> None:
        progress_pc_: float
        if do_report(step):
            progress_pc_ = int(100*(step+1)/(n_steps_))
            print(f"   {progress_pc_}%:  step={self.i_step}"
                + f"  t={np.round(self.t_total, rnd)}")
        if ((step+1) % n_record_steps)==0:
            self.record_T_slice(
                self.i_step if self.i_step is not None else 1, 
                self.t_total, 
            )

    #######################################################################
    if self.t_total is None or do_initial:
        # Initial do-nothing step...
        self.record_T_slice(0, 0,)
        self.i_step = 0
        self.update(0, do_timing=do_timing,)
    else:
        # ...or a full run
        for step_ in range(n_steps_):
            self.i_step = self.i_step+1 if self.i_step is not None else 1
            self.update(Δt, Δt_reinitialize,)
            record_report(step_)
        self.update(Δt, do_finalize=True, do_timing=do_timing,)
    self.n_total = self.i_step

export_data 

export_data(dir: str, filename: str = 'T_slices') -> None

Tool to export time T slice data to HDF5 file.

Parameters

dir: str Path to directory to write HDF5 file into. filename: str="T_slices" Name of HDF file.

Source code in erosionfront/levelset/export.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def export_data(
        self,
        dir: str, 
        filename: str="T_slices",
    ) -> None:
    """
    Tool to export time T slice data to HDF5 file.

    Parameters
    ----------
    dir: str
        Path to directory to write HDF5 file into.
    filename: str="T_slices"
        Name of HDF file.
    """
    file_path: str = os.path.join(dir, f"{filename}.hdf5",)
    file_: hd5File
    with h5py.File(file_path, "w",) as file_:
        i_: int
        for i_ in self.T_slices:
            T_group_ = file_.create_group(f"{self.T_slices[i_]['T']}")
            T_group_.create_dataset("x", data=self.T_slices[i_]["x"])
            T_group_.create_dataset("z", data=self.T_slices[i_]["z"])
            T_group_.create_dataset("β", data=self.T_slices[i_]["β"])

export_metadata 

export_metadata(dir: str, filename: str = 'metadata') -> dict

Tool to simulation metadata to JSON file.

Parameters

dir: str Path to directory to write JSON file into. filename: str="metadata" Name of JSON file.

Source code in erosionfront/levelset/export.py 
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
def export_metadata(
        self, 
        dir: str, 
        filename: str="metadata",
    ) -> dict:
    """
    Tool to simulation metadata to JSON file.

    Parameters
    ----------
    dir: str
        Path to directory to write JSON file into.
    filename: str="metadata"
        Name of JSON file.
    """
    file_path: str = os.path.join(dir, f"{filename}.json",)
    metadata: dict = {}
    for group_ in (".", "domain", "surface", "substrate", "model", "gma",):
        if group_==".":
            object_ = self
            group_ = "sim"
        else:
            object_ = getattr(self, group_)
        metadata[group_] = {}
        for attr_ in object_.__dict__:
            value_ = object_.__dict__[attr_]
            type_ = type(value_)
            if is_serializable(value_):
                if type_ is tuple and not isinstance(value_[0], np.ndarray):
                    metadata[group_].update({
                        attr_: [
                            make_serializable(element_)
                            for element_ in value_
                        ]
                    })
                else:
                    metadata[group_].update({
                        attr_: make_serializable(value_)
                    })
            else:
                if type_ is np.ndarray and value_.size<100:
                    metadata[group_].update({
                        attr_: make_serializable(value_)
                    })
    file_: TextIOWrapper
    with open(file_path, "w",) as file_:
        json.dump(
            metadata, 
            file_, 
            indent=4,
            ensure_ascii=False,
        )
    return metadata

find_shock 

find_shock(
    i_step: int,
    method: ChangepointMethods = BINSEG,
    i_max: int = 300,
    n_filterwidth: int = 50,
    which_changepoint: int = 1,
    do_plot: bool = False,
) -> tuple[float, float]

Locate main break-in-slope by smoothing & changepoint detection.

Parameters

i_step: int

method: ChangepointMethods=ChangepointMethods.BINSEG

i_max: int=300

n_filterwidth: int=50

which_changepoint: int=1

do_plot: bool=False

Returns

tuple[float, float]: x, z coordinate of break-in-slope shock changepoint.

Source code in erosionfront/levelset/shock.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def find_shock(
        self, 
        i_step: int, 
        method: ChangepointMethods=ChangepointMethods.BINSEG,
        i_max: int=300,
        n_filterwidth: int=50,
        which_changepoint: int=1,
        do_plot: bool=False, 
    ) -> tuple[float, float]:
    """
    Locate main break-in-slope by smoothing & changepoint detection.

    Parameters
    ----------
    i_step: int

    method: ChangepointMethods=ChangepointMethods.BINSEG

    i_max: int=300

    n_filterwidth: int=50

    which_changepoint: int=1

    do_plot: bool=False

    Returns
    -------
    tuple[float, float]:
        x, z coordinate of break-in-slope shock changepoint.
    """
    changepoint_fn: ruptures.Pelt | ruptures.Window | ruptures.Binseg
    changepoints: list
    β_max: float = np.pi/2 * 0.95
    β_: NDArray = median_filter(self.T_slices[i_step]["β"], n_filterwidth,)
    i_β_cliff: NDArray = np.where(β_>β_max)[0]
    i_end: int = max(
        i_β_cliff[0] if i_β_cliff.shape[0]>1
        else len(β_), i_max
    )
    β_ = β_[:i_end]
    match method:
        case ChangepointMethods.PELT: 
            changepoint_fn = ruptures.Pelt(model="l2", min_size=30,)
            changepoint_fn.fit(β_)
            changepoints = changepoint_fn.predict(pen=3,)
        case ChangepointMethods.WINDOW:
            changepoint_fn = ruptures.Window(model="l2", width=10,)
            changepoint_fn.fit(β_)
            changepoints = changepoint_fn.predict(n_bkps=2,)
        case ChangepointMethods.BINSEG:
            changepoint_fn = ruptures.Binseg(model="l2", min_size=10,)
            changepoint_fn.fit(β_)
            changepoints = changepoint_fn.predict(n_bkps=2,)
    i_shock: int = changepoints[which_changepoint]
    if do_plot:
        ruptures.display(β_, [], changepoints)
    return (
        self.T_slices[i_step]["x"][i_shock], 
        self.T_slices[i_step]["z"][i_shock],
    )

find_shocks 

find_shocks(
    method: ChangepointMethods = BINSEG,
    subset: slice = s_[100:1000:50],
    i_max: int = 1000,
    n_filterwidth: int = 50,
    which_changepoint=0,
    i_plot: int = 200,
    i_plot_offset: int = 50,
) -> None

Find breaks-in-slope in the set of T-sliced surfaces.

Parameters

method: ChangepointMethods=ChangepointMethods.BINSEG

subset: slice=np.s_[100:1000:50]

i_max: int=1000

n_filterwidth: int=50

which_changepoint=0,

i_plot: int=200

i_plot_offset: int=50

Attributes

Modifies shocks_xz with estimated locations of breaks-in-slope.

Source code in erosionfront/levelset/shock.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
def find_shocks(
        self,
        method: ChangepointMethods=ChangepointMethods.BINSEG,
        subset: slice=np.s_[100:1000:50],
        i_max: int=1000,
        n_filterwidth: int=50,
        which_changepoint=0,
        i_plot: int=200,
        i_plot_offset: int=50,
    ) -> None:
    """
    Find breaks-in-slope in the set of T-sliced surfaces.

    Parameters
    ----------
    method: ChangepointMethods=ChangepointMethods.BINSEG

    subset: slice=np.s_[100:1000:50]

    i_max: int=1000

    n_filterwidth: int=50

    which_changepoint=0,

    i_plot: int=200

    i_plot_offset: int=50

    Attributes
    ----------
    Modifies `shocks_xz` with estimated locations of breaks-in-slope.
    """
    slices: list = list(self.T_slices.keys())[subset]
    self.shocks_xz: NDArray = np.array([
        self.find_shock(
            i_step_, 
            method=method,
            i_max=i_max,
            n_filterwidth=n_filterwidth,
            which_changepoint=which_changepoint,
            do_plot=(
                (i_step_==slices[-1]) | (i_step_-i_plot_offset)%i_plot==0
            ), 
        )
        for i_step_ in slices
    ]).T

first_update 

first_update() -> None

Set up time-stepping solution with requisite initializations.

Source code in erosionfront/levelset/solver.py 
118
119
120
121
122
123
124
125
126
127
128
def first_update(self) -> None:
    """Set up time-stepping solution with requisite initializations."""
    Timer.start("fast-march travel time on signed distance grid",)
    self.β = np.zeros_like(self.φ)
    self.ξ = np.ones_like(self.φ)
    if self.model is None:
        raise ValueError("Erosion model not specified")
    self.φ_next = None
    self.φ_aux = None
    self.T_slices = {}
    Timer.stop()

fm_distance staticmethod 

fm_distance(
    f: NDArray, band_width: float | None, resolution: float, order: int | None
) -> MaskedArray

Use fast marching to compute signed distance from a ref grid (zeroset).

The reference grid (typically travel time, but may be rasterized line offsets) provides a zero contour from which the signed distance is computed by fast marching.

Parameters

f: NDArray Reference grid to provide a zero contour. band_width: float | None, Width of the narrow band. resolution: float Grid pixel width (and height). order: int Fast-marching order of computational stencil (1 or 2).

Returns

MaskedArray Signed distance grid.

Source code in erosionfront/levelset/grid.py 
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
@staticmethod
def fm_distance(
        f: NDArray, 
        band_width: float | None,
        resolution: float,
        order: int | None,
    ) -> MaskedArray:
    """
    Use fast marching to compute signed distance from a ref grid (zeroset).

    The reference grid (typically travel time, but may be rasterized line
    offsets) provides a zero contour from which the signed distance is 
    computed by fast marching.

    Parameters
    ----------
    f: NDArray
        Reference grid to provide a zero contour.
    band_width: float | None,
        Width of the narrow band. 
    resolution: float
        Grid pixel width (and height).
    order: int
        Fast-marching order of computational stencil (1 or 2).

    Returns
    -------
    MaskedArray
        Signed distance grid.
    """
    φ: NDArray | MaskedArray \
          = skfmm.distance(
                f, 
                dx=resolution, 
                order=order, 
                narrow=band_width,
            )
    return (
        φ if type(φ)==MaskedArray
        else np.ma.masked_array(φ, np.zeros_like(φ, dtype=np.bool,))
    )

fm_distance_everywhere 

fm_distance_everywhere(
    φ: MaskedArray, φ_everywhere: NDArray, resolution: float, order: int | None
) -> NDArray

FM compute a narrowband signed distance and then extend everywhere.

The previous 'everywhere' grid of signed distance is used to supply values beyond the narrowband.

Parameters

φ: MaskedArray Reference grid to provide a zero contour. φ_everywhere: NDArray Previous grid of signed distance everywhere. resolution: float Grid pixel width (and height). order: int Fast-marching order of computational stencil (1 or [2]).

Attributes

Uses the fm_distance method.

Returns

MaskedArray Grid with signed distance values everywhere, not just on the nb.

Source code in erosionfront/levelset/grid.py 
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
def fm_distance_everywhere(
        self, 
        φ: MaskedArray, 
        φ_everywhere: NDArray,
        resolution: float,
        order: int | None,
    ) -> NDArray:
    """
    FM compute a narrowband signed distance and then extend everywhere.

    The previous 'everywhere' grid of signed distance is used to supply
    values beyond the narrowband.

    Parameters
    ----------
    φ: MaskedArray
        Reference grid to provide a zero contour.
    φ_everywhere: NDArray
        Previous grid of signed distance everywhere.
    resolution: float
        Grid pixel width (and height).
    order: int
        Fast-marching order of computational stencil (1 or [2]).

    Attributes
    ----------
    Uses the fm_distance method.

    Returns
    -------
    MaskedArray
        Grid with signed distance values everywhere, not just on the nb.
    """
    previous_φ_everywhere: NDArray \
        = φ_everywhere
    previous_distant_pixels: NDArray = φ.mask
    φ_everywhere_: NDArray \
        = self.fm_distance(φ, 0, resolution, order,).data
    φ_everywhere_[previous_distant_pixels] \
        = previous_φ_everywhere[previous_distant_pixels]
    return φ_everywhere_

fwd_back_differences staticmethod 

fwd_back_differences(
    φ: MaskedArray, Δx: float
) -> tuple[MaskedArray, MaskedArray, MaskedArray, MaskedArray]

Compute forward, backward, mean & double finite differences.

Could perhaps use np.roll to make this process more efficient?

Parameters

φ: MaskedArray Signed distance grid. Δx: float Grid spacing in both x and z directions (assumed equal).

Returns

(Δφ_x, Δφ_z, Δ2φ_x2, Δ2φ_z2,): tuple[MaskedArray, MaskedArray, MaskedArray, MaskedArray] 1st-order (center-weighted) and 2nd-order φ finite differences.

Source code in erosionfront/levelset/gradient.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
@staticmethod
def fwd_back_differences(
        φ: MaskedArray,
        Δx: float,
    ) -> tuple[MaskedArray, MaskedArray, MaskedArray, MaskedArray]:
    """
    Compute forward, backward, mean & double finite differences.

    Could perhaps use np.roll to make this process more efficient?

    Parameters
    ----------
    φ: MaskedArray
        Signed distance grid.
    Δx: float
        Grid spacing in both x and z directions (assumed equal).

    Returns
    -------
    (Δφ_x, Δφ_z, Δ2φ_x2, Δ2φ_z2,): tuple[MaskedArray, MaskedArray, \
        MaskedArray, MaskedArray]
        1st-order (center-weighted) and 2nd-order φ finite differences.
    """
    Δz: float = Δx
    φ_xp: MaskedArray = φ.copy()
    φ_xm: MaskedArray = φ.copy()
    φ_zp: MaskedArray = φ.copy()
    φ_zm: MaskedArray = φ.copy()

    # Generate shifted grids to make differencing cleaner
    # Remember arrays are [rows, cols] and that rows=z, cols=x
    φ_xp[:,:-1] = φ[:,1:]
    φ_xm[:,1:]  = φ[:,:-1]
    φ_zp[:-1,:] = φ[1:,:]
    φ_zm[1:,:]  = φ[:-1,:]

    # Forward and backward differences in x and z
    Δφ_xp: MaskedArray = φ_xp - φ
    Δφ_xm: MaskedArray = φ - φ_xm
    Δφ_zp: MaskedArray = φ_zp - φ
    Δφ_zm: MaskedArray = φ - φ_zm

    # Mean differences ~dφ/dx∇∇𝛻𝚺∇
    Δφ_x: MaskedArray = ((Δφ_xp + Δφ_xm)/2) / Δx
    Δφ_z: MaskedArray = ((Δφ_zp + Δφ_zm)/2) / Δz

    # Double differences ~d2φ/dx2
    Δ2φ_x2: MaskedArray = ((Δφ_xp - Δφ_xm)/2) / Δx
    Δ2φ_z2: MaskedArray = ((Δφ_zp - Δφ_zm)/2) / Δz

    return (Δφ_x, Δφ_z, Δ2φ_x2, Δ2φ_z2,)

get_extent 

get_extent() -> NDArray

Get domain (x,z) extent.

Returns

NDArray: Domain bounds as [x_min, x_max, z_max, z_min].

Source code in erosionfront/levelset/elementary.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
def get_extent(self) -> NDArray:
    """
    Get domain (x,z) extent.

    Returns
    -------
    NDArray:
        Domain bounds as [x_min, x_max, z_max, z_min].

    """
    if self.domain is None or self.domain.extent is None:
        raise ValueError("Domain/extent not specified")
    return self.domain.extent

get_front_points 

get_front_points(f: NDArray | MaskedArray) -> NDArray

Fetch points along erosion front surface aka zeroset of f.

Parameters

f: NDArray | MaskedArray Presumably a signed distance grid (but not necessarily).

Attributes

Uses compute_contour method.

Returns

NDArray Coordinate sequence (x,z)*N along surface as (2,N) array.

Source code in erosionfront/levelset/grid.py 
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
def get_front_points(
        self, 
        f: NDArray | MaskedArray,
    ) -> NDArray:
    """
    Fetch points along erosion front surface aka zeroset of f.

    Parameters
    ----------
    f: NDArray | MaskedArray
        Presumably a signed distance grid (but not necessarily).

    Attributes
    ----------
    Uses compute_contour method.

    Returns
    -------
    NDArray
        Coordinate sequence (x,z)*N along surface as (2,N) array.
    """
    return self.compute_contour(f, level=0,)

get_line_points 

get_line_points() -> NDArray

Convert Shapely geometry Line into (x,z) points.

Attributes

Uses line attribute.

Returns

NDArray Line as N coordinates in (2,N) array format, i.e., as two rows of N-length sequences.

Source code in erosionfront/levelset/grid.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def get_line_points(self) -> NDArray:
    """
    Convert Shapely geometry Line into (x,z) points.

    Attributes
    ----------
    Uses line attribute.

    Returns
    -------
    NDArray
        Line as N coordinates in (2,N) array format, 
        i.e., as two rows of N-length sequences.
    """
    return (shapely.get_coordinates(self.line)).T

get_points_xz 

get_points_xz() -> NDArray

Get surface (x,z) coordinates.

Returns

NDArray Points along surface as (x,z) pairs.

Source code in erosionfront/levelset/elementary.py 
134
135
136
137
138
139
140
141
142
143
def get_points_xz(self) -> NDArray:
    """
    Get surface (x,z) coordinates.

    Returns
    -------
    NDArray
        Points along surface as (x,z) pairs.
    """
    return self.points_xz.T

get_substrate_ξ0 

get_substrate_ξ0(φ: MaskedArray) -> MaskedArray

Generate relative η grid (effective speed ξ0).

Parameters

φ: MaskedArray Any grid of the size being used in the sim.

Attributes

substrate: OneLayerSubstrate | TwoLayerSubstrate | MultiLayerSubstrate Instance of class containing substrate info. surface.height: float Top height of surface above z=0 base. η_upperlayer: float Relative erodiblity of strong layer (weak layer assumed = 1). do_hardtop: bool Whether to fake a highly resistant, very thin, upper "caprock". do_hardbase: bool Whether to fake a highly resistant, very thin, lower "bedrock". η_hardtopbase: float Relative erodiblity of base/top layers (weak layer assumed = 1).

Returns

MaskedArray Modulated erodbility (ref speed) grid ξ0 across the narrow band.

Source code in erosionfront/levelset/speed.py 
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def get_substrate_ξ0(
        self, 
        φ: MaskedArray,
    ) -> MaskedArray:
    """
    Generate relative η grid (effective speed ξ0).

    Parameters
    ----------
    φ: MaskedArray
        Any grid of the size being used in the sim.

    Attributes
    ----------
    substrate: OneLayerSubstrate | TwoLayerSubstrate | MultiLayerSubstrate
        Instance of class containing substrate info.
    surface.height: float
        Top height of surface above z=0 base.
    η_upperlayer: float
        Relative erodiblity of strong layer (weak layer assumed = 1).
    do_hardtop: bool
        Whether to fake a highly resistant, very thin, upper "caprock".
    do_hardbase: bool
        Whether to fake a highly resistant, very thin, lower "bedrock".
    η_hardtopbase: float
        Relative erodiblity of base/top layers (weak layer assumed = 1).

    Returns
    -------
    MaskedArray
        Modulated erodbility (ref speed) grid ξ0 across the narrow band.
    """
    substrate: OneLayerSubstrate | TwoLayerSubstrate | MultiLayerSubstrate
    if self.substrate is not None:
        substrate = self.substrate 
    else:
        raise ValueError("Substrate not specified")
    surface: (
        StraightLineSurface | 
        ErrorFunctionSurface | 
        SemiCircularArcSurface | 
        QuarterCircularArcSurface
    )
    if self.surface is not None:
        surface = self.surface 
    else:
        raise ValueError("Surface not specified")
    z_max: float = surface.height
    z_change: Callable = lambda z: int(self.xz_to_ij((0, z,))[1])
    ξ0: MaskedArray = np.ones_like(φ)
    # Prevent erosion of top surface
    if type(substrate)==TwoLayerSubstrate:
        ξ0[z_change(0.5):,:] = substrate.η_upperlayer
    if substrate.do_hardtop:
        # Make the hard top either 5% of height, or 10 pixels high,
        #    - whichever is smaller
        ξ0[z_change(z_max-min(self.Δx*substrate.nz_hardtop, 
                              0.001*substrate.nz_hardtop)):,:] \
            = substrate.η_hardtopbase
    # Prevent erosion of bottom surface and thus motion to left
    if substrate.do_hardbase:
        #BUG: this doesn't work if the initial surface isn't z \in [0,1]
        ξ0[:(z_change(0)),:] = substrate.η_hardtopbase
    # Smear out a bit to stabilize
    ξ0_ = scipy.ndimage.uniform_filter(
        ξ0, 
        size=substrate.nz_hardtop,
    )
    ξ0 = np.ma.masked_array(ξ0_, mask=ξ0.mask,)
    return ξ0

get_x_limits 

get_x_limits() -> tuple

Get domain x min, max.

Returns

tuple: x_min, x_max pair.

Source code in erosionfront/levelset/elementary.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def get_x_limits(self) -> tuple:
    """
    Get domain x min, max.

    Returns
    -------
    tuple:
        x_min, x_max pair.
    """
    if self.domain is None or self.domain.extent is None:
        raise ValueError("Domain/extent not specified")
    return tuple(self.domain.extent[:2])

get_z_limits 

get_z_limits() -> tuple

Get domain z min, max.

Returns

tuple: z_min, z_max pair.

Source code in erosionfront/levelset/elementary.py 
107
108
109
110
111
112
113
114
115
116
117
118
def get_z_limits(self) -> tuple:
    """
    Get domain z min, max.

    Returns
    -------
    tuple:
        z_min, z_max pair.
    """
    if self.domain is None or self.domain.extent is None:
        raise ValueError("Domain/extent not specified")
    return tuple(self.domain.extent[2:][::-1])

get_φ_next 

get_φ_next() -> MaskedArray

Return next signed-distance grid φ_next or just φ if not available.

Returns

MaskedArray: Next signed-distance φ grid or just current φ if next unavailable.

Source code in erosionfront/levelset/speed.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
def get_φ_next(self) -> MaskedArray:
    """
    Return next signed-distance grid φ_next or just φ if not available.

    Returns
    -------
    MaskedArray:
        Next signed-distance φ grid or just current φ if next unavailable.
    """
    return (
        self.φ if self.φ_next is None
        else self.φ_next
    )

ij_to_xz 

ij_to_xz(ij: NDArray | tuple) -> NDArray

Convert pixel indices [i,j] into pixel center coords (x,z).

Parameters

ij: NDArray | tuple Pixel index pairs [i,j] as numpy array or tuple.

Returns

NDArray: Pixel (x,z) coordinate pairs as float64 numpy array.

Source code in erosionfront/levelset/elementary.py 
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def ij_to_xz(self, ij: NDArray | tuple,) -> NDArray:
    """
    Convert pixel indices [i,j] into pixel *center* coords (x,z).

    Parameters
    ----------
    ij: NDArray | tuple
        Pixel index pairs [i,j] as numpy array or tuple.

    Returns
    ------- 
    NDArray:
        Pixel (x,z) coordinate pairs as float64 numpy array.
    """
    if self.domain is None or (origin:=self.domain.origin()) is None:
        raise ValueError("Domain/origin not specified")
    return (
        self.Δx * np.array(ij).astype(np.float64) 
        + origin 
        + np.array([1,1])*self.Δx/2
    )

initialize 

initialize(do_timing: bool = False) -> None

Initialize level-set solver.

Parameters

do_timing: bool=False Perform timing estimates for each computational stage.

Attributes

Many class attributes are modified.

Source code in erosionfront/levelset/base.py 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
def initialize(self, do_timing: bool=False,) -> None:
    """
    Initialize level-set solver.

    Parameters
    ----------
    do_timing: bool=False
        Perform timing estimates for each computational stage.

    Attributes
    ----------
    Many class attributes are modified.
    """
    Timer(do_timing, "Initializing",)

    Timer.start("create initial surface line; generate subsurface polygon",)
    self.domain = copy(self.raw_domain)
    if self.surface is None:
        raise ValueError("Surface not specified")
    self.points_xz = self.surface.get_points_xz()
    self.line = self.make_line(self.points_xz)
    if self.n_pixels_xz is None:
        raise ValueError("Grid size not specified")
    self.in_rock_polygon \
        = self.make_polygon(
            self.points_xz, 
            self.Δx,
            self.n_pixels_xz, 
        )
    Timer.stop()

    Timer.start("rasterize line",)
    self.ρ = self.rasterize_line(self.points_xz, self.do_extend_line,)
    Timer.stop()
    Timer.start("pad the grid vertically",)
    if self.raw_domain is None:
        raise ValueError("Domain not specified")
    if self.n_pad_pixels is None:
        raise ValueError("Number of pad pixels not specified")
    self.domain = self.pad_domain(
        self.raw_domain, 
        self.Δx,
        self.n_pad_pixels,
    )
    self.n_pixels_xz \
        = self.compute_pixel_dimensions(self.domain, self.Δx,)
    self.n_pixels_xz0 \
        = self.compute_pixel_origin(self.domain, self.Δx,)
    if self.n_pad_pixels is None:
        raise ValueError("Number of padding pixels not specified")
    self.ρ_padded \
        = self.pad_grid(
            grid=self.ρ, 
            n_pad_pixels=self.n_pad_pixels,
        )
    Timer.stop()

    Timer.start("measure pixel offsets from line",)
    self.ρ_offsets \
        = self.compute_line_pixel_offsets(
            line=self.line,
            ρ=self.ρ_padded,
            do_φ_everywhere=self.do_φ_everywhere,
        )
    Timer.stop()

    Timer.start("fast-march signed distance from pixel offsets",)
    self.φ \
        = self.fm_distance(
            f=self.ρ_offsets, 
            band_width=self.band_width,
            resolution=self.Δx,
            order=self.fm_order,
        )
    Timer.stop()

    Timer.start("fill whole signed-distance grid with very-far values",)
    self.φ_everywhere \
        = self.redo_φ_everywhere(
            φ_everywhere=self.ρ_offsets, 
            φ=self.φ,
        )
    Timer.stop()

ji_to_xz 

ji_to_xz(ji: NDArray | tuple) -> NDArray

Convert reversed pixel indices [j,i] into pixel center coords (x,z).

Parameters

ji: NDArray | tuple Pixel (reversed) index pairs [j,i] as numpy array or tuple.

Returns

NDArray: Pixel (x,z) coordinate pairs as float64 numpy array.

Source code in erosionfront/levelset/elementary.py 
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def ji_to_xz(self, ji: NDArray | tuple,) -> NDArray:
    """
    Convert reversed pixel indices [j,i] into pixel *center* coords (x,z).

    Parameters
    ----------
    ji: NDArray | tuple
        Pixel (reversed) index pairs [j,i] as numpy array or tuple.

    Returns
    ------- 
    NDArray:
        Pixel (x,z) coordinate pairs as float64 numpy array.

    """
    return self.ij_to_xz(np.array(ji)[::-1])

make_line staticmethod 

make_line(points_xz: NDArray) -> LineString

Convert surface (x,z) coordinates into a Shapely Line geometry.

Parameters

points_xz: NDArray Stacked vectors of x and z surface point coordinates.

Returns

LineString: Shapely LineString geometry of surface line.

Source code in erosionfront/levelset/elementary.py 
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
@staticmethod
def make_line(points_xz: NDArray,) -> LineString:
    """
    Convert surface (x,z) coordinates into a Shapely Line geometry.

    Parameters
    ----------
    points_xz: NDArray
        Stacked vectors of x and z surface point coordinates.

    Returns
    -------
    LineString:
        Shapely LineString geometry of surface line.
    """
    return LineString(points_xz)

make_polygon staticmethod 

make_polygon(
    points_xz: NDArray, resolution: float, n_pixels_xz: tuple | NDArray
) -> Polygon

Create a 'bounding' polygon from surface x,z points.

This polygon is used to decide whether a point lies within the bedrock or outside it.

BUG: Crudely done. Should be much tighter geometry close to domain.

Parameters

points_xz: NDArray Surface coordinates n_pixels_xz: tuple | NDArray Grid dimensions resolution: float Grid pixel size.

Returns

Polygon Shapely geometry of bounding polygon.

Source code in erosionfront/levelset/elementary.py 
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
@staticmethod
def make_polygon(
        points_xz: NDArray,
        resolution: float,
        n_pixels_xz: tuple | NDArray,
    ) -> Polygon:
    """
    Create a 'bounding' polygon from surface x,z points.

    This polygon is used to decide whether a point lies within the bedrock
    or outside it.

    BUG: Crudely done. Should be much tighter geometry close to domain.

    Parameters
    ----------
    points_xz: NDArray
        Surface coordinates
    n_pixels_xz: tuple | NDArray
        Grid dimensions
    resolution: float
        Grid pixel size.

    Returns
    -------
    Polygon
        Shapely geometry of bounding polygon.
    """
    if n_pixels_xz is None:
        raise ValueError("Grid size not specified")
    x_offset: float = n_pixels_xz[0]*resolution*10
    z_offset: float = n_pixels_xz[1]*resolution*10
    left_anchor: NDArray \
        = np.array([points_xz[0][0]-x_offset, 
                    points_xz[0][1]-z_offset])
    right_anchor: NDArray \
        = np.array([points_xz[-1][0]+x_offset, 
                    points_xz[-1][1]-z_offset])
    poly_points: NDArray \
        = np.vstack((
            left_anchor, 
            points_xz[0]-np.array([resolution*10,0]),
            points_xz, 
            points_xz[-1]+np.array([resolution*10,0]),
            right_anchor,
        ))
    return Polygon(poly_points)

measure_normal_distance staticmethod 

measure_normal_distance(
    line: tuple | NDArray | LineString, point: tuple | NDArray | Point
) -> float

Compute distance between line and point using Shapely.

Parameters

line: tuple | NDArray | LineString Line represented by various possible types. point: tuple | NDArray | Point Coordinate (x,z) represented by various possible types.

Returns

float Perpendicular distance from line curve to point.

Source code in erosionfront/levelset/grid.py 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
@staticmethod
def measure_normal_distance(
        line: tuple | NDArray | LineString, 
        point: tuple | NDArray | Point,
    ) -> float:
    """
    Compute distance between line and point using Shapely.

    Parameters
    ----------
    line: tuple | NDArray | LineString
        Line represented by various possible types.
    point: tuple | NDArray | Point
        Coordinate (x,z) represented by various possible types.

    Returns
    -------
    float
        Perpendicular distance from line curve to point.
    """
    return float(shapely.distance(
        line if type(line) is LineString else LineString(line), 
        point if type(point) is Point else Point(point)
    ))

measure_φ_ji 

measure_φ_ji(line: LineString, ji: tuple[int, int] | NDArray) -> float

Compute signed distance between line and i,j pixel using Shapely.

The sign is -1 for "in the air" and +1 for "in the rock".

Parameters

line: LineString Surface line. ji: tuple[int,int] | NDArray Pixel indexes in reversed format.

Attributes

Uses measure_φ_xz method to test which side of the surface the pixel lies.

Returns

float Signed perpendicular distance from the line to the pixel.

Source code in erosionfront/levelset/grid.py 
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
def measure_φ_ji(
        self,
        line: LineString,
        ji: tuple[int,int] | NDArray,
    ) -> float:
    """
    Compute signed distance between line and i,j pixel using Shapely.

    The sign is -1 for "in the air" and +1 for "in the rock".

    Parameters
    ----------
    line: LineString
        Surface line.
    ji: tuple[int,int] | NDArray
        Pixel indexes in reversed format.

    Attributes
    ----------
    Uses measure_φ_xz method to test which side of the 
    surface the pixel lies.

    Returns
    -------
    float
        Signed perpendicular distance from the line to the pixel.
    """
    return self.measure_φ_xz(line, self.ji_to_xz(ji),)

measure_φ_xz 

measure_φ_xz(line: LineString, xz: tuple[float, float] | NDArray) -> float

Compute signed distance between line and xz point using Shapely.

The sign is -1 for "in the air" and +1 for "in the rock".

Parameters

line: LineString Surface line. xz: tuple[float,float] | NDArray Point coordinates.

Attributes

Uses above_or_below_xz method to test which side of the surface the point lies.

Returns

float Signed perpendicular distance from the line to the point.

Source code in erosionfront/levelset/grid.py 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
def measure_φ_xz(
        self,
        line: LineString,
        xz: tuple[float,float] | NDArray,
    ) -> float:
    """
    Compute signed distance between line and xz point using Shapely.

    The sign is -1 for "in the air" and +1 for "in the rock".

    Parameters
    ----------
    line: LineString
        Surface line.
    xz: tuple[float,float] | NDArray
        Point coordinates.

    Attributes
    ----------
    Uses above_or_below_xz method to test which side of the surface the 
    point lies.

    Returns
    -------
    float
        Signed perpendicular distance from the line to the point.
    """
    return (
        self.measure_normal_distance(line, xz)
        *self.above_or_below_xz(xz, self.in_rock_polygon,)
    )

pad_domain staticmethod 

pad_domain(domain: Domain, resolution: float, n_pad_pixels: int) -> Domain

Expand the domain vertically by pixel width amounts.

Parameters

domain: Domain (Initial) domain to be padded. resolution: float Grid pixel size. n_pad_pixels: int Number of pixel widths to pad by below the current domain.

Returns

New domain instance padded vertically.

Source code in erosionfront/levelset/elementary.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
@staticmethod
def pad_domain(
        domain: Domain, 
        resolution: float,
        n_pad_pixels: int,
    ) -> Domain:
    """
    Expand the domain vertically by pixel width amounts.

    Parameters
    ----------
    domain: Domain
        (Initial) domain to be padded.
    resolution: float
        Grid pixel size.
    n_pad_pixels: int
        Number of pixel widths to pad by *below* the current domain.

    Returns
    -------
    New domain instance padded vertically.
    """
    dz_bot: float = (n_pad_pixels)*resolution
    dz_top: float = (n_pad_pixels)*resolution
    # dz_top: float = resolution*5 if n_pad_pixels>0 else 0
    n_digits: int = int(round(-np.log10(resolution)+0.5))+1
    return Domain(
        round(domain.x_min, n_digits),
        round(domain.x_max, n_digits),
        round(domain.z_min - dz_bot, n_digits),
        round(domain.z_max + dz_top, n_digits),
    )

pad_grid staticmethod 

pad_grid(grid: NDArray, n_pad_pixels: int) -> NDArray

Pad grid with specified number of pixels below and above.

Parameters

grid: NDArray Grid to be padded. n_pad_pixels: int Number of pixels to pad below.

Returns

NDArray Padded grid.

Source code in erosionfront/levelset/elementary.py 
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
@staticmethod
def pad_grid(
        grid: NDArray, 
        n_pad_pixels: int, 
    ) -> NDArray:
    """
    Pad grid with specified number of pixels below and above.

    Parameters
    ----------
    grid: NDArray
        Grid to be padded.
    n_pad_pixels: int
        Number of pixels to pad below.

    Returns
    -------
    NDArray
        Padded grid.
    """
    pad_ztop: NDArray \
        = np.zeros((n_pad_pixels,grid.shape[1],), dtype=np.int64,)
    pad_zbot: NDArray \
        = np.zeros((n_pad_pixels,grid.shape[1],), dtype=np.int64,)
    expanded_grid: NDArray 
    if n_pad_pixels>0:
        expanded_grid =  np.vstack((pad_zbot, grid, pad_ztop,))
    else:
        expanded_grid = grid
    return expanded_grid

rasterize_T_slices 

rasterize_T_slices(n_subset: int = 1) -> None

Rasterize time T slices onto grids.

Parameters

n_subset: int=1 Rate of subsetting of T slices (n=1 means all of them).

Source code in erosionfront/levelset/base.py 
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
def rasterize_T_slices(self, n_subset: int=1,) -> None:
    """
    Rasterize time T slices onto grids.

    Parameters
    ----------
    n_subset: int=1
        Rate of subsetting of T slices (n=1 means all of them).
    """
    print("Gridding T slices...", end="",)
    if self.domain is None:
        raise ValueError("Domain not defined")
    bounds: NDArray | None= self.domain.bounds()
    if bounds is None:
        raise ValueError("Domain bounds not defined")
    if self.n_pixels_xz is None:
        raise ValueError("Number of x,z pixels not defined")
    n_x: int = self.n_pixels_xz[0]
    n_z: int = self.n_pixels_xz[1]
    x_bounds: NDArray = (bounds.T)[0]
    z_bounds: NDArray = (bounds.T)[1]
    x_pts: NDArray = np.linspace(x_bounds[0], x_bounds[1], n_x,) + self.Δx/2
    z_pts: NDArray = np.linspace(z_bounds[0], z_bounds[1], n_z,) + self.Δx/2
    assert x_pts.shape[0]==n_x
    assert z_pts.shape[0]==n_z
    x_mg: NDArray
    z_mg: NDArray
    (z_mg, x_mg,) = np.meshgrid(z_pts, x_pts, indexing="ij",)
    x_points: NDArray = np.concatenate(
        [slice_["x"] for slice_ in self.T_slices.values()]
    )[::n_subset]
    z_points: NDArray = np.concatenate(
        [slice_["z"] for slice_ in self.T_slices.values()]
    )[::n_subset]
    T_values: NDArray = np.concatenate(
        [slice_["x"]*0+slice_["T"] for slice_ in self.T_slices.values()]
    )[::n_subset]
    zx_points: NDArray = np.vstack((z_points, x_points,)).T
    assert type(zx_points) is np.ndarray
    assert type(T_values) is np.ndarray
    T_arrival: NDArray \
        = griddata(zx_points, T_values, (z_mg, x_mg), method="linear",)
    not_reached: NDArray = np.zeros_like(self.φ_everywhere).astype(np.bool)
    not_reached[self.φ_everywhere<0] = True
    self.T_arrival = np.ma.masked_array(T_arrival, mask=not_reached,)
    print(" done")

rasterize_line 

rasterize_line(points_xz: NDArray, do_extend_line: bool = True) -> NDArray

Rasterize surface line as x,z points into grid pixels using Rasterio.

The x,z bounds of the raster grid are precisely set by the bounding box of the line.

Parameters

points_xz: NDArray Surface line as a series of x,z points. do_extend_line: bool=True Artificially extend the surface line so it reaches the boundary.

Attributes

Uses n_pixels_xz and xz_to_ij methods.

Returns

NDArray Grid of rasterized line.

Source code in erosionfront/levelset/grid.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def rasterize_line(
        self, 
        points_xz: NDArray,
        do_extend_line: bool=True,
    ) -> NDArray:
    """
    Rasterize surface line as x,z points into grid pixels using Rasterio.

    The x,z bounds of the raster grid are precisely set by the
    bounding box of the line.

    Parameters
    ----------
    points_xz: NDArray
        Surface line as a series of x,z points.
    do_extend_line: bool=True
        Artificially extend the surface line so it reaches the boundary.

    Attributes
    ----------
    Uses n_pixels_xz and xz_to_ij methods.

    Returns
    -------
    NDArray
        Grid of rasterized line.
    """
    #HACK: artificially extend line horizontally to ensure
    #      we get correct bounds
    points_xz_extended: NDArray
    if do_extend_line:
        points_xz_extended  = np.vstack((
            [2*points_xz[0][0],0],
            points_xz,
            [2*points_xz[-1][0],1],
        ))
    else:
        points_xz_extended = points_xz.copy()
    pts_ij = self.xz_to_ij(points_xz_extended)
    line: LineString = LineString(pts_ij)
    if self.n_pixels_xz is None:
        raise ValueError("Grid size not specified")
    # Specific caching actually slows rasterization down!
    #   with rasterio.Env(GDAL_CACHEMAX=4_096_000_000) as env:
    ρ: NDArray \
        = rasterio.features.rasterize(
            [line],
            out_shape=self.n_pixels_xz[::-1], 
            all_touched=True, 
            fill=0,
        )
    # assert ρ.shape==self.n_pixels_xz[::-1]
    return ρ #, line

record_T_slice 

record_T_slice(step: int, t_total: float) -> None

Record selected sim data for a single time step.

Recorded data are placed in a 'T_slice' dictionary.

Parameters

step: int Simulation step number. t_total: float Total time in simulation so far.

Attributes

Data for this sim time step are appended to the time-slicing dict.

Source code in erosionfront/levelset/base.py 
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
def record_T_slice(
        self, 
        step: int, 
        t_total: float, 
    ) -> None:
    """
    Record selected sim data for a single time step.

    Recorded data are placed in a 'T_slice' dictionary.

    Parameters
    ----------
    step: int
        Simulation step number.
    t_total: float
        Total time in simulation so far.

    Attributes
    ----------
    Data for this sim time step are appended to the time-slicing dict.
    """
    if self.φ_next is None:
        d = self.φ/np.max(self.φ)
    else:
        d = self.φ_next/np.max(self.φ_next)
    # d: NDArray | MaskedArray = self.φ_everywhere/np.max(self.φ_everywhere)
    xz_profile: NDArray = self.get_front_points(d)
    x: NDArray = xz_profile[0]
    z: NDArray = xz_profile[1]
    β: NDArray = np.arctan2(np.gradient(z), np.gradient(x))
    self.T_slices[step] = {
        "x" : x,
        "z" : z,
        "φ": self.φ,
        "β" : β,
        "T" : float(t_total),
    }

redo_φ_everywhere staticmethod 

redo_φ_everywhere(φ_everywhere: NDArray, φ: MaskedArray) -> NDArray

Reassign signed distance values to pixels beyond the narrow band.

Assumes that pixels beyond the narrow band have signed-distance signs that are currently correct: given these signs, each beyond-nb pixel is assigned either the min or max value of the nb signed-distance values.

Parameters

φ_everywhere: NDArray Beyond narrow-band signed-distance grid to be recomputed. φ: MaskedArray Narrow-band signed-distance grid.

Returns

NDArray Recomputed signed-distances beyond narrow band.

Source code in erosionfront/levelset/grid.py 
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
@staticmethod
def redo_φ_everywhere(
        φ_everywhere: NDArray, 
        φ: MaskedArray,
    ) -> NDArray:
    """
    Reassign signed distance values to pixels beyond the narrow band.

    Assumes that pixels beyond the narrow band have signed-distance signs
    that are currently correct: given these signs, each beyond-nb pixel is 
    assigned either the min or max value of the nb signed-distance values.

    Parameters
    ----------
    φ_everywhere: NDArray
        Beyond narrow-band signed-distance grid to be recomputed.
    φ: MaskedArray
        Narrow-band signed-distance grid.

    Returns
    -------
    NDArray
        Recomputed signed-distances beyond narrow band.
    """
    mask: NDArray = np.array(φ.mask)
    far_recomputed: NDArray = φ_everywhere.copy()
    far_recomputed[φ_everywhere<0] = np.min(φ)
    far_recomputed[φ_everywhere>0] = np.max(φ)
    far_recomputed[~mask] = φ[~mask]
    return far_recomputed

update 

update(
    Δt: float,
    Δt_reinitialize: float | None = None,
    do_finalize: bool = False,
    do_timing: bool = False,
) -> None

Update grids for one time step of duration Δt.

The basic cycle needs to be:

set t=0, set up grids

  • loop while t<=T_limit
    • sometimes reinitialize signed distance φ(t) by fast marching
    • back & forward gradients (Δφ_xp, Δφ_xm), (Δφ_zp, Δφ_zm)
    • mean grad Δφ_x, Δφ_z = (Δφ_xp+Δφ_xm)/2, (Δφ_zp+Δφ_zm)/2
    • dbl-diff Δ2φ_x2, Δ2φ_z2 = (Δφ_xp-Δφ_xm)/2, (Δφ_zp-Δφ_zm)/2
    • slope angle from mean grad β = arctan2(Δφ_z, Δφ_x)
    • front-normal speed from slope angle ξ(β)
    • level-set Hamiltonian lsH(Δφ_x,Δφ_z) = ξ(β).|(Δφ_x, Δφ_z)|
    • Hamiltonian derivatives lsH_{∂φ/∂x}, lsH_{∂φ/∂z} wrt ∂φ/∂x, ∂φ/∂z
    • numerical viscosities κx, κz = max(lsH_{∂φ/∂x}), max(lsH_{∂φ/∂z})
    • update rate dφ/dt = - lsH(Δφ_x, Δφ_z) - κx.Δ2φ_x2 - κz.Δ2φ_z2
    • signed distance change Δφ = (∂φ/∂t).Δt
    • signed distance update φ(t+Δt) = φ(t) + Δφ
    • front T(x,z) from zero level-set {x,z} | φ(t)=0
    • update time t -> t+Δt
  • repeat

Parameters

Δt: float Time step. Δt_reinitialize: float | None=None Periodicity of reinitialization of φ signed-distance grid. do_finalize: bool=False Just finish up, don't do another loop of computation. do_timing: bool=False Measure and report compute times taken at steps during the loop.

Attributes

Most of the sim grids are modified.

Source code in erosionfront/levelset/solver.py 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
def update(
        self, 
        Δt: float, 
        Δt_reinitialize: float | None=None,
        do_finalize: bool=False,
        do_timing: bool=False,
    ) -> None:
    """
    Update grids for one time step of duration Δt.

    The basic cycle needs to be:

    set t=0, set up grids

    - **loop while t<=T_limit**
        - _sometimes_ reinitialize signed distance φ(t) by fast marching
        - back & forward gradients (Δφ_xp, Δφ_xm), (Δφ_zp, Δφ_zm)
        - mean grad  Δφ_x, Δφ_z = (Δφ_xp+Δφ_xm)/2, (Δφ_zp+Δφ_zm)/2
        - dbl-diff  Δ2φ_x2, Δ2φ_z2 = (Δφ_xp-Δφ_xm)/2, (Δφ_zp-Δφ_zm)/2
        - slope angle from mean grad  β = arctan2(Δφ_z, Δφ_x)
        - front-normal speed from slope angle ξ(β)
        - level-set Hamiltonian  lsH(Δφ_x,Δφ_z) = ξ(β).|(Δφ_x, Δφ_z)|
        - Hamiltonian derivatives lsH_{∂φ/∂x}, lsH_{∂φ/∂z} wrt ∂φ/∂x, ∂φ/∂z
        - numerical viscosities  κx, κz = max(lsH_{∂φ/∂x}), max(lsH_{∂φ/∂z})
        - update rate  dφ/dt = - lsH(Δφ_x, Δφ_z) - κx.Δ2φ_x2 - κz.Δ2φ_z2
        - signed distance change  Δφ = (∂φ/∂t).Δt
        - signed distance update  φ(t+Δt) = φ(t) + Δφ
        - front T(x,z) from zero level-set  {x,z} | φ(t)=0
        - update time  t -> t+Δt
    - **repeat** 

    Parameters
    ----------
    Δt: float
        Time step.
    Δt_reinitialize: float | None=None
        Periodicity of reinitialization of φ signed-distance grid.
    do_finalize: bool=False
        Just finish up, don't do another loop of computation.
    do_timing: bool=False
        Measure and report compute times taken at steps during the loop.

    Attributes
    ----------
    Most of the sim grids are modified.
    """
    if self.model is None:
        raise ValueError("Model not specified")

    ################### propagate front by Δt #####################
    Timer(do_timing, "Evolving")

    ################### get ready #####################
    Timer.start("getting ready",)
    n_reinitialize: int | None = (
        int(np.round(Δt_reinitialize/Δt))
        if Δt>0 and Δt_reinitialize is not None
        else None
    )
    if not do_finalize:
        self.φ = self.get_φ_next()
    Timer.stop()

    ################### 1st & 2nd order gradients of φ ###############
    Timer.start("compute gradients  ∂φ/∂x, ∂φ/∂z, ∂2φ/∂x2, ∂2φ/∂z2",)
    (self.Δφ_x, self.Δφ_z, self.Δ2φ_x2, self.Δ2φ_z2,) = \
        self.fwd_back_differences(self.φ, self.Δx,)
    assert type(self.Δφ_x) is MaskedArray
    assert type(self.Δφ_z) is MaskedArray
    assert type(self.Δ2φ_x2) is MaskedArray
    assert type(self.Δ2φ_z2) is MaskedArray
    Timer.stop()

    ################ convert gradient φ into surface angle β ############
    Timer.start("gradient ∂φ/∂x, ∂φ/∂z -> angle β",)
    self.β = -np.arctan2(self.Δφ_x, self.Δφ_z,)  #type: ignore
    self.β[~np.isfinite(self.β)] = 0
    assert type(self.β) is MaskedArray
    Timer.stop()

    if do_finalize:
        return

    ################ convert angle β into front speed ξ ############
    Timer.start("angle β -> model ξ",)
    ξ_model: float | NDArray | MaskedArray = self.model.ξ_fn_β(self.β,)
    if type(ξ_model) is not np.ma.masked_array:
        raise ValueError("ξ_model is not a masked array")
    else:
       self.ξ_model = ξ_model
    # assert type(self.ξ_model) is MaskedArray
    Timer.stop()

    ################ get ξ0 grid, compute ξ, combine ############
    Timer.start("get ξ0 grid, compute ξ, combine",)
    self.ξ0 = self.get_substrate_ξ0(self.φ)
    self.ξ = self.ξ0 * self.ξ_model
    assert type(self.ξ0) is MaskedArray
    assert type(self.ξ) is MaskedArray
    Timer.stop()

    ################### level-set Hamiltonian #####################
    Timer.start("model ξ -> level-set Hamiltonian")
    # Remove the mask!
    self.H_ls = self.ξ.data.copy()
    assert type(self.H_ls) is np.ndarray
    Timer.stop()

    ################### ls-Hamiltonian derivatives ###################
    Timer.start("level-set Hamiltonian derivatives")
    # Need to scale by substrate speeds ξ0
    self.dHlsdφx = self.ξ0*self.model.dHlsdφ_fn_Δφx(self.Δφ_x, self.Δφ_z,)
    self.dHlsdφz = self.ξ0*self.model.dHlsdφ_fn_Δφz(self.Δφ_x, self.Δφ_z,)
    assert type(self.dHlsdφx) is MaskedArray
    assert type(self.dHlsdφz) is MaskedArray
    Timer.stop()

    ################### numerical viscosities ###################
    Timer.start("numerical viscosities κx, κz")
    self.κx = np.max(self.dHlsdφx)
    self.κz = np.max(self.dHlsdφz)
    Timer.stop()

    ################### compute dφdt, next φ ######################
    Timer.start("compute dφdt, next φ")
    self.dφdt = (
        self.H_ls + 
        (self.κx*self.Δ2φ_x2 + self.κz*self.Δ2φ_z2)*self.κ_boost
        * (
            1 +
              np.heaviside(self.Δ2φ_x2 + self.Δ2φ_z2, 0)
            * np.heaviside(np.abs(self.β)-(np.pi/2)*0.5, 0)
            * 0
        )
    )
    self.φ_next = self.φ + self.dφdt*self.Δt
    assert type(self.dφdt) is MaskedArray
    assert type(self.φ_next) is MaskedArray
    Timer.stop()

    ################ signed distance φ everywhere #################
    Timer.start("compute signed distance φ everywhere",)
    self.φ_everywhere = self.fm_distance_everywhere(
        self.φ_next,
        self.φ_everywhere,
        self.Δx,
        self.fm_order if self.fm_order is not None else 2,
    )
    assert type(self.φ_everywhere) is np.ndarray
    Timer.stop()

    ################# φ thresholding, slab failure ###################
    if (
        self.model.do_slabfailure
        and self.i_step is not None
        and self.i_step>0 
        and self.n_slabfailure is not None 
        and (self.i_step % self.n_slabfailure==0)
    ):
        self.φ_aux = self.φ_everywhere.copy()
        self.φ_aux[self.Δφ_z>0] = 0
        # Need to pick a minimum distance from which to trigger slab failure
        #   - the 3Δx here is a heuristic minimum
        #   - depends critically on κ_boost, here assumed=1
        #   - the maximum is arbitrary
        φ_thresholds: tuple[float,float] = (self.Δx*3, 1,)
        self.φ_aux[(self.φ<φ_thresholds[0]) | (self.φ>=φ_thresholds[1])] = 0
        self.φ_aux[self.φ_aux!=0] = 1
        i_φ_break = np.argmax(self.φ_aux, axis=0)
        for j_, i_ in enumerate(i_φ_break):
            if i_>0:
                self.φ_everywhere[i_:,j_] = φ_thresholds[0]

    ################# possibly reinitialize φ_next ###################
    if (
        n_reinitialize is not None 
        and self.i_step is not None
        and (self.i_step % n_reinitialize)==0
    ):
        Timer.start("reinitializing φ_next",)
        # print(f"Reinitializing: t={self.i_step}  t={self.t_total}")
        self.φ_next = self.fm_distance(
            self.φ_everywhere, 
            self.band_width if self.band_width is not None else 50,
            self.Δx,
            self.fm_order if self.fm_order is not None else 2,
        )
        Timer.stop()
    elif self.i_step is not None and self.i_step>0:
        self.φ_next = np.ma.masked_array(self.φ_next, mask=self.φ.mask,)
    assert type(self.φ_next) is MaskedArray

    ################# advance ###################
    self.t_total += Δt

xz_to_ij 

xz_to_ij(xz: NDArray | tuple) -> NDArray

Convert pixel center coords (x,z) into pixel indices [i,j].

Parameters

xz: NDArray | tuple Pixel (x,z) coordinate pairs as numpy array or tuple.

Returns

NDArray: Pixel index pairs [i,j] as int64 numpy array

Source code in erosionfront/levelset/elementary.py 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
def xz_to_ij(self, xz: NDArray | tuple,) -> NDArray:
    """
    Convert pixel *center* coords (x,z) into pixel indices [i,j].

    Parameters
    ----------
    xz: NDArray | tuple
        Pixel (x,z) coordinate pairs as numpy array or tuple.

    Returns
    ------- 
    NDArray:
        Pixel index pairs [i,j] as int64 numpy array
    """
    if self.domain is None or (origin:=self.domain.origin()) is None:
        raise ValueError("Domain/origin not specified")
    return np.int64(  #type: ignore
        np.array(xz)/self.Δx 
        - origin/self.Δx
    )