Tool for Weertman to Regularized Coulomb basal friction conversion - #752
Draft
matthewhoffman wants to merge 55 commits into
Draft
Tool for Weertman to Regularized Coulomb basal friction conversion#752matthewhoffman wants to merge 55 commits into
matthewhoffman wants to merge 55 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>
Add a new map panel 'impliedC' showing the per-cell implied Coulomb Friction Coefficient (tau_b_weertman / N_albany), masked to the fast-flowing/full-Coulomb fit region (fit_mask) used to fit the single global scalar C. Plotted on a log color scale (turbo colormap) since values can span orders of magnitude; the title also reports the fitted scalar C value for reference. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add fit_mask parameter to plot_transects(); along-transect runs falling inside the fast-flowing/full-Coulomb fit region (used to fit the single global scalar C) are now shaded orange on every panel, with a corresponding legend entry, so the fit region is visible alongside the existing grounded/floating/ice-free shading. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add a new 'implied C' panel to plot_transects() showing local_C (Weertman Tau_b / N), restricted to the Coulomb C-fit region (fit_mask) -- the line is NaN'd outside that region so it only draws where meaningful. Overlay a dashed reference line at the fitted global scalar C, restricted to the same along-transect fit-region spans, for direct visual comparison with the local values. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
New coulomb_N_sensitivity.py: a lightweight, tangential diagnostic that plots the effective pressure N implied by assuming pure Coulomb sliding (Tau_b = C * N) at a range of candidate scalar Coulomb Friction Coefficients C, restricted to the fast-flowing region (same grounded + speed > critical-velocity test friction_law_conversion.py uses to fit its own C), along the same glacier transects. Reuses friction_law_conversion.py's plot_transects() directly via import (including its existing multi-line-panel support, used here to overlay one N(C) curve per candidate C) rather than duplicating the transect-loading/projection/plotting machinery; only the handful of lines needed to read mu/thickness/bed/velocity and compute Tau_b_weertman = mu * speed^qW are duplicated. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
…ript Extend coulomb_N_sensitivity.py to also derive, for each candidate scalar C, the floatation fraction (Pw/Pice) and Shreve hydraulic potential implied by the pure-Coulomb N(C), using the same formulas as friction_law_conversion.py. Panel order (elevation, N, floatation fraction, hydropotential) now matches the main script's --plot-transects panel ordering; the elevation panel was already included automatically via plot_transects(). Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add optional per-panel ylim support to plot_transects() (4th tuple element in the fields dict, applied via ax.set_ylim()), and use it in coulomb_N_sensitivity.py to fix the floatation-fraction(C) panel to [0, 1]. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add 'Pine_Island' (matching the geometric_features landice/transect subdirectory name) to the default --transect-names list in both friction_law_conversion.py and coulomb_N_sensitivity.py. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Change plot_transects()'s optional per-panel ylim from a hard set_ylim() to a clamp on the autoscaled range: get the natural autoscaled y-limits, then bound them to at most the given (ymin, ymax). A hard set_ylim(0, 1) was squashing the floatation-fraction(C) panel's genuinely narrow-but-valid data (e.g. 0.95-1.0) into an invisible sliver hugging the top axis border, making the curves appear to not render at all. Clamping instead preserves the real autoscaled range while still preventing runaway/outlier-driven scaling outside [0, 1]. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add a 'Tau_b (Weertman)' panel to coulomb_N_sensitivity.py, showing the Weertman basal shear stress (mu * speed^qW) restricted to the same fast-flowing region as the N(C)/floatation-fraction(C)/ hydropotential(C) panels below it, so the shear stress driving those pure-Coulomb-implied quantities is visible alongside them. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
MALI's muFriction field's correct physical units, per a correction to MPAS-Tools' Registry.xml, are kPa * (yr/m)^qW (not Pa * (yr/m)^qW). mu * speed^qW therefore evaluates to kPa, not Pa. friction_law_conversion.py needed no functional change: it only ever uses mu * speed^qW (Tau_b_weertman) as a ratio against N_albany (N / ALBANY_EFFECTIVE_PRESSURE_PA_PER_UNIT, already kPa-scaled), so the ratio-based C-fit and Lambda-solve computations there are dimensionally correct regardless of mu's absolute unit -- added clarifying comments/docstrings there, and fixed the mislabeled muFriction map-panel colorbar units (was "Pa (m yr-1)^-q", now "kPa (m yr-1)^-q"). coulomb_N_sensitivity.py, however, uses Tau_b_weertman directly as an absolute quantity (N(C), floatation fraction, hydropotential), so it needed a real fix: multiply by ALBANY_EFFECTIVE_PRESSURE_PA_PER_UNIT (1000) to convert kPa -> Pa before using it. This was the source of the reported ~1000x-too-small N(C) values: N(C) is now directly comparable in magnitude to the ice overburden pressure, so the floatation-fraction(C) and hydropotential(C) panels now show physically meaningful variation across candidate C values instead of appearing nearly flat. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Pass aa=False to mosaic.polypcolor() in plot_maps(). Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
The marine (below-sea-level) term in downs_johnson_effective_pressure() had a sign flip introduced in an earlier commit: it computed max(rho_w*bed, 0.0) instead of Albany's max(-rho_w*bed, 0.0) (LandIce_BasalFrictionCoefficient_Def.hpp). Since bed is negative below sea level, the buggy version always evaluated to 0 there, causing N to collapse to the full inland overburden floor and making floatation fraction incorrectly read ~0 (and hydropotential correspondingly too negative) wherever bed < 0. Also removed a leftover debug print statement. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Albany reads its "Bed Roughness Field Name" field (bedRoughnessRC) verbatim, with no internal unit conversion, but its hardcoded scaling factor (scaling = secsInYr * pow(1000, n+1) in LandIce_BasalFrictionCoefficient_Def.hpp; comment: "//bedRoughness in km") is built assuming that stored value is already expressed in km. The exact-solve derivation for Lambda correctly produces a value in physical meters, but the script wrote that meters-based value directly to the output field without the required division by 1000, making the stored Lambda 1000x too large. This inflated basal resistance across the exact-solve (slow-flowing) region, consistent with the near-zero/no-slip velocities reported downstream in Albany runs. Added ALBANY_LAMBDA_METERS_PER_KM (1000.0) and applied it when computing Lambda; updated related comments/docstrings and the mislabeled "Pa (m yr-1)^-1/3" map units string (now "km"). Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
- Fixed the suggested-YAML "Effective Pressure Type" string, which incorrectly printed "Hydrostatic At Nodes" (not a valid Albany parameter value; Albany's validated string is "Hydrostatic Computed At Nodes", which would throw an exception at runtime otherwise). - Added a new --effective-pressure-type=ocean-connection option: same Albany "Hydrostatic Computed At Nodes" Effective Pressure Type as "downs-johnson", but with "Use Pressurized Bed Above Sea Level: false" -- i.e. full ocean (hydrostatic) water pressure wherever bed is below sea level, with no inland tapering/floor, and (since Albany does not even read them in this mode) no dependence on --min-fraction-overburden/--pressure-length-scale. - Added ocean_connection_effective_pressure(), matching Albany's formula with use_pressurized_bed=false (f_p == 0 identically): N = max(rho_i*g*H - max(-rho_w*g*bed, 0), 0). - Made --min-fraction-overburden/--pressure-length-scale optional when --effective-pressure-type=ocean-connection (not required/used). Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
The script was dividing the solved-for Lambda (in meters) by 1000
before writing it to the bedRoughnessRC field, under the assumption
that Albany reads bedRoughnessRC as-is and expects it pre-converted
to km.
That assumption is incorrect. Tracing the actual MALI/Albany
coupling:
- MALI's Registry.xml declares bedRoughnessRC with units="m".
- mode_forward/mpas_li_velocity_external.F passes the field straight
through to the C++ interface with no conversion.
- mode_forward/Interface_velocity_solver.cpp already divides it by
unit_length (1000) before Albany ever sees it:
bedRoughnessData[index] = bedRoughnessRC_F[iCell] / unit_length;
exactly as it does for bedTopography and thickness.
So MALI's own coupling layer performs the m -> km conversion
required by Albany's hardcoded "scaling" factor in
LandIce_BasalFrictionCoefficient_Def.hpp. The script's additional
/1000 was therefore double-converting the value, making the Lambda
Albany actually uses at runtime ~1,000,000x too small (1000x from
the script, times another 1000x from MALI's interface).
Fix: write Lambda directly in meters, with no further scaling, and
update the surrounding docstring/comments/NetCDF attrs to document
the verified unit pipeline instead of the previous incorrect
assumption.
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Lambda/bedRoughnessRC is now solved for and stored in meters (see previous commit removing the erroneous double m->km conversion), but the diagnostic map-plot colorbar label for this field was left as "km". Update it to "m" to match the actual stored/plotted units. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Previously the input Weertman muFriction field was dropped from the output NetCDF since it is not used directly by the Regularized Coulomb law. Per updated project plan, the fitted, spatially uniform RC coefficient C (Albany's "Mu") is now written into that same field (--mu-field, default muFriction) instead, so an Albany YAML using "Mu Type: Field" can read it directly from the IC file rather than needing a hardcoded "Mu: <value>" constant. The YAML template/output is intentionally left unchanged here; the YAML update will be made separately. 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
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
1 participant
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