API Reference¶

Primary solvers¶

Both solvers follow the same three-step workflow:

flex.initialize()
flex.run()
flex.finalize()

After finalize(), the deflection is available as flex.w. Call output() to save results to file or display plots.

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. The deflection is available as flex.w after finalize().

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

PlateSolutionType¶

FD stencil variant: 'vWC1994' (van Wees & Cloetingh 1994, recommended) or 'G2009' (Govers et al. 2009).

Type:: str

Solver¶

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

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

Te¶

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_W, BC_E, BC_N, BC_S

Boundary conditions on the west, east, north, and south edges. FD options: '0Displacement0Slope', '0Slope0Shear', '0Moment0Shear', 'Mirror', 'Periodic'. SAS option: 'NoOutsideLoads' (the default when unset).

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

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.PlateSolutionType = 'vWC1994'
flex.Solver = 'direct'
flex.g = 9.8
flex.E = 65e9
flex.nu = 0.25
flex.rho_m = 3300.
flex.rho_fill = 0.
flex.Te = 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_W = flex.BC_E = flex.BC_S = flex.BC_N = '0Moment0Shear'
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 configuration file (INI or YAML). 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.

finalize()[source]¶

Clean up after the solver.

Restores self.Te to its pre-run value if gFlex padded it internally, so the object can be reused cleanly in a model-coupling loop. Clears the cached coefficient matrix.

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. The deflection is available as flex.w after finalize().

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) or 'iterative'.

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

Te¶

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_W, BC_E

Boundary conditions on the west (left) and east (right) ends. FD options: '0Displacement0Slope', '0Slope0Shear', '0Moment0Shear', 'Mirror', 'Periodic'. SAS option: 'NoOutsideLoads' (the default when unset).

Type:: str

sigma_xx¶

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

Type:: float, optional

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.Te = 30e3
flex.qs = np.zeros(300)
flex.qs[100:200] = 1e6      # 100-cell load
flex.dx = 4000.             # 4 km grid
flex.BC_W = '0Displacement0Slope'
flex.BC_E = '0Moment0Shear'
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 configuration file (INI or YAML). 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.

finalize()[source]¶

Clean up after the solver.

Restores self.Te to its pre-run value if gFlex padded it internally, so the object can be reused cleanly in a model-coupling loop. Clears the cached coefficient matrix.

Output¶

Flexure.output()[source]¶

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

Does nothing if neither wOutFile nor plotChoice has been set. Set wOutFile to a path ending in '.npy' for a binary NumPy array, or any other extension for an ASCII grid. Set plotChoice to '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
`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 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() and smooth_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 F2D with the padded grids:

flex = gflex.F2D()
flex.Te = 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 F2D with 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.Te = 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:

int

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() and smooth_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 F1D with the padded arrays:

flex = gflex.F1D()
flex.Te = 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 F1D with 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.Te = 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:

int

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:

'0Moment0Shear' — assumes a free broken plate end (zero moment and shear force). This is physically appropriate for a rifted margin or spreading ridge where the plate really is broken, but is often applied uncritically elsewhere in the literature.
'0Slope0Shear' — requires the plate to be simultaneously horizontal and shear-free at the boundary. No clear geological analog is known for this combination in a non-trivial deflection setting.

Proximity warnings fire for '0Displacement0Slope' 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=".*0Moment0Shear.*")

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()

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:

Te (float) – Elastic thickness [m].
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

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: 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.

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 configuration file (INI or YAML).

update() → None[source]¶

Compute flexural deflection for the current load.

Writes the current lithosphere__load_pressure 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[tuple[Any, ...], dtype[Any]]) → ndarray[tuple[Any, ...], dtype[Any]][source]¶: Copy the flattened values of variable name into dest and return it.

set_value(name: str, src: ndarray[tuple[Any, ...], dtype[Any]]) → None[source]¶: Overwrite the entire array for variable name with values from src.