Repository navigation
Tool for Weertman to Regularized Coulomb basal friction conversion - #752
Draft
matthewhoffman wants to merge 77 commits into
Draft
matthewhoffman wants to merge 77 commits into
matthewhoffman wants to merge 77 commits into
Conversation
- Add SECONDS_PER_YEAR constant and use it in the Lambda formula (Lambda = uc / (SECONDS_PER_YEAR * A * N^n)); it was previously missing, making Lambda ~3.15e7x too large (excess friction/no-slip). - Scale N by ALBANY_EFFECTIVE_PRESSURE_PA_PER_UNIT (1000) when deriving the Coulomb Friction Coefficient C, to match Albany's internal km/kPa-scaled effective pressure representation; C was previously 1000x too small relative to Albany's runtime N. - Fix suggested YAML: Effective Pressure Type corrected from the invalid "From Surface" to "Hydrostatic". - Add "Use Pressurized Bed Above Sea Level: true" to the suggested YAML, without which Albany ignores Minimum Fraction Overburden Pressure / Length Scale Factor entirely. - Convert the printed Length Scale Factor to km to match Albany's internally km-scaled bed field. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
- Replace the previous critical-velocity-based Lambda formula with an exact per-cell solve: Lambda is chosen so that Albany's Regularized Coulomb law reproduces the same basal shear stress the Weertman law would produce at the cell's actual current sliding speed (from uReconstructX/uReconstructY, last nVertInterfaces level, converted m/s -> m/yr). - Weertman basal shear stress is computed as Tau_b = mu * u^qW, with NO effective-pressure term, matching MALI's actual Weertman sliding law (muFriction was calibrated for this convention, not Albany's generic N-dependent Power Law). - Cells where the Weertman stress at the current speed already meets or exceeds the Coulomb limit C*N have no valid non-negative Lambda; these are set to 0 (maximal Coulomb sliding) and counted/reported. - Generalized the vertical-level extraction helper to also match nVertInterfaces (not just nVertLevels), and added --velocity-x-field/--velocity-y-field CLI options. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Replace area_weighted_optimal_C() with fit_coulomb_C_fast_region(): - Identify the fast-flowing region as grounded cells with actual current sliding speed (from velocity-x/y-field) greater than --critical-velocity. - Assume that region is already in the fully-plastic Coulomb regime of the RC law (Tau_b = C * N), and fit C as the area-weighted least-squares best match of C * N to the Weertman law's basal shear stress (mu * speed^qW, no effective-pressure term) in that region. - Lambda is then solved exactly per cell (unchanged logic) so the RC law reproduces the Weertman shear stress at the cell's actual current speed, using this new C. - Cells where the Weertman shear stress at the current speed already meets or exceeds the fitted Coulomb limit C*N have no valid Lambda solution; Lambda is set to 0 there, which is not an arbitrary filler but the physically-correct full-Coulomb limit (Tau_b_RC saturates at C*N as Lambda -> 0, independent of speed). - Updated diagnostics/local_C to reflect the new fast-region fit. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…egion; add diagnostic masks - fit_coulomb_C_fast_region(): replace N^2-weighted least-squares fit (which was dominated by a small number of high-N, high-leverage cells) with a straight area-weighted mean of the per-cell implied C = Tau_b_Weertman / N over the fast-flowing region. - Lambda/bedRoughnessRC: fast-flowing cells (speed > critical velocity, the same region used to fit C) are now always assumed to be in the full-Coulomb regime and forced to a reference value (new --lambda-reference-value, default 0.0), rather than attempting (and failing) an exact per-cell solve there. The exact per-cell solve is now only attempted for the slow-flowing grounded region. - Added three new int8 mask fields to the output NetCDF for evaluation: maskGrounded, maskFastFlowing, and maskValidBedRoughnessRC (marks cells where Lambda was solved exactly vs. set to the full-Coulomb reference value). - Fixed a diagnostics bug (local_C/basal sliding speed range) that arose from separating the fast-flowing and slow-flowing regions. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
- Fixed broken syntax (missing commas in signature, stray duplicate
closing paren/return-type annotation) and a real bug (q_start =
min(q_start, q_inland) used Python's builtin min() on numpy arrays,
which does not broadcast correctly; replaced with np.minimum()) in
the user-supplied effective_pressure4() function.
- Guarded against division by zero at ice-free cells (H == 0).
- Exposed h_ocean as a new parameter/CLI flag (--transition-h-ocean,
default 0.025) instead of a hardcoded magic number; changed the
function to return only N (matching
downs_johnson_effective_pressure()'s convention), dropping the
previously-unused N_inland/N_ocean tuple outputs.
- Added --effective-pressure-type {downs-johnson,transition} to
select between Albany's own internal Hydrostatic-At-Nodes N formula
(default, unchanged) and this new near-ocean/inland-transition
parameterization; added --transition-alpha/--transition-length-
scale/--transition-h-ocean CLI args (required only for
'transition'), with mutual-exclusivity validation against the
downs-johnson-specific args.
- Since Albany does not yet support consuming a precomputed
effective-pressure Field for the Regularized Coulomb law, the
suggested YAML for 'transition' prints a placeholder
'Effective Pressure Type: TRANSITION OPTION TO BE ADDED' line and a
NOTE explaining this is offline-evaluation-only pending upstream
Albany support.
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
The default h_ocean value of 0.025 was originally specified in km (25 m), but effective_pressure4() compares it directly against height_above_flotation, which is computed in meters (consistent with thickness/bed elsewhere in the script). Changed the default to 25.0 m so it is self-consistent with the rest of the script's unit conventions. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…ydropotential fields
- New --diagnostics/--no-diagnostics mutually-exclusive flag
(default: --diagnostics, i.e. enabled) controls whether diagnostic
fields (effectivePressure, flowRateA, floatationFraction,
hydropotential, maskGrounded, maskFastFlowing,
maskValidBedRoughnessRC) are written to the output NetCDF. The
required bedRoughnessRC field is always written regardless of this
flag; muFriction is always removed regardless of this flag.
- Added two new diagnostic fields:
- floatationFraction = Pw / Pice, with Pw = Pice - N (basal water
pressure) and Pice = rho_i * gravity * thickness (ice overburden
pressure).
- hydropotential: Shreve hydraulic potential,
phi = rho_water * gravity * bedTopography + Pw.
- Added --floatation-fraction-field/--hydropotential-field to name
these new output fields (defaults: floatationFraction,
hydropotential).
- effectivePressure's long_name now reflects which
--effective-pressure-type was actually used (Downs-Johnson vs. the
near-ocean/inland-transition parameterization), rather than always
saying 'Downs-Johnson'.
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…potential - New load_transect()/project_transect()/plot_transects() functions: reads a geometric_features-style transect.geojson directly (no dependency on the geometric_features python package, since these flowline transects are not yet exposed by it), projects the lon/lat LineString into the MALI mesh's planar CRS (EPSG:3031 for Antarctica, matching the convention used elsewhere in MPAS-Tools' ismip7_postprocessing/grid_and_mapping.py), samples effectivePressure/floatationFraction/hydropotential at the nearest MALI cell center (scipy cKDTree) along the transect, and saves a stacked-panel PNG plot vs. along-transect distance. - New CLI args: --transects-dir, --transect-names (default: Thwaites, Totten, Lambert, Foundation, Bindschadler), --plot-dir, and a --plot-transects/--no-plot-transects flag pair (default: disabled). - --plot-transects requires --diagnostics (the fields plotted are diagnostic-only) and --transects-dir; validated via parser.error(). - pyproj/scipy/matplotlib are only imported lazily inside plot_transects(), so the rest of the script has no new hard dependencies when --plot-transects is not used. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
- --min-fraction-overburden and --pressure-length-scale are now shared between --effective-pressure-type=downs-johnson and --effective-pressure-type=transition, instead of duplicating them as --transition-alpha/--transition-length-scale. This makes it convenient to switch --effective-pressure-type without having to re-specify a second, differently-named set of equivalent args. - effective_pressure4()'s alpha parameter (a deficit fraction with no particular meaning beyond a quick prototype) is replaced with min_fraction_overburden (a retained fraction), matching downs_johnson_effective_pressure()'s convention exactly: N_inland = rho_i * gravity * H * min_fraction_overburden. - --transition-h-ocean remains transition-only (no analog exists in the downs-johnson formula). - Simplified the effective-pressure-type CLI validation and output NetCDF global-attribute writing accordingly, since both types now share the same two required args. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
- plot_transects() now takes thickness/bed/rho_i/rho_w and classifies each MALI cell as grounded ice (unshaded), floating ice, ice-free ocean, or ice-free land, using the same grounded-ice test used elsewhere in the script (rho_i * H + rho_w * bed > 0). - Background shading (axvspan) is drawn on every panel for contiguous along-transect runs of each non-grounded class, with a shared legend at the top of the figure. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…panel - plot_transects() now accepts min_fraction_overburden and draws a horizontal dashed reference line at 1 - min_fraction_overburden on the floatation-fraction panel: the maximum floatation fraction (Pw/Pice) --min-fraction-overburden implies inland, since floatation_fraction = 1 - N/Pice and N/Pice is bounded below by min_fraction_overburden. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Two bugs found by inspecting real output: 1. effective_pressure4() and the transect-plot bound line used min_fraction_overburden as the *retained* overburden fraction, but Albany's actual HYDROSTATIC formula treats it as the *subtracted* fraction, so the floatation-fraction floor attained far inland is min_fraction_overburden itself, not its complement. This was introduced when --min-fraction-overburden/--transition-alpha were merged: effective_pressure4()'s q_inland is now 1.0 - min_fraction_overburden (was min_fraction_overburden), and the plotted bound line is now min_fraction_overburden directly (was 1 - min_fraction_overburden). Docstrings/help text corrected to match. 2. Shreve hydraulic potential reused --rho-water (ocean/seawater density, 1028 kg/m^3) for the bed-elevation term, but subglacial water is fresh, not salty. Added a new --rho-freshwater CLI arg (default 1000.0 kg/m^3) used only in the hydropotential formula; --rho-water is now documented as ocean/seawater-only. Also documented (via a NetCDF attribute on effectivePressure) that some floating cells can retain nonzero N/floatationFraction != 1 near the grounding line: this is expected Albany HYDROSTATIC behavior (f_p depends on bed elevation only, not the true flotation criterion), not a bug. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
plot_transects() now always shows a first panel with bed elevation and ice-surface elevation (computed via the flotation criterion for floating ice), ahead of the caller-supplied N/floatation-fraction/ hydropotential panels. Useful for interpreting the other panels (e.g. seeing bed troughs that explain elevated floatation fractions inland). Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
New plot_maps() function uses the 'mosaic' package (available in the e3sm-unified conda environment) to render unstructured-mesh maps directly in the MALI mesh's native planar polar-stereographic coordinates -- no reprojection needed. Automatically produced alongside the transect plots (same --plot-transects trigger, no new CLI flag) in a 'maps' subdirectory of --plot-dir, covering bedRoughnessRC (Lambda), all diagnostic fields (effectivePressure, floatationFraction, hydropotential, flowRateA), and all diagnostic masks (maskGrounded, maskFastFlowing, maskValidBedRoughnessRC). Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
plot_maps() now supports multi-panel figures with an optional log-scale colormap per panel (non-positive values masked out). bedRoughnessRC is now plotted log-scale, side by side with the original input muFriction field (also log-scale) for comparison. All diagnostic plots (transects and maps) are now written directly into a single --plot-dir (default: diagnostic_plots), with no maps/ subdirectory, per user request for one flat CLI-controlled directory. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
plot_maps() panels now accept an optional per-panel colormap name.
bedRoughnessRC and the comparison muFriction panel now use 'turbo'
(more perceptually distinct hues than viridis, revealing more spatial
detail), with muFriction's colormap reversed ('turbo_r') so that
similar colors indicate the physically inverse relationship between
the two fields (low mu ~ high Lambda and vice versa), making them
easier to visually compare.
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
plot_maps() panel specs are now dicts (values/units/title/log/cmap/ vmin/vmax/mask) instead of positional tuples, adding support for clipping the colormap to a fixed range and masking out cells (set to NaN, not displayed) based on a boolean array. floatationFraction's map now clips the colorbar to [0, 1] (a handful of near-zero-thickness outlier cells previously produced ratios up to ~1500, which is only physically meaningful for grounded ice with non-negligible overburden) and masks out non-grounded cells entirely, since Pw/Pice is not a meaningful/valid quantity there. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
N is only physically meaningful for grounded ice; hide floating/ ice-free cells on the map (same mask approach as floatationFraction). Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Same grounded-ice masking as effectivePressure/floatationFraction, for consistency across the diagnostic maps. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
plot_maps() now accepts an optional transects=[(name, x, y), ...] list (already projected into the mesh's planar CRS) and draws each as a labeled red line on every panel of every map. Enabled by default (--plot-transects-on-maps / --no-plot-transects-on-maps), reusing the existing --transect-names/--transects-dir plumbing and the load_transect()/project_transect() helpers already used for the transect line plots -- no new transect-specific inputs needed. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Adds interpolated bed=0 crossing distances (a rough grounding- line/coastline proxy) as vertical gray dotted lines on every panel of each transect plot, to help visually correlate feature changes (e.g. in N, floatationFraction, hydropotential) with topographic context. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Adds an alternative --method transition-velocity for deriving
bedRoughnessRC (Lambda) and the Regularized Coulomb coefficient
(--mu-field), alongside the existing stress-match-fit method:
- Choose a fixed transition velocity u0 (--transition-velocity) and
compute, in closed form for every grounded cell:
Lambda = u0 / (SECONDS_PER_YEAR * A * N^n)
C = tau_b * (ub + u0)^(1/3) / (N * ub^(1/3))
so the RC law exactly reproduces the Weertman law's basal shear
stress at each cell's current sliding speed, with no fast/slow
region split (see solve_transition_velocity()).
- C is spatially uniform for stress-match-fit (unchanged) but
spatially varying for transition-velocity.
- New --transition-velocity/--u0 and --mu-reference-value CLI args;
--critical-velocity is now required only for stress-match-fit.
- Diagnostics/masks/plots/output attrs adapted per method:
maskFastFlowing means speed > u0 for the new method,
maskValidBedRoughnessRC covers every valid closed-form solve, and
the "implied C" diagnostic/plot is stress-match-fit only.
- Added a map panel showing the actual output --mu-field (C) values
alongside Lambda and the original input muFriction, for both
methods.
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Move the output --mu-field (calculated Regularized Coulomb C) panel out of the bedRoughnessRC map figure and into its own dedicated map file (map_<mu-field>.png, e.g. map_muFriction.png), alongside the original input muFriction, so it's easy to find for both methods. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add --extrapolate-terminus-cells (default disabled), which discards the computed bedRoughnessRC (Lambda), and muFriction (C) for --method=transition-velocity, at the outermost row of grounded cells adjacent to the grounding line or a grounded marine terminus, plus every non-grounded cell (floating ice, ice-free ocean, ice-free land). Those discarded values are then creep-filled from the remaining grounded interior across MPAS mesh connectivity (cellsOnCell/nEdgesOnCell), using either inverse-distance weighting or a minimum-value rule (--creep-fill-method), adapted from MPAS-Tools' conversion_exodus_init_to_mpasli_mesh.py extrapolation loop. Adds two new helpers, cell_has_neighbor_where() and creep_fill_extrapolate(), a maskTerminusExtrapolated diagnostic output field, and updated field descriptions/console summaries. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add a regimeRatio Antarctic-wide map panel showing u / (SECONDS_PER_YEAR * Lambda * A * N^n), where u is the basal sliding speed computed from the input file's velocity fields. This ratio is the input speed divided by the Regularized Coulomb law's implied critical/transition velocity at each cell, so it directly indicates how close each cell is to the fully-plastic Coulomb regime (ratio >> 1) versus the power-law/Weertman-like regime (ratio << 1). Plotted on a log scale from 0.1 to 10 with a diverging colormap (RdBu_r) centered at 1.0 (the transition itself), masked to grounded ice. Also extends plot_maps() so panel vmin/vmax are honored for log-scale panels (previously ignored, auto-ranged from data only), which this new panel and any future fixed-range log panel needs. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
matthewhoffman
marked this pull request as draft
August 31, 2026 15:12
Previously the script assumed the input MALI friction law was always
classic Weertman (no effective-pressure term: Tau_b = mu * speed^qW).
Generalize this to a power law with an explicit effective pressure
N_source (Tau_b_source = mu * N_source * speed^qW), of which classic
Weertman is the special case N_source == 1 everywhere.
- Add --source-effective-pressure-type {constant,downs-johnson,
ocean-connection,transition,field}, independent of
--effective-pressure-type used for the RC law's own N. Default
"constant" with --source-effective-pressure=1.0 exactly reproduces
the old Weertman-only behavior. The other types compute N_source via
the same downs_johnson/ocean_connection/transition parameterizations
(with independently-configurable parameters:
--source-min-fraction-overburden, --source-pressure-length-scale,
--source-transition-h-ocean), or read a precomputed field
(--source-effective-pressure-field), giving a genuine Budd-type law.
The exponent on N_source is fixed at 1 (not user-configurable).
- Add compute_named_effective_pressure() to share the
downs-johnson/ocean-connection/transition dispatch between the RC
law's N and the source law's N_source.
- Rename tau_b_weertman -> tau_b_source throughout; refactor
fit_coulomb_C_fast_region()/solve_transition_velocity() to take
tau_b_source directly instead of deriving it internally from mu and
speed.
- Add an "effectivePressureSource" diagnostic output field when
--diagnostics is set, and new regularizedCoulomb_source* global
NetCDF attributes recording the source effective-pressure
configuration used.
- Update module docstring, CLI help text, printed summary, and map/
diagnostic titles/units accordingly (mu's physical units now depend
on --source-effective-pressure-type: kPa*(yr/m)^qW when constant,
plain (yr/m)^qW otherwise).
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…tion_law_conversion.py MALI extends the input mesh by one cell around its boundary when building the FEM mesh handed to (standalone) Albany, assigning a fixed minimum thickness to that extended ring, so the ice thickness/bed elevation Albany actually uses can differ from the plain MALI initial-condition file at/near the domain margin. friction_law_conversion.py now always overrides thickness/bed with values read from MALI's config_write_albany_ascii_mesh ascii output (thickness.ascii, bed_topography.ascii, mpas_cellID.ascii, required in --ascii-mesh-dir, default '.') for every cell covered by that output, before computing effective pressure, grounded/floating classification, grounding-line/terminus identification, diagnostics, and plots. Cells outside the ascii files' coverage keep their netCDF thickness/bedTopography values. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
friction_law_conversion.py's transect_fields (used for the per-transect line plots) previously omitted the Lambda and mu (input/output) panels entirely -- those were only plotted on the Antarctic-wide maps (map_fields). Add them unconditionally to transect_fields so every transect plot always shows Lambda and mu (input vs. output C), regardless of --method or other options. Also add a panel for the input source-law basal shear stress (tau_b_source = mu * N_source * speed^qW). Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…o transect plots friction_law_conversion.py's transect Lambda and mu (input/output) panels now zero values outside grounded ice, since these fields are only meaningful there, making the grounded-ice variation easier to read. Also add background shading to those same panels indicating the Regularized Coulomb regime ratio (u / (Lambda * A * N^n), the same diagnostic used by the regimeRatio map): a strong shade marks fully Coulomb (ratio >= 10) or fully power-law (ratio <= 0.1) regions, and a lighter shade marks runs trending toward either regime (ratio >= 2.0 or <= 0.5, respectively). Unlike the existing region-class/fit-region shading, this is drawn only on the Lambda/mu panels, not every panel. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add --effective-pressure-type=field, paired with a new --effective-pressure-input-field option, so the target Regularized Coulomb effective pressure N can be read directly from a precomputed field in the input file (e.g. an Albany-diagnosed effectivePressure field from a previous run) instead of being recomputed here from thickness/bed. The field is copied through unchanged to the output (--effective-pressure-field), written unconditionally when --effective-pressure-type=field since it is required (not just diagnostic) for Albany's own "Effective Pressure Type: Field" at runtime, and the suggested Albany YAML block is updated accordingly. Both --effective-pressure-input-field and the existing --source-effective-pressure-field now default to "effectivePressure", matching the default --effective-pressure-field output name. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add --rc-power-exponent (default: 1/3, matching the prior fixed project convention) so the Albany Regularized Coulomb law's Power Exponent qR can be overridden for diagnostic experiments -- e.g. setting qR=1 to make the slow-flowing branch linear in speed instead of a steep power law, to test whether excessive converted velocities are driven by sensitivity to effective-pressure (N) errors near the default qR. Threaded through solve_transition_velocity() (--method= transition-velocity) and the per-cell Lambda solve (--method= stress-match-fit), plus the printed summary, output attrs, and suggested Albany YAML block. Not recommended for production runs unless Albany's own YAML "Power Exponent" is changed to match. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Reconstructs the FE triangulation implied by the MPAS mesh (nodes = cell centers, elements = the 3 cells surrounding each MPAS vertex, via cellsOnVertex), linearly interpolates N, muFriction (C), bedRoughnessRC (Lambda), flow-rate A, and sliding speed from nodal (cell-center) values to quadrature points (--fe-quadrature-rule: centroid or edge-midpoint), and only then nonlinearly combines them into the Regularized Coulomb stress formula -- mirroring Albany's actual 'interpolate-then-combine' FEM evaluation order instead of this script's own 'combine-at-nodes' closed-form values. The resulting per-quadrature-point residual against the source law's stress is aggregated back onto MPAS cells and written as feStressResidual/feStressResidualMax/feStressResidualRelative/ feStressResidualRelativeMax (enabled via new --fe-diagnostics flag, requires --diagnostics), for direct comparison against observed velocity errors. Diagnostic only -- does not modify Lambda/C. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…ure-point residual Adds solve_fe_corrected_mu(), which corrects the closed-form muFriction (C) field so that Albany's actual FE-quadrature-point evaluation of the Regularized Coulomb stress (as reconstructed by --fe-diagnostics) better matches the source law's stress, for the fixed N/bedRoughnessRC/flow-rate/speed fields already computed. For fixed N/Lambda/A/speed, the RC stress is exactly linear in the quadrature-point-interpolated C, so this is solved as a genuine Tikhonov-regularized (ridge) linear least-squares problem (via scipy.sparse.linalg.lsqr) rather than an ad hoc iterative patch, anchored to the existing closed-form C as both the starting point and the regularization prior -- cells with little/no FE data support are left unchanged, while well-constrained cells are corrected towards better FE-level agreement. A --fe-correction- bound-factor safety net resets any corrected value that strays too far from (or becomes non-finite/non-positive relative to) its prior. Only C is corrected (not Lambda), since C dominates the RC stress in the fast-flowing/thin-margin (deep Coulomb-plateau) regime this is meant to help with. This makes muFriction spatially varying even for --method=stress-match-fit. New CLI: --fe-correct-mu, --fe-correction-regularization, --fe-correction-bound-factor. The --fe-diagnostics block is refactored to share the FE interpolation/residual logic with the new correction and to report before/after residual stats when both are used together. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
The previous implementation of solve_fe_corrected_mu() used a single global scalar regularization weight (sqrt(regularization) * median(k_qp)), but k_qp (essentially proportional to Albany's effective pressure N) spans orders of magnitude across a typical ice-sheet domain: near zero at thin margins and very large in the thick interior. A single global scale made the ridge term far too weak in the (large, dominant-by-area) thick interior relative to its own local data scale -- letting C drift with little justification there -- and far too strong at thin margins, exactly where the correction is actually needed. There was also no gating: every grounded cell touched by any valid FE triangle participated as a free unknown, even cells whose pre-existing FE residual was already negligible, so noise/bias could spread across the whole grounded interior. Together these explain the reported symptom of --fe-correct-mu making velocities faster almost everywhere instead of concentrating the fix in genuinely bad regions. This commit fixes both issues: 1. Local regularization scaling: the ridge weight for each cell is now scaled by that cell's own local characteristic k (computed by aggregating only the valid FE quadrature-point equations that touch it via aggregate_quadrature_to_cells()), falling back to the global median only for cells with no touching valid equations (a value that has no effect on their solution, which is pinned to C_prior in that case anyway). 2. Residual-based gating: cells whose pre-correction FE quadrature-point relative stress residual (computed once already, before correction, for reporting) is below a new --fe-correction-residual-threshold (default 0.05) are pinned near C_prior by multiplying their local ridge weight by a new --fe-correction-pin-factor (default 1e4). NaN residuals (cells without a meaningful diagnostic) are treated as below-threshold. solve_fe_corrected_mu() now takes residual_relative, residual_threshold, and pin_factor parameters, and returns an additional n_pinned count. The call site in main() passes the already-computed pre-correction relative residual and reports the number of cells pinned versus reset by the bound-factor safety net. Verified end-to-end on a synthetic Delaunay-triangulated dataset for both --method options and both --fe-quadrature-rule choices: the correction is now concentrated on a small fraction of cells (in one test, 13/400 vs. all 400 previously), with comparable or better overall residual reduction. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
The Albany ascii mesh override message compared overridden thickness/bedTopography values to the input netCDF using exact floating-point equality (!=). Since the ascii files round-trip these values through decimal text in km (double -> text -> double), nearly every cell picks up harmless last-bit floating-point noise (~1e-13 m) even when nothing physically changed, causing the reported 'differ' count to vastly overstate how many cells were actually affected by MALI's mesh extension/fixed-margin thickness. Verified against a real ~500k-cell Antarctic mesh: exact equality reported 352,996 of 354,506 ascii-covered cells as differing, while a 1 mm tolerance -- which is well below any meaningful geometric difference -- found only 12,416 actually changed, consistent with a one-cell-ish boundary/fixed-margin ring rather than 60% of the mesh. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
… solve The previous fix for --fe-correct-mu's regularization scaling kept every cell (pinned or not) as a free unknown in the sparse least-squares system, and instead multiplied pinned cells' ridge weight by a large pin_factor (1e4) to push their solved value toward C_prior. On a small test mesh (400 cells) this was fast enough to go unnoticed, but on a real ~1e6-triangle mesh it made scipy.sparse.linalg.lsqr converge extremely slowly (observed as an indefinite hang): the ridge diagonal now spanned many more orders of magnitude than the data-row weights (already highly variable via k_qp), which is a badly-scaled system for a Krylov solver like LSQR, and the unknown/row count never shrank (iter_lim=5*n_cells stayed at millions of iterations for a full-size mesh). This commit removes pinned cells from the unknown vector entirely, instead of down-weighting them: - Each data equation's pinned-node contributions (bary * k * C_prior) are folded into the right-hand side rather than kept as matrix entries, leaving only genuinely free nodes as unknowns. - Equations that reference zero free nodes (all 3 triangle vertices pinned) are dropped outright, since they cannot affect the solution. - The linear system solved by lsqr is now sized by the number of free cells/equations only, not the full mesh, so iter_lim (5x the unknown count) shrinks along with it. On a real-world-scale synthetic test (60k cells, ~120k triangles, ~92% of grounded cells pinned), the corrected system solves in ~2 seconds; the previous approach kept the full ~54k-unknown system with an ill-scaled diagonal for every run regardless of how many cells were pinned. Verified the small (400-cell) synthetic test produces the same residual-reduction numbers as before for both --method options and both --fe-quadrature-rule choices, confirming the reformulation is mathematically equivalent to the original weighted approach for cells that are pinned, while being far cheaper and better conditioned to solve. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…gress output Eliminating pinned cells from the unknown vector (previous commit) is enough when only a small fraction of cells are free, but on a real-world dataset where a much larger share of grounded cells sit above --fe-correction-residual-threshold, the reduced system can still have many tens of thousands of free unknowns. Even restricted to just the free cells, k_qp (and hence the per-cell data/ridge coefficients) can vary by orders of magnitude within that free set alone -- e.g. a nearly-grounding-line cell next to a merely somewhat-mismatched interior cell -- leaving a badly-scaled linear system that scipy.sparse.linalg.lsqr converges on very slowly, which is indistinguishable from a hang with no progress output and a 5x-unknown-count iteration budget. This commit: - Column-equilibrates the system before solving (each free cell's column is scaled to unit norm, and the solution is rescaled back afterwards) -- a standard, cheap Krylov-solver preconditioner that removes the magnitude disparity between differently-scaled free cells without changing the solution. - Caps the lsqr iteration limit at min(5 * n_free, 20000) instead of letting it scale unbounded with problem size; any cells lsqr hasn't resolved well within that budget still fall back to C_prior via the existing bound_factor safety net. - Prints the actual problem size (equations/free cells/pinned cells) before solving, and the iteration count/stop reason after, so a long-running solve is visible and diagnosable instead of appearing to hang silently. Verified on a synthetic 80k-cell mesh deliberately made rough enough that ~22% of grounded cells were left free (17,785 unknowns): lsqr now converges (istop=2) in 84 iterations. Also re-verified the 400-cell regression test reproduces identical residual-reduction numbers as before for all --method/--fe-quadrature-rule combinations, confirming the preconditioning does not change the solution. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…t script
Adds a standalone script implementing an alternative to
friction_law_conversion.py's --fe-correct-mu correction: rather than
approximating how Albany's FEM evaluates the Regularized Coulomb
(RC) friction law at quadrature points, this script refines C using
Albany's own actual simulated velocity from a completed MALI/Albany
run.
Workflow: run friction_law_conversion.py once to get an initial
Lambda/C conversion, run Albany with it, then run this script with
(1) the original MALI input (for mu, N_source, and the original
target velocity/Tau_b_target, held fixed across iterations), (2) the
previous conversion output (for C_old and Lambda), and (3) the
Albany run's output (for u_model, the actual simulated velocity) to
produce a refined C. Since the RC law is exactly linear in C at a
fixed N/Lambda/A, this reduces to a closed-form per-cell Picard
update:
C_new = C_old * (Tau_b_target / (C_old * k(u_model))) ** relaxation
with relaxation=1.0 (full step) by default, an optional
--bound-factor safety clip, and no update applied to cells with
undefined/invalid data (kept at C_old). The user reruns Albany with
the refined C and repeats as needed, guided by a printed
velocity-mismatch convergence summary and two diagnostic output
fields (muFrictionRefinementRatio, velocityMismatchRelative).
Verified against a synthetic three-file dataset: refined C exactly
reproduces the target stress at the (perturbed) model velocity
(relative error ~1e-16), and --relaxation/--bound-factor behave as
expected.
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…ard update Tau_b_target was being evaluated at the *original* input file's velocity while k(u) was evaluated at Albany's actual simulated velocity u_model. Since k(u) is monotonically increasing in u (saturating at N_albany as u -> infinity), this mismatch made the update structurally backwards: whenever u_model exceeded the original velocity (i.e. exactly the too-fast regions this tool is meant to fix), C_new came out *smaller* than C_old, reducing friction and making already-too-fast regions even faster -- and, symmetrically, needlessly perturbing regions where u_model was already close to the original velocity, since the two evaluation points still differed there too. Fix: evaluate Tau_b_target using the ORIGINAL input file's mu/ N_source/qW, but AT u_model (Albany's own simulated velocity), not at the original file's velocity. This exactly reduces to friction_law_conversion.py's own solve_transition_velocity() C formula, re-solved with u_model substituted for the original conversion's speed -- a standard Picard fixed point (recalibrate C so the RC law agrees with the trusted source law at whatever velocity Albany is actually producing right now). This is self-limiting (near-identity where u_model ~= target) and pushes C in the correct direction elsewhere. The original input file's velocity is now used only for the printed convergence/mismatch diagnostic, not for the update itself. Verified on synthetic data mimicking --method=transition-velocity (u_c held to a fixed u0 for all cells, as solve_transition_velocity() actually produces): too-fast synthetic margin cells now get C increased (ratio 1.12-1.28, more friction), while already-good cells are left nearly unchanged (ratio ~0.995-1.004), matching the intended behavior. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…b_mu.py Adds a printed self-consistency check: using this script's own N/ N_source/u_c formulas, recompute what C_old *should* have been at the ORIGINAL target velocity u_target (i.e. re-run solve_transition_velocity()'s C formula with u_target in place of u_model) and compare to the actual C_old read from previous_conversion_output. Restricted to cells where the printed velocity mismatch is already small (<5%), since a mismatch there should be ~0 to floating-point precision if every physics option passed to this script (--weertman-q, --glen-n, --rc-power-exponent, --effective-pressure-type/--source-effective- pressure-type and their thresholds, --flow-rate-type/--flow-rate, density/gravity values, --ascii-mesh-dir) matches what was used for the original friction_law_conversion.py run that produced previous_conversion_output. Motivated by a user report that already-good regions were still being perturbed (and made faster) after the previous sign fix. A mismatch confined to already-good cells is a strong signal of a physics-option mismatch between the two script invocations (e.g. a different --rc-power-exponent, effective-pressure parameters, or flow-rate), not a flaw in the Picard update itself -- verified on synthetic data: residual is ~1e-16 with matched options, and jumps to tens-of-percent when --rc-power-exponent is deliberately mismatched, with the reported ratio in already-good cells shifting well away from 1.0 in that case, exactly matching the reported symptom. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
friction_law_conversion.py now writes a <output>.config.json sidecar alongside its converted IC, capturing every physics option (weertman-q, effective-pressure-type, flow-rate, etc.) that refine_regularized_coulomb_mu.py needs to reproduce the same physics. refine_regularized_coulomb_mu.py auto-loads that config (derived from previous_conversion_output, or via --config) and fills in any physics option not given explicitly on its own command line, falling back to friction_law_conversion.py's own hardcoded defaults only if no config is found. Explicit CLI options still take precedence, with a warning printed if they differ from the config-file value. This eliminates the need to manually keep physics options in sync between the two scripts' invocations, which was a real, hard-to-spot source of error in practice (confirmed via the self-consistency diagnostic showing large spurious residuals from option mismatches). Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
refine_regularized_coulomb_mu.py's output previously had no <output>.config.json sidecar, so a second refinement iteration using that output as previous_conversion_output would find no config file, silently fall back to hardcoded physics defaults, and reintroduce the option-mismatch bug the config-file mechanism was built to prevent. Now each refine pass writes its own resolved config forward, so repeated iterations stay physics-consistent automatically. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Previously every grounded cell with valid inputs received a Picard update, even when the correction was negligible -- since an exact ratio of 1 is rare, essentially the entire grounded domain was perturbed by a small amount each pass. With Albany/MALI's globally coupled membrane-stress force balance, thousands of small, spatially incoherent perturbations can ripple into a widespread velocity change well beyond the handful of cells that actually needed correcting, which is consistent with a user report of many regions getting visibly worse (on a log scale) after a refinement pass despite "good" aggregate Picard-ratio statistics. Add --min-taub-mismatch (default 0.1): cells whose Tau_b stress mismatch at u_model, |Tau_b_target(u_model) / Tau_b_RC(u_model, C_old) - 1| (the unrelaxed Picard ratio minus 1), is below this threshold are now left at C_new = C_old exactly. This is gated on the stress ratio rather than relative velocity mismatch because velocity mismatch is a much noisier metric wherever the target velocity is small (small denominators inflate relative error even when the underlying stress mismatch is modest), and is only an indirect proxy for the local friction-law error since velocity also reflects nonlocal coupling to neighboring cells. Use --min-taub-mismatch 0 to recover the old behavior. Also updates the printed refinement summary to report how many eligible cells were left unchanged by the new threshold, and clarifies muFrictionRefinementRatio's attrs accordingly. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
At terrestrial-margin cells with a thin minimum ice thickness (e.g. 1 m) and zero water pressure, the overburden-based effective pressure N is tiny purely from the thin ice, driving Lambda = u0 / (SECONDS_PER_YEAR * A * N^n) to a huge value that is extremely sensitive to small drifts in N during a simulation. Add --thin-ice-thickness-threshold and --thin-ice-u0-factor (default 0.1) to shrink the transition velocity u0 at cells below the thickness threshold before the closed-form Lambda/C solve in solve_transition_velocity(), bounding Lambda's magnitude at these cells. solve_transition_velocity() now accepts a per-cell u0 array in addition to a scalar. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
matthewhoffman
force-pushed
the
landice/weertman-to-rc-conversion
branch
from
September 29, 2026 19:25
b03dcc6 to
43eea2a
Compare
This ensures N will be continuous where ocean connection extends to locations above sea level.
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This PR adds a script that converts an existing MALI initial condition using a Weertman basal friction law to regularized Coulomb friction law.
WIP