gflex 2.0.0b1

One- and two-dimensional plate bending, designed for Earth's lithosphere

gFlex

Multiple methods to solve elastic plate flexure, designed for applications to Earth's lithosphere.

These instructions are meant to take a user familiar with computers but new to (or a beginner with) Python through the basics of how to get gFlex to work. The Python scripting part towards the end should be pretty straightforward as well, insofar as information is provided on how to get and set the chosen values inside gFlex. Please leave a message if you have trouble working with gFlex; your comments could assist both you and the more general improvement of this documentation.

When you use gFlex, please cite:

Wickert, A. D. (2016), Open-source modular solutions for flexural isostasy: gFlex v1.0, Geosci. Model Dev., 9(3), 997–1017, doi:10.5194/gmd-9-997-2016.

If you additionally want an up-to-date citation for the latest source-code release, please see that given in CITATION.cff.

Documentation

Full documentation, including a configuration file parameter reference, accuracy benchmarks, and the full API, is available at gflex.readthedocs.io.

Installation

pip install gflex

gFlex requires Python ≥ 3.11. Dependencies (numpy, scipy, matplotlib, pyyaml) are installed automatically.

For a development install from source:

git clone https://github.com/awickert/gFlex.git
cd gFlex
pip install -e .

Running

Once gFlex is installed, it is possible to run it in four ways:

1. With a configuration file
2. Within a Python script
3. Within GRASS GIS
4. As part of the Landlab Earth-surface modeling framework, with a CSDMS Basic Model Interface (BMI)

For options 1 and 2, there are pre-built methods that can be selected along the way to visualize results. These use Python's Matplotlib plotting library. For option 3, GRASS GIS is used for visualization. In option 4, output from Landlab can be visualized with Matplotlib.

With configuration file

A configuration file can be generated to run gFlex; see examples in the input/ directory. To run gFlex using this file, one simply opens a terminal window and types:

gflex <path-to-configuration-file>

This can be run from any directory, as the installation of gFlex adds the program "gflex" to the system path.

Configuration files must be YAML (.yaml / .yml extension). See input/input_f1d.yaml and input/input_f2d.yaml for complete examples.

A minimal 2-D YAML configuration file looks like:

mode:
  dimension: 2
  method: fd
parameter:
  youngs_modulus: 6.5e10
  poissons_ratio: 0.25
  gravitational_acceleration: 9.8
  mantle_density: 3300
  infill_material_density: 0
input:
  loads: path/to/loads.txt
  elastic_thickness: path/to/Te.txt
output:
  plot: both
numerical:
  grid_spacing_x: 4000
  boundary_condition_west: zero_moment_zero_shear
  boundary_condition_east: zero_displacement_zero_slope
numerical2D:
  grid_spacing_y: 4000
  boundary_condition_north: mirror
  boundary_condition_south: zero_slope_zero_shear

For a full parameter reference, see the Configuration Files page on ReadTheDocs.

Finite-difference boundary conditions

Six boundary conditions are available for FD solutions (see also Table 1 in Wickert, 2016):

Name	Condition	Physical interpretation
zero_displacement_zero_slope	w = 0, dw/dx = 0	Clamped end: zero deflection and zero slope (no rotation)
zero_displacement_zero_moment	w = 0, d²w/dx² = 0	Simply supported (pinned) end: zero deflection, free to rotate
zero_moment_zero_shear	d²w/dx² = d³w/dx³ = 0	Broken plate: free end with no moment or shear
zero_slope_zero_shear	dw/dx = d³w/dx³ = 0	Plate is level at the boundary but free to deflect there; no shear transmitted
mirror	w(b − x) = w(b + x)	Mirror-symmetry plane — model only half of a symmetric system
periodic	w(0) = w(L)	Wrap-around: the domain tiles infinitely in both directions

For SAS and SAS_NG, no_outside_loads (or a blank entry) is used instead; the plate is assumed undeflected at infinity.

FD boundary-condition warnings: when running F1D or F2D with the finite-difference solver, gFlex issues UserWarning messages for 'zero_moment_zero_shear' (free broken plate end — verify a rifted margin is intended), 'zero_slope_zero_shear' (no clear geological analog), and when the nearest loaded cell is within one flexural wavelength of a 'zero_displacement_zero_slope' boundary (the forebulge would be suppressed). See the API reference for how to suppress or re-enable these warnings.

A note on zero_slope_zero_shear: the label in Wickert (2016) is "free displacement of a horizontally clamped boundary." The plate is forced to be exactly level at the boundary (dw/dx = 0) — as if it were clamped against rotation — while its vertical position is unconstrained and no shear force is transmitted (d³w/dx³ = 0). This is superficially similar to mirror (both enforce zero slope at the boundary), but zero_slope_zero_shear uses a different finite-difference stencil and the two produce noticeably different solutions even far from the boundary. For symmetry problems — e.g., modelling half of a symmetric mountain range or ice sheet — mirror is the more accurate choice.

Within a Python script (with or without a configuration file)

You may run gFlex from other Python programs. When you install it (above), this also produces a Python module that you may import to access it while scripting.

With no configuration file (recommended)

input/run_in_script_2D.py, reproduced below, is a good example of how to set the variables and run the model. This method requires no input file, as all of the values are set inside the Python script that imports gflex. This is essentially how the GRASS GIS interface was written, and is a way to embed the abilities of gFlex into another model. A one-dimensional example, input/run_in_script_1D.py, is also available.

#! /usr/bin/env python

import gflex
import numpy as np

flex = gflex.F2D()

flex.quiet = False

flex.method = 'fd' # Solution method: * fd (finite difference)
                   #                  * fft (spectral)
                   #                  * sas (superposition of analytical solutions)
                   #                  * sas_ng (ungridded SAS)

flex.g = 9.8 # acceleration due to gravity
flex.E = 65E9 # Young's Modulus
flex.nu = 0.25 # Poisson's Ratio
flex.rho_m = 3300. # mantle_density
flex.rho_fill = 0. # infill_material_density

flex.T_e = 35000.*np.ones((50, 50)) # Elastic thickness [m] -- scalar but may be an array
flex.T_e[:,-3:] = 0.
flex.qs = np.zeros((50, 50)) # Template array for surface load stresses
flex.qs[10:40, 10:40] += 1E6 # Populating this template
flex.dx = 5000. # grid cell size, x-oriented [m]
flex.dy = 5000. # grid cell size, y-oriented [m]
# Boundary conditions can be:
# (FD): zero_displacement_zero_slope, zero_displacement_zero_moment, zero_moment_zero_shear, zero_slope_zero_shear, mirror, or periodic
# For SAS or SAS_NG, no_outside_loads is valid, and no entry defaults to this
flex.bc_west = 'zero_displacement_zero_slope' # west boundary condition
flex.bc_east = 'zero_moment_zero_shear' # east boundary condition
flex.bc_south = 'zero_displacement_zero_slope' # south boundary condition
flex.bc_north = 'zero_displacement_zero_slope' # north boundary condition

# latitude/longitude solutions are exact for SAS, approximate otherwise
#latlon = # true/false: flag to enable lat/lon input. Defaults False.
#planetary_radius = # radius of planet [m], for lat/lon solutions

# Optional: in-plane stresses [Pa] (supported by FD and FFT solvers)
#flex.sigma_xx = 0.  # east–west compression/tension
#flex.sigma_yy = 0.  # north–south compression/tension
#flex.sigma_xy = 0.  # shear — couples x and y deflection

flex.initialize()
flex.run()

# Read deflection BEFORE calling finalize (finalize clears w and qs)
deflection = flex.w  # assign to a variable for later use

# If you want to plot or save output, do it before finalize:
flex.plot_choice='both'
# An output file for deflections could also be defined here
# flex.w_out_file =
flex.output() # Plots and/or saves output, or does nothing, depending on
              # whether flex.plot_choice and/or flex.w_out_file have been set

flex.finalize()

With a configuration file

If you would like to use a Python script with a configuration file, this is also possible.

import gflex

# Pass the YAML file to the constructor; F1D for 1-D, F2D for 2-D problems.
filename = 'input/input_f1d.yaml'
flex = gflex.F1D(filename)

flex.initialize()
flex.run()

# Read deflection before finalize (finalize clears w and qs)
deflection = flex.w

# Standalone plotting output if you so desire
flex.plot_choice = 'w'
flex.output()

flex.finalize()

Within GRASS GIS

To run gFlex inside of GRASS GIS 8, install the addons from within a GRASS GIS session:

g.extension r.flexure
g.extension v.flexure

r.flexure is used for raster grids by either finite difference or analytical methods. v.flexure takes advantage of the ungridded analytical method to solve for flexure at an arbitrary set of load points, albeit limited to cases with constant elastic thickness. The source code and manual pages are available in the GRASS GIS addons repository.

When running r.flexure, it is important to ensure that the elastic thickness map is at or properly interpolated to the computational region (g.region) resolution before solving. A nearest-neighbor interpolated Te map will cause perceived gradients in elastic thickness to be very sharp, and this will strongly affect (and misdirect) the flexural solutions.

As part of Landlab and CSDMS

Landlab is an Earth-surface modeling framework built to facilitate easy integration of geomorphic, ecological, hydrological, geological, and other Earth-surface models. gFlex can be used as a Landlab component; see the Landlab repository for details.

gFlex also implements the CSDMS Basic Model Interface (BMI), enabling it to be coupled with other models in the CSDMS framework. The BMI wrapper is available as gflex.BmiGflex and requires the optional bmipy dependency:

pip install gflex[bmi]

Plotting

There are four plot choices, defined via self.plot_choice:

• 'q': plots the load in mantle-density-equivalent units of length
• 'w': plots the deflection in units of length
• 'both': plots both deflection and loads in separate panels of a 2-subplot figure
• 'combo': (1D only): plots lithospheric deflections and the deflected mantle-density-equivalent load atop it.
  • Note that the load does not affect the area above/below the datum filled when rho_fill != 0. This affects the buoyant balance associated with the motion of the plate, with no additional considerations for topography. If you would like to include topography, an iterative approach (e.g., finding areas below sea level, filling them, flexing, finding new areas below sea level, and so on) is recommended.

Utilities

The input/ directory contains several example scripts. The public API also provides standalone utility functions importable from gflex:

• flexural_wavelengths(Te, ...) — computes the flexural parameter α, first zero-crossing, and flexural wavelength for a given elastic thickness; useful for choosing grid spacing and domain size.

Domain-padding utilities (reduce spurious boundary effects when using variable Te):

2-D (F2D):

• recommended_pad_width(Te, dx, ...) — returns the recommended padding width (in cells).
• smooth_pad_Te(Te, pad_width, ...) — extends a 2-D variable-Te array with a smooth linear taper.
• pad_domain(Te, qs, dx, ...) — pads both the 2-D elastic thickness and load arrays and returns the padding width.

1-D (F1D):

• recommended_pad_width_1d(Te, dx, ...) — returns the recommended 1-D padding width (in cells).
• smooth_pad_Te_1d(Te, pad_width, ...) — extends a 1-D variable-Te array with a smooth linear taper.
• pad_domain_1d(Te, qs, dx, ...) — pads both the 1-D elastic thickness and load arrays and returns the padding width.

See the API reference for full documentation of these functions.

Credits

Development Lead

• Andrew D. Wickert

Contributors

• Eric W. H. Hutton

Release Notes

[2.0.0b1] - 2026-06-04

Breaking changes

• finalize() now clears all model state — w, qs, and the coefficient matrix are deleted. Read w and call output() before finalize(). The CLI (gflex.py) has been fixed to follow this order.

• INI configuration file format removed — only YAML (.yaml / .yml) configuration files are supported. Passing a .ini file now raises ValueError. Convert any remaining INI configs to YAML format.

• Iterative FD solver removed — F1D now always uses a direct sparse LU factorization. The ConvergenceTolerance configuration key is ignored.

• G2009 plate solution removed — the finite-difference solver now always uses the vWC1994 stencil. Setting PlateSolutionType = 'G2009' previously selected a different stencil; it is now a no-op attribute (ignored silently).

• PlateSolutionType removed from the public interface — the attribute is no longer documented or read by any solver. Existing scripts that set it will not raise an error, but the value has no effect.

• te attribute renamed to T_e — the elastic-thickness attribute on F1D and F2D is now T_e (PEP 8 subscript notation; avoids ambiguity with a time variable t and a generic elastic parameter e). flex.te now raises AttributeError. Update all occurrences: flex.te = … → flex.T_e = ….

• "mirror" is now an alias for "zero_slope_zero_shear" — the canonical BC name is "zero_slope_zero_shear" (consistent with all other long-form names); "mirror" remains a permanent short alias (like "clamped" and "free") and does not trigger any warning.

• solver attribute now raises ValueError for unsupported values — only "direct" is supported in this release. Passing any other string (e.g. "iterative") previously silently fell back to the direct solver; it now raises ValueError immediately. An iterative solver may be added in a future version, at which point additional values will become valid.

• BC string case normalised; v1.x PascalCase strings removed — all boundary-condition strings are now lowercase ("zero_displacement_zero_slope", "zero_moment_zero_shear", "zero_slope_zero_shear", "periodic"). The old v1.x PascalCase names ("0Displacement0Slope", "0Moment0Shear", "Mirror", etc.) now raise ValueError. See the boundary_conditions page for the full mapping.

New features

• Inhomogeneous (prescribed-value) boundary conditions in 1-D and 2-D FD solver. Any edge can carry a combination of prescribed displacement, slope, moment, and/or shear by passing a dict instead of a string:

  flex.bc_west = {"displacement": w_arr, "slope": dwdx_arr}
  flex.bc_east = {"moment": M0, "shear": V0}
  

  The primary use case is nested-domain (coarse-to-fine) modelling: extract w and np.gradient-based slopes at the boundary of a fine sub-domain from a coarse regional solve, impose them as BCs, and the fine solver inherits the full far-field loading through those four boundary arrays. Also supports broken-plate problems with applied edge loads. Both 1-D and 2-D; verified by a round-trip unit test (TestNestedModelGradientRoundTrip).

• gflex.VALID_BC_STRINGS_1D and gflex.VALID_BC_STRINGS_2D — public frozenset constants listing every accepted BC string (canonical names and aliases). Wrappers can validate user input against these sets instead of maintaining a parallel copy that drifts with new releases.

• "clamped", "free", and "mirror" BC aliases — concise alternatives to the full canonical names: "clamped" normalises to "zero_displacement_zero_slope"; "free" normalises to "zero_moment_zero_shear"; "mirror" normalises to "zero_slope_zero_shear". All three produce bit-identical results to their canonical names.

Performance

• LU factorization cache — set flex.cache_factorization = True to cache the sparse-LU factorization of the FD coefficient matrix across run() calls. The cache is reused automatically when the matrix is unchanged (Te, grid, BCs, and physical parameters all fixed) and invalidated by a hash comparison when any of those inputs change. Set flex.cache_factorization = "no_check" to skip the hash and reuse the factorization unconditionally — maximum performance for coupling workflows where the user guarantees matrix stability. Default False preserves existing behaviour.
• Multithreaded FFT — scipy.fft transforms now use all available CPU threads (workers=-1). Benefit is marginal at typical grid sizes; meaningful at very large grids (≳ 2000 × 2000 cells).

Improvements

• twoSurfplots() now shows three panels (load, elastic thickness, deflection) when Te is a 2-D array, using an asymmetric diverging colormap for the deflection.
• export_for_blender() writes load-cylinder geometry (radius from loaded area, height proportional to load magnitude) and per-vertex colour weights to the mesh file, enabling the companion Blender script to render a physically scaled load cylinder.
• Blender scene script: cylinder base deforms with the plate surface; camera frames the full scene vertically.

Bug fixes

• Fixed SyntaxWarning for invalid escape sequences (\sigma) in the F2D class docstring (Python 3.12+).
• Fixed a longstanding ghost-node error in the 1-D zero_displacement_zero_slope FD boundary condition. The original implementation (present since the first version) dropped ghost-node stencil terms rather than eliminating them via even reflection, and left the boundary row coupled to interior nodes rather than decoupling it as a Dirichlet constraint. The practical effect was that the boundary deflection was not strictly zero and the zero-slope condition was not enforced; results were accurate only when the boundary was far enough from any load that deflection was naturally negligible there (the recommended use case, and the intent of the existing proximity warning). The corrected implementation decouples the boundary row (enforcing w = 0 exactly) and folds the even-reflected ghost into the adjacent interior node (enforcing dw/dx = 0).
• Fixed the same two ghost-node bugs in the 2-D zero_displacement_zero_slope boundary condition (boundary rows/columns not decoupled as Dirichlet constraints; even-reflection ghost absent at the first interior row/column on all four edges). Convergence recovers from first order (O(Δx^0.92)) to second order (O(Δx^1.99)).
• Fixed a ghost-node inconsistency in the zero_moment_zero_shear FD boundary condition (1-D and 2-D). The first interior node's stencil used a shear ghost evaluated at a staggered location one cell inward, inconsistent with the moment ghost used at the boundary node itself. The corrected implementation uses the moment condition consistently at both rows, recovering second-order convergence (O(Δx^2.00) in 1-D, O(Δx^2.01) in 2-D) from first order.

Documentation

• New boundary_conditions page: comprehensive reference covering all BC types with a summary table (structural-mechanics and geophysical names), SVG diagrams for each condition, physical guidance on when each is appropriate, and a deprecation table mapping v1.x PascalCase strings to their v2.0 lowercase equivalents.
• New greenland_example page: a realistic Greenland ice-sheet flexure model (BedMachine v6 ice thickness; Steffen et al. spatially variable elastic thickness) and a hypothetical nested-domain seamount scenario that illustrates how inhomogeneous BCs couple a coarse regional model to a fine local one.
• Accuracy page extended with MMS convergence tables and figures for the 2-D zero_displacement_zero_slope ghost-node fix and the 1-D/2-D zero_moment_zero_shear ghost-node fix.

Tests

• 335 tests passing across 1-D and 2-D FD, FFT, SAS/SAS_NG solvers, all BC types, inhomogeneous BCs, domain padding, warnings, and BC aliases.
• TestNestedModelGradientRoundTrip: verifies that np.gradient-extracted slopes used as inhomogeneous Dirichlet BCs reproduce the full-domain interior to within 2 % of peak deflection, confirming the slope sign convention for all four edges in the nested-domain workflow.
• Parametrised alias tests confirm "clamped" and "free" produce bit-identical results to their canonical names in both 1-D and 2-D.

1.4.0 - 2026-05-29

New features

• FFT spectral solver for 1-D and 2-D problems. Exact for periodic boundaries; other boundary conditions are handled with 4α zero-padding (NoOutsideLoads approximation). Requires scalar (uniform) elastic thickness. In-plane stresses (sigma_xx, sigma_yy, sigma_xy) are fully supported.
• YAML configuration file format alongside the legacy INI format. Pass a .yaml / .yml file to the CLI or to the F1D / F2D constructor.
• Domain-padding utilities to reduce spurious boundary effects when using spatially variable elastic thickness:
  • 2-D: pad_domain, smooth_pad_Te, recommended_pad_width
  • 1-D: pad_domain_1d, smooth_pad_Te_1d, recommended_pad_width_1d
• FD boundary-condition warnings — F1D and F2D now issue UserWarning messages for '0Moment0Shear' (free broken end — verify a rifted margin is intended), '0Slope0Shear' (no clear geological analog), and when the nearest loaded cell is within one flexural wavelength of a '0Displacement0Slope' boundary (forebulge suppression).

Improvements

• Iterative FD solver (Solver = iterative) upgraded from Jacobi to ILU preconditioning; falls back to a direct solve if LGMRES does not converge. The ConvergenceTolerance parameter is now passed as rtol (relative residual) to SciPy's LGMRES; default remains 1e-3.
• In-plane membrane stresses (sigma_xx, sigma_yy, sigma_xy) now supported by all FD and FFT solvers in both 1-D and 2-D.

Bug fixes

• Fixed 2-D FD constant-Te stencil: corrected dx/dy swap, σ swap, and missing /dx² factors that produced incorrect results for asymmetric grids.
• Fixed 2-D FFT sigma_xy assembly: double-wrapped corner entries and wrong coefficient in the Periodic right-roll buffer.
• Fixed 0Slope0Shear description in configuration reference and README.

Documentation

• New Sphinx / Read the Docs site with a theory page (governing equations, solution methods, in-plane stresses), an accuracy page, a configuration reference, API reference, and changelog.
• Theory page covers the full governing PDE with in-plane stress terms for both 1-D and 2-D, the variable-D FD expansion, and the FFT transfer function.

Tests

• Comprehensive new test suites for 1-D SAS/SAS_NG, 2-D SAS/SAS_NG, 1-D FFT, 2-D FFT (including sigma_xy), 1-D FD boundary conditions with analytical cross-validation, 2-D FD, FD boundary-condition warnings, and all domain-padding utilities.

1.3.0 - 2026-05-27

• New: BmiGflex — CSDMS Basic Model Interface (BMI) v2 implementation; bmipy is an optional dependency and gFlex works without it.
• New: F1D, F2D, and BmiGflex are now accessible directly from the gflex package namespace (e.g. import gflex; gflex.F2D()); previously this raised AttributeError.
• Fix: updated scipy.sparse.linalg import for compatibility with modern scipy.
• Fix: latent typo self.T_e → self.Te in the scalar-Te branch of get_coeff_values (introduced 2021-06-26); would produce incorrect results when tectonic stress terms (sigma_xx, sigma_yy, sigma_xy) are non-zero.
• Removed outdated root-level gflex_bmi.py.

1.2.0 - 2024-01-08

• Python 3 modernisation: black formatting, isort, updated imports, argparse CLI, pyproject.toml replacing setup.py.
• GitHub Actions CI replacing Travis CI.
• Eric W. H. Hutton added as author.

1.1.1 - 2021-06-26

• Updated PyPI support: twine upload, README.md
• Throw meaningful error if a nonuniform self.Te grid is used with the analytical solution

1.1.0 - 2018-05-28

• Support for both Python 2 and Python 3
• Main code
• Examples
• README.md updates
• PATH updates in
• Code tests included
• Code testing on commit by Travis
• Updated on PyPI

1.0.1 - 2018-05-25

• Final Python 2 (only) release
• Additional documentation added

1.0.0 - 2017-05-10

This is the update is what is available now from pypi.

• Minor updates to the main gFlex codes
• Addition of a missing "12" in the flexural wavelength calculator in the utilities.

1.0 - 2016-01-28

• First full release of gFlex in association with the now-accepted GMD paper, "Open-source modular solutions for flexural isostasy: gFlex v1.0", by A. D. Wickert.

0.9 - 2015-03-05

• Release submitted to GMD and that appears on PyPI; this is the same as v1.0a on umn-earth-surface.

0.8.1 - 2015-03-04

• Fixed error in PyPI integration from v0.8 and updated README.md to include PyPI integration.

0.8 - 2015-03-04

• First fully-functional and error-checked release, and therefore the first true non-"pre-release".

0.7 - 2015-01-21

This release may be short-lived, but is the turning point at which work on gFlex will go from fixing deficiencies to adding new capabilities. As such, it is receiving version 0.7 (the last version of its parent project, "Flexure", was 0.6 back in 2012). It now seems to be functional and stable, and thus this tagged release will be the basis for this stage of the model to be pulled to @csdms-contrib and @umn-earth-surface.