API Reference¶

Primary solvers¶

Both solvers follow the same lifecycle:

flex.initialize()
flex.run()
w = flex.w          # read deflection before finalize clears it
flex.output()       # optional: save to file or display plots
flex.finalize()     # releases w, qs, and the coefficient matrix

Warning

finalize() deletes w, qs, and the cached coefficient matrix. Read w (and call output() if needed) before calling finalize. Accessing w afterwards raises AttributeError.

class gflex.F2D(filename=None)[source]¶

Bases: Flexure

Two-dimensional lithospheric flexure solver.

Computes the deflection w(x, y) of a thin elastic plate overlying an inviscid fluid (mantle) given a surface load stress qs. Supports spatially variable elastic thickness Te.

Set instance attributes, then call initialize(), run(), and finalize() in sequence. Read flex.w before calling finalize(); finalize clears all model state including w.

method¶

Solution method. 'fd' (finite difference, supports variable Te), 'fft' (spectral, requires scalar Te; 2-D only), 'sas' (superposition of analytical solutions, constant Te only), or 'sas_ng' (SAS on an ungridded point cloud).

Type:: str

solver¶

Linear solver: 'direct' (sparse LU, default).

Type:: str

g¶

Gravitational acceleration [m s⁻²].

Type:: float

E¶

Young’s modulus [Pa].

Type:: float

nu¶

Poisson’s ratio.

Type:: float

rho_m¶

Mantle density [kg m⁻³].

Type:: float

rho_fill¶

Infill material density [kg m⁻³] (0 for air, ~1000 for water, ~2700 for rock).

Type:: float

T_e¶

Elastic thickness [m]. A scalar is broadcast to the full grid.

Type:: float or ndarray of shape (M, N)

qs¶

Surface load stress [Pa].

Type:: ndarray of shape (M, N)

dx¶

Grid spacing in the x (column) direction [m].

Type:: float

dy¶

Grid spacing in the y (row) direction [m].

Type:: float

bc_west, bc_east, bc_north, bc_south

Boundary conditions on the west, east, north, and south edges. FD options: 'zero_displacement_zero_slope' (alias 'clamped'), 'zero_displacement_zero_moment' (alias 'pinned'), 'zero_moment_zero_shear' (alias 'free'), 'zero_slope_zero_shear' (alias 'mirror'), 'periodic', 'no_outside_loads' (auto-pad by one flexural wavelength and apply 'zero_displacement_zero_slope' at the new outer edges; self.w is trimmed to the original domain). SAS option: 'no_outside_loads' (the default when unset). FFT: each opposite pair (west/east, north/south) is treated independently — set both sides of a pair to 'periodic' for exact periodicity along that axis, or leave them unset (or set to 'no_outside_loads') to zero-pad that axis. Mixing periodic and non-periodic axes is valid. Setting only one side of a pair to 'periodic' raises a UserWarning and falls back to zero-padding for that axis.

Type:: str

sigma_xx¶

Normal in-plane stress in the x-direction \(\sigma_{xx}\) [Pa]. Supported by fd and fft. Default 0.

Type:: float, optional

sigma_yy¶

Normal in-plane stress in the y-direction \(\sigma_{yy}\) [Pa]. Supported by fd and fft. Default 0.

Type:: float, optional

sigma_xy¶

In-plane shear stress \(\sigma_{xy}\) [Pa]. Supported by fd and fft. Default 0.

Type:: float, optional

fft_pad_n_alpha¶

Number of 2-D flexural-parameter units (α₂D = (D/Δρg)^0.25) to zero-pad on each side for non-periodic FFT runs. Periodic images of the load are separated by 2 × fft_pad_n_alpha × α₂D. Default 4 (8α₂D total separation). Ignored when method != 'fft' or when both sides of the relevant axis are 'periodic'.

Type:: int or float

cache_factorization¶

Controls LU factorisation caching for the FD 'direct' solver. False (default) — re-factorises on every run() call. True — caches the LU factorisation and reuses it (the coefficient matrix is freed once the factors are built). Reuse is safe because smart invalidation clears the cache automatically when T_e, dx, dy, boundary conditions, or physical parameters are reassigned, and array inputs are read-only (in-place edits raise). "no_check" is a deprecated alias for True. Ignored when method != 'fd'.

Type:: bool

quiet¶

Suppress timing output. Default False.

Type:: bool

verbose¶

Print progress messages. Default True.

Type:: bool

Examples

Minimal finite-difference run:

import numpy as np
from gflex import F2D

flex = F2D()
flex.quiet = True
flex.method = 'fd'
flex.solver = 'direct'
flex.g = 9.8
flex.E = 65e9
flex.nu = 0.25
flex.rho_m = 3300.
flex.rho_fill = 0.
flex.T_e = 30e3 * np.ones((50, 50))
flex.qs = np.zeros((50, 50))
flex.qs[20:30, 20:30] = 1e6   # 10 × 10 cell load
flex.dx = flex.dy = 5000.     # 5 km grid
flex.bc_west = flex.bc_east = flex.bc_south = flex.bc_north = 'zero_moment_zero_shear'
flex.initialize()
flex.run()
flex.finalize()
deflection = flex.w           # (50, 50) array, negative downward

initialize(filename=None)[source]¶

Validate inputs and prepare the solver.

Must be called once before run(). If a configuration-file path was passed to the constructor (or to this method), parameters are read from that file; otherwise they are taken from the instance attributes set by the caller.

Parameters:: filename (str, optional) – Path to a gFlex YAML configuration file. Overrides any filename supplied to the constructor.

run()[source]¶

Execute the flexural solution.

Selects and runs the method specified by self.method. The deflection array is stored in self.w on return. Call finalize() afterwards to restore any internally modified state.

For repeated solves (e.g. a coupling loop), set cache_factorization = True before initialize() to reuse the LU factorisation when only qs changes.

finalize()[source]¶

Release all model state.

Calls the base finalize, which deletes w, qs, and the cached coefficient matrix. Read w before calling this method. self.T_e is never modified by a solve, so nothing to restore.

class gflex.F1D(filename=None)[source]¶

Bases: Flexure

One-dimensional lithospheric flexure solver.

Computes the deflection w(x) of a thin elastic beam overlying an inviscid fluid (mantle) given a surface load stress qs. Supports spatially variable elastic thickness Te.

Set instance attributes, then call initialize(), run(), and finalize() in sequence. Read flex.w before calling finalize(); finalize clears all model state including w.

method¶

Solution method. 'fd' (finite difference, supports variable Te), 'fft' (spectral, requires scalar Te), 'sas' (superposition of analytical solutions, constant Te only), or 'sas_ng' (SAS on an ungridded point array).

Type:: str

solver¶

Linear solver: 'direct' (sparse LU, default).

Type:: str

g¶

Gravitational acceleration [m s⁻²].

Type:: float

E¶

Young’s modulus [Pa].

Type:: float

nu¶

Poisson’s ratio.

Type:: float

rho_m¶

Mantle density [kg m⁻³].

Type:: float

rho_fill¶

Infill material density [kg m⁻³] (0 for air, ~1000 for water, ~2700 for rock).

Type:: float

T_e¶

Elastic thickness [m]. A scalar is broadcast to the full grid.

Type:: float or ndarray of shape (N,)

qs¶

Surface load stress [Pa].

Type:: ndarray of shape (N,)

dx¶

Grid spacing [m].

Type:: float

bc_west, bc_east

Boundary conditions on the west (left) and east (right) ends. FD options: 'zero_displacement_zero_slope' (alias 'clamped'), 'zero_displacement_zero_moment' (alias 'pinned'), 'zero_moment_zero_shear' (alias 'free'), 'zero_slope_zero_shear' (alias 'mirror'), 'periodic', 'no_outside_loads' (auto-pad by one flexural wavelength and apply 'zero_displacement_zero_slope' at the new outer edge; self.w is trimmed to the original domain), 'sandbox'. SAS option: 'no_outside_loads' (the default when unset). FFT: set both to 'periodic' for exact periodic behavior; any other value (including unset) uses zero-padding to approximate 'no_outside_loads'. Setting only one to 'periodic' raises a UserWarning and falls back to zero-padding.

Type:: str

sigma_xx¶

Normal stress applied at the plate ends [Pa]. FD only.

Type:: float, optional

fft_pad_n_alpha¶

Number of 1-D flexural-parameter units (α₁D = (4D/Δρg)^0.25) to zero-pad on each side for non-periodic FFT runs. Periodic images of the load are separated by 2 × fft_pad_n_alpha × α₁D. Default 4 (8α₁D total separation). Ignored when method != 'fft' or when all BCs are 'periodic'.

Type:: int or float

cache_factorization¶

Controls LU factorisation caching for the FD 'direct' solver. False (default) — re-factorises on every run() call. True — caches the LU factorisation and reuses it (the coefficient matrix is freed once the factors are built). Reuse is safe because smart invalidation clears the cache automatically when T_e, dx, boundary conditions, or physical parameters are reassigned, and array inputs are read-only (in-place edits raise). "no_check" is a deprecated alias for True. Ignored when method != 'fd'.

Type:: bool

quiet¶

Suppress timing output. Default False.

Type:: bool

verbose¶

Print progress messages. Default True.

Type:: bool

Examples

Minimal finite-difference run:

import numpy as np
from gflex import F1D

flex = F1D()
flex.quiet = True
flex.method = 'fd'
flex.solver = 'direct'
flex.g = 9.8
flex.E = 65e9
flex.nu = 0.25
flex.rho_m = 3300.
flex.rho_fill = 1000.
flex.T_e = 30e3
flex.qs = np.zeros(300)
flex.qs[100:200] = 1e6      # 100-cell load
flex.dx = 4000.             # 4 km grid
flex.bc_west = 'zero_displacement_zero_slope'
flex.bc_east = 'zero_moment_zero_shear'
flex.initialize()
flex.run()
flex.finalize()
deflection = flex.w         # (300,) array, negative downward

initialize(filename=None)[source]¶

Validate inputs and prepare the solver.

Parameters:: filename (str, optional) – Path to a gFlex YAML configuration file. Overrides any filename supplied to the constructor.

run()[source]¶

Execute the flexural solution.

Selects and runs the method specified by self.method. The deflection array is stored in self.w on return. Call finalize() afterwards to restore any internally modified state.

For repeated solves (e.g. a coupling loop), set cache_factorization = True before initialize() to reuse the LU factorisation when only qs changes.

finalize()[source]¶

Release all model state.

Calls the base finalize, which deletes w, qs, and the cached coefficient matrix. Read w before calling this method. self.T_e is never modified by a solve, so nothing to restore.

Boundary-condition properties¶

Added in version 2.0.0.

The boundary conditions are set as properties on the solver instance before calling initialize(). Each accepts a canonical BC string, a short alias, or a dict for inhomogeneous (prescribed-value) conditions:

flex.bc_west  = "zero_displacement_zero_slope"   # or "clamped"
flex.bc_east  = "zero_moment_zero_shear"          # or "free"
flex.bc_south = "zero_slope_zero_shear"           # or "mirror"
flex.bc_north = {"displacement": w_arr, "slope": dw_arr}

All valid strings are listed in gflex.VALID_BC_STRINGS_1D and gflex.VALID_BC_STRINGS_2D.

gflex.VALID_BC_STRINGS_1D¶: Added in version 2.0.0.

frozenset of every accepted BC string for F1D (canonical names and aliases). Use this to validate user input without maintaining a parallel copy that may drift with new releases.

gflex.VALID_BC_STRINGS_2D¶: Added in version 2.0.0.

frozenset of every accepted BC string for F2D.

Output¶

Flexure.output()[source]¶

Save deflection to file and/or plot, based on optional attributes.

Does nothing if neither w_out_file nor plot_choice has been set. Set w_out_file to a path ending in '.npy' for a binary NumPy array, or any other extension for an ASCII grid. Set plot_choice to 'q', 'w', 'both', or (1D only) 'combo' to display plots.

In-plane stresses¶

Added in version 1.4.0.

In-plane stresses are set as attributes directly on the solver instance before calling initialize(). They are not available as configuration file keys.

Attribute	Solvers	Description
`sigma_xx`	FD, FFT (1-D and 2-D)	Normal stress in the x-direction \(\sigma_{xx}\) [Pa]. Default `0`.
`sigma_yy`	FD, FFT (2-D only)	Normal stress in the y-direction \(\sigma_{yy}\) [Pa]. Default `0`.
`sigma_xy`	FD, FFT (2-D only)	Shear stress \(\sigma_{xy}\) [Pa]. Default `0`.

All three default to zero if not assigned; setting any of them with SAS or SAS_NG raises a RuntimeWarning and has no effect. See Theory and Numerics for the governing equations that include these terms.

Domain-padding utilities¶

These functions help when running F1D or F2D with a spatially variable elastic thickness grid. A smooth padding zone reduces spurious deflections at the domain boundary caused by sharp rigidity gradients, and ensures that the flexural forebulge can develop freely before reaching the boundary.

All-in-one helper

pad_domain() handles both 1-D and 2-D grids and both scalar and array elastic thickness. It calls the lower-level helpers below and is the recommended starting point.

Added in version 1.4.0: 1-D support (dispatching on qs.ndim) added in 2.0.0.

gflex.pad_domain(Te, qs, dx, dy=None, n_wavelengths=1.0, Te_out=None, E=65000000000.0, nu=0.25, rho_m=3300.0, rho_fill=0.0, g=9.8)[source]¶

Pad a flexure domain for use with F1D or F2D.

Dispatches to a 1-D or 2-D implementation based on the shape of qs. For array Te, the padding region is tapered from the inner-domain edge values toward Te_out to avoid an abrupt rigidity step. For scalar Te, only the load array is zero-padded; Te is returned unchanged as a float.

The returned pad width p can be used to trim the deflection output after the run:

# 2-D
w_inner = flex.w[p:-p, p:-p]
# 1-D
w_inner = flex.w[p:-p]

Parameters:

Te (scalar or array) – Elastic thickness [m]. A scalar is broadcast to the full padded grid internally by F1D / F2D; no tapering is applied. A 1-D array is expected when qs is 1-D; a 2-D array when qs is 2-D.
qs (1-D or 2-D array) – Surface load [Pa] for the inner domain. Shape determines whether 1-D or 2-D padding is applied.
dx (float) – Grid cell size [m]. For 2-D grids, this is the x-direction spacing.
dy (float, optional) – Grid cell size in the y-direction [m]. 2-D only; ignored for 1-D. Defaults to dx.
n_wavelengths (float, optional) – Padding width expressed as a number of flexural wavelengths. Default 1.0; use 0.5 for a less conservative (narrower) padding.
Te_out (float, optional) – Te value at the outer edge of the padding region. Only used for array Te; defaults to Te.mean().
E (float, optional) – Young’s modulus [Pa]. Default 65 GPa.
nu (float, optional) – Poisson’s ratio. Default 0.25.
rho_m (float, optional) – Mantle density [kg m^-3]. Default 3300.
rho_fill (float, optional) – Infill density [kg m^-3]. Default 0 (air).
g (float, optional) – Gravitational acceleration [m s^-2]. Default 9.8.

Returns:

Te_padded (float or array) – Elastic thickness for the padded domain. Float when Te is scalar; array of shape (len(qs) + 2p,) (1-D) or (M+2p, N+2p) (2-D) when Te is an array.
qs_padded (array) – Surface load zero-padded to the padded domain shape.
p (int) – Pad width in grid cells (same on both ends / all four sides).

Examples

2-D scalar Te:

>>> import numpy as np
>>> from gflex import pad_domain
>>> qs = np.zeros((5, 5))
>>> Te_pad, qs_pad, p = pad_domain(10e3, qs, dx=10000.,
...     E=65e9, nu=0.25, rho_m=3300., rho_fill=0., g=9.8)
>>> p
13
>>> qs_pad.shape
(31, 31)
>>> Te_pad   # scalar returned unchanged
10000.0

2-D array Te:

>>> Te = 10e3 * np.ones((5, 5))
>>> Te_pad, qs_pad, p = pad_domain(Te, qs, dx=10000.,
...     E=65e9, nu=0.25, rho_m=3300., rho_fill=0., g=9.8)
>>> Te_pad.shape
(31, 31)

1-D:

>>> qs_1d = np.zeros(5)
>>> Te_pad, qs_pad, p = pad_domain(10e3, qs_1d, dx=10000.)
>>> qs_pad.shape
(43,)

Lower-level building blocks

These two pairs of functions are the building blocks used by pad_domain(). Use them directly when you need finer control — for example, to compute the pad width once and apply it to multiple arrays, or to inspect the tapered Te grid before running the solver.

Recommended pad width (number of cells):

Added in version 1.4.0.

gflex.recommended_pad_width(Te, dx, E=65000000000.0, nu=0.25, rho_m=3300.0, rho_fill=0.0, g=9.8, n_wavelengths=1.0)[source]¶

Return the recommended padding width in grid cells for a variable-Te run.

The padded domain boundary should be at least one flexural wavelength from the load so that the plate’s response to the load is negligible at the boundary. Half a wavelength is often sufficient in practice.

The 2-D flexural wavelength is computed from the maximum Te value via flexural_wavelengths(), giving the most conservative (widest) padding estimate.

Parameters:

Te (scalar or 2-D array) – Elastic thickness [m]. The maximum value is used.
dx (float) – Grid cell size [m]. Use the smaller of dx and dy if they differ.
E (float, optional) – Young’s modulus [Pa]. Default 65 GPa.
nu (float, optional) – Poisson’s ratio. Default 0.25.
rho_m (float, optional) – Mantle density [kg m^-3]. Default 3300.
rho_fill (float, optional) – Infill density [kg m^-3]. Default 0 (air).
g (float, optional) – Gravitational acceleration [m s^-2]. Default 9.8.
n_wavelengths (float, optional) – Number of flexural wavelengths to use as the padding width. Default 1.0. Use 0.5 for a less conservative estimate.

Returns:

Recommended padding width in grid cells.

Return type:

int

gflex.recommended_pad_width_1d(Te, dx, E=65000000000.0, nu=0.25, rho_m=3300.0, rho_fill=0.0, g=9.8, n_wavelengths=1.0)[source]¶

Return the recommended padding width in grid cells for a 1-D variable-Te run.

The padded domain boundary should be at least one 1-D flexural wavelength from the load so that the plate’s response is negligible at the boundary.

The 1-D flexural wavelength is computed from the maximum Te value, giving the most conservative (widest) padding estimate.

Parameters:

Te (scalar or 1-D array) – Elastic thickness [m]. The maximum value is used.
dx (float) – Grid cell size [m].
E (float, optional) – Young’s modulus [Pa]. Default 65 GPa.
nu (float, optional) – Poisson’s ratio. Default 0.25.
rho_m (float, optional) – Mantle density [kg m^-3]. Default 3300.
rho_fill (float, optional) – Infill density [kg m^-3]. Default 0 (air).
g (float, optional) – Gravitational acceleration [m s^-2]. Default 9.8.
n_wavelengths (float, optional) – Number of flexural wavelengths to use as the padding width. Default 1.0. Use 0.5 for a less conservative estimate.

Returns:

Recommended padding width in grid cells.

Return type:

int

Examples

>>> recommended_pad_width_1d(Te=35e3, dx=5000.)
94

Smooth Te taper (extends a Te array into the padding zone):

Added in version 1.4.0.

gflex.smooth_pad_Te(Te, pad_width, Te_out=None)[source]¶

Pad a 2-D elastic thickness array with a smooth linear taper.

When a spatially variable Te grid is padded with a constant value before being passed to F2D, the abrupt step in flexural rigidity D at the inner/outer boundary drives spurious deflections via the D-derivative terms in the vWC1994 stencil (issue #45).

This function eliminates that step by linearly blending the inner-domain edge values toward Te_out across the padding ring, reducing the rigidity gradient at the inner/outer boundary by a factor of ~pad_width compared with an abrupt step.

The corresponding surface load array should be padded with zeros, e.g.:

qs_padded = numpy.pad(qs, pad_width, mode='constant')

Parameters:

Te ((M, N) array) – Elastic thickness [m] for the inner domain.
pad_width (int) – Width of the padding ring in grid cells.
Te_out (float, optional) – Te value at the outer edge of the padding ring. Defaults to Te.mean().

Returns:

Te_padded – Padded elastic thickness with a smooth linear taper.

Return type:

(M + 2*pad_width, N + 2*pad_width) array

gflex.smooth_pad_Te_1d(Te, pad_width, Te_out=None)[source]¶

Pad a 1-D elastic thickness array with a smooth linear taper.

When a spatially variable Te array is padded with a constant value before being passed to F1D, the abrupt step in flexural rigidity D at the inner/outer boundary drives spurious deflections via the D-derivative terms in the FD stencil.

This function eliminates that step by linearly blending the inner-domain edge values toward Te_out across the padding region.

The corresponding surface load array should be padded with zeros, e.g.:

qs_padded = numpy.pad(qs, pad_width, mode='constant')

Parameters:

Te (1-D array) – Elastic thickness [m] for the inner domain.
pad_width (int) – Width of the padding on each end in grid cells.
Te_out (float, optional) – Te value at the outer edge of the padding. Defaults to Te.mean().

Returns:

Te_padded – Padded elastic thickness with a smooth linear taper.

Return type:

1-D array of length len(Te) + 2 * pad_width

FD boundary-condition warnings¶

When running F1D or F2D with the finite-difference solver, gFlex issues UserWarning messages for two categories of potentially problematic boundary conditions.

BC-type warnings fire whenever a side carries a BC whose physical interpretation deserves verification:

'zero_moment_zero_shear' (alias 'free') — assumes a free broken plate end (zero moment and shear force). Physically appropriate for rifted or passive continental margins, subduction trenches with an applied edge load (slab pull), and broken-plate flexure (Turcotte & Schubert). Often applied uncritically elsewhere in the literature — verify that one of these settings applies.

Proximity warnings fire for 'zero_displacement_zero_slope' (alias 'clamped') boundaries when the nearest loaded cell is within one flexural wavelength (\(\lambda = 2\pi\alpha\), where \(\alpha = (4D / \Delta\rho g)^{1/4}\)) of the boundary. Within this distance the flexural forebulge — which peaks at \(\approx \pi\alpha\) from the load — will be suppressed by the zero-displacement condition, contaminating the solution. The warning message reports the distance as a fraction of the local flexural wavelength and directs you to the domain-padding utilities.

Warning deduplication in model-coupling loops

Python’s default warning filter shows each unique warning once per call site per interpreter session. In a time-stepping or iterative-coupling loop such as:

for load in loads:
    flex.qs = load
    flex.run()          # warning fires on first iteration, silenced thereafter

the proximity warning fires on the first iteration and is not repeated, even if the load subsequently moves closer to the boundary. To re-enable the warning on every call:

import warnings
warnings.filterwarnings("always", category=UserWarning, module="gflex")

Suppressing warnings you have verified

Once you have confirmed that a boundary condition is appropriate for your setup, suppress the corresponding warning by message text:

import warnings
warnings.filterwarnings("ignore", message=".*zero_moment_zero_shear.*")

To suppress all gFlex warnings:

warnings.filterwarnings("ignore", module="gflex")

Or use a context manager for a single run:

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    flex.run()

LU factorization cache¶

Added in version 2.0.0.

For coupling workflows that call run() (or run()) repeatedly with the same grid, elastic thickness, and boundary conditions, the sparse-LU factorization of the coefficient matrix can be cached to avoid re-factorizing on every call. Set the attribute before calling initialize():

Value	Behaviour
`False` (default)	No caching. The matrix is factorized on every `run()` call.
`True`	Cache the LU factorization and reuse it on every call. The coefficient matrix is freed from memory immediately after factorization; only the LU factors are retained. Reuse is safe because smart invalidation (see below) clears the cache when a matrix-determining input is reassigned, and array inputs are read-only, so the cached factorization can never silently desynchronise.
`"no_check"`	Deprecated alias for `True` (emits a `DeprecationWarning`). The two were distinct before the cache redesign, when `True` validated a per-call matrix hash; that check is no longer needed.

Smart invalidation

Reassigning any matrix-determining attribute — T_e, E, nu, g, rho_m, rho_fill, dx, dy, boundary conditions, or in-plane stresses — automatically clears the cached coefficient matrix and LU factorization. No explicit cache management is needed between solves when only qs changes.

Note

Smart invalidation is triggered by assignment (flex.T_e = new_array). Because the cache is keyed to these values, the matrix-determining array inputs (T_e, sigma_xx, sigma_yy, sigma_xy) are stored read-only, so an in-place edit raises ValueError rather than silently leaving the cache stale. To change a value, reassign the whole array:

te = flex.T_e.copy()
te[5] = 40e3
flex.T_e = te        # triggers invalidation

Example (coupling loop):

flex.cache_factorization = True
flex.initialize()

for load in load_sequence:
    flex.qs = load
    flex.run()
    w = flex.w
    # ... process w ...

flex.finalize()

The cache is cleared by finalize().

Unpaired periodic boundaries¶

Added in version 2.0.0.

For the finite-difference solver, a one-sided periodic boundary — 'periodic' on only one side of an opposite pair — is not well-posed, and run() raises a ValueError. Set allow_unpaired_periodic = True to override the guard and solve anyway; enabling it emits a one-time UserWarning that the safety check is disabled. Defaults to False.

Timing attributes¶

Added in version 2.0.0.

After each run() call, the following read-only attributes report wall-clock times measured with time.perf_counter():

Attribute	Solvers	Description
`time_to_solve`	All	Total solve wall time [s] from the start of `run()` to the end.
`coeff_creation_time`	FD only	Time [s] to construct the sparse coefficient matrix. Not set for SAS or FFT solvers. In a coupling loop with `cache_factorization` enabled, this is zero on cache-hit calls.
`linear_solve_time`	FD only	Time [s] for the LU backsolve (triangular solve). Not set for SAS or FFT solvers.

coeff_creation_time and linear_solve_time together account for most of time_to_solve; the remainder is boundary-condition setup and array housekeeping. Both are most useful in coupling loops: coeff_creation_time confirms the coefficient matrix was reused (value ≈ 0) and linear_solve_time shows the marginal cost per run() call.

Flexural wavelengths¶

Added in version 1.4.0.

gflex.flexural_wavelengths(Te, E, nu, rho_m, rho_fill, g)[source]¶

Compute flexural parameters and wavelengths for a thin elastic plate.

Parameters:

Te (float or ndarray of float) – Elastic thickness [m]. Arrays are supported and return arrays of the same shape for every output key (used, e.g., by the per-cell load-proximity guard).
E (float) – Young’s modulus [Pa].
nu (float) – Poisson’s ratio.
rho_m (float) – Mantle density [kg m^-3].
rho_fill (float) – Infill density [kg m^-3] (e.g. 0 for air, 1000 for water).
g (float) – Gravitational acceleration [m s^-2].

Returns:

Keys alpha_1D, lambda_1D, zero_crossing_1D, alpha_2D, lambda_2D, zero_crossing_2D — all in metres.

Return type:

dict

Examples

>>> from gflex import flexural_wavelengths
>>> r = flexural_wavelengths(Te=30e3, E=65e9, nu=0.25,
...                          rho_m=3300., rho_fill=0., g=9.8)
>>> round(r["alpha_2D"] / 1e3, 1)  # km
46.9
>>> round(r["lambda_1D"] / r["lambda_2D"], 6)
1.414214

Coupling guide¶

gFlex is material-agnostic: it receives a surface-normal stress [Pa] and returns a deflection [m], regardless of whether the load comes from ice, water, sediment, volcanic edifices, or any combination. The caller is responsible for converting source-specific quantities into Pa before passing them to gFlex.

Load conversion

# Glacial isostasy
qs = rho_ice * g * ice_thickness        # kg m⁻³ × m s⁻² × m → Pa

# Sediment or volcanic load
qs = rho_sediment * g * sediment_thickness

# Multiple sources: sum them
qs = rho_ice * g * h_ice + rho_sed * g * h_sed + rho_water * g * h_water

The rho_fill parameter should be set to the density of the material that replaces the load inside the flexural depression — rho_fill=0 for subaerially exposed basins, rho_fill=1030 for submarine basins, rho_fill=2000–2700 for sediment-filled basins.

Applying deflection to topography

gFlex returns the total instantaneous deflection, not an increment. In a time-stepping loop, apply only the change to topography:

flex.qs = qs
flex.run()
w_new = flex.w
topo += w_new - w_prev
w_prev = w_new.copy()

QGIS Processing provider

For GIS-based workflows where load and elastic-thickness data are already rasters, processing_gflex exposes gFlex as a no-code algorithm in the QGIS Processing Toolbox and Graphical Modeler. It supports all 2-D solution methods (FD, FFT, SAS), variable or scalar \(T_e\), all boundary conditions, and in-plane stresses. Installation:

pip install "gflex>=2.0.0"   # auto-installed by the plugin on first use

Requires QGIS ≥ 3.16. Usable from the Toolbox, Graphical Modeler, or headlessly via processing.run().

BMI interface¶

BmiGflex exposes the CSDMS Basic Model Interface, enabling gFlex to be coupled with other models in the CSDMS framework. It requires the optional bmipy dependency (pip install gflex[bmi]).

BMI variables

Grid 0 — spatial flexure grid:

Name	Direction	Units	Description
`load__normal_component_of_stress`	input	Pa	Surface-normal load stress \(q_s\). Material-agnostic: convert ice, water, sediment, etc. to Pa before calling `set_value()`.
`lithosphere__elastic_thickness`	input	m	Elastic thickness \(T_e\). Usually set once at initialisation; updating it between `update()` calls invalidates the cached LU factorisation so the next solve rebuilds the stiffness matrix.
`lithosphere__vertical_displacement`	output	m	Deflection \(w\) (downward negative).

Grid 1 — scalar physical constants, exposed as the BMI 'scalar' grid type (rank 0) and read or written as single-element arrays:

Name	Direction	Units	Description
`lithosphere__young_modulus`	input	Pa	Young’s modulus \(E\).
`lithosphere__poisson_ratio`	input	1	Poisson’s ratio \(\nu\) (dimensionless).
`mantle__mass-per-volume_density`	input	kg m⁻³	Mantle density \(\rho_m\).
`infill_material__mass-per-volume_density`	input	kg m⁻³	Infill material density \(\rho_\text{fill}\). The only constant with a runtime-update use case: a basin transitioning from subaerial (\(\rho_\text{fill}=0\)) to subaqueous (\(\rho_\text{fill}=1030\)) during a simulation.
`planet_surface__gravitational_acceleration`	input	m s⁻²	Gravitational acceleration \(g\).

Updating any scalar constant via set_value() propagates to the solver immediately and invalidates the cached LU factorisation, so the next update() rebuilds the stiffness matrix automatically.

Coupling example

from gflex import BmiGflex
import numpy as np

bmi = BmiGflex()
bmi.initialize("my_config.yaml")

n = bmi.get_grid_size(0)
load = np.zeros(n)

for step in range(n_steps):
    load[:] = rho_ice * g * ice_thickness.ravel()
    bmi.set_value("load__normal_component_of_stress", load)
    bmi.update()
    w = np.empty(n)
    bmi.get_value("lithosphere__vertical_displacement", w)
    # … apply w to topography …

bmi.finalize()

class gflex.BmiGflex[source]¶

Bases: object

BMI wrapper for gFlex lithospheric flexure.

Implements the CSDMS Basic Model Interface v2 specification. Supports 1-D and 2-D gridded flexure solutions (FD, FFT, SAS methods). The SAS_NG point-load method is not suited to the BMI grid model.

Grids¶

Grid 0 — the spatial flexure grid (uniform rectilinear).: Shape (nrows,) in 1-D or (nrows, ncols) in 2-D, with spacing (dy, dx) and origin at (0, 0).
Grid 1 — scalar parameter grid (uniform rectilinear, shape (1,)).: Holds the five physical constants below. These are spatially uniform by assumption; exposing them as single-element arrays supports introspection and ensemble initialisation via the BMI.

Time¶

gFlex solves instantaneous elastic equilibrium. Time is therefore nominal: start=0, step=1, end=inf. Each call to update() applies the current load and computes deflection.

Variables — grid 0¶

Input: load__normal_component_of_stress [Pa]: Surface-normal load stress q_s = ρ g h. Material-agnostic.
Input: lithosphere__elastic_thickness [m]: Elastic thickness T_e. Updating it invalidates the cached LU factorisation; the next update() rebuilds the stiffness matrix.
Output: lithosphere__vertical_displacement [m]: Lithospheric deflection w (downward negative).

Variables — grid 1 (scalar constants)¶

Input: lithosphere__young_modulus [Pa] Input: lithosphere__poisson_ratio [1] Input: mantle__mass-per-volume_density [kg m-3] Input: infill_material__mass-per-volume_density [kg m-3] Input: planet_surface__gravitational_acceleration [m s-2]

Changes to any scalar constant via set_value() are pushed to the solver immediately and invalidate the cached LU factorisation, so the next update() rebuilds the stiffness matrix automatically.

initialize(config_file: str) → None[source]¶

Initialize gFlex from a configuration file.

Parameters:: config_file (str) – Path to a gFlex YAML configuration file.

update() → None[source]¶

Compute flexural deflection for the current load.

Writes the current load__normal_component_of_stress array into the model’s internal qs field, runs the solver, then copies the result into lithosphere__vertical_displacement.

finalize() → None[source]¶: Tear down the model and release resources.

get_value(name: str, dest: NDArray[Any]) → NDArray[Any][source]¶: Copy the flattened values of variable name into dest and return it.

set_value(name: str, src: NDArray[Any]) → None[source]¶

Overwrite the entire array for variable name with values from src.

For lithosphere__elastic_thickness and the five scalar physical constants, the new values are pushed to the solver immediately, invalidating the cached coefficient matrix so that the next update() uses the updated parameters.

Landlab component¶

The Landlab component (landlab.components.gFlex) exposes gFlex within the Landlab Earth-surface modelling framework. It uses the same underlying solver and the same CSDMS Standard Name fields as the BMI, but follows Landlab conventions: construction via __init__ and time-stepping via run_one_step().

Comparison with the BMI

	gFlex BMI	Landlab component
Lifecycle	`initialize` / `update` / `finalize`	`__init__` / `run_one_step`
Load input	`set_value("load__normal_component_of_stress", q)`	`grid.at_node["load__normal_component_of_stress"][:] = q`
\(T_e\) update	`set_value("lithosphere__elastic_thickness", te)`	`grid.at_node["lithosphere__elastic_thickness"][:] = te`
Deflection output	`get_value("lithosphere__vertical_displacement", w)`	`grid.at_node["lithosphere__vertical_displacement"]`

Field names

CSDMS Standard Name	Direction	Units	Notes
`load__normal_component_of_stress`	input	Pa	Required. Maps to gFlex internal \(q_s\).
`lithosphere__elastic_thickness`	input	m	Optional. Re-read on every `run_one_step()` call if present, enabling runtime \(T_e\) updates without re-initialisation.
`lithosphere__vertical_displacement`	output	m	Total deflection \(w\) (downward negative).

rho_fill defaults to 0.0 (air; no infill) in the Landlab component. Set it explicitly at construction for marine (rho_fill=1030) or sediment-filled (rho_fill=2000–2700) basins.

Coupling example

import numpy as np
from landlab import RasterModelGrid
from landlab.components import gFlex

mg = RasterModelGrid((100, 100), xy_spacing=5000.0)
mg.add_zeros("load__normal_component_of_stress", at="node")
mg.add_zeros("topographic__elevation", at="node")

gf = gFlex(mg, Youngs_modulus=65e9, Poissons_ratio=0.25,
           rho_mantle=3300, rho_fill=0, elastic_thickness=35e3)

w_prev = np.zeros(mg.number_of_nodes)

for step in range(n_steps):
    mg.at_node["load__normal_component_of_stress"][:] = (
        rho_ice * g * ice_thickness.ravel()
    )
    gf.run_one_step()
    w = mg.at_node["lithosphere__vertical_displacement"]
    mg.at_node["topographic__elevation"] += w - w_prev
    w_prev = w.copy()

Installation

pip install landlab

Requires a Landlab release that includes the gFlex v2 component (landlab/landlab#2420).