geodynamics · baagaard-usgs · Oct 10, 2025 · Sep 6, 2025 · Oct 6, 2025 · Aug 19, 2025
diff --git a/configure.ac b/configure.ac
@@ -114,7 +114,7 @@ AC_CHECK_HEADER([mpi.h], [], [AC_MSG_ERROR([header 'mpi.h' not found])])
 
 dnl PETSC
 AC_LANG(C)
-CIT_PATH_PETSC([3.23.5])
+CIT_PATH_PETSC([3.24.0])
 CIT_HEADER_PETSC
 CIT_CHECK_LIB_PETSC
 
@@ -205,6 +205,7 @@ AC_CONFIG_FILES([Makefile
 		libsrc/pylith/problems/Makefile
 		libsrc/pylith/topology/Makefile
 		libsrc/pylith/utils/Makefile
+		libsrc/pylith/scales/Makefile
         libsrc/pylith/testing/Makefile
 		modulesrc/Makefile
 		modulesrc/include/Makefile
@@ -216,6 +217,7 @@ AC_CONFIG_FILES([Makefile
 		modulesrc/meshio/Makefile
 		modulesrc/mpi/Makefile
 		modulesrc/problems/Makefile
+		modulesrc/scales/Makefile
 		modulesrc/topology/Makefile
 		modulesrc/utils/Makefile
 		tests/Makefile
@@ -233,6 +235,7 @@ AC_CONFIG_FILES([Makefile
 		tests/libtests/problems/Makefile
 		tests/libtests/problems/data/Makefile
 		tests/libtests/meshio/data/Makefile
+		tests/libtests/scales/Makefile
 		tests/libtests/topology/Makefile
 		tests/libtests/topology/data/Makefile
 		tests/libtests/testing/Makefile

diff --git a/developer/uncrustify.cfg b/developer/uncrustify.cfg
@@ -51,8 +51,8 @@ sp_endif_cmt = Add
 sp_after_new = Add
 sp_paren_brace = Add
 sp_fparen_brace = Add
-sp_before_tr_emb_cmt = Force
-sp_num_before_tr_emb_cmt = 1
+sp_before_tr_cmt = Force
+sp_num_before_tr_cmt = 1
 
 # Align
 align_left_shift = True
@@ -98,7 +98,7 @@ nl_if_leave_one_liners = True
 nl_after_brace_open = False
 nl_after_brace_close = False
 nl_after_brace_open_cmt = False
-nl_max = 2
+nl_max = 3
 nl_after_func_proto = 2
 nl_after_func_body = 3
 nl_after_func_body_class = 2

diff --git a/docs/Makefile.am b/docs/Makefile.am
@@ -679,6 +679,7 @@ dist_noinst_DATA = \
 	user/problems/figs/hdf5layout.svg \
 	user/problems/figs/hdf5layout.tex \
 	user/problems/index.md \
+	user/problems/nondimensionalization.md \
 	user/problems/output.md \
 	user/problems/problems.md \
 	user/run-pylith/figs/cells2d.pdf \

diff --git a/docs/developer/implementation/figs/classdiagram_problem.tex b/docs/developer/implementation/figs/classdiagram_problem.tex
@@ -19,7 +19,7 @@
   \node (pyre-component) [abstract-class, anchor=west] at ($(pythia-component.east)+(2em,0)$) {\umlemptyclass{PyreComponent}};
 
   \node (python-problem) [abstract-class] at ($(petsc-component.south)-(0,4em)$) {\umlclass{Python Problem}{%
-      normalizer\\
+      scales\\
       materials\\
       bc\\
       interfaces\\
@@ -28,7 +28,7 @@
       observers
 }};
   \node (cxx-problem) [abstract-class, anchor=north] at ($(python-problem.north)+(12em,0)$) {\umlclass{C++ Problem}{%
-      normalizer\\
+      scales\\
       materials\\
       bc\\
       interfaces\\
@@ -53,8 +53,8 @@
 }};
 
 
-  \node (normalizer) [abstract-class, anchor=west] at ($(cxx-problem.east)+(12em,12em)$) {\umlemptyclass{spatialdata::units::Nondimensional}};
-  \node (material) [abstract-class] at ($(normalizer.south)-(0,1em)$) {\umlemptyclass{Material}};
+  \node (scales) [abstract-class, anchor=west] at ($(cxx-problem.east)+(12em,12em)$) {\umlemptyclass{pylith::scales::Scales}};
+  \node (material) [abstract-class] at ($(scales.south)-(0,1em)$) {\umlemptyclass{Material}};
   \node (bc) [abstract-class] at ($(material.south)-(0,1em)$) {\umlemptyclass{BoundaryCondition}};
   \node (interface) [abstract-class] at ($(bc.south)-(0,1em)$) {\umlemptyclass{FaultCohesive}};
   \node (gravity-field) [concrete-class] at ($(interface.south)-(0,1em)$) {\umlemptyclass{spatialdata::spatialdb::GravityField}};
@@ -76,7 +76,7 @@
   \draw[inherit] (cxx-problem) -- (pyre-component);
   \draw[inherit] (cxx-time-dependent) -- (cxx-problem);
 
-  \draw[aggregate] ($(cxx-problem.east)+(0,8ex)$) -- (normalizer.west);
+  \draw[aggregate] ($(cxx-problem.east)+(0,8ex)$) -- (scales.west);
   \draw[aggregate] ($(cxx-problem.east)+(0,5.5ex)$) -- (material.west);
   \draw[aggregate] ($(cxx-problem.east)+(0,3.0ex)$) -- (bc.west);
   \draw[aggregate] ($(cxx-problem.east)+(0,0.5ex)$) -- (interface.west);

diff --git a/docs/developer/testing/libtests.md b/docs/developer/testing/libtests.md
@@ -60,7 +60,7 @@ namespace pylith {
 class pylith::problems::TestPhysics : public CppUnit::TestFixture {
     CPPUNIT_TEST_SUITE(TestPhysics); // CppUnit macro used to define the `TestPhysics` test suite.
 
-    CPPUNIT_TEST(testSetNormalizer); // CppUnit macro to add test implemented by class method.
+    CPPUNIT_TEST(testSetScales); // CppUnit macro to add test implemented by class method.
     CPPUNIT_TEST(testSetAuxiliaryFieldDB);
     CPPUNIT_TEST(testSetAuxiliarySubfieldDiscretization);
     CPPUNIT_TEST(testObservers);
@@ -81,7 +81,7 @@ public:
 
     // Methods that implement tests. The naming convention is testMethodName.
     // All test methods must return void and not have any arguments.
-    void testSetNormalizer(void);
+    void testSetScales(void);
     void testSetAuxiliaryFieldDB(void);
     void testSetAuxiliarySubfieldDiscretization(void);
     void testObservers(void);
@@ -129,20 +129,20 @@ pylith::problems::TestPhysics::tearDown(void) {
 
 
 void
-pylith::problems::TestPhysics::testSetNormalizer(void) {
+pylith::problems::TestPhysics::testSetScales(void) {
     PYLITH_METHOD_BEGIN;
 
-    spatialdata::units::Nondimensional normalizer;
+    pylith::scales::Scales scales;
     const PylithReal lengthScale = 3.0;
-    normalizer.setLengthScale(lengthScale);
+    scales.setLengthScale(lengthScale);
 
     CPPUNIT_ASSERT(_physics);
-    _physics->setNormalizer(normalizer);
+    _physics->setScales(scales);
 
-    CPPUNIT_ASSERT_DOUBLES_EQUAL(lengthScale, _physics->_normalizer->getLengthScale(), 1.0e-6);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(lengthScale, _physics->_scales->getLengthScale(), 1.0e-6);
 
     PYLITH_METHOD_END;
-} // testSetNormalizer
+} // testSetScales
 
 
 void

diff --git a/docs/user/components/problems/GreensFns.md b/docs/user/components/problems/GreensFns.md
@@ -25,9 +25,9 @@ Implements `Problem`.
 * `materials`: Materials in problem.
   - **current value**: 'homogeneous', from {default}
   - **configurable as**: homogeneous, materials
-* `normalizer`: Nondimensionalizer for problem.
+* `scales`: Nondimensionalizer for problem.
   - **current value**: 'nondimelasticquasistatic', from {default}
-  - **configurable as**: nondimelasticquasistatic, normalizer
+  - **configurable as**: nondimelasticquasistatic, scales
 * `petsc_defaults`: Flags controlling which default PETSc options to use.
   - **current value**: 'petscdefaults', from {default}
   - **configurable as**: petscdefaults, petsc_defaults

diff --git a/docs/user/components/problems/Problem.md b/docs/user/components/problems/Problem.md
@@ -29,9 +29,9 @@ If the nonlinear (SNES) solver requires multiple iterations to converge for thes
 * `materials`: Materials in problem.
   - **current value**: 'homogeneous', from {default}
   - **configurable as**: homogeneous, materials
-* `normalizer`: Nondimensionalizer for problem.
+* `scales`: Nondimensionalizer for problem.
   - **current value**: 'nondimelasticquasistatic', from {default}
-  - **configurable as**: nondimelasticquasistatic, normalizer
+  - **configurable as**: nondimelasticquasistatic, scales
 * `petsc_defaults`: Flags controlling which default PETSc options to use.
   - **current value**: 'petscdefaults', from {default}
   - **configurable as**: petscdefaults, petsc_defaults

diff --git a/docs/user/components/problems/TimeDependent.md b/docs/user/components/problems/TimeDependent.md
@@ -28,9 +28,9 @@ Implements `Problem`.
 * `materials`: Materials in problem.
   - **current value**: 'homogeneous', from {default}
   - **configurable as**: homogeneous, materials
-* `normalizer`: Nondimensionalizer for problem.
+* `scales`: Nondimensionalizer for problem.
   - **current value**: 'nondimelasticquasistatic', from {default}
-  - **configurable as**: nondimelasticquasistatic, normalizer
+  - **configurable as**: nondimelasticquasistatic, scales
 * `petsc_defaults`: Flags controlling which default PETSc options to use.
   - **current value**: 'petscdefaults', from {default}
   - **configurable as**: petscdefaults, petsc_defaults
@@ -89,8 +89,8 @@ ic = [domain]
 # Turn on gravitational body forces
 gravity_field = spatialdata.spatialdb.GravityField
 
-# Set the normalizer for nondimensionalizing the problem
-normalizer = spatialdata.units.NondimElasticQuasistatic
+# Set the scales for nondimensionalizing the problem
+scales = pylith.scales.NondimElasticQuasistatic
 
 # Set the subfields in the solution
 solution = = pylith.problems.SolnDispLagrange

diff --git a/docs/user/examples/box-2d/step05-sheardisptractrate.md b/docs/user/examples/box-2d/step05-sheardisptractrate.md
@@ -34,7 +34,7 @@ start_time = -1.0*year
 end_time = 5.0*year
 initial_dt = 1.0*year
 
-[pylithapp.problem.normalizer]
+[pylithapp.problem.scales]
 relaxation_time = 10.0*year
 ```
 

diff --git a/docs/user/examples/box-3d/step05-sheardisptractrate.md b/docs/user/examples/box-3d/step05-sheardisptractrate.md
@@ -34,7 +34,7 @@ initial_dt = 1.0*year
 start_time = -1.0*year
 end_time = 5.0*year
 
-[pylithapp.problem.normalizer]
+[pylithapp.problem.scales]
 relaxation_time = 10.0*year
 ```
 

diff --git a/docs/user/examples/magma-2d/common-information.md b/docs/user/examples/magma-2d/common-information.md
@@ -38,10 +38,10 @@ pressure.basis_order = 1
 trace_strain.basis_order = 1
 
 [pylithapp.problem]
-normalizer = spatialdata.units.NondimElasticQuasistatic
-normalizer.length_scale = 100.0*m
-normalizer.relaxation_time = 0.2*year
-normalizer.shear_modulus = 10.0*GPa
+scales = pylith.scales.NondimElasticQuasistatic
+scales.length_scale = 100.0*m
+scales.relaxation_time = 0.2*year
+scales.shear_modulus = 10.0*GPa
 ```
 
 We use the material properties in all of the simulations in this directory, so we specify them in `pylithapp.cfg` to avoid repeating the information in the file with parameters for each simulation.

diff --git a/docs/user/examples/poroelastic-outerrise-2d/common-information.md b/docs/user/examples/poroelastic-outerrise-2d/common-information.md
@@ -50,10 +50,10 @@ auxiliary_subfields.isotropic_permeability.basis_order = 2
 caption: Nondimensionalization parameters for the 2D outer-rise examples with poroelasticity.
 ---
 [pylithapp.problem]
-normalizer = spatialdata.units.NondimElasticQuasistatic
-normalizer.length_scale = 100.0*m
-normalizer.relaxation_time = 1*year
-normalizer.shear_modulus = 10.0*GPa
+scales = pylith.scales.NondimElasticQuasistatic
+scales.length_scale = 100.0*m
+scales.relaxation_time = 1*year
+scales.shear_modulus = 10.0*GPa
 ```
 
 ```{code-block} cfg

diff --git a/docs/user/examples/reverse-2d/step07-twofaults-maxwell.md b/docs/user/examples/reverse-2d/step07-twofaults-maxwell.md
@@ -22,7 +22,7 @@ initial_dt = 4.0*year
 start_time = -4.0*year
 end_time = 100.0*year
 
-normalizer.relaxation_time = 20.0*year
+scales.relaxation_time = 20.0*year
 ```
 
 ```{code-block} cfg

diff --git a/docs/user/governingeqns/elasticity/index.md b/docs/user/governingeqns/elasticity/index.md
@@ -3,6 +3,7 @@
 
 :::{toctree}
 derivation.md
+nondimensionalization.md
 infinitesimal-strain.md
 prescribed-slip.md
 bulk-rheologies/index.md

diff --git a/docs/user/governingeqns/elasticity/nondimensionalization.md b/docs/user/governingeqns/elasticity/nondimensionalization.md
@@ -0,0 +1,125 @@
+# Nondimensionalization
+
+Starting with the elasticity equation,
+%
+\begin{gather}
+\rho(\vec{x})\frac{\partial^{2}\vec{u}}{\partial t^{2}}-\vec{f}(\vec{x},t)-\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}=\vec{0}\text{ in }\Omega,\\
+\boldsymbol{\sigma}(\vec{u})\cdot\vec{n}=\vec{\tau}(\vec{x},t)\text{ on }\Gamma_{\tau}\text{,}\\
+\vec{u}^{+}-\vec{u}^{-}=\vec{d}\text{ on }\Gamma_{f}.
+\end{gather}
+%
+we define nondimensional values:
+%
+\begin{align}
+\vec{x}^* &= \frac{\vec{x}}{x_o}, \\
+\vec{u}^* &= \frac{\vec{u}}{u_o}, \\
+\rho^* &= \frac{\rho}{\rho_o}, \\
+\vec{f}^* &= \frac{\vec{f}}{f_o}, \\
+\boldsymbol{\sigma}^* &= \frac{\boldsymbol{\sigma}}{\sigma_o}, \\
+\vec{\tau}^* &= \frac{\vec{\tau}}{\sigma_o}.
+\end{align}
+%
+We also recognize that
+%
+\begin{equation}
+\boldsymbol{\nabla}^* = x_o \boldsymbol{\nabla}.
+\end{equation}
+
+Substituting into the equations, we have
+%
+\begin{gather}
+\rho_o  \rho^*(\vec{x}^*) \frac{u_o}{t_o^2} \frac{\partial^{2}\vec{u}^*}{\partial {t^*}^{2}} - f_o \vec{f}^*(\vec{x}^*,t^*) - \frac{\sigma_o}{x_o} \boldsymbol{\nabla}^*\cdot\boldsymbol{\sigma}^* = \vec{0}\text{ in }\Omega,\\
+\sigma_o \boldsymbol{\sigma}^*(\vec{u^*}) \cdot \vec{n} = \sigma_o \vec{\tau}^*(\vec{x}^*,t^*)\text{ on }\Gamma_{\tau}\text{,}\\
+u_o \left(\vec{u}^{*^{+}}-\vec{u}^{*^{-}}\right) = u_o \vec{d}^* \text{ on }\Gamma_{f}.
+\end{gather}
+
+For the second two equations, the nondimensional scales in for the terms in each equation are consistent, so we will limit our discussion to the first equation.
+Grouping terms and multiplying by $\frac{x_o}{\sigma_o}$, we have
+%
+\begin{equation}
+\left( \rho_o \frac{u_o}{t_o^2}\frac{x_o}{\sigma_o}\right) \rho^*(\vec{x}^*) \frac{\partial^{2}\vec{u}^*}{\partial {t^*}^{2}} - \left(f_o \frac{x_o}{\sigma_o}\right) \vec{f}^*(\vec{x}^*,t^*) -  \boldsymbol{\nabla}^*\cdot\boldsymbol{\sigma}^* = \vec{0}\text{ in }\Omega.
+\end{equation}
+%
+All terms should be nondimensional, which implies
+%
+\begin{align}
+f_o &= \frac{\sigma_o}{x_o}, \\
+\rho_o &= \frac{t_o^2 \sigma_o}{u_o x_o}.
+\end{align}
+
+We want to determine the stress scale, $\sigma_o$.
+Considering isotropic, linear elasticity we have
+%
+\begin{equation}
+\boldsymbol{\sigma} = \boldsymbol{C} : \boldsymbol{\epsilon} = \boldsymbol{C} : \frac{1}{2}\left(\boldsymbol{\nabla} + \boldsymbol{\nabla}^T \right) \vec{u}.
+\end{equation}
+%
+Substituting in our nondimensional values yields
+%
+\begin{equation}
+\sigma_o \boldsymbol{\sigma}^* = \mu_o \boldsymbol{C}^* : \frac{u_o}{x_o} \frac{1}{2}\left(\boldsymbol{\nabla}^* + \boldsymbol{\nabla}^{*^T}\right) \vec{u}^*.
+\end{equation}
+%
+We recognize that for the equation to be nondimensional,
+\begin{equation}
+\sigma_o = \mu_o \frac{u_o}{x_o}.
+\end{equation}
+
+Returning to the expression for $\rho_o$ and substituting in the expression for $\sigma_o$, we have
+\begin{align}
+\rho_o &= \frac{t_o^2 \sigma_o}{u_o x_o}, \\
+\rho_o &= \mu_o \frac{t_o^2}{x_o^2}.
+\end{align}
+
+## Inertia
+
+The scale of the inertial term is $\Pi_\mathit{inertia} = \frac{\rho_o x_o^2}{\mu_o t_o^2}$.
+If this term is O(1), then we should include inertia and solve the dynamic form of the elasticity equation.
+If this term is very small, then we can neglect inertia and solve the quasistatic form of the elasticity equation.
+
+**If the time scale equals the time it takes the shear wave ($v_o^2 = \frac{\mu_o}{\rho_o}$) to propagate over the length scale**, then
+\begin{align}
+t_o &= \frac{x_o}{v_o}, \\
+t_o^2 &= \frac{x_o^2}{v_o^2}, \\
+t_o^2 &= \frac{x_o^2 \rho_o}{\mu_o}.
+\end{align}
+Substituting into the expression for $\Pi_\mathit{inertia}$, we have
+\begin{align}
+\Pi_\mathit{inertia} &= \frac{\rho_o x_o^2}{\mu_o t_o^2}, \\
+\Pi_\mathit{inertia} &= \frac{\rho_o x_o^2}{\mu_o} \frac{\mu_o}{\rho_o x_o^2}, \\
+\Pi_\mathit{inertia} &= 1.
+\end{align}
+In this case, the inertial term is important, and we have the dynamic case with propagating seismic waves.
+
+**If the time scale is much larger than the time it takes the shear wave to propagate over the length scale**, then
+\begin{align}
+t_o &\gg \frac{x_o}{v_o}, \\
+t_o^2 &\gg \frac{x_o^2 \rho_o}{\mu_o}.
+\end{align}
+Substituting into the expression for $\Pi_\mathit{inertia}$, we have
+\begin{align}
+\Pi_\mathit{inertia} &= \frac{\rho_o x_o^2}{\mu_o t_o^2}, \\
+\Pi_\mathit{inertia} &\ll \frac{\rho_o x_o^2}{\mu_o} \frac{\mu_o}{\rho_o x_o^2}, \\
+\Pi_\mathit{inertia} &\ll 1.
+\end{align}
+In this case, the inertial term is negligible, and we have quasistatic elasticity.
+:::
+
+:::{note}
+In PyLith v4 and earlier, we used a displacement scale equal to the length scale.
+Using separate length and displacement scales facilitates accurately solving for displacements many orders of magnitude smaller than the length scale.
+It also separates the stress scale used to nondimensionalize stress, tractions, and fluid pressure from the rigidity scale.
+:::
+
+```{table} Nondimensionalization of elasticity.
+:name: tab:elasticity:scales
+| **Quantity** | **Scale**          |
+| :---------: | :------------------ |
+| $x_o$       | Length scale        |
+| $u_o$       | Displacement scale  |
+| $\mu_o$     | Rigidity scale      |
+| $t_o$       | Time scale          |
+| $\sigma_o = \mu_o \frac{u_o}{x_o}$  | Stress scale        |
+| $f_o = \frac{\sigma_o}{x_o}$    | Body force scale       |
+| $\rho_o = \mu_o \frac{t_o^2}{x_o^2}$    | Density scale       |
+```
diff --git a/docs/user/governingeqns/incompressible-elasticity/index.md b/docs/user/governingeqns/incompressible-elasticity/index.md
@@ -3,5 +3,6 @@
 
 :::{toctree}
 infinitesimal-strain.md
+nondimensionalization.md
 bulk-rheologies/index.md
 :::
diff --git a/docs/user/governingeqns/incompressible-elasticity/infinitesimal-strain.md b/docs/user/governingeqns/incompressible-elasticity/infinitesimal-strain.md
@@ -15,7 +15,7 @@ The strong form is
   % Solution
   \vec{s}^T = \left( \vec{u} \quad \ p \right)^T, \\
   % Elasticity
-\rho \frac{\partial^2\vec{u}}{\partial t^2} - \vec{f}(\vec{x},t) - \left(\boldsymbol{\sigma}^\mathit{dev}(\vec{u}) - p\boldsymbol{I}\right) = \vec{0} \text{ in }\Omega,
+\rho \frac{\partial^2\vec{u}}{\partial t^2} - \vec{f}(\vec{x},t) - \left(\boldsymbol{\sigma}^\mathit{dev}(\vec{u}) - p\boldsymbol{I}\right) = \vec{0} \text{ in }\Omega, \\
   % Pressure
   \vec{\nabla} \cdot \vec{u} + \frac{p}{K} = 0 \text{ in }\Omega, \\
   % Neumann