Enforce exact fraction closure for fractional fields - #39
Open
matthewhoffman wants to merge 14 commits into
Open
Enforce exact fraction closure for fractional fields#39matthewhoffman wants to merge 14 commits into
matthewhoffman wants to merge 14 commits into
Conversation
Author
|
@jsta and @rfiorella , I found these changes were necessary to get ELM to run without error using a surface file created by dapper. That said, I'm not 100% sure if the changes here are more extensive than necessary, or possibly not comprehensive for certain ELM configurations (e.g. urban areas). I assigned you both as reviewers, but if you'd rather not, feel free to remove yourself. |
rfiorella
approved these changes
Jul 8, 2026
rfiorella
left a comment
rfiorella
left a comment
Collaborator
There was a problem hiding this comment.
@matthewhoffman seems reasonable to me from the ELM side, I'll let @jsta or Jon comment on how it fits in w/ the package python conventions
jsta
approved these changes
Jul 8, 2026
Collaborator
|
Even after merging this PR, I see:
|
Collaborator
|
maybe this is missing some topounit awareness |
Author
|
@jsta , I had only tested it with a basic surface file without topounits, so it's possible it needs more work. Feel free to try to fix that if you have a working failure case, or I can look at it later this week. |
Collaborator
|
I think something has to be done to modify how |
Collaborator
|
still doesn't have correct topounit pft indexing afaict: # reference topounits-enabled surface file
PCT_NAT_PFT (natpft, topounit, gridcell) float64 1kB 4.0 4.0 ... 0.0
# dapper produced topounits-enabled surface file with this PR
PCT_NAT_PFT (lsmlat, natpft, lsmlon) float64 136B 7.517 ... 0.0 |
Collaborator
|
I hacked on dapper and got past the PFT errors but now I'm seeing that dapper is reporting |
ELM's surfrdUtilsMod.F90 performs a strict equality check on the sum of
natural-patch weights (wt_nat_patch) and landunit weights (wt_lunit).
Dapper's zonal aggregation produces each fraction band independently via
area-weighted mean, so small per-band floating-point residuals accumulate
across all bands. Combined with float32 serialization, the written sums
could deviate from 1.0 by enough to trigger ELM's fatal error:
surfrd_veg_all ERROR: sum of wt_nat_patch not 1.0 ...
surfrd_get_data ERROR: sum of wt_lunit not 1.0 ...
Changes
-------
sfile.py
- Add _scale_to_target_sum(): scales a DataArray along a given dimension
so its sum matches a target, working in float64.
- Add _snap_partition_sum(): after scaling, assigns the exact residual to
the last band so that sum is exactly equal to the target in float64
arithmetic. This eliminates any remaining machine-epsilon drift.
- Add _normalize_surface_fraction_closure(): orchestrates normalization
for all fraction groups before write:
* Landunit scalars (PCT_NATVEG, PCT_CROP, PCT_WETLAND, PCT_LAKE,
PCT_GLACIER, PCT_URBAN) are proportionally scaled to 100 where the
existing total is already near 100 (gating prevents distorting
partial datasets), then an exact residual snap ensures the final
sum is exactly 100 in float64.
* PCT_NAT_PFT is normalized then snapped to exactly 100 (ELM treats
these as natural-patch weights that must sum to 1.0 after
dividing by 100).
* PCT_CFT and PCT_GLC_MEC are scaled and snapped to their parent
percent fields (PCT_CROP and PCT_GLACIER respectively).
* FSURF and FGRD (irrigation split) are normalized to sum to 1.
- Closure-critical fraction variables are now encoded as float64 in the
written NetCDF file instead of float32 to prevent quantization
reintroducing the drift that was just corrected in memory. The
closure-critical set is:
PCT_NAT_PFT, PCT_CFT, PCT_GLC_MEC, PCT_NATVEG, PCT_CROP,
PCT_GLACIER, PCT_WETLAND, PCT_LAKE, PCT_URBAN.
- _normalize_surface_fraction_closure() is called inside write_surface_nc()
so the fix applies to every code path that writes a surface file,
regardless of sampling method.
landuse.py
- Add _scale_to_target_sum() and _normalize_landuse_fraction_closure()
mirroring the surface logic, covering the same fraction groups.
- _normalize_landuse_fraction_closure() is called immediately before
_write_nc() in the zonal sampling path of sample_landuse_timeseries()
so exported landuse files satisfy the same closure constraints.
- Note: the landuse write path does not yet have the float64 encoding
override; it relies on normalization-in-memory being sufficient. If
ELM-level landuse checks exhibit a similar issue the encoding should
be hardened there as well.
validate.py
- V-105 consistency check now verifies sum(PCT_NAT_PFT) against 100
(previously compared against PCT_NATVEG, which is the wrong target
for ELM's natural-patch weight convention).
- Tightened tolerance from 1e-3 to 1e-6 (max absolute difference,
not mean).
- Diagnostic message now reports both mean and max absolute residual.
test_fraction_closure.py (new file)
- test_write_surface_nc_enforces_fraction_closure: verifies PCT_NAT_PFT
closes to 100, PCT_CFT closes to PCT_CROP, landunit total closes to
100, and FSURF+FGRD closes to 1 after a write/read roundtrip.
- test_landuse_normalizer_enforces_fraction_closure: unit test for
_normalize_landuse_fraction_closure() in isolation.
- test_write_surface_nc_natural_patch_weights_near_exact_one: regression
for the first ELM failure; uses a 17-PFT distribution with deliberate
sub-ULP perturbations and asserts sum(PCT_NAT_PFT)/100 == 1.0 within
1e-12 after write/read.
- test_write_surface_nc_landunit_weights_near_exact_one: regression for
the second ELM failure; verifies that the sum of all landunit percent
fields divided by 100 equals 1.0 within 1e-12 after write/read.
This commit refactors the previous commit by introducing a shared module, fraction_closure.py, that handles closure operations and is called from multiple places (surface and landuse variables). It also moves the closure requirements into the schema/specs, rather than having them ad hoc. It also extends the fractional closure calculations to topounits. Tests are added for the topounit closure.
…imal fraction - Update fraction_closure to handle only TopounitFracArea with decimal target (1.0) - Remove PCT_TOPUNIT from closure_critical_variables set - Update validate.py to check TopounitFracArea closure (target=1.0, not 100) - Consolidate separate validation checks for TopounitFracArea - Update test to verify TopounitFracArea closes to 1.0 (not 100) This fixes the unit mismatch where ELM expects TopounitFracArea as a decimal fraction (0-1) but dapper was producing percent units (0-100).
After add_topounits_from_domain(), all SURFACE_VAR_SPECS variables with 'topounit' in their dims are now expanded to include the topounit dimension if they lack it. This uses uniform broadcast (all topounits inherit the parent cell's distribution). Fixes the issue where PCT_NAT_PFT and other topounit-dependent variables were written without a topounit dimension, causing ELM to fail with dimension mismatch errors. Also updates add_topounits_from_domain() to use TopounitFracArea instead of PCT_TOPUNIT. Includes comprehensive test to verify: - PCT_NAT_PFT gains topounit dimension with correct ordering - PFT fractions still sum to 100 per topounit - all topounit-indexed vars are properly expanded
matthewhoffman
force-pushed
the
roundoff-fix-fractional-vars
branch
from
July 16, 2026 20:36
6ea185b to
da03322
Compare
Author
|
@jsta , I force-pushed this branch after adding some commits and rebasing on main (which replaces the need for your merge of main into this branch). There are two new commits that I believe should fix the most recent issues you've identified. I ran a modified version of your domain/surface generation script that you shared with me on Teams, and I confirmed that the resulting surf file has the required dimensions: Try the current head of this branch with your workflow, and let me know if you run into any other issues. |
…ce variables table
- Add topounit dimension to phosphorus pools (APATITE_P, LABILE_P, OCCLUDED_P, SECONDARY_P) with corrected units (gP/m2) from ncdump metadata - Add topounit dimension to soil properties (ORGANIC, PCT_SAND, PCT_CLAY, PCT_GRVL) with pattern: nlevsoi,topounit,lsmlat,lsmlon - Reorder class dimensions before topounit for consistency: PCT_URBAN, PCT_CFT, NFERT, PFERT, PCT_NAT_PFT, PCT_GLC_MEC, TOPO_GLC_MEC (class dimension now precedes topounit, not vice versa) - Add lsmpft and topounit dimensions to monthly vegetation variables (MONTHLY_LAI, MONTHLY_SAI, MONTHLY_HEIGHT_TOP, MONTHLY_HEIGHT_BOT) with dimension order: time,lsmpft,topounit,lsmlat,lsmlon - Update MaxTopounitElv and TOPO2 documentation to clarify topounit relationships - Add topounits context tag to all affected variables Ensures variables align with actual ELM surface file structure (topounit-based discretization).
Add 8 new urban albedo variables (diffuse and direct components) for impervious/pervious roads, roofs, and walls across spectral bands. These variables control solar radiation absorption in urban land model. Variables added: - ALB_IMPROAD_DIF, ALB_IMPROAD_DIR - ALB_PERROAD_DIF, ALB_PERROAD_DIR - ALB_ROOF_DIF, ALB_ROOF_DIR - ALB_WALL_DIF, ALB_WALL_DIR All have dimensions: numrad,numurbl,topounit,lsmlat,lsmlon
Add 16 new variables defining urban material properties and geometry: Thermal properties: - CV_* (3): Volumetric heat capacity for impervious road, roof, wall - TK_* (3): Thermal conductivity for impervious road, roof, wall - EM_* (4): Emissivity for impervious/pervious road, roof, wall Structural/geometric properties: - HT_ROOF, THICK_ROOF, THICK_WALL: Roof/wall height and thickness - NLEV_IMPROAD: Number of impervious road layers - T_BUILDING_MAX, T_BUILDING_MIN: Temperature constraints Urban geometry: - CANYON_HWR: Canyon height-to-width ratio - WIND_HGT_CANYON: Wind reference height in canyon - WTLUNIT_ROOF, WTROAD_PERV: Urban landunit weight fractions
Add 8 new variables controlling water table dynamics, infiltration, and surface runoff generation in the Variable Infiltration Capacity (VIC) model: VIC ARNO parameters: - Ds: Fractional saturated area infiltration parameter - Dsmax: Maximum infiltration rate (mm/day) - binfl: Infiltration nonlinearity parameter Saturation/inundation parameters: - F0: Maximum fractional inundated area - FMAX: Maximum fractional saturated area - ZWT0: Water table decay factor Runoff parameters: - P3: Surface runoff lag coefficient - LAKEDEPTH: Average water depth for lake landunit All have dimensions: topounit,lsmlat,lsmlon
Add 4 new variables characterizing soil properties and topographic heterogeneity: Soil classification: - SOIL_COLOR: Categorical soil color affecting albedo - SOIL_ORDER: Categorical soil order for pedogenic classification Topographic properties: - SLP_P10: Slope quantile distribution (nlevslp,topounit,lsmlat,lsmlon) captures topographic variability at multiple percentiles - aveDTB: Average depth to bedrock (topounit,lsmlat,lsmlon) critical for groundwater and subsurface hydrology Enhances subgrid topographic and pedologic characterization for improved process representation in topographic discretization.
Add 6 new variables defining basal emission factors for isoprene (EF1) across plant functional types, used in volatile organic compound (VOC) biogeochemistry model: Vegetation-specific emission factors: - EF1_BTR: Broadleaf tree isoprene emission factor - EF1_CRP: Crop isoprene emission factor - EF1_FDT: Deciduous forest isoprene emission factor - EF1_FET: Evergreen forest isoprene emission factor - EF1_GRS: Grass isoprene emission factor - EF1_SHR: Shrub isoprene emission factor All have dimensions: topounit,lsmlat,lsmlon Enables subgrid-resolved VOC and atmospheric chemistry calculations.
Add 8 new variables spanning hydrologic, disturbance, and erosion processes: Hydrologic VIC parameter: - Ws: VIC maximum soil moisture storage parameter Agricultural/disturbance indicators: - abm: Agricultural fire peak month - peatf: Peatland fraction - gdp: Economic activity indicator Hillslope erosion parameters: - parEro_c1: Rainfall-driven erosion coefficient - parEro_c2: Runoff-driven erosion coefficient - parEro_c3: Sediment transport capacity parameter All have dimensions: topounit,lsmlat,lsmlon Completes topounit discretization for terrestrial and disturbance processes.
Enable flexible filename patterns in export_met() by supporting {var}
placeholder in the filename parameter. This allows generating E3SM PR #25
compatible filenames (ERA5_{var}_1950-2025_z01.nc) directly without
post-processing file renaming.
Changes:
- Update _nc_filename() to detect {var} placeholder and use .format()
- Maintain backward compatibility with simple prefix format
- Update docstrings in Domain.export_met() and Exporter.run()
Example usage:
domain.export_met(..., filename='ERA5_{var}_1950-2025_z01')
# Generates: ERA5_TBOT_1950-2025_z01.nc, ERA5_FSDS_1950-2025_z01.nc, etc.
Legacy format still works:
domain.export_met(..., filename='ERA5')
# Generates: ERA5_TBOT.nc, ERA5_FSDS.nc, etc.
Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
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
3 participants
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.
ELM's surfrdUtilsMod.F90 performs a strict equality check on the sum of natural-patch weights (wt_nat_patch) and landunit weights (wt_lunit). Dapper's zonal aggregation produces each fraction band independently via area-weighted mean, so small per-band floating-point residuals accumulate across all bands. Combined with float32 serialization, the written sums could deviate from 1.0 by enough to trigger ELM's fatal error:
This PR fixes this be enforcing roundoff level closure of fields that ELM expects to add to 1.0. This eliminates the errors in ELM.