docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.py on development builds its gmsh geometry in kilometres and then hands the plex to uw.discretisation.Mesh unscaled. UW3 reads raw DM coordinates as already non-dimensional, so the model it actually solves is 30x too large in every direction.
The script's own comment says the domain is 4H x H = 120 km x 30 km. It builds 3600 x 900 km.
Where
H = uw.quantity(30, "km")
model.set_reference_quantities(length = H, ...)
S = H.magnitude # 30.0 - scale factor: non-dimensional -> km
... # every gmsh point is multiplied by S
mesh1 = uw.discretisation.Mesh(plex[1], boundaries=boundaries, ...) # no rescale
Measured
Loading the mesh both ways and printing what the solver sees:
| build |
non-dim extent |
dimensional |
| as written |
x = -60..60, y = -30..0 |
3600 x 900 km |
| DM coordinates divided by H first |
x = -2..2, y = -1..0 |
120 x 30 km |
eta_bg non-dimensionalises to exactly 1 in both cases, so viscosity is not the tell — only the mesh extent is.
Why it is fatal rather than cosmetic
Non-dimensionally the momentum equation has one group, rho*g*L^2/(eta*V), and the yield law adds C*L/(eta*V). Inflating L by 30 multiplies the first by 900 and the second by 30, so the body force overwhelms the imposed convergence and the cohesion/friction balance is shifted to a point no correctly-scaled version of this benchmark reaches.
Concretely, at the benchmark's own default parameters (eta_bg = 1e24, V = 2.5 mm/yr), taking the domain-average viscous stress 2*eta*edot with edot = 2V/width:
| build |
tau_y / viscous stress |
| as written |
136.5 — nothing yields anywhere |
| corrected |
0.188 — partially yielded |
So as shipped, the plasticity benchmark runs with no active plasticity at its own reference parameters.
Fix
Divide the DM coordinates by H before constructing the Mesh, keeping length = 30 km as the reference:
plex = uw.discretisation._from_gmsh(msh, useMultipleTags=True, useRegions=True)
_cv = plex[1].getCoordinates()
_cv.scale(1.0 / H.magnitude)
plex[1].setCoordinates(_cv)
mesh1 = uw.discretisation.Mesh(plex[1], ...)
Equivalently, build the gmsh geometry in non-dimensional units in the first place (drop the * S), which is simpler and removes the trap at source.
Note on the misleading comment
The script says "The Model must be set up before mesh creation so that the mesh inherits the correct coordinate units (km)". That is true but does not do what the comment implies: patch_coordinate_units only labels the coordinate symbols with units, it never rescales the DM. Worth correcting the comment wherever it appears, since it is the reason the bug looks safe on reading.
This was found while auditing a research harness that had copied the same mesh-construction block; the same defect is in drivers/convergence.py and drivers/fmg_contraction_probe.py in that harness, so it has propagated.
docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.pyondevelopmentbuilds its gmsh geometry in kilometres and then hands the plex touw.discretisation.Meshunscaled. UW3 reads raw DM coordinates as already non-dimensional, so the model it actually solves is 30x too large in every direction.The script's own comment says the domain is
4H x H = 120 km x 30 km. It builds 3600 x 900 km.Where
Measured
Loading the mesh both ways and printing what the solver sees:
eta_bgnon-dimensionalises to exactly 1 in both cases, so viscosity is not the tell — only the mesh extent is.Why it is fatal rather than cosmetic
Non-dimensionally the momentum equation has one group,
rho*g*L^2/(eta*V), and the yield law addsC*L/(eta*V). Inflating L by 30 multiplies the first by 900 and the second by 30, so the body force overwhelms the imposed convergence and the cohesion/friction balance is shifted to a point no correctly-scaled version of this benchmark reaches.Concretely, at the benchmark's own default parameters (
eta_bg = 1e24,V = 2.5 mm/yr), taking the domain-average viscous stress2*eta*edotwithedot = 2V/width:So as shipped, the plasticity benchmark runs with no active plasticity at its own reference parameters.
Fix
Divide the DM coordinates by H before constructing the Mesh, keeping
length = 30 kmas the reference:Equivalently, build the gmsh geometry in non-dimensional units in the first place (drop the
* S), which is simpler and removes the trap at source.Note on the misleading comment
The script says "The Model must be set up before mesh creation so that the mesh inherits the correct coordinate units (km)". That is true but does not do what the comment implies:
patch_coordinate_unitsonly labels the coordinate symbols with units, it never rescales the DM. Worth correcting the comment wherever it appears, since it is the reason the bug looks safe on reading.This was found while auditing a research harness that had copied the same mesh-construction block; the same defect is in
drivers/convergence.pyanddrivers/fmg_contraction_probe.pyin that harness, so it has propagated.