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:
FlexureTwo-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(), andfinalize()in sequence. Readflex.wbefore callingfinalize(); finalize clears all model state includingw.- 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:
- rho_fill¶
Infill material density [kg m⁻³] (0 for air, ~1000 for water, ~2700 for rock).
- Type:
- 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)
- 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','zero_moment_zero_shear'(alias'free'),'zero_slope_zero_shear'(alias'mirror'),'periodic'. SAS option:'no_outside_loads'(the default when unset).- Type:
- sigma_xx¶
Normal in-plane stress in the x-direction \(\sigma_{xx}\) [Pa]. Supported by
fdandfft. Default0.- Type:
float, optional
- sigma_yy¶
Normal in-plane stress in the y-direction \(\sigma_{yy}\) [Pa]. Supported by
fdandfft. Default0.- Type:
float, optional
- sigma_xy¶
In-plane shear stress \(\sigma_{xy}\) [Pa]. Supported by
fdandfft. Default0.- Type:
float, optional
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 inself.won return. Callfinalize()afterwards to restore any internally modified state.For repeated solves (e.g. a coupling loop), set
cache_factorization = Truebeforeinitialize()to reuse the LU factorisation when onlyqschanges; use"no_check"for maximum throughput when the coefficient matrix is guaranteed stable.
- class gflex.F1D(filename=None)[source]¶
Bases:
FlexureOne-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(), andfinalize()in sequence. Readflex.wbefore callingfinalize(); finalize clears all model state includingw.- 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:
- rho_fill¶
Infill material density [kg m⁻³] (0 for air, ~1000 for water, ~2700 for rock).
- Type:
- 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,)
- 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','zero_moment_zero_shear'(alias'free'),'zero_slope_zero_shear'(alias'mirror'),'periodic','sandbox'. SAS option:'no_outside_loads'(the default when unset).- Type:
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.
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 inself.won return. Callfinalize()afterwards to restore any internally modified state.For repeated solves (e.g. a coupling loop), set
cache_factorization = Truebeforeinitialize()to reuse the LU factorisation when onlyqschanges; use"no_check"for maximum throughput when the coefficient matrix is guaranteed stable.
Output¶
- Flexure.output()[source]¶
Save deflection to file and/or plot, based on optional attributes.
Does nothing if neither
w_out_filenorplot_choicehas been set. Setw_out_fileto a path ending in'.npy'for a binary NumPy array, or any other extension for an ASCII grid. Setplot_choiceto'q','w','both', or (1D only)'combo'to display plots.
In-plane stresses¶
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 |
|---|---|---|
|
FD, FFT (1-D and 2-D) |
Normal stress in the x-direction \(\sigma_{xx}\) [Pa].
Default |
|
FD, FFT (2-D only) |
Normal stress in the y-direction \(\sigma_{yy}\) [Pa].
Default |
|
FD, FFT (2-D only) |
Shear stress \(\sigma_{xy}\) [Pa]. Default |
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 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.
2-D (F2D)
- 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 both the elastic thickness and surface load arrays for use with F2D.
Combines
recommended_pad_width()andsmooth_pad_Te()into a single call, and zero-pads qs to match. The returned pad width p can be used to trim the deflection output after the run:w_inner = flex.w[p:-p, p:-p]
- Parameters:
Te ((M, N) array) – Elastic thickness [m] for the inner domain.
qs ((M, N) array) – Surface load [Pa] for the inner domain.
dx (float) – Grid cell size in the x-direction [m].
dy (float, optional) – Grid cell size in the y-direction [m]. 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 ring. 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 ((M + 2p, N + 2p) array) – Smoothly tapered elastic thickness.
qs_padded ((M + 2p, N + 2p) array) – Surface load zero-padded to match.
p (int) – Pad width in grid cells (same on all four sides).
Examples
>>> import numpy as np >>> from gflex import pad_domain >>> Te = 10e3 * np.ones((5, 5)) >>> qs = np.zeros((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) >>> p 13 >>> Te_pad.shape (31, 31)
To run
F2Dwith the padded grids:flex = gflex.F2D() flex.T_e = Te_pad flex.qs = qs_pad # ... set other parameters and run ... w_inner = flex.w[p:-p, p:-p] # trim padding from output
- 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. Use
recommended_pad_width()to obtain a suitable value.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
Examples
>>> import numpy as np >>> from gflex import smooth_pad_Te >>> Te = 35e3 * np.ones((10, 10)) >>> Te_pad = smooth_pad_Te(Te, pad_width=4) >>> Te_pad.shape (18, 18)
To run
F2Dwith the padded grid:import numpy as np from gflex import F2D, smooth_pad_Te, recommended_pad_width p = recommended_pad_width(Te, dx=5000.) Te_pad = smooth_pad_Te(Te, p) qs_pad = np.pad(qs, p, mode='constant') flex = F2D() flex.T_e = Te_pad flex.qs = qs_pad # ... set other parameters and run ... w_inner = flex.w[p:-p, p:-p] # trim padding from output
- 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:
1-D (F1D)
- gflex.pad_domain_1d(Te, qs, dx, 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 both the elastic thickness and surface load arrays for use with F1D.
Combines
recommended_pad_width_1d()andsmooth_pad_Te_1d()into a single call, and zero-pads qs to match. The returned pad width p can be used to trim the deflection output after the run:w_inner = flex.w[p:-p]
- Parameters:
Te (1-D array) – Elastic thickness [m] for the inner domain.
qs (1-D array) – Surface load [Pa] for the inner domain.
dx (float) – Grid cell size [m].
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. 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 (1-D array of length
len(Te) + 2p) – Smoothly tapered elastic thickness.qs_padded (1-D array of length
len(qs) + 2p) – Surface load zero-padded to match.p (int) – Pad width in grid cells (same on both ends).
Examples
>>> import numpy as np >>> Te = 10e3 * np.ones(5) >>> qs = np.zeros(5) >>> Te_pad, qs_pad, p = pad_domain_1d(Te, qs, dx=10000., ... E=65e9, nu=0.25, rho_m=3300., rho_fill=0., g=9.8) >>> p 19 >>> Te_pad.shape (43,)
To run
F1Dwith the padded arrays:flex = gflex.F1D() flex.T_e = Te_pad flex.qs = qs_pad # ... set other parameters and run ... w_inner = flex.w[p:-p] # trim padding from output
- 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. Use
recommended_pad_width_1d()to obtain a suitable value.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
Examples
>>> import numpy as np >>> Te = 35e3 * np.ones(10) >>> Te_pad = smooth_pad_Te_1d(Te, pad_width=4) >>> Te_pad.shape (18,)
To run
F1Dwith the padded arrays:import numpy as np from gflex import F1D, smooth_pad_Te_1d, recommended_pad_width_1d p = recommended_pad_width_1d(Te, dx=5000.) Te_pad = smooth_pad_Te_1d(Te, p) qs_pad = np.pad(qs, p, mode='constant') flex = F1D() flex.T_e = Te_pad flex.qs = qs_pad # ... set other parameters and run ... w_inner = flex.w[p:-p] # trim padding from output
- 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:
Examples
>>> recommended_pad_width_1d(Te=35e3, dx=5000.) 94
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¶
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 |
|---|---|
|
No caching. The matrix is factorized on every |
|
Cache with hash check. The factorization is reused when a hash of the coefficient matrix matches the stored hash; it is recomputed when any of Te, dx/dy, BCs, or physical parameters change. |
|
Cache without hash check. The stored factorization is reused on every call without computing a hash. Gives the maximum performance benefit. The coefficient matrix is freed from memory immediately after factorization; only the LU factorization is retained. Smart invalidation (see below) still applies: reassigning a matrix-determining input clears the cache and triggers a rebuild on the next call. |
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),
not by in-place mutation of a NumPy array (flex.T_e[5] = 40e3). If
you mutate an array in place, reassign it afterwards to ensure the cache
is correctly invalidated:
flex.T_e[5] = 40e3
flex.T_e = flex.T_e # trigger 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().
Flexural wavelength¶
- gflex.flexural_wavelengths(Te, E, nu, rho_m, rho_fill, g)[source]¶
Compute flexural parameters and wavelengths for a thin elastic plate.
- Parameters:
- Returns:
Keys
alpha_1D,lambda_1D,zero_crossing_1D,alpha_2D,lambda_2D,zero_crossing_2D— all in metres.- Return type:
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
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]).
- class gflex.BmiGflex[source]¶
Bases:
objectBMI 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.
Grid¶
All variables share grid identifier 0, a uniform rectilinear grid. Grid shape follows NumPy / C order: (nrows,) in 1-D or (nrows, ncols) in 2-D, with spacing (dy, dx) and origin at (0, 0).
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¶
Input:
lithosphere__load_pressure[Pa] — surface normal stress q0 Output:lithosphere__vertical_displacement[m] — downward deflection w- 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
lithosphere__load_pressurearray into the model’s internalqsfield, runs the solver, then copies the result intolithosphere__vertical_displacement.