Surface energy balance on immersed-boundary building faces (erf.ibseb), with immersed-boundary fixes for the forcing and the MRF/YSUNew schemes - #3960
Conversation
…ance Eight phases on faces of resolved buildings: face storage, shortwave with ray-cast shadowing, longwave through sky/ground/building view fractions, sensible heat through the existing immersed-forcing wall model, slab conduction, the prognostic balance, and two canonical cases. Radiation reaches the balance through a provider interface (prescribed, or the two-stream column when that branch is merged), so the branch stays based on development and carries only its own commits. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…balance erf.ibseb.enable builds, on every level, a compact list of the wall faces of the resolved buildings: a face between a fluid cell (blanking < 0.5) and a solid cell (>= 0.5), stored once on the rank that owns the fluid cell as device arrays (cell, direction, side, building id, material id, area, skin temperature, slab temperatures, view fractions, fluxes), contiguous per local box. Building ids label the solid columns of the blanking. Nothing evolves yet. Output through cell-centred fields on demand: ibseb_nfaces and ibseb_tskin in the plotfile, IBSEBState in the checkpoint, refilled into the rebuilt list on restart so restarts do not depend on the rank count. A [IBSEB] summary line and a per-building CSV every erf.ibseb.csv_int steps. Regtest Exec/CanonicalTests/SEB/Storage: the ImmersedForcingTest skyscraper; face counts per direction checked against the plotfile mask, identical on one and four ranks (2056 faces, 51400 m2), checkpoint at step 2 and restart to step 4 reproducing the CSV row. Sphinx: inputs table and a theory page. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
plot_storage.py draws, with yt, a horizontal slice of the mask and of the faces-per-cell field through the building and vertical slices of the mean skin temperature through its centroid, with the mask outline, in absolute metres (ERF plotfiles carry no unit metadata) and with the y-normal slice swapped to (x, z). The four sources of the balance get file, class and function documentation in the Doxygen style of the existing code, spelling out the face convention (fluid cell, direction, side, outward normal), the ownership and per-fab layout of the list, the building labelling, the six-slot checkpoint field, and which phase fills which array. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
erf.ibseb.debug now prints, in the fire module's [FIRE DEBUG] style: the inputs and the global face counts at build, each building's footprint columns and index box, every rank's face count with its fab ranges and device memory, a summary with per-building rows every step, and the checkpoint save and load of the face state. The regtests move to Exec/CanonicalTests/SEB/Phase<n>_<name>, starting with Phase1_Storage, so the phase order is visible in the tree; the plan, the theory page and the README follow. The y-normal slice of plot_storage.py gets its labels after yt's axis swap, so they read x and z. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
A prescribed radiation provider (erf.ibseb.radiation = prescribed) gives the direct-normal and horizontal diffuse irradiances and the sun vector, either fixed (sun_mode = fixed, for analytic tests) or from the site and time (sun_mode = solar: Spencer solar geometry, Bird direct beam, a diffuse share of the attenuated beam), through a single copied header so the branch has no dependency on ERF-Radiation. Per face and per step: direct = DNI max(0, n.s) unless the ray from the face centre toward the sun hits a building, decided by a 2D walk over the column tops of the level (a small array replicated on every rank; the ray only rises, so a column blocks it when the entry height is below its top); diffuse = f_sky D + f_ground albedo_ground (DNI cos z + D) with the placeholder view fractions (roof: sky 1; wall: sky and ground 0.5 each); absorbed = (1 - albedo) times the sum. Called at the start of every step. New plot fields ibseb_sw_abs and ibseb_shadow, shortwave columns in the per-building CSV, a per-rank face dump (erf.ibseb.dump_faces_file), and the sun in the [IBSEB DEBUG] output. Regtest Exec/CanonicalTests/SEB/Phase2_Shortwave: a short box 40 m east of a tall one; the shadow flag of every face matches an independent Python ray cast, the fluxes match the formulas on every face, the tall core roof and west wall are unshadowed with the exact incidence, the shadow on the short box's core west wall stops at H - gap tan(elevation), one and four ranks agree, and the solar mode gives the solstice-noon zenith at Boulder. The embedded-boundary reader steps each edge over one cell, which the test reads from the dump rather than assuming. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
… balance Once at initialisation every face samples a cosine-weighted hemisphere around its outward normal (erf.ibseb.view_n_az x view_n_el rays, stratified so the counts are view factors) through the column walk of the shadow test made direction-aware: a rising ray is blocked by a column whose top is above its entry height, a falling ray by a column it descends into, and reaches the ground otherwise. The fractions of rays ending on the sky, the ground and a building replace the phase 2 placeholders in the diffuse shortwave. Longwave per step: LW_in = f_sky LW_sky + f_ground eps_g sigma T_g^4 + f_bldg sigma T_skin^4, the sky term fixed (erf.ibseb.lw_down) or gray (sky_emissivity sigma T_air^4 with the air temperature of the face's fluid cell, through the equation of state), the building term the isothermal- surroundings approximation; LW_net = eps (LW_in - sigma T_skin^4). No face-to-face view factors, no radiosity. New plot fields ibseb_lw_net and ibseb_f_sky, longwave in the CSV and the face dump, the sampling summary in the debug output. Regtest Exec/CanonicalTests/SEB/Phase3_Longwave on the two-box deck: the three fractions of every face sum to one and equal an independent Python hemisphere sampling; the tall core roof sees only sky, no roof sees the ground, the tall west wall sees exactly half sky with its rim ledge filling a quarter of the view just above it and little 50 m up, the short box's core west wall sees the tall box; the longwave formulas hold on every face for both sky modes; one and four ranks agree. Phase 2 still passes with the sampled fractions. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Per face and per step: the tangential wind of the fluid cell at half a cell from the wall gives u* with the roughness z0_wall, the skin-to-air potential-temperature difference (the skin temperature converted with the cell's Exner function) gives theta* with z0h_wall, and H = rho c_p u* theta*, positive out of the face. Walls are neutral; erf.ibseb.stability_correction applies the surface layer's similarity functions on roofs. The latent flux stays zero with its argument in place. The flux enters the atmosphere as an explicit source: every face deposits H A / (c_p V Pi) into the rho-theta equation of its fluid cell, added after make_sources at every slow stage with atomic adds; erf.ibseb.couple_heat = false diagnoses without applying. The immersed forcing's own surface- temperature inputs are refused when the balance is on, since it now owns the temperature condition at the buildings. New plot field ibseb_H, the flux in the CSV, the dump and the debug output, and H_total_W per step. Regtest Exec/CanonicalTests/SEB/Phase4_Sensible: a 40 m cube held at 320 K in an 8 m/s wind at 300 K. u* and H match the formulas on every face of the dump, H is positive and largest on the windward wall, one and four ranks agree, and with the flux applied the extra internal energy of the air against the diagnostic run at the same step matches the summed face flux (the domain is a rigid closed box, so it is c_v, not c_p, that closes). A mass-inflow, pressure-outflow variant (as the Askervein canonical) runs with a wake warmer than the inflow. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Every face carries a slab of erf.ibseb.n_slab_layers uniform layers between its skin and the building interior at T_interior, solved implicitly with the Thomas algorithm in the form of the SLUCM branch's slab solver, with the skin temperature as the top boundary instead of a flux (ERF_IBSEBSlab.H); G = 2k/dz (T_skin - T_0) is the conduction into the wall, positive inward. Unconditionally stable, up to 256 layers, advanced every step. Materials: erf.ibseb.material_file names a CSV in the SLUCM schema (rank-0 read, broadcast), material_default and material_by_building assign them by building id; the conductivity, heat capacity, thickness, albedo and emissivity become per-face arrays, and the shortwave and longwave now use the per-face optical properties. Without a file the uniform k_therm, rho_cp, thickness, albedo and emissivity apply. New plot field ibseb_G, G in the summary and the dump with the material columns and the top and bottom slab layers. Regtest Exec/CanonicalTests/SEB/Phase5_Ground: a 200 mm slab in 1 mm layers follows the semi-infinite erfc solution for a boundary step (bottom layer 300.925 K against 300.921 K expected at 50 s) with no flux through the skin; a light 20 mm slab reaches G = k dT / L exactly with a linear profile; two buildings carry the concrete and timber of the CSV; the slab and its flux restart exactly through a checkpoint. The atmosphere-derived columns differ by about 1e-4 after a restart, which is the immersed-forcing atmosphere of development, not the balance, and is noted in the plan. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Every face now solves SW_abs + eps Q_ext + LW_net - H - LE - G = 0 for its skin temperature at the end of each step (ERF_IBSEBBalance.H, lifted from the SLUCM facet solver). The conduction term is the flux the implicit slab step will take, G = a T - b from two trial slab steps, so the balance and the slab agree to rounding and the slab energy changes by exactly dt (G - G_bottom). The wall term of the incoming longwave folds into the emission as (1 - f_bldg); the wall-function coefficient is frozen at the wind of the step. Q_ext is an incident external flux absorbed with the emissivity, the hook for a fire's radiation (erf.ibseb.Q_ext_uniform for tests). The bounds and the step cap are inputs since fire exposure exceeds the urban canopy model's 380 K. erf.ibseb.prognostic defaults to true; the phase 2-5 decks pin it false to keep checking each term on its own. Fixes the solar azimuth of phase 2, which was mirrored east-west (sign of the sine term); the noon check could not see it, the sunrise deck did. Regtest SEB/Phase6_Prognostic on a cube: residual below 1e-8 W/m2 on every face at every step, every stored flux consistent with the skin temperature, slab energy exact per step, closure over the run within the summed residual, an independent Python model (own Newton, dense implicit slab) driven by the per-step face dumps reproducing the skin temperature to 1e-9 K, an external flux run past 380 K, a checkpoint restart, and the sun rising over the cube at Boulder with the east wall warming first. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…l_cells) A building from a nodal height map has sliver cells with 1-20 percent solid at the corners of the reader's one-cell ramp. The drag on a partial cell fades with its fraction, but the wall law was applied at full rate to any cell above 0.005 solid whose normal neighbour was exactly zero, so the top cell of a sliver column was forced like a wall on top of an almost free cell. In a neutral 3 m/s run on a 40 m cube at 10 m and 0.5 s a vertical two-cell checkerboard in theta grew there over hours (10 K, 5 m/s spikes) until the floating-point trap fired after 1.8-2.5 h, with the surface energy balance off, a no-slip ground, the forcing outside the substeps and a coarser small-cell threshold as well. With erf.if_snap_partial_cells the six forcing functions place the wall law and the surface temperature, flux and Obukhov conditions on cells at least half solid whose normal neighbour is less than half solid, the rule the face balance uses; cells below half solid keep the fraction-weighted drag only. Faces passing several wall tests at once have their relaxations averaged so the explicit substep forcing stays within its limit. Default false: with the switch off every result is bit-identical (checked on the phase 6 closure deck). Regtest Exec/RegTests/ImmersedForcingTest/PartialCells: the height-map cube for 2.64 h with the switch on stays neutral to 1e-3 K with a wake; --reproduce runs the original selection, which traps. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Exec/CanonicalTests/SEB/Phase7_IsolatedBuilding runs a 40 m concrete cube at Boulder on the June solstice from midnight for 24 hours in a light westerly with the prescribed clear-sky provider and a gray sky. The cube is an exact box (eb2.geometry = box) so the immersed forcing has no partial cells; the per-building CSV gains the sun (zenith, azimuth, DNI, diffuse) so the day can be plotted without debug output. The day: every face radiates below the air at night, the roof coldest at dawn; the east wall is the first face to rise above the air after sunrise, forty minutes before the roof; the east wall peaks at 09:39, the roof at 13:55 at 340 K two hours after the sun, the west wall at 16:45; the south wall is 9 K warmer than the north at 13:00; conduction turns around at 18:10 and the slab's stored 16.5 MJ/m2 keeps the roof above the air past 23:00. The checker asserts that sequence, the balance residual over the day (4e-8 W/m2), the absorbed shortwave on the roof against the clear-sky formulas integrated independently in Python (24.882 MJ/m2 both) and the slab energy against the integrated conduction (0.01 %). The plot script draws the skin temperature by orientation with the air, the roof budget, the sun path, a slab Hovmoller and yt slices. Records under findings that the sensible flux off the 340 K roof is only 10-30 W/m2 with the neutral wall function, which the next PR addresses with the stability functions and a convective velocity scale, both behind switches; and that the restart non-exactness of the immersed forcing persists after erf-model erf-model#3956. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
The neutral log law on the tangential wind cannot shed heat from a hot face in calm air, which is why the day canonical's roof reached 340 K. Two switches, both off by default so every earlier result is unchanged: erf.ibseb.convective_velocity = deardorff puts a convective velocity scale into the wind the wall function sees, U_eff^2 = U_tan^2 + (beta w*)^2 (Beljaars), with w* = (g/theta H/(rho c_p) depth)^(1/3) from the previous step's flux out of the face. The depth is the mixed layer above a roof (z_i - z_face, floored at the building height) and the building height for a wall; z_i follows the day by the bulk Richardson diagnostic on the horizontal-mean profile (erf.ibseb.z_i_mode = bulk_ri), from the surface layer's pblh at the column (pblh), or fixed (z_i). erf.ibseb.stability_correction now iterates Dyer's similarity functions on the roofs to convergence on the face's own Obukhov length, under- relaxed as the surface layer's iteration is in erf-model erf-model#3486, seeded from the ground surface layer's 2D field at the column (erf.ibseb.obukhov_seed). Walls stay on the log law. The face's L stays its own because a roof in a separation zone or a sunlit wall can be in the opposite regime from the ground. The face dump gains w_star, Olen, z_i and h_bld; the summary gains w_star_max. Regtest SEB/Phase8_WallFunction: the cube in calm air under a strong sun; the neutral law sheds 0.2 W/m2 from a 342 K roof and the scale 270 W/m2; w*, the depth, u* and H follow the formulas to 1e-9; the roofs' L is negative and consistent with u* and theta* and the corrected log law to 1e-7; the bulk Richardson depth on a capped sounding is the first cell above the inversion. The default path is bit-identical to the phase 6 reference. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…dary aware erf.if_snap_partial_cells now reads the cell and face blanking snapped to solid or fluid at half, so a height-map building becomes the staircase of whole cells an exact box is (wall faces fully blanked, interiors damped, nothing on the sliver cells), and it selects the point-implicit drag: the staircase puts full-strength forcing on every wall face, whose explicit rate at rims and corners drives the density negative within seconds on a thin slab (the first version, which only moved the selection threshold, survived the cube for 2.6 h but killed a 3-cell slab in two minutes). With the switch off every result is bit-identical. The PartialCells regtest runs the height-map cube 2.6 h clean. erf.pbl_ib_aware (MRF and YSUNew only; MYNN, MYJ, YSU and SHOC untouched) makes each column's surface the first fluid cell above the immersed solid: the bulk Richardson heights, the boundary-layer depth and the K profile are measured from it, the diffusivities vanish inside the solid, and the surface scales of a column with solid cells are a neutral log law at its top with erf.pbl_ib_z0, since the ground surface layer evaluates u*, theta* and L on cells inside the building. Without immersed cells the results are bit-identical; not supported with terrain-fitted coordinates. Regtest ImmersedForcingTest/PBL_IBAware: with the switch the schemes are finite everywhere and zero inside the cube, and the profile over the roof has the shape of the ground's against the local height (MRF peaks 15 m above each); without it MRF fills the domain with NaN and YSUNew drives the density negative at its second step, reported not asserted. The balance prints a cost line (per-step time of the slowest rank, faces per rank, initialisation time) for estimating city-scale cases. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Exec/CanonicalTests/SEB/Phase8_BuildingSet: four buildings from a nodal height map on a 480 m periodic domain at 10 m (a 60 m concrete slab, a 40 m brick cube east of it, two 20 m timber blocks), three materials by building, a 3 m/s westerly, Boulder on the solstice from 05:00 for six hours with the prescribed provider, the convective velocity scale and the stability functions on, and the immersed forcing snapped to whole cells. Checked: the four buildings and their materials; the residual (3e-8 W/m2 over 1444 rows); the mutual shadowing (at sunrise the cube's shadow lies on the slab's east wall, 17 percent of its faces, clearing to 6 percent; the 20 m blocks free of shadow by late morning); the building view fractions of the facing walls (28 and 48 percent) against the far block's (18); the timber roofs ending at 333 K against the concrete slab's 320 K; w* on every sunlit face and 93 percent of the roofs unstable; the cost line (0.4 ms per step for 157 faces per rank). The plot script draws the per-building temperatures, the shadow fractions, a face map from above and a yt slice. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…3 through a restart The report row gained net-longwave and sensible-flux columns after phase 1 was written; they read the air temperature and wind, which the immersed forcing does not restart bit-for-bit (2e-5 relative, see Phase5_Ground). The geometry, skin and slab columns stay exact. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Single-precision builds (Apple Clang and SYCL, Mesh SINGLE) refused two brace-initialised arrays in ERF_IBFaceSet.cpp that narrowed double to float: the face-centre coordinates and the cell-centred velocity. Both are now built from Real expressions. The IB-aware YSUNew kernels read the surface height into a local that seven of them never used; those locals are removed (no change in behaviour). codespell flagged the Fourier-number variable "Fo" in the slab solver and its Python re-implementation, and a "tha" alias in the phase 6 check; renamed to Fourier and th_air. Phase 6 and the PBL IB-awareness regtests pass on the rebuilt binary. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
plot_buildingset.py gains three figures from the existing run output: every face coloured by the incident direct beam at 06:00, 08:00 and 10:00 solar time (a view from the south-east and a top-down map with roofs as squares and wall columns as bars), and the horizontal wind in the first two cells above the ground with a vertical slice of u through the slab and the cube from the plotfile. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Conflicts in the MRF and YSUNew kernels, where development's PBLH work (erf-model#3486) meets this branch's immersed-boundary awareness (local surface ksurf and effective surface scales us_eff/ts_eff/qs_eff/ol_eff): - MRF passes 2 and 4: keep development's convective/shear w* blend, evaluated on the column's local surface (theta_v at ksurf) with the effective u* and theta*, as the branch does for every other surface term. - MRF K-profile: development's ghost-cell guard and notes, with the branch's k < ksurf (no diffusivity inside a solid) condition. - YSUNew K-profile: development's rho guard, plus the branch's ksurf and effective Obukhov length. - VH96 shear correction (new on development, off by default): use the effective u*. On columns without immersed cells it equals u*, so runs without erf.pbl_ib_aware are unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ents - M_PI is undefined on MSVC (C2065 in ERF_IBSEBSolar.H, which failed the WIN32 jobs of the previous head); the solar geometry and the reports use PI from ERF_Constants.H. - The reduced face counts were `long`, which is 32-bit on Windows and has no ReduceLongSum overload; they are amrex::Long. - init_ibseb() aborts on a non-uniform vertical grid (the face areas, heights and the ray cast take the level's constant cell sizes) and on erf.regrid_int > 0 (the face list is built once from the blanking); ibseb_report() bounds its level loop by the face sets it holds. - Inputs.rst states the grid and regrid requirements. - The comments describe the routines instead of the development phases. Verified: a Debug build with MPI, no FFT, all warnings, assertions and bound checking is warning-free and runs the 471 unit tests; the Phase 6 closure, external-flux and restart checks and the PBL_IBAware regtest pass on 2 ranks; the six touched translation units compile in single precision. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
The face detection takes every solid column of the blanking for a building, so with erf.terrain_type = ImmersedForcing the terrain would be put under the balance as well, which the documentation says it is not. init_ibseb() now requires erf.buildings_type = ImmersedForcing and aborts on terrain by immersed forcing; Inputs.rst says so. IBSEB_Cube runs the Phase6_Prognostic cube (32x32x16, 40 steps, fixed sun, prognostic skin, slab, heat flux into the air) and compares the plotfile with its face diagnostics against a gold; nothing in CI exercised the balance before. Passes in a Debug build with MPI and no FFT in 57 s on 2 ranks; a run with a different albedo fails the comparison. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Resolves the MRF conflict with erf-model#3972 (PBL passes on the tile work box): the immersed-boundary work arrays of MRF and YSUNew now cover the same grown tile box as the PBLH passes, and the blanking's halo is part of the smoothing halo check. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Development merged erf-model#3972 and erf-model#3992, which put the PBL passes on a grown tile work box; the IB-aware per-column arrays now live on that box too (merge commit), and this registers PBL_IBAware_MRF_Tiling and PBL_IBAware_YSUNew_Tiling in the tiling-parity harness: the 40 m cube under a heated surface layer with a capped mixed layer, so the PBL height differs between the roof columns and open ground. Tiled and untiled runs agree bitwise for both schemes. The solar header takes PI from ERF_NumericalConstants.H, as the constants split of erf-model#3994 asks. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…at they test Source/ImmersedBoundarySEB/IBSEB_DEVELOPMENT.md is a development log; it is removed from the branch and listed in .gitignore so a local copy stays private, and nothing in the code or the documentation refers to it. The cases under Exec/CanonicalTests/SEB are renamed after what they test (FaceStorage, Shortwave, Longwave, SensibleHeat, SlabConduction, PrognosticSkin, IsolatedBuilding, WallFunction, BuildingSet) and their READMEs, decks, scripts and the theory page no longer describe the work as numbered phases. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
asalmgren
left a comment
asalmgren
left a comment
There was a problem hiding this comment.
Automated correctness review of the diff (8 inline findings), posted as comments with no verdict attached.
Three findings concern erf.if_snap_partial_cells silently disabling existing immersed-forcing wall/roof treatments, one is a pblh convention mismatch between MRF and YSUNew under erf.pbl_ib_aware, and four are in the new erf.ibseb module.
Things that checked out: the slab tridiagonal coefficients and the affine slab_skin_response / Newton Jacobian are mutually consistent (the balance closes exactly against the implicit slab step); the solar geometry is self-consistent with the documented azimuth convention; the new ibseb_* plot variables line up between ERF.H::derived_names and the mf_comp sequence in Write3DPlotFile; m_fab_start[mfi.LocalIndex()] is valid in every consumer; and the default-off paths look bit-identical (snapb returns b unchanged, and with pbl_ib_aware = false the *_eff work arrays are plain copies with ksrf = klo, zib = 0). All six modified/added .cpp files compile clean against the existing CMake configuration.
| myhalf * (t_blank_arr(i, j-1, k ) + t_blank_arr(i-1, j-1, k )); | ||
| Real t_blank_south_raw = (t_blank_xface_arr) ? snapb(t_blank_xface_arr(i, j-1, k )) : | ||
| myhalf * (snapb(t_blank_arr(i, j-1, k )) + snapb(t_blank_arr(i-1, j-1, k ))); | ||
| const Real t_blank_south = (t_blank_south_raw < small_volfrac) ? zero : t_blank_south_raw; |
There was a problem hiding this comment.
erf.if_snap_partial_cells makes the MOST roof wall function unreachable.
(The issue is on line 461 just below, and identically on line 636 in ImmersedForcingBuildings_Ymom; anchored here because 461 is not inside a diff hunk.)
const Real roof_mask = (t_blank > zero && t_blank < t_blank_below && t_blank_above == zero && l_use_most) ? one : zero;Once snapb() maps the blanking to {0, 1}, the top solid cell of a building has t_blank == 1 and t_blank_below == 1, so t_blank < t_blank_below can never be true. roof_mask is therefore identically zero everywhere whenever erf.if_snap_partial_cells = true.
Concrete scenario: a height-map building run with erf.buildings_type = ImmersedForcing and erf.if_use_most = true currently gets a log-law target velocity on its roof cells (the top cell is partially blanked, e.g. t_blank = 0.6 < 1, so roof_mask == 1 and the compute_if_most_target_vel branch fires). Turning on the new switch silently drops that branch: the roof cell no longer matches roof_mask, falls through to interior_mask == 1, and instead gets the full interior drag rho*CdM*ux*windspeed toward zero. No warning is printed, and the same happens for y-momentum at line 636.
Suggest either relaxing the roof test to t_blank <= t_blank_below (matching the <= already used by the side-wall masks on lines 462-463) or warning/aborting when if_snap_partial_cells and if_use_most are combined.
There was a problem hiding this comment.
Fixed in 41c458d. With the snap the roof is the top solid cell: the test is t_blank <= t_blank_below when erf.if_snap_partial_cells is on and stays < otherwise, so the raw-fraction path is unchanged and roof_mask fires on a snapped staircase (x- and y-momentum). PartialCells reran for 19000 steps with the roof law active: theta within 0.000 K of neutral, max |w| 0.48 m/s (was 0.60), the wake unchanged.
| myhalf * (snapb(t_blank_arr(i, j-1, k )) + snapb(t_blank_arr(i-1, j-1, k ))); | ||
| const Real t_blank_south = (t_blank_south_raw < small_volfrac) ? zero : t_blank_south_raw; | ||
|
|
||
| const Real dx_z = (z_cc_arr) ? (z_cc_arr(i,j,k) - z_cc_arr(i,j,k-1)) : dx_arr[2]; |
There was a problem hiding this comment.
Under the snap, a wall cell gets the MOST wall target and the full interior drag at the same time.
(Refers to lines 462-467 just below, and 637-642 in ImmersedForcingBuildings_Ymom; anchored here because those lines are not inside a diff hunk.)
The side-wall masks use <=, while the interior mask tests for exactly 1:
const Real south_mask = (t_blank > zero && t_blank <= t_blank_north && t_blank_south == zero && l_use_most) ? one : zero;
const Real interior_mask = (t_blank == 1.0) ? one : zero;With snapb() forcing t_blank to exactly 0 or 1, the southern-most solid cell of a building satisfies south_mask == 1 and interior_mask == 1. Before the snap the two were mutually exclusive: MOST fired only on partially blanked cells (0 < t_blank < 1) and interior_mask only on fully blanked ones.
Concrete scenario: on a snapped building wall the code first adds -(u_target - ux) * rho * CdM * U_s to pull the face toward the log-law target from compute_if_most_target_vel, then immediately adds interior_mask * rho * CdM * ux * windspeed pulling the same face toward zero. Both are summed into drag (and into lambda for the point-implicit form), so the wall law is cancelled by the interior damping and erf.if_use_most degenerates into a no-slip staircase. If that is intended it should be stated, since the wall-law inputs (z0, the stability correction, the theta_surf reconstruction) are still evaluated at cost with no effect.
There was a problem hiding this comment.
Fixed in 41c458d. Under the snap interior_mask excludes any cell that carries a wall or roof law (t_blank == 1 && !(l_snap && most_mask > 0), in all three momentum components, with the roof drag counted in z), so a boundary cell gets the wall law alone and the interior cells the drag; the point-implicit lambda follows the same masks. With the snap off the masks are as before, so the exact-box path (where a boundary cell does carry both) is bit-identical; the DataStruct comment, Inputs.rst and the PartialCells README now say so.
| const Real t_blank_north = snapb(t_blank_arr(i , j+1, k)); | ||
| const Real t_blank_south = snapb(t_blank_arr(i , j-1, k)); | ||
| const Real t_blank_east = snapb(t_blank_arr(i+1, j , k)); | ||
| const Real t_blank_west = snapb(t_blank_arr(i-1, j , k)); |
There was a problem hiding this comment.
Under the snap, the building-wall thermal boundary conditions become unreachable while roofs keep working.
(Refers to the else if on lines 1122-1125 below, and the identical wall test in the SURFACE HEAT FLUX block around line 1195+.)
} else if (((t_blank > zero && t_blank < t_blank_west && t_blank_east == zero) || ...Each disjunct requires t_blank > 0 and t_blank < t_blank_<neighbour>, i.e. a strictly partially blanked cell. With snapb() above forcing t_blank into {0, 1}, t_blank < t_blank_west implies t_blank == 0, which fails t_blank > zero. No cell can ever satisfy the wall test.
The roof branch on line 1118 (t_blank > 0 && t_blank_above == zero && t_blank_below == one) is not affected and still fires for the top solid cell.
Concrete scenario: a deck that sets erf.if_snap_partial_cells = true together with erf.if_init_surf_temp (or erf.if_surf_temp_flux) on a height-map building silently heats only the roofs. The walls stop receiving any surface-temperature relaxation, with no diagnostic, and the run looks superficially fine. Worth either aborting when the two are combined, or reformulating the wall test in terms of the neighbour blanking alone (e.g. t_blank > zero && t_blank_west == one && t_blank_east == zero) so that it survives the snap.
There was a problem hiding this comment.
Fixed in 41c458d. The four wall faces are named once per cell (south_face etc.): with the snap they read the neighbour blanking (t_blank == 1 && t_blank_north == 1 && t_blank_south == 0 for the south face), without it the original t_blank < neighbour forms, and both the surface-temperature and the surface-heat-flux blocks use them. The roof test was left as it was.
| zval = (use_terrain_fitted_coords) | ||
| ? Compute_Zrel_AtCellCenter(i, j, kpbl, z_nd_arr) | ||
| : (kpbl + myhalf) * gdata.CellSize(2); | ||
| : ((kpbl + myhalf) * gdata.CellSize(2) - zib); |
There was a problem hiding this comment.
MRF and YSUNew adopt opposite pblh conventions under erf.pbl_ib_aware, yet both feed the same SurfLayer->set_pblh().
Here MRF subtracts zib from every height:
zval = ... : ((kpbl + myhalf) * gdata.CellSize(2) - zib);and clamps with pblh_emp = myhalf * gdata.CellSize(2) (line 434), so pblh_pred_arr / pblh_corr_arr end up as a depth measured from the top of the immersed column.
Source/PBL/ERF_ComputeDiffusivityYSUNew.cpp does the opposite: zval stays absolute and the surface height is carried as z_sfc = zib (line 666), with pblh_min = z_sfc + 0.5*dz_terrain (line 681), so its pblh_corr_arr is a height from the domain bottom.
Each is internally consistent, but both then do pblh_out(i,j,k) = pblh_corr_arr(i,j,0) (MRF:688, YSU:1116) and SurfLayer->set_pblh(level, pblh_mf) (MRF:1213, YSU:2025), and every downstream consumer -- the Beljaars correction, the pblh plot variable, erf.ibseb.z_i_mode = pblh -- reads that single field with one interpretation.
Concrete scenario: a 40 m building with a 500 m mixed layer above its roof. With erf.pbl_type = MRF the stored pblh for those columns is 500; with erf.pbl_type = YSUNew and the identical state it is 540. One of the two must be wrong for any consumer of get_pblh(). Please pick a single convention (absolute height looks safer, since that is what the non-IB columns already report) and make both schemes match.
There was a problem hiding this comment.
Fixed in 41c458d, absolute height in both schemes as you suggest. MRF keeps its internal arrays relative to the immersed column (that is what its K profile and bounds want) and adds zib when it copies out: pblh_out = pblh_corr + zib and Turb_lengthscale = pblh_corr + zib, so get_pblh(), the pblh plot variable and Lturb carry one convention. zib is zero without erf.pbl_ib_aware, so nothing else moves. The TurbStruct comment and Inputs.rst state the convention, and the new PBL_IBAware_MRF_Tiling test checks the stored field varies and is tiling-invariant.
| if (d == 2) { | ||
| Real zi = use_pblh ? pb(i, j, 0) : zi_fallback; | ||
| if (!(zi > 0.0 && zi < 1.0e5)) { zi = zi_fallback; } | ||
| depth = amrex::max(zi - (pzf[f] - z_lo), phb[f]); |
There was a problem hiding this comment.
erf.ibseb.z_i_mode = "pblh" double-subtracts the roof height when combined with MRF + erf.pbl_ib_aware.
depth = amrex::max(zi - (pzf[f] - z_lo), phb[f]);zi comes from pb(i,j,0), i.e. m_SurfaceLayer->get_pblh(lev), and this line treats it as a height above the ground by subtracting the roof-face height pzf[f] - z_lo. But with erf.pbl_ib_aware on and the MRF scheme, ComputeDiffusivityMRF already stores pblh relative to the top of the immersed column (see ERF_ComputeDiffusivityMRF.cpp:388, where zib is subtracted from every zval), so the roof height is removed twice.
Concrete scenario: a 200 m building with a 300 m mixed layer above the roof. MRF stores pblh = 300; this line computes depth = max(300 - 200, 200) = 200 m instead of 300 m, so w* = cbrt(g/theta * H/(rho c_p) * depth) is about 13% low on every roof face, propagating into Ut_eff, u*, H_coeff and the skin temperature. With YSUNew (which keeps pblh absolute) the same input gives the intended 300 m, so the answer also depends on which PBL scheme is selected. Settling the MRF/YSUNew convention would fix this too; otherwise this line needs to know which convention it is handed.
There was a problem hiding this comment.
Resolved by the convention fix (41c458d): get_pblh() is the absolute height for both schemes, so subtracting the roof-face height here is now right with MRF as well as YSUNew.
| const Real U = std::sqrt(uu[k] * uu[k] + vv[k] * vv[k]); | ||
| const Real dU = U - U1; | ||
| const Real rib = CONST_GRAV * (k * dz) * (th - th1) / (th1 * (dU * dU + ustar_floor2)); | ||
| if (rib > ibseb_params.ri_crit) { z_i = geom[lev].ProbLo(2) + (k + 0.5) * dz - geom[lev].ProbLo(2); break; } |
There was a problem hiding this comment.
Mixed height origins in ibseb_bulk_richardson_height.
if (rib > ibseb_params.ri_crit) { z_i = geom[lev].ProbLo(2) + (k + 0.5) * dz - geom[lev].ProbLo(2); break; }The + ProbLo(2) ... - ProbLo(2) cancels, so the crossing branch returns (k + 0.5) * dz, a depth measured above ProbLo(2). The fallback initialised on line 205, Real z_i = z_top; with z_top = geom[lev].ProbHi(2) (line 201), is an absolute coordinate. The two disagree whenever ProbLo(2) != 0.
Concrete scenario: a domain with geometry.prob_lo = 0 0 500 and geometry.prob_hi = ... 2500. On a neutral profile (rib never exceeds ri_crit) the function returns z_i = 2500 m instead of the 2000 m domain depth, 25% too deep. That value is handed to compute_sensible() as z_i_bulk, becomes zi_fallback, and inflates w* on every face. Suggest dropping the dead arithmetic and using z_top = geom[lev].ProbHi(2) - geom[lev].ProbLo(2) so that both branches are depths above the domain bottom.
There was a problem hiding this comment.
Fixed in 41c458d: z_top = ProbHi(2) - ProbLo(2) and the crossing branch returns (k + 0.5) * dz, so both branches are depths above the domain bottom; the dead + ProbLo - ProbLo is gone and the doc comment says which origin is used.
| fill(d_ustar, nf, Real(0.0)); | ||
| fill(d_SW_abs, nf, Real(0.0)); | ||
| fill(d_LW_net, nf, Real(0.0)); | ||
| fill(d_H, nf, Real(0.0)); |
There was a problem hiding this comment.
The previous step's d_H is not checkpointed, so erf.ibseb.convective_velocity = deardorff does not restart consistently.
save_state() writes only d_T_skin and d_T_slab into the IBSEBState field, and build() re-fills d_H with 0 here. compute_sensible() reads d_H as the previous step's flux to form the convective velocity scale:
if (conv && pH[f] > 0.0) { ... wstar = std::cbrt(CONST_GRAV / th * pH[f] / (rho * c_p) * depth); }On restart, ERF::init_ibseb() does call compute_sensible() once, but with olen_ground = nullptr, pblh_ground = nullptr and z_i_bulk = 0 (so zi_fallback = m_params.z_i), and against the t_new atmosphere -- not the value solve_balance() left at the end of the pre-checkpoint step, which was computed with that step's old atmosphere and the converged skin temperature.
Concrete scenario: run 200 steps with erf.ibseb.convective_velocity = deardorff, checkpoint, restart. On the first step after the restart every face sees a different w* than the continuous run, hence a different Ut_eff, u*, H_coeff and skin temperature, and the two trajectories separate. The PR description states "the balance itself restarts exactly", which holds only for the default convective_velocity = none. Adding d_H (one more slot per face) to save_state() / load_state() would close this.
There was a problem hiding this comment.
Fixed in 41c458d. The sensible flux is the second component of each slot of IBSEBState (six slots of 2 + n_layers), and init_ibseb() reloads the state after its initial diagnostics so the checkpointed flux, not the diagnostic one, is what the next step's convective scale reads. The WallFunction case now runs the Deardorff deck through a checkpoint at step 300 and compares the step-599 face dumps: identical to 0.0 with the fix; on the previous binary the same check fails with the skin 4e-3 K off, so it discriminates. The PR text about restarting exactly now holds for the convective scale too.
| if (N < 1 || N > SLAB_MAX_LAYERS) return 0.0; | ||
| const amrex::Real Fourier = k * dt / (rho_cp * dz * dz); // interior Fourier number | ||
| // Coefficients of a l T_{l-1} + b_l T_l + c_l T_{l+1} = d_l (all-plus form). | ||
| amrex::Real alpha[SLAB_MAX_LAYERS], gamma[SLAB_MAX_LAYERS]; |
There was a problem hiding this comment.
Fixed 256-element stack arrays in a per-face GPU kernel, regardless of n_slab_layers.
amrex::Real alpha[SLAB_MAX_LAYERS], gamma[SLAB_MAX_LAYERS];plus amrex::Real W[SLAB_MAX_LAYERS]; in slab_skin_response (line 107), on top of Real T[ibseb::SLAB_MAX_LAYERS]; in IBFaceSet::compute_ground (ERF_IBFaceSet.cpp:451) and IBFaceSet::solve_balance (ERF_IBFaceSet.cpp:485). With SLAB_MAX_LAYERS = 256 and an 8-byte Real that is roughly 8 KB of thread-local storage per face, independent of erf.ibseb.n_slab_layers -- the default of 4 layers pays exactly the same as 256.
Both call sites are amrex::ParallelFor bodies over device vectors, so on a GPU build every thread spills about 8 KB to local memory. Concrete scenario: solve_balance launched over 1e4 faces on an A100 -- each thread needs ~8 KB of local memory, so occupancy collapses to a couple of warps per SM and the per-step cost of the balance is dominated by local-memory traffic rather than by the Newton solve. The PR notes GPU is untested; this is the part most likely to bite there.
Suggest sizing the arrays from a much smaller compile-time bound (16 or 32 layers covers every realistic wall) and tightening the n_slab_layers check in IBSEBParams::init_params() to match -- no cost on CPU, and it makes these kernels usable on GPU.
There was a problem hiding this comment.
Fixed in 41c458d: SLAB_MAX_LAYERS = 32, erf.ibseb.n_slab_layers is checked against it in init_params() and Inputs.rst says [1, 32]. That is 3 x 32 x 8 bytes per thread across the two kernels instead of 8 kB.
|
Re the plotfile documentation finding in the review body: fixed in 9510a3e. Also in that commit, from the same review: the GNU Make wiring, the YSUNew mixing length above the column's surface, the removal of the dead |
asalmgren
left a comment
asalmgren
left a comment
There was a problem hiding this comment.
Correctness review of the diff (code only; data/inputs/docs checked for name consistency). Five findings below. The first is a hard build break that I reproduced locally — it blocks everything else, including CI, so the rest have not been exercised.
| #include <AMReX_REAL.H> | ||
| #include <AMReX_GpuQualifiers.H> | ||
| #include <ERF_NumericalConstants.H> | ||
| #include <ERF_SolarGeometry.H> |
There was a problem hiding this comment.
The branch does not compile — ERF_SolarGeometry.H does not exist.
This header is not present anywhere in the branch (or in origin/development), and the five functions the wrappers below call — compute_solar_declination, compute_solar_hour_angle, compute_solar_zenith_angle, compute_solar_azimuth_angle, compute_earth_sun_distance_factor — are defined nowhere either.
Scenario: any build at all. ERF_IBFaceSet.cpp includes this header, and it is added unconditionally by both CMake/BuildERFExe.cmake and Source/ImmersedBoundarySEB/Make.package, so every configuration fails — not just runs with erf.ibseb.enable = true. Reproduced locally with this repo's own compile flags:
ERF_IBSEBSolar.H:7:10: fatal error: ERF_SolarGeometry.H: No such file or directory
Probably the shared solar-geometry header lived on a branch that was renamed/reworked before merging. With a stub header in place every other file in this PR compiles clean (double and single precision, -Wall -Wextra -Wshadow), so this is the only build blocker.
There was a problem hiding this comment.
Fixed in 8d20391 (reference outputs in 27073f4). The five solar functions (Spencer declination and equation of time, hour angle from UTC and longitude, zenith, clockwise-from-north azimuth, Earth-Sun distance) now live in ERF_IBSEBSolar.H itself, in the ibseb namespace; the header includes what it uses (AMReX_Algorithm.H for amrex::min, which a standalone compile of the header caught). They were the wrappers' targets before #4004 replaced the shared header with the CESM orbital formulation of ERF_OrbCosZenith.H, which has no azimuth, so the balance keeps its own. Verified bitwise against the removed ERF_SolarGeometry.H over one million random (day, latitude, longitude, time) samples and the four case sites (0 mismatches), so the sunrise case and the IBSEB_Cube gold are unchanged. The lesson taken: a new module includes only headers that are on development; a header from an open PR is copied, not included.
| // With the snap every solid cell has t_blank = 1, so a roof cell is | ||
| // the top solid cell (t_blank <= t_blank_below) and a cell carrying a | ||
| // wall law is not also an interior cell. | ||
| const Real roof_mask = (t_blank > zero && (l_snap ? t_blank <= t_blank_below : t_blank < t_blank_below) && t_blank_above == zero && l_use_most) ? one : zero; // roof cell |
There was a problem hiding this comment.
erf.if_snap_partial_cells does not actually give a 0/1 face blanking in the default build, which is what these masks assume.
snapb() snaps the cell-centred blanking, but the momentum kernels then average two cells onto the face:
myhalf * (snapb(t_blank_arr(i,j,k)) + snapb(t_blank_arr(i-1,j,k)))The face-centred arrays that would give a true 0/1 are only allocated under USE_FC_FACTORY, which defaults to 0 (Source/EB/ERF_EB.H:98), so t_blank_xface_arr is an empty Array4 and the averaging path is what actually runs.
Scenario: a default build with erf.if_snap_partial_cells = true. Every wall face gets t_blank == 0.5, not 1, so:
interior_mask(t_blank == 1.0) never fires on a wall face — there is no full blanking and no no-penetration condition, contrary to the comment above (“wall faces fully blanked”, “every solid cell has t_blank = 1”);- the
<→<=relaxation on this line fires on vertical east/west wall x-faces at the top row of a building (t_blank == t_blank_below == 0.5,t_blank_above == 0), handing a side face the horizontal-roof wall law it did not get before.
Same in ImmersedForcingBuildings_Ymom (line 651). Either snap the face value (and the cell-average fallback) to {0,1}, or restate the masks in terms of the cell blanking, which ImmersedForcingBuildings_Scalar already does correctly.
There was a problem hiding this comment.
Fixed in 8d20391 (reference outputs in 27073f4), and you are right about the default build: snapb() snapped the two cells and the average of 1 and 0 is 0.5 on every wall face, so interior_mask never fired there and the <= relaxation handed side faces the roof law. The momentum kernels now take a face's blanking, with the snap on, from the two cells it joins (solid when either is), in both builds, so the staircase does not depend on USE_FC_FACTORY; a face between a solid and a fluid cell is wall-normal and gets the interior drag toward zero (no penetration) whatever row it lies in, and a face between two solid cells carries the roof or wall law of its row or the interior drag, never both. With the snap off the expressions reduce to the previous ones exactly, so the raw-fraction path is bit-identical. Same in y and z. While at it: the terrain kernels (ImmersedForcingTerrain_*) had also been reading snapb(), and their wall law is weighted by the fluid fraction (1 - t_blank), which a snapped face has not, so the top faces of terrain got no forcing at all under the snap; the switch now applies to the buildings kernels only and the doc rows say so. The comment above the kernel, SolverChoice::if_snap_partial_cells and the Inputs.rst row describe the face rule. PartialCells rerun with the new faces (4 ranks, 19000 steps): theta within 0.000 K of neutral, max |w| 0.49 m/s (0.48 before), the wake's minimum u -0.47 m/s (-0.33 before, the wall-normal faces now damped to zero rather than half-weighted). BuildingSet (height map, 4 ranks, 6 h) passes all eight checks; the roof means at 11:00 move by about -2 K (slab 318.6, north block 329.5, cube 321.8, far block 329.4 K, from 320.2/332.2/323.7/332.1), w* 0.13-0.75 m/s and every roof unstable by then; the README numbers are refreshed.
| // SurfaceLayer::get_pblh(). zib is zero without erf.pbl_ib_aware. | ||
| ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { | ||
| pblh_out(i, j, k) = pblh_corr_arr(i, j, 0); | ||
| pblh_out(i, j, k) = pblh_corr_arr(i, j, 0) + zib_arr(i, j, 0); |
There was a problem hiding this comment.
(Anchored here because the line the problem is really on — ApplyPBLHSmoothing(...) at line 673 — is unchanged context and not in a diff hunk. This line is the other half of the same issue.)
PBLH smoothing and erf.pbl_ib_aware operate in different vertical frames.
With pbl_ib_aware, pbl_height_corrector now holds a height above each column's own immersed surface, and zib is only added here. ApplyPBLHSmoothing runs in between, at line 673.
Scenario: erf.pbl_ib_aware = true together with erf.enable_pblh_smoothing = true, over a building. The smoothing stencil averages relative heights across neighbouring columns whose surfaces differ by the full building height, mixing a roof-relative value with a ground-relative one, and only then is the per-column offset restored — so the stored PBLH near a building edge is wrong by roughly the building height times the stencil weight.
Either smooth after adding zib, or exclude IB columns from the stencil. Note the new PBL_IBAware_*_Tiling tests run with smoothing off (the default), so CI will not catch this.
There was a problem hiding this comment.
Fixed in 8d20391 (reference outputs in 27073f4). MRF now adds zib on the work box before ApplyPBLHSmoothing and takes it off after, so the stencil mixes absolute heights (the stored convention, as YSUNew already did), and both schemes keep the smoothed value at or above the column's own floor (half a cell above its surface, at least 10 m), which the smoothing could otherwise violate beside lower neighbours. Measured on the tiling deck with erf.enable_pblh_smoothing = true (2 ranks, 10 steps, old binary 9510a3e vs this one), stored PBLH along the row through the 40 m cube: at the roof edge 98.3 m before, 93.3 m now; on the ground column beside it 80.1 m before (it had been averaged with the roof's roof-relative value), 85.1 m now; the largest change is 10.0 m, the building height times the stencil weight. The floor mattered too: on a strongly stable sounding (ground columns at the 10 m floor) YSUNew stored 44.8 m on the roof-edge columns and 39.7 m on the corner ones, below the roof's own first cell at 45 m, and now stores 45.0; MRF's roof-edge column goes 53.5 -> 50.0 (its 10 m floor above the roof). zib is zero without erf.pbl_ib_aware, and the flat case is bit-identical (the tiling parity tests pass unchanged). New regression test PBL_IBAware_MRF_Smoothing (gold plotfile, the same deck with the smoothing on): the previous binary fails it with Lturb 10 percent off, so CI now covers the pair.
| * non-positive interval disables both. | ||
| */ | ||
| void | ||
| ERF::ibseb_report (int nstep, Real time) |
There was a problem hiding this comment.
The [IBSEB] summary and CSV step number are off by one, and step 0 is written twice.
ERF::Evolve calls post_timestep(step, ...) with the 0-based loop index, while the step just completed is step + 1; the other reporters in post_timestep compensate (e.g. cloud_chamber_budget->report(nstep + 1, ...)). ibseb_report passes nstep straight through to IBFaceSet::report(time, nstep, ...).
Scenario: a fresh run with erf.ibseb.csv_int = 100. init_ibseb() writes a row labelled step=0 for the initial state; after the first step this function writes a second row also labelled step=0 (with a different time_s), and every later row is one step earlier than the matching plotfile number. csv_now = (nstep % csv_int == 0) inherits the same shift, so rows land after steps 1, 101, 201 rather than 100, 200, 300.
There was a problem hiding this comment.
Fixed in 8d20391 (reference outputs in 27073f4): post_timestep passes nstep + 1, as its other reporters do, so the summary and the CSV rows carry the number of completed steps (the plotfiles' numbering); the initial report of init_ibseb() is step 0 and is written once, the rows land after steps csv_int, 2 csv_int, ... and the face dumps are tagged the same way. The case scripts that had adapted to the old labels (the sunrise and building-set checkers used (step + 1) dt as the time of a dump, the closure and Deardorff restart checks compared step000199 / step000599) now use the dump's own step, and the shortwave-energy check of the isolated building takes the shortwave of a dump as the one computed at the start of its last step. The READMEs' reference lines are relabelled.
| if (!restart_chkfile.empty()) { | ||
| const std::string name = MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "IBSEBState"); | ||
| if (FileExists(name + "_H")) { | ||
| restored = std::make_unique<MultiFab>(grids[lev], dmap[lev], m_ibseb[lev]->state_ncomp(), 0); |
There was a problem hiding this comment.
Restarting with a different erf.ibseb.n_slab_layers reads a mismatched checkpoint field.
restored is sized state_ncomp() = 6 * (2 + n_slab_layers) from the current inputs and handed straight to VisMF::Read on the next line, with no check against what the checkpoint holds.
Scenario: run with erf.ibseb.n_slab_layers = 4, checkpoint, then restart with n_slab_layers = 6. The on-disk IBSEBState has 36 components and the target has 48; the read either aborts inside VisMF or, worse, load_state() then indexes with the new slot stride (dir*2 + (side>0)) * (2 + nl) and silently loads skin temperatures and slab profiles into the wrong slots, so the restart continues from garbage face state.
Writing n_slab_layers alongside the field (or comparing against the header ncomp before reading) would turn this into a clear error.
There was a problem hiding this comment.
Fixed in 8d20391 (reference outputs in 27073f4): init_ibseb() reads the on-disk width with VisMF(name).nComp() before the VisMF::Read and aborts when it differs from 6 x (2 + n_slab_layers) of the deck, naming the checkpoint's layer count and the deck's ("IBSEBState in chk00002 has 36 components, written with erf.ibseb.n_slab_layers = 4; the deck sets erf.ibseb.n_slab_layers = 6. Restart with the checkpoint's value."). The FaceStorage case now restarts its 4-layer checkpoint with n_slab_layers = 6 and requires that abort (README updated); a restart with the same count is unaffected.
…ed face blanking, smoothing frame, report steps, restart width - The five solar functions live in ERF_IBSEBSolar.H again (erf-model#4004 replaced the shared header with the CESM orbital form, which has no azimuth); bitwise identical to the removed header over 1e6 samples; the header includes what it uses. - With erf.if_snap_partial_cells the momentum kernels take a face's blanking from the two cells it joins (solid when either is), in both builds; a face between a solid and a fluid cell is wall-normal and gets the drag toward zero, a face between two solid cells the roof or wall law of its row or the interior drag, never both; the snap off is bit-identical; the terrain kernels no longer read the switch (their wall law is weighted by the fluid fraction). - MRF smooths the PBL height in the absolute frame (zib on before ApplyPBLHSmoothing, off after) and both MRF and YSUNew keep the smoothed height at or above the column's own floor; new gold CTest PBL_IBAware_MRF_Smoothing (the previous binary fails it, Lturb 10 percent off). - ibseb_report takes the number of completed steps (post_timestep's nstep + 1), so the step-0 report is written once and the rows land after csv_int, 2 csv_int, ...; the case scripts and READMEs follow. - A restart whose erf.ibseb.n_slab_layers differs from the checkpoint's aborts naming both (VisMF header width); FaceStorage runs that negative case. - Shortwave/Longwave rank-independence checks: geometry, view, shadow and shortwave exact, atmosphere columns to 1e-9 (the dump moved from step 1 to step 2). Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…k slab within the 32-layer bound, reference outputs rerun - SlabConduction's thick-slab deck used 200 layers, outside the 32-layer bound set in the first review round; it is now 32 mm in 1 mm layers over 20 s (the wave from the interior travels 4 mm and the skin side stays untouched), README and checker text follow. - PartialCells (max |w| 0.49 m/s, wake -0.47 m/s), BuildingSet (roof means about 2 K lower with the wall-normal faces damped, every roof unstable), PrognosticSkin (step labels, 201 dumps) and FaceStorage reference outputs rerun with the new binary. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
asalmgren
left a comment
asalmgren
left a comment
There was a problem hiding this comment.
Review of the erf.ibseb surface energy balance and the two immersed-boundary switches. I read every source hunk (~5,100 lines under Source/) and opened the surrounding code for context; the rest of the diff is docs, decks and gold files.
Five findings, left inline. The first is the one I would hold the merge on -- it is a silent capability loss in a feature the PR advertises, and it changes what the PartialCells regtest is actually testing. The other four are robustness and scaling concerns rather than correctness bugs in the CI configurations.
Things I checked and found correct:
- Backward compatibility of the
snapb/fbrefactor in all four buildings kernels: withl_snap == falsethe expressions reproduce the previous ones exactly, including thek == 0and face-centred-blanking branches. if_snap_partial_cellsforcingif_implicit_dragis ordered after theif_implicit_dragqueryAdd, so it does take effect.- Plotfile component ordering:
setPlotVariablesrebuildsplot_var_namesinderived_namesorder, and the newibseb_*write block sits betweenterrain_IB_maskandvolfrac, matchingERF.H. m_fab_start[mfi.LocalIndex()]is safe -- everyMFIterused sharesgrids[lev]/dmap[lev]withterrain_blanking[lev], and none is tiled or inside an OpenMP region.terrain_blankingis rebuilt andFillBoundary-ed on both the fresh-start and restart paths beforeInitData_post()callsinit_ibseb().- The slab tridiagonal coefficients,
slab_skin_responselinearity, the Newton Jacobian sign, theLW_net/e_effalgebra, and theH*A/(c_p*V*Pi)source conversion all check out dimensionally and algebraically. - Solar geometry (declination, hour angle, zenith, clockwise-from-north azimuth, sun vector) is self-consistent, and the DDA nudge
t0*salways steps into the fluid column for both wall sides. - MRF's relative-height convention with the
+/-zibaroundApplyPBLHSmoothing(which does write intogrow(xybx_tile, passes-1)) and YSUNew's absolute-height convention with only the floor clamp are each internally consistent; theng_availcheck now includesterrain_blank's halo. Inputs.rstandPlotfile3DReference.rstcover all the newerf.*inputs and all eightibseb_*plot variables, per the project convention.
Thank you for the thorough test coverage and documentation -- the per-term verification cases made this much easier to review.
| // With the snap every solid face has t_blank = 1: a roof face lies in | ||
| // the top solid row (t_blank <= t_blank_below, the face above fluid), | ||
| // and a face carrying a wall law is not also an interior face. | ||
| const Real roof_mask = (!normal_face && t_blank > zero && (l_snap ? t_blank <= t_blank_below : t_blank < t_blank_below) && t_blank_above == zero && l_use_most) ? one : zero; // roof cell |
There was a problem hiding this comment.
Under erf.if_snap_partial_cells, the wall/roof law here degenerates to no-slip.
compute_if_most_target_vel() (line 62, unchanged by this PR) returns
Real uTarget = (1 - t_blank) * ustar / KAPPA * (std::log(myhalf * delta / z0) - psi_m);With the snap on, fb() returns max(snapb(c1), snapb(c2)), which is exactly 0 or 1. Every face that can set roof_mask / south_mask / north_mask requires both t_blank > zero and !normal_face (i.e. the two cells snap alike), so t_blank == 1 on every such face and uTarget is identically zero. bc_forcing_x then collapses to ux, and the three MOST branches below just relax the momentum toward zero at rate CdM * U_s (with U_s = 1). erf.if_z0, erf.if_stability_correction and the log profile have no effect at all on a snapped building.
The same (1 - t_blank) weighting also makes wall_mask and east_west_mask unreachable under the snap, since both require t_blank < one.
This contradicts the stated intent in the PR description ("a face between two solid cells [gets] the roof or wall law of its row"), and it means Exec/RegTests/ImmersedForcingTest/PartialCells is validating a no-slip staircase rather than a wall-modelled one. The (1 - t_blank) factor is a partial-cell weighting that no longer makes sense once the blanking is snapped -- it needs to be bypassed (or replaced by a snap-appropriate weight) on the snap path.
Same issue at the Ymom and Zmom mask blocks (lines 636 and 826).
There was a problem hiding this comment.
Fixed in 665b641, and thank you: with the snap t_blank was 1 on every face that carries a law, so (1 - t_blank) made the target zero and the PartialCells deck was indeed running a no-slip staircase. The snap path now passes a zero blanking to compute_if_most_target_vel() (t_blank_law), so a snapped boundary face takes the full log-law target u* / kappa (ln(delta/2 z0) - psi_m); the raw-fraction path keeps its partial-cell weight and is unchanged. The kernel comment, SolverChoice::if_snap_partial_cells and the Inputs row say so, and note that the partial-cell branches (wall_mask, east_west_mask) do not arise under the snap. run_partial_cells.sh now proves the law is live: two 400-step runs with erf.if_z0 0.01 and 0.1 must differ in x_velocity (max |du| 0.20 m/s with this commit; the previous binary gives exactly 0). The 19000-step run and BuildingSet rerun with the law live: PartialCells (19000 steps): theta within 0.000 K of neutral, max |w| 0.57 m/s (0.49 with the no-slip faces), the wake's minimum u -0.41 m/s (-0.47). BuildingSet (height map, 6 h): all eight checks pass; the roof means at 11:00 rise by about 1 K with the log-law walls (slab 319.6, north block 331.1, cube 323.0, far block 331.0 K), w* 0.12-0.74 m/s, 96 percent of the roofs unstable; the README numbers are refreshed.
| ERF::ibseb_write_checkpoint (const std::string& checkpointname, int lev) const | ||
| { | ||
| if (!ibseb_params.enable || lev >= static_cast<int>(m_ibseb.size()) || !m_ibseb[lev]) { return; } | ||
| MultiFab state(grids[lev], dmap[lev], m_ibseb[lev]->state_ncomp(), 0); |
There was a problem hiding this comment.
The face state is checkpointed through a full 3-D MultiFab over the whole level: state_ncomp() is 6 * (2 + n_slab_layers), i.e. 36 components at every cell with the default n_slab_layers = 4, and 204 at the allowed maximum of 32 -- to carry data that only exists on a thin shell of building faces.
On the 32x32x16 CI case this is invisible. On a realistic urban deck (say 1024x1024x100) it is on the order of 30 GB of extra checkpoint per level, plus the same amount of transient allocation inside WriteCheckpointFile, dwarfing the conserved state and likely OOM-ing the run at checkpoint time.
A face-indexed side file, or a sparse layout restricted to the building bounding boxes, would avoid this without changing the restart semantics.
There was a problem hiding this comment.
Fixed in 665b641. IBSEBState no longer lives on the level's grids: IBFaceSet::build() now also derives, from a reduced block-height map, a BoxArray of 8 x 8 column blocks that reach up to the highest face-owning cell of each block (state_boxarray()), and the field is defined on those boxes, so it scales with the built volume rather than the domain (on the 1024 x 1024 x 100 example with 30-cell buildings over a fifth of the columns that is a few GB against 30, and nothing where there are no buildings). A transfer layer keeps the restart semantics: each rank's faces are written into a global BoxArray of its own grids cut by the state boxes (AllGatherBoxes, owner = the rank the grid belongs to), ParallelCopy moves it onto the field's distribution, and the reverse on restart, so a checkpoint restarts on any rank count. init_ibseb() compares the on-disk BoxArray with the one the blanking gives and aborts naming both if they differ (a checkpoint of other buildings), next to the layer-count check; a level without faces writes nothing. FaceStorage restarts the 4-rank checkpoint on 1 rank and requires the same last CSV row as the 4-rank restart (both exact against the straight run); the debug summary prints the field's boxes and cells against the level's (on the FaceStorage deck, a 128^3 level with one skyscraper: 12 boxes, 15,616 of the level's 2,097,152 cells, 2.1 MB on disk against 288 MB before). The transient in WriteCheckpointFile is the same two sparse fields. On the four-building BuildingSet, the checkpoint written at the end of the six hours restarts on 4 and on 2 ranks with identical skin and slab temperatures on all 432 faces after two further steps (fluxes to 1e-10).
| // Largest blanking in each column of this rank; reduced below so every | ||
| // rank labels the same solid columns. The highest solid cell of each | ||
| // column gives the column top for the ray cast. | ||
| std::vector<Real> colmax(static_cast<size_t>(nx) * ny, 0.0); |
There was a problem hiding this comment.
build() allocates four full-domain host arrays of size nx*ny (colmax, colk, h_col_top, label), does two nx*ny all-reduces (lines 158-159), and uploads an nx*ny device copy (d_col_top, line 168) on every rank.
At city scale (4096x4096 columns) that is roughly 130 MB of host temporaries plus 33 MB of device memory per rank, and two 33 MB all-reduces. Multiplied by ranks-per-node this is a real OOM risk, and it is not covered by the performance caveats in the PR description, which mention only the ray-sampling cost.
Distributing the column arrays over the box array, or restricting them to the buildings' bounding box, would bound this.
There was a problem hiding this comment.
Fixed in 665b641 as far as a replicated column map allows: build() first reduces the bounding box of the solid columns and allocates colmax, colk, the column tops, the labels and the device copy over that box only (8 bytes per built column per rank), and the two ray walks look a column up through ibseb::column_top(), which returns open ground outside the box; the faces' building ids and heights index the box. On an isolated building or a district that is a small fraction of the domain, which is what the cases run; a city that fills the domain still carries an nx x ny map per rank, which the ray cast needs because a ray can cross the domain. That is now stated in the PR's performance row, with a distributed column map as the follow-up if the replicated one ever binds. Shortwave and Longwave rerun: every view-fraction, shadow and shortwave column of the face dumps is identical to the previous run's (max difference 0.0 over 2616 faces), and IBSEB_Cube's gold is unchanged.
| for (int k = 1; k < nz; ++k) { | ||
| const Real th = rth[k] / rho[k]; | ||
| const Real U = std::sqrt(uu[k] * uu[k] + vv[k] * vv[k]); | ||
| const Real dU = U - U1; |
There was a problem hiding this comment.
dU is the difference of horizontal-mean speeds (U = sqrt(uu[k]^2 + vv[k]^2)), but the function's own doc comment specifies |U(z) - U_1|^2 -- the magnitude of the vector shear.
For a veering profile (an Ekman spiral, or any case where the wind turns while keeping its magnitude) the true shear can be large while dU is near zero. The denominator then falls back to the 100*ustar^2 = 1 m^2/s^2 floor, Ri_b is overestimated, and z_i comes back far too shallow -- which feeds straight into w* and the Deardorff wall function.
Suggest const Real dU2 = (uu[k]-uu[0])*(uu[k]-uu[0]) + (vv[k]-vv[0])*(vv[k]-vv[0]); and using that directly.
There was a problem hiding this comment.
Fixed in 665b641: dU2 = (u(z) - u_1)^2 + (v(z) - v_1)^2 of the horizontal-mean wind vector, as the doc comment specified, so a veering wind of constant speed counts as shear. The WallFunction bulk-Richardson case (a unidirectional profile) gives the same depth as before.
| amrex::Real dt, int N, amrex::Real& a, amrex::Real& b) | ||
| { | ||
| amrex::Real W[SLAB_MAX_LAYERS]; | ||
| for (int l = 0; l < N; ++l) { W[l] = T[l]; } |
There was a problem hiding this comment.
advance_slab_dirichlet guards its stack arrays with if (N < 1 || N > SLAB_MAX_LAYERS) return 0.0; (line 55); slab_skin_response has no such check and copies N entries into a Real W[SLAB_MAX_LAYERS].
This is safe today only because IBSEBParams::init_params() rejects n_slab_layers > 32. If that input bound is ever raised, or a caller passes N from somewhere else, this is a silent stack overflow instead of the graceful return 0.0 the sibling routine gives. Worth adding the same guard for symmetry.
There was a problem hiding this comment.
Fixed in 665b641: slab_skin_response() returns a = b = 0 for N < 1 or N > SLAB_MAX_LAYERS, as advance_slab_dirichlet() does, so the two agree if the input bound is ever raised or another caller passes a larger N.
… snap, checkpoint field on column blocks, column map over the built box, vector shear, slab guard - With erf.if_snap_partial_cells the boundary solid faces take the full log-law target (the partial-cell weight 1 - t_blank had made them no-slip); PartialCells now requires that erf.if_z0 changes the flow (0.20 m/s after 400 steps; the previous binary gave exactly 0). - IBSEBState lives on 8 x 8 column blocks up to the highest face-owning cell around the buildings, with a transfer layer (this rank's grids cut by the blocks, AllGatherBoxes, ParallelCopy) so a restart works on any rank count; init_ibseb aborts on a checkpoint of other buildings; FaceStorage restarts the 4-rank checkpoint on 1 rank and rejects a rotated height map. On the 128^3 FaceStorage level the field covers 15616 of 2097152 cells (2.1 MB against 288 MB). - The ray cast's column map, labels and building heights cover the bounding box of the built columns; the two ray walks look a column up through ibseb::column_top(), open ground outside; the shortwave and longwave dumps are unchanged. - The bulk Richardson depth uses |U(z) - U_1|^2 of the wind vector; slab_skin_response guards N like advance_slab_dirichlet. - Reference outputs rerun: PartialCells (w 0.57, wake -0.41), BuildingSet (roof means 319.6/331.1/323.0/331.0 K), PrognosticSkin closure lines. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Follows the fourth review of erf-model#3960, which left the checkpoint scaling with the built volume and the ray cast's column map replicated at full width on every rank. - IBSEBState goes on 4 x 4 column blocks clipped to the k-range that owns faces, instead of 8 x 8 blocks running from the ground to the roof. Face- owning cells are the fluid cells against solid, a one-cell shell, so a block over the interior of a wide footprint now keeps only the layer above the roof while a block on a wall still spans its height. - The field carries one slot per face of a cell, numbered in the (dir, side) order build() already walks, rather than a fixed six; its width is the largest face count on any cell of the level, which build() reproduces from the same blanking on restart. Three nodal fields would have been the obvious way to drop the slot index, but nodal block boxes share their faces at block boundaries and the ParallelCopy of save_state() would have two destinations for one node with only one of them written. On FaceStorage the field is now 32 boxes over 7872 cells at 18 components against 15616 cells at 36, a quarter of the size; n_slots is 3 there because three faces meet on a cell at the rim, and it is never worse than the six the layout used to assume. - The slot numbering is relative, so it would shift if a cell's face list depended on the decomposition. It cannot -- the blanking carries ComputeGhostCells() + 2 ghost cells -- and build() now asserts the ghost layer it rests on. - The column map of the ray cast drops from three Real arrays plus a mask to one int array of the highest solid cell of each column: colmax was redundant (a column is solid exactly when its top is not -1), h_col_top held the same fact again as a height, and ibseb::column_top() now builds the height from the index. That is 4 bytes per built column per rank instead of 24, and one all-reduce instead of two. The flood-fill label and stack are scoped so they free before the per-face uploads. - Docs and the FaceStorage README follow the new layout, and the theory notes record what still scales with the built area. FaceStorage, Shortwave and Longwave pass on 1 and 4 ranks: the checkpoint round trip and its redistribution onto another rank count, both restart aborts, 0/2616 shadow mismatches against the checkers' independent ray cast at both zenith angles with the recorded 738 and 644 shadowed faces, and 0/2616 view-fraction mismatches against the independent hemisphere sampling. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
I have opened two follow-up issues for the parts of my earlier review that this PR narrows but does not close, so they are not lost when it merges:
Neither is a correctness problem and neither needs to hold up this PR; both are about behaviour at city scale. I am also pushing a commit to this branch that takes the cheap, answer-preserving part of each: the checkpoint field goes onto 4 x 4 blocks clipped to the k-range that owns faces and carries one slot per face of a cell instead of a fixed six (a quarter of the size on FaceStorage), and the column map collapses from three |
…f-model#3960 The 24 h day passes unchanged; the table now starts from the initial state (step 0) and two last digits move with the relabelled dumps. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
|
On 3d95fc7 (the 4 x 4 blocks clipped to the face-owning k-range, one slot per face of a cell, the integer column map): rebuilt and rerun here on top of the fourth-round fixes. Release and Debug (assertions + bound check, all warnings) build clean; 556 gtests; the four IBSEB / IB-aware CTests pass with the golds unchanged; FaceStorage, Shortwave, Longwave, SlabConduction, WallFunction and PrognosticSkin pass, the shadow / view-fraction / shortwave columns of the dumps bitwise equal to the previous layout's; a fresh checkpoint of the four-building BuildingSet (27 boxes, 1520 of 36864 cells at 3 slots) restarts on 4 and on 2 ranks with identical skin and slab temperatures on all 432 faces, and the checkpoint written by the 8 x 8 layout is rejected by the width check as it should be. GNU Make links. One thing the CI's single-precision jobs would have caught: |
Merge upstream/development (IBSEB, erf-model#3960 and erf-model#4022) into ERF-Fire
…IBSEB PR #3860 turned ERF::m_SurfaceLayer from a single unique_ptr into one entry per domain face, and PR #3960 added ibseb_advance, which reads the surface layer for the wall function beyond neutral. The two merged independently and git merged them cleanly, so development did not compile: ERF_IBSEB.cpp still dereferenced m_SurfaceLayer as a pointer. Index zlo. What the wall function needs is the ground beneath the buildings; a surface layer on a lateral or upper wall says nothing about the stability of the column, and this is the same face that the other zlo-only consumers select (see ERF.cpp:876). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…IBSEB PR #3860 turned ERF::m_SurfaceLayer from a single unique_ptr into one entry per domain face, and PR #3960 added ibseb_advance, which reads the surface layer for the wall function beyond neutral. The two merged independently and git merged them cleanly, so development did not compile: ERF_IBSEB.cpp still dereferenced m_SurfaceLayer as a pointer. Index zlo. What the wall function needs is the ground beneath the buildings; a surface layer on a lateral or upper wall says nothing about the stability of the column, and this is the same face that the other zlo-only consumers select (see ERF.cpp:876). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
* Add split/merge AMR API stubs + flag + Redistribute hook * Implement SplitMergeAtLevelBoundary: 3-part split-new/merge-departees/repopulate * Implement SplitParticlesForRefinement (regrid-time): split coarse particles in newly-refined cells * Implement MergeParticlesAfterDerefining (regrid-time): masked merge in de-refined cells plus escapee absorption * Wire SplitParticlesForRefinement / MergeParticlesAfterDerefining into regrid + FillPatch S_new before microphysics * Use SuperDropletsRealIdx/IntIdx (dev names) instead of stale SoA aliases * Skip refRatio scaling of coalescence bins when split/merge AMR is on (splitting keeps per-cell SD density uniform across levels) * Normalize active flag to 0/1 when used as a count weight in deposition + diagnostics (active=2 split tag must not double-count) * Generalize SDM split-merge tag to level+1 (Phase 1 of multi-level redesign) * Cumulative finest-level split in SplitParticlesForRefinement (Phase 2) * Fold MergeParticlesAfterDerefining into per-level tag-normalizing sweep (Phase 3) * Resync particle container dummy MFs at start of SplitParticlesForRefinement to fix MFIter assert after regrid creates/remakes a level * Redistribute particles before SDM regrid-time split; iterate leaf level with per-tag cumulative factor so the m_particles state stays consistent with the post-regrid BoxArrays * Restrict SDM Part 3 empty-cell repopulation to same coarse parent to prevent SD leakage and runaway L1 growth under count-based refinement * Add ventilation factor (Bayley et al. 2025) to SDM condensation/evaporation phase change * Revert per-level Part 2 no-host branch to retag-only to fix evaporation bias in long moist runs (matches per-step behaviour of 215b705) * Add MergeParticlesAtDerefinement regrid-time pass to reduce np_bin to np_bin/merge_factor in cells that lost fine-level coverage, capping SD count growth from oscillating qt_refine in moist runs * WIP: gradual SDM ice melting (Seifert-Beheng 2006 rate) with mixed ice-water super-droplets; replaces instantaneous melt, untested Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * Add SDM_Box3D_IceMelt regression test seeding pure ice in a warm box that melts partially over 10 steps, exercising the gradual melt model [run-ci] Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * Restrict SDM AMR merges to same log-radius bin (per-cell adaptive range, N_R = target SDs/cell) so only same-size super-droplets combine, preserving the rain-seeding large-droplet tail * Add multiplicity floor to per-step SplitMergeAtLevelBoundary Part 1 split: entrants with multiplicity < split_factor are retagged native instead of split, preventing sub-unity ghost super-droplets * Integer multiplicity split at both SDM split sites: each piece gets floor(mult/k) droplets plus one each of the remainder, keeping multiplicities integer >=1 and conserved so coalescence cannot produce sub-unity ghost super-droplets * Fix missing PlaneAverage include in ERF_AdvanceDycore.cpp * using std::ceil for density scaling of multiplicities to avoid non-integer and <1 values * Implemented Gunn-Kinzer clamp on vterm calculation * Bound wet/mixed-ice terminal velocity: blend Boehm dry-ice speed toward the equivalent raindrop speed by liquid water fraction (Frick et al. 2013), fixing the IceBohm divergence for near-melted ice cores Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> * updated gold files * Update SDM theory docs: gradual melting and mixed ice-water particles, plus terminal-velocity equations (raindrop 4 mm clamp, Bohm ice, Frick mixed-phase blend) with citations * Add terminal-velocity unit tests (Gunn-Kinzer clamp, Bohm ice vs Boehm 1989 Fig.3, Frick melt blend) on host and device; factor the Frick blend into TerminalVelocity helpers * Allow seeded/injected ice particles to take a prescribed apparent density (low-density snowflakes) via ice_apparent_density; default keeps the solid-ice sphere * Add ice deposition-growth unit test (Chen-Lamb habit vs Takahashi 1991) on host and device: grows a 2um crystal 10 min at constant T + water saturation via the MassChange_SV routines, reproducing the plate/column habit transitions * TerminalVelocity gtest: add EvalFromCsv mode (CSV-driven, phase-routed fall speed replicating the SD advection step) for Monte-Carlo diagnosis of anomalous terminal velocities * SDM sublimation: release a fully-sublimated ice core back to a dry aerosol (m_ice=0, a=c=rime=nmono=0) instead of pinning a 1 nm ice sphere -- fixes the spurious IceBohm settling of sub-resolution remnants (residual large terminal velocities / Dane FPE); mirrors the full-melt path * SDM freezing: if the freezing water is below the minimum ice mass (an essentially dry aerosol carrying only trace water), leave it as a dry aerosol instead of seeding a sub-floor ice crystal -- stops the spurious 1 nm ice remnants at the source (complements the sublimation-side release) * SDM IceFrzDep test: seed resolvable droplet water so freezing yields ice above the minimum-ice floor The trace default water (4.19e-42 kg) now freezes to a dry aerosol rather than ice, leaving the test with no ice to grow. Seed 1.0e-14 kg of water per super-droplet so the frozen crystal is well above the minimum ice mass, and regenerate the gold plotfile. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * SDM MultiSpecies gold: regenerate after the trace-water-to-aerosol change The freezing/sublimation fixes release trace default water as a dry aerosol, which shifts the noise-level species fields in this test. Update the gold. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * SDM SV habit table: function-local static constexpr for host+device Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * Fixed precision-related errors and warnings * Bugfix: visc coeff needs to be converted to SI units * TerminalVelocity test: recenter graupel band, add Stokes-regime check * TerminalVelocity test: dump L-H habit + ice-sphere(+Stokes) curves for the size figure * Updated ice benchmarks * Add SDM microphysics unit tests: condensation, melting/sublimation, INAS freezing, coalescence * Factor riming helpers into header; add riming unit tests * Add fractional super-droplet injection with SDM_Bubble2D_FracInjection ctest [run-ci] * Clean up branch-unique SDM inputs (DYCOMS/GATE/TempSourceCold); sync DYCOMS with run-repo Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> * SDM: add prescribed_w steady-updraft advection option Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NiSV5oMFPUeWiTBnLuQU7K * Add SDM ice-process ctests (sublimation/melting/freezing/aggregation/coalescence/riming) Drop the redundant 0D box IceSub/IceMelt tests; EXTRA collision/freezing tests gated and CPU-only. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NiSV5oMFPUeWiTBnLuQU7K * Add SDM ice-microphysics example inputs and READMEs in Exec Sublimation/melting/freezing/riming columns (TemperatureSourceSpatial_cold) and the 0D ice-aggregation box (Bubble), with prescribed soundings. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NiSV5oMFPUeWiTBnLuQU7K * Removing sublimation shaft test for AMD GPUs * SDM: per-box high-multiplicity injection [run-ci] Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_016y46YhAwKTUV94hTCkfwEy * Add ice_aggregation_efficiency input for SDM ice-ice aggregation (default 0.1) * switch kernel velocity terms * Fix fill boundary issue * Rebuild FFT solvers on restart * Fixed GATE readme * SDM riming: at T>=0C route the collected droplet to meltwater instead of freezing it as rime -- guards the supercooled Heymsfield-Pflaum rime density (erf_esati=0 above 0C gives inf/inf in iceSurfaceTemperature; a CPU FPE, silently masked on GPU by std::min(-0.01,NaN)=-0.01) * SDM: allow zero-mass constant aerosols; require at least one non-zero aerosol per mode * SDM: add erf.random_seed input to vary the RNG (init sampling + collisions) per run Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NiSV5oMFPUeWiTBnLuQU7K * SDM/particles: fix OOB indexing, div-by-zero, half-step drift, ghost fill, leak, mean free path [run-ci] * SDM/particles: fix tiled in-place saturation ratio, RNG race, zeta convention, subsidence reduction, precip wiring [run-ci] * Add OpenMP, ASan, terrain-SDM, tiled and single-rank test coverage for particle/SDM paths [run-ci] * Narrow recycle terrain clamp to the lookup, clamp tail index instead of truncating, drop rank-variant test [run-ci] * Add in-repo gold files for SDM_Bubble2D_WoA and SDM_Bubble2D_Tiled [run-ci] * Use ThreadSanitizer for the threaded job; assert the sampled tail index is in range [run-ci] * Update RadiationSimple to current IRadiation interface (double time/dt, updated_lsm) Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Restore lsf_data/nudge_data sizing and the RadiationSimple branch lost when upstream relocated the constructors Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Advertise ice, snow and graupel plot variables for the super-droplet moisture model Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Drop input keys upstream no longer reads; correct recycle-bounds names and terrain text in the SDM docs Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Rename TemperatureSourceSpatial_cold to SDM_Congestus3D_cold to match the warm case Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Reject a species that is already modeled; drop the duplicate ice entry from the aggregation decks Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Remove local agent config and notes accidentally committed under .claude Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Point LC CI back at the main gold-file ref Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Restore SDM_Congestus3D and SDM_RICO3D golds as whole plotfiles; the merge had mixed two box layouts Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Strip trailing whitespace from new decks and profile output Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Regenerate SDM_SublimationShaft, SDM_Congestus3D and SDM_RICO3D golds after the upstream merge Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Drop the ThreadSanitizer and AddressSanitizer CI jobs Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Run the gated SDM cases as assertion-checked smoke tests in the GNU build Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Fix single-precision build: keep the injected super-droplet count in Real [run-ci] Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Leave the coalescence box tests out of the smoke pass; the kernels have unit tests Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Update the moisture capability truth table for the super-droplet ice variables Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Fix single-precision build in the large-scale forcing update; drop an unused variable [run-ci] Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Temporarily point LC CI at the SDM_w_ice gold-file ref Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Revert "Temporarily point LC CI at the SDM_w_ice gold-file ref" This reverts commit d3f78f3. * Rename fo to nbr_off and ofs, and fix a hyphenated word, to satisfy codespell Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Build the super-droplet microphysics unit tests only when particles are enabled Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Silence unused-parameter, unused-variable and shadow warnings in the cold-process code Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Rename the inner reduction locals so they no longer shadow the ReduceData object Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Review fixes: drop duplicate theta nudging, index LSF/nudge data by level, write valid tiles only, guard moisture paths, correct riming pressure units, conserve meltwater on sub-freezing riming Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Review fixes: correct the momentum subsidence gradient and k-range, set lsf.start_time, parenthesize NVAR_max and restore the moist-state comparison * Review fixes: guard log10 in the riming fit, gate the ghost refresh on Lagrangian moisture, drop redundant syncs and dead params, assert on malformed LSF input, document the LSF/nudging/seed options * Compare against the signed Long returned by amrex::Vector::size() * Temporarily point LC CI at the SDM_w_ice gold-file ref * Revert "Temporarily point LC CI at the SDM_w_ice gold-file ref" This reverts commit 63623d5. * Regenerate in-repo golds for the gated SDM tests and drop the particle data from them * Guard the ice-category active tag, convert erf_esati to Pa for the riming heat balance, call snowAccumulation and emit the accumulation fields in canonical order * Regenerate the riming gold and keep the plotfile-selection semantic check on graupel accumulation * Clean up old code * Update unit test for new diffusion call signature * Fix warning * Update goldfile branch * Cleanup diffusion * Remove unused variables * Update checkpoint names used for single face case * Update PBL MOST average indices * Add missing test * Turn off FFT option for CI * Fix warnings * Add gold file for MOST walls test * Skip test if not built with FFT enabled * Add option for FFT to LC CI and temporarily enable * Style * Add doc section for surface layer on sides * Remove goldfile treated as text * Readd goldfile as binary * Move LC gitlab goldfile ref back to main * Fix stress mapping on x face * Fix review issues * Minor cleanup for plotfile variables * Update gold files for new changes * Fix unused warnings * Fix extra ghost cells for z face in compute_fluxes * Fix more unused parameter warnings * Add fixes from review * Generalize tau transpose writes for face corners * Add unit tests for SurfaceLayer and MOSTAverage on side walls * Update gold file * Cleanup MOSTAverage k index functions and add tests * Fix codespell error * Split out unit tests into serial and parallel, add tiling checks * Fix tiling issues for planar average/qsurf * Updates for grid check and lateral boundary fill * Add extra asserts for unsupported options * Add additional asserts on layouts, update tiling and no-terrain paths * Remove unused w rot field and add additional checks * Update k indices for non terrain to be wall-relative * Add check for regional average difference * Update gold file * Update new get_ng_fill and extrap function for lateral faces * Cleanup redundant runtime assertions * Fix warnings and CUDA/HIP builds * Fix CUDA extended lambda in unit tests * Fix another CUDA extended lambda issue * Fixes from review * Update tests to work with new constructor change * Fix issue for SST test * Fix zref interval for face orientation * Update zref handling for explicit specified value * Fix surface layer check to be face aware after incoming changes * Fixes from review * Reduce test max_step * More fixes from review * Fixes from review * Generalize PlanarBoundary for any face * Update regional average test with PlanarBoundary behavior * Fix bounds issue for rans * Fix error from new upstream changes to use zlo SurfaceLayer * Fix RRTMGP SurfaceLayer temperature conversion * Polish RRTMGP surface temperature regression * Fixes from review * Update qsurf tests to still check intended lateral behavior * Fix the build: take the ground face of the per-face surface layer in IBSEB PR #3860 turned ERF::m_SurfaceLayer from a single unique_ptr into one entry per domain face, and PR #3960 added ibseb_advance, which reads the surface layer for the wall function beyond neutral. The two merged independently and git merged them cleanly, so development did not compile: ERF_IBSEB.cpp still dereferenced m_SurfaceLayer as a pointer. Index zlo. What the wall function needs is the ground beneath the buildings; a surface layer on a lateral or upper wall says nothing about the stability of the column, and this is the same face that the other zlo-only consumers select (see ERF.cpp:876). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Detect legacy surface-layer inputs from what the user wrote, not the live table With more than one surface-layer face, ERF decides whether to keep reading the historical unqualified erf.most.* / erf.surface_layer.* keys for zlo by asking whether any such key exists. It asked ParmParse directly, which cannot answer the question: ERF_InputSoundingData.H queries most.surf_temp and most.surf_moist with queryAdd while reading the sounding, and queryAdd inserts the default it was handed whenever the key is absent. By the time InitData_post runs, both keys are in the table holding their negative sentinels whether or not the user ever typed them. Every input_sounding run with more than one surface-layer face therefore looked as though it carried legacy inputs. A deck that also wrote the face-qualified spelling -- Tests/test_files/ABL_MOST_Cloudchamber, for one -- hit the "both spellings present" abort and could not start at all; a deck without it silently took the legacy path. Snapshot the inputs table at the top of ReadParameters, which is the last point at which presence in the table still means the user wrote it, and scan the snapshot. This stays correct if another erf.most.* queryAdd is added elsewhere later, which a fix targeting the two known keys would not. ABL_MOST_Cloudchamber is the regression test: it aborts before this change and initializes cleanly after it. Also document the prefixing rule in SurfaceLayer.rst -- unqualified for a lone zlo face, face-qualified once a second face is declared, unqualified applied to zlo with a warning, both spellings an error -- along with the configurations that now abort away from zlo, which the page previously covered only as "not all existing options are supported". Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * fix gmake. (#4028) * more fixes for gmake. (#4029) * SurfaceLayer::fill_qsurf_with_qsat: intersect with the cell box of the nodal heights The loop box is cell-centred and z_phys_nd is nodal, so `gtbx &= z_phys_nd->fabbox()` mixes index types, and with AMReX_ASSERTIONS every run that reaches this function aborts at initialisation (sameType in AMReX_Box.H). The kernel reads nodes i..i+1, j..j+1, k..k+1 through Compute_Zrel_AtCellCenter, so the right limit is the nodal box converted to cells. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Set up case for noah + rad, remove unused inputs, fix wps ml inputs. * clean up inputs. * Address comments. --------- Co-authored-by: Debojyoti Ghosh <debojyoti.ghosh@gmail.com> Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com> Co-authored-by: Cole Kendrick <kendrick6@llnl.gov> Co-authored-by: Emma Catherine Ware <ware10@llnl.gov> Co-authored-by: Debojyoti Ghosh <ghosh5@llnl.gov> Co-authored-by: Kyle <kyle.pressel@pnnl.gov> Co-authored-by: Ann Almgren <asalmgren@lbl.gov> Co-authored-by: Harish <harish.gopalan@nrel.gov>
* Add split/merge AMR API stubs + flag + Redistribute hook * Implement SplitMergeAtLevelBoundary: 3-part split-new/merge-departees/repopulate * Implement SplitParticlesForRefinement (regrid-time): split coarse particles in newly-refined cells * Implement MergeParticlesAfterDerefining (regrid-time): masked merge in de-refined cells plus escapee absorption * Wire SplitParticlesForRefinement / MergeParticlesAfterDerefining into regrid + FillPatch S_new before microphysics * Use SuperDropletsRealIdx/IntIdx (dev names) instead of stale SoA aliases * Skip refRatio scaling of coalescence bins when split/merge AMR is on (splitting keeps per-cell SD density uniform across levels) * Normalize active flag to 0/1 when used as a count weight in deposition + diagnostics (active=2 split tag must not double-count) * Generalize SDM split-merge tag to level+1 (Phase 1 of multi-level redesign) * Cumulative finest-level split in SplitParticlesForRefinement (Phase 2) * Fold MergeParticlesAfterDerefining into per-level tag-normalizing sweep (Phase 3) * Resync particle container dummy MFs at start of SplitParticlesForRefinement to fix MFIter assert after regrid creates/remakes a level * Redistribute particles before SDM regrid-time split; iterate leaf level with per-tag cumulative factor so the m_particles state stays consistent with the post-regrid BoxArrays * Restrict SDM Part 3 empty-cell repopulation to same coarse parent to prevent SD leakage and runaway L1 growth under count-based refinement * Add ventilation factor (Bayley et al. 2025) to SDM condensation/evaporation phase change * Revert per-level Part 2 no-host branch to retag-only to fix evaporation bias in long moist runs (matches per-step behaviour of 215b705) * Add MergeParticlesAtDerefinement regrid-time pass to reduce np_bin to np_bin/merge_factor in cells that lost fine-level coverage, capping SD count growth from oscillating qt_refine in moist runs * WIP: gradual SDM ice melting (Seifert-Beheng 2006 rate) with mixed ice-water super-droplets; replaces instantaneous melt, untested Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * Add SDM_Box3D_IceMelt regression test seeding pure ice in a warm box that melts partially over 10 steps, exercising the gradual melt model [run-ci] Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * Restrict SDM AMR merges to same log-radius bin (per-cell adaptive range, N_R = target SDs/cell) so only same-size super-droplets combine, preserving the rain-seeding large-droplet tail * Add multiplicity floor to per-step SplitMergeAtLevelBoundary Part 1 split: entrants with multiplicity < split_factor are retagged native instead of split, preventing sub-unity ghost super-droplets * Integer multiplicity split at both SDM split sites: each piece gets floor(mult/k) droplets plus one each of the remainder, keeping multiplicities integer >=1 and conserved so coalescence cannot produce sub-unity ghost super-droplets * Fix missing PlaneAverage include in ERF_AdvanceDycore.cpp * using std::ceil for density scaling of multiplicities to avoid non-integer and <1 values * Implemented Gunn-Kinzer clamp on vterm calculation * Bound wet/mixed-ice terminal velocity: blend Boehm dry-ice speed toward the equivalent raindrop speed by liquid water fraction (Frick et al. 2013), fixing the IceBohm divergence for near-melted ice cores Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> * updated gold files * Update SDM theory docs: gradual melting and mixed ice-water particles, plus terminal-velocity equations (raindrop 4 mm clamp, Bohm ice, Frick mixed-phase blend) with citations * Add terminal-velocity unit tests (Gunn-Kinzer clamp, Bohm ice vs Boehm 1989 Fig.3, Frick melt blend) on host and device; factor the Frick blend into TerminalVelocity helpers * Allow seeded/injected ice particles to take a prescribed apparent density (low-density snowflakes) via ice_apparent_density; default keeps the solid-ice sphere * Add ice deposition-growth unit test (Chen-Lamb habit vs Takahashi 1991) on host and device: grows a 2um crystal 10 min at constant T + water saturation via the MassChange_SV routines, reproducing the plate/column habit transitions * TerminalVelocity gtest: add EvalFromCsv mode (CSV-driven, phase-routed fall speed replicating the SD advection step) for Monte-Carlo diagnosis of anomalous terminal velocities * SDM sublimation: release a fully-sublimated ice core back to a dry aerosol (m_ice=0, a=c=rime=nmono=0) instead of pinning a 1 nm ice sphere -- fixes the spurious IceBohm settling of sub-resolution remnants (residual large terminal velocities / Dane FPE); mirrors the full-melt path * SDM freezing: if the freezing water is below the minimum ice mass (an essentially dry aerosol carrying only trace water), leave it as a dry aerosol instead of seeding a sub-floor ice crystal -- stops the spurious 1 nm ice remnants at the source (complements the sublimation-side release) * SDM IceFrzDep test: seed resolvable droplet water so freezing yields ice above the minimum-ice floor The trace default water (4.19e-42 kg) now freezes to a dry aerosol rather than ice, leaving the test with no ice to grow. Seed 1.0e-14 kg of water per super-droplet so the frozen crystal is well above the minimum ice mass, and regenerate the gold plotfile. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * SDM MultiSpecies gold: regenerate after the trace-water-to-aerosol change The freezing/sublimation fixes release trace default water as a dry aerosol, which shifts the noise-level species fields in this test. Update the gold. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * SDM SV habit table: function-local static constexpr for host+device Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> * Fixed precision-related errors and warnings * Bugfix: visc coeff needs to be converted to SI units * TerminalVelocity test: recenter graupel band, add Stokes-regime check * TerminalVelocity test: dump L-H habit + ice-sphere(+Stokes) curves for the size figure * Updated ice benchmarks * Add SDM microphysics unit tests: condensation, melting/sublimation, INAS freezing, coalescence * Factor riming helpers into header; add riming unit tests * Add fractional super-droplet injection with SDM_Bubble2D_FracInjection ctest [run-ci] * Clean up branch-unique SDM inputs (DYCOMS/GATE/TempSourceCold); sync DYCOMS with run-repo Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> * SDM: add prescribed_w steady-updraft advection option Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NiSV5oMFPUeWiTBnLuQU7K * Add SDM ice-process ctests (sublimation/melting/freezing/aggregation/coalescence/riming) Drop the redundant 0D box IceSub/IceMelt tests; EXTRA collision/freezing tests gated and CPU-only. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NiSV5oMFPUeWiTBnLuQU7K * Add SDM ice-microphysics example inputs and READMEs in Exec Sublimation/melting/freezing/riming columns (TemperatureSourceSpatial_cold) and the 0D ice-aggregation box (Bubble), with prescribed soundings. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NiSV5oMFPUeWiTBnLuQU7K * Removing sublimation shaft test for AMD GPUs * SDM: per-box high-multiplicity injection [run-ci] Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_016y46YhAwKTUV94hTCkfwEy * Add ice_aggregation_efficiency input for SDM ice-ice aggregation (default 0.1) * switch kernel velocity terms * Fix fill boundary issue * Rebuild FFT solvers on restart * Fixed GATE readme * SDM riming: at T>=0C route the collected droplet to meltwater instead of freezing it as rime -- guards the supercooled Heymsfield-Pflaum rime density (erf_esati=0 above 0C gives inf/inf in iceSurfaceTemperature; a CPU FPE, silently masked on GPU by std::min(-0.01,NaN)=-0.01) * SDM: allow zero-mass constant aerosols; require at least one non-zero aerosol per mode * SDM: add erf.random_seed input to vary the RNG (init sampling + collisions) per run Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NiSV5oMFPUeWiTBnLuQU7K * SDM/particles: fix OOB indexing, div-by-zero, half-step drift, ghost fill, leak, mean free path [run-ci] * SDM/particles: fix tiled in-place saturation ratio, RNG race, zeta convention, subsidence reduction, precip wiring [run-ci] * Add OpenMP, ASan, terrain-SDM, tiled and single-rank test coverage for particle/SDM paths [run-ci] * Narrow recycle terrain clamp to the lookup, clamp tail index instead of truncating, drop rank-variant test [run-ci] * Add in-repo gold files for SDM_Bubble2D_WoA and SDM_Bubble2D_Tiled [run-ci] * Use ThreadSanitizer for the threaded job; assert the sampled tail index is in range [run-ci] * Update RadiationSimple to current IRadiation interface (double time/dt, updated_lsm) Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Restore lsf_data/nudge_data sizing and the RadiationSimple branch lost when upstream relocated the constructors Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Advertise ice, snow and graupel plot variables for the super-droplet moisture model Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Drop input keys upstream no longer reads; correct recycle-bounds names and terrain text in the SDM docs Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Rename TemperatureSourceSpatial_cold to SDM_Congestus3D_cold to match the warm case Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Reject a species that is already modeled; drop the duplicate ice entry from the aggregation decks Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Remove local agent config and notes accidentally committed under .claude Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Point LC CI back at the main gold-file ref Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Restore SDM_Congestus3D and SDM_RICO3D golds as whole plotfiles; the merge had mixed two box layouts Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Strip trailing whitespace from new decks and profile output Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Regenerate SDM_SublimationShaft, SDM_Congestus3D and SDM_RICO3D golds after the upstream merge Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Drop the ThreadSanitizer and AddressSanitizer CI jobs Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Run the gated SDM cases as assertion-checked smoke tests in the GNU build Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Fix single-precision build: keep the injected super-droplet count in Real [run-ci] Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Leave the coalescence box tests out of the smoke pass; the kernels have unit tests Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Update the moisture capability truth table for the super-droplet ice variables Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Fix single-precision build in the large-scale forcing update; drop an unused variable [run-ci] Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Temporarily point LC CI at the SDM_w_ice gold-file ref Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Revert "Temporarily point LC CI at the SDM_w_ice gold-file ref" This reverts commit d3f78f3. * Rename fo to nbr_off and ofs, and fix a hyphenated word, to satisfy codespell Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Build the super-droplet microphysics unit tests only when particles are enabled Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Silence unused-parameter, unused-variable and shadow warnings in the cold-process code Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Rename the inner reduction locals so they no longer shadow the ReduceData object Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Review fixes: drop duplicate theta nudging, index LSF/nudge data by level, write valid tiles only, guard moisture paths, correct riming pressure units, conserve meltwater on sub-freezing riming Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Fd1XcC2jtby94LspffcTMb * Review fixes: correct the momentum subsidence gradient and k-range, set lsf.start_time, parenthesize NVAR_max and restore the moist-state comparison * Review fixes: guard log10 in the riming fit, gate the ghost refresh on Lagrangian moisture, drop redundant syncs and dead params, assert on malformed LSF input, document the LSF/nudging/seed options * Compare against the signed Long returned by amrex::Vector::size() * Temporarily point LC CI at the SDM_w_ice gold-file ref * Revert "Temporarily point LC CI at the SDM_w_ice gold-file ref" This reverts commit 63623d5. * Regenerate in-repo golds for the gated SDM tests and drop the particle data from them * Guard the ice-category active tag, convert erf_esati to Pa for the riming heat balance, call snowAccumulation and emit the accumulation fields in canonical order * Regenerate the riming gold and keep the plotfile-selection semantic check on graupel accumulation * Clean up old code * Update unit test for new diffusion call signature * Fix warning * Update goldfile branch * Cleanup diffusion * Remove unused variables * Update checkpoint names used for single face case * Update PBL MOST average indices * Add missing test * Turn off FFT option for CI * Fix warnings * Add gold file for MOST walls test * Skip test if not built with FFT enabled * Add option for FFT to LC CI and temporarily enable * Style * Add doc section for surface layer on sides * Remove goldfile treated as text * Readd goldfile as binary * Move LC gitlab goldfile ref back to main * Fix stress mapping on x face * Fix review issues * Minor cleanup for plotfile variables * Update gold files for new changes * Fix unused warnings * Fix extra ghost cells for z face in compute_fluxes * Fix more unused parameter warnings * Add fixes from review * Generalize tau transpose writes for face corners * Add unit tests for SurfaceLayer and MOSTAverage on side walls * Update gold file * Cleanup MOSTAverage k index functions and add tests * Fix codespell error * Split out unit tests into serial and parallel, add tiling checks * Fix tiling issues for planar average/qsurf * Updates for grid check and lateral boundary fill * Add extra asserts for unsupported options * Add additional asserts on layouts, update tiling and no-terrain paths * Remove unused w rot field and add additional checks * Update k indices for non terrain to be wall-relative * Add check for regional average difference * Update gold file * Update new get_ng_fill and extrap function for lateral faces * Cleanup redundant runtime assertions * Fix warnings and CUDA/HIP builds * Fix CUDA extended lambda in unit tests * Fix another CUDA extended lambda issue * Fixes from review * Update tests to work with new constructor change * Fix issue for SST test * Fix zref interval for face orientation * Update zref handling for explicit specified value * Fix surface layer check to be face aware after incoming changes * Fixes from review * Reduce test max_step * More fixes from review * Fixes from review * Generalize PlanarBoundary for any face * Update regional average test with PlanarBoundary behavior * Fix bounds issue for rans * Fix error from new upstream changes to use zlo SurfaceLayer * Fix RRTMGP SurfaceLayer temperature conversion * Polish RRTMGP surface temperature regression * Fixes from review * Update qsurf tests to still check intended lateral behavior * Fix the build: take the ground face of the per-face surface layer in IBSEB PR #3860 turned ERF::m_SurfaceLayer from a single unique_ptr into one entry per domain face, and PR #3960 added ibseb_advance, which reads the surface layer for the wall function beyond neutral. The two merged independently and git merged them cleanly, so development did not compile: ERF_IBSEB.cpp still dereferenced m_SurfaceLayer as a pointer. Index zlo. What the wall function needs is the ground beneath the buildings; a surface layer on a lateral or upper wall says nothing about the stability of the column, and this is the same face that the other zlo-only consumers select (see ERF.cpp:876). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Detect legacy surface-layer inputs from what the user wrote, not the live table With more than one surface-layer face, ERF decides whether to keep reading the historical unqualified erf.most.* / erf.surface_layer.* keys for zlo by asking whether any such key exists. It asked ParmParse directly, which cannot answer the question: ERF_InputSoundingData.H queries most.surf_temp and most.surf_moist with queryAdd while reading the sounding, and queryAdd inserts the default it was handed whenever the key is absent. By the time InitData_post runs, both keys are in the table holding their negative sentinels whether or not the user ever typed them. Every input_sounding run with more than one surface-layer face therefore looked as though it carried legacy inputs. A deck that also wrote the face-qualified spelling -- Tests/test_files/ABL_MOST_Cloudchamber, for one -- hit the "both spellings present" abort and could not start at all; a deck without it silently took the legacy path. Snapshot the inputs table at the top of ReadParameters, which is the last point at which presence in the table still means the user wrote it, and scan the snapshot. This stays correct if another erf.most.* queryAdd is added elsewhere later, which a fix targeting the two known keys would not. ABL_MOST_Cloudchamber is the regression test: it aborts before this change and initializes cleanly after it. Also document the prefixing rule in SurfaceLayer.rst -- unqualified for a lone zlo face, face-qualified once a second face is declared, unqualified applied to zlo with a warning, both spellings an error -- along with the configurations that now abort away from zlo, which the page previously covered only as "not all existing options are supported". Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * fix gmake. (#4028) * more fixes for gmake. (#4029) * SurfaceLayer::fill_qsurf_with_qsat: intersect with the cell box of the nodal heights The loop box is cell-centred and z_phys_nd is nodal, so `gtbx &= z_phys_nd->fabbox()` mixes index types, and with AMReX_ASSERTIONS every run that reaches this function aborts at initialisation (sameType in AMReX_Box.H). The kernel reads nodes i..i+1, j..j+1, k..k+1 through Compute_Zrel_AtCellCenter, so the right limit is the nodal box converted to cells. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Canonical_RANS checks: a range check must fail outside its band Report.check(kind="range") compared the distance outside [lo, hi] against half the band width instead of zero, so a band [lo, hi] accepted [lo - w/2, hi + w/2] with w = hi - lo. The midpoint and half width are still computed for the printed table; only the verdict changes. test_rans_checks.py states what each kind of comparison must decide for values inside and just outside the tolerance or band, 14 cases, and runs as the RANS_Checks_SelfTest entry (labels rans and unit; pure Python, no ERF run). On the old logic it reports 11 of 14 and exits non-zero, disagreeing exactly on the three range cases. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Tests: guard the RANS self-test on a found interpreter, unpin a band edge Two follow-ups to the range-check fix. RANS_Checks_SelfTest is the only test carrying the "unit" label that is not a built binary, and "ctest -L unit" runs in the gcc, macos, ci and windows workflows. Its command is ERF_RANS_PYTHON, which falls back to the bare name "python3" when find_package(Python3 ... QUIET) comes up empty, so a configuration without that executable on PATH -- the MSVC job being the likely one -- would fail the whole unit stage on a test that runs no ERF code. Register it only under Python3_Interpreter_FOUND. check_stable.py's "theta(k=0) - imposed surface theta" band had its lower edge at exactly 0.0, a hard physical boundary. The old half-width slack absorbed 0.75 K below it; with the band now enforced exactly, any transient or round-off putting the first cell a hair below the imposed surface value fails the check. Move the edge to -0.05 K, small enough that a real unstable surface layer is still caught. Note the other --physics range bands (check_stable.py 41, 55, 58; check_convective.py 47, 57; check_hill.py 144, 149; check_hill3d.py 127, 130) have not been rerun under the strict comparison; only the band whose edge sat on a physical zero is changed here. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> --------- Co-authored-by: Debojyoti Ghosh <debojyoti.ghosh@gmail.com> Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com> Co-authored-by: Cole Kendrick <kendrick6@llnl.gov> Co-authored-by: Emma Catherine Ware <ware10@llnl.gov> Co-authored-by: Debojyoti Ghosh <ghosh5@llnl.gov> Co-authored-by: Kyle <kyle.pressel@pnnl.gov> Co-authored-by: Ann Almgren <asalmgren@lbl.gov> Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com>
2 participants
Surface energy balance on immersed-boundary building faces
Adds a prognostic surface energy balance on the faces of buildings represented by immersed forcing (
erf.ibseb.*, off by default). Every wall and roof face carries a skin temperature and a conduction slab and closeseach step, with the sensible flux deposited into the adjacent fluid cell. Each term has a verification case under
Exec/CanonicalTests/SEB/(one and four ranks, checked against independent Python implementations or analytic solutions), documented inDocs/sphinx_doc/theory/ImmersedBoundarySEB.rstandDocs/sphinx_doc/Inputs.rst. CI runs the balance through the CTestIBSEB_Cube(the PrognosticSkin cube for 40 steps against a gold plotfile) and the immersed-boundary-aware PBL schemes through the tiling-parity testsPBL_IBAware_MRF_TilingandPBL_IBAware_YSUNew_Tiling. The balance needserf.buildings_type = ImmersedForcingon a uniform vertical grid without regridding; terrain by immersed forcing is rejected.Q_exthook for external radiation; bounds as inputsDemo
The isolated building over a solstice day (IsolatedBuilding): the roof lags the sun by two hours on 30 cm of concrete, conduction reverses at 18:10, and every face radiates below the air at night.
The building set through a morning (BuildingSet): four height-map buildings, three materials; the two identical timber blocks track each other to 0.2 K, the brick cube and the concrete slab lag by their heavier walls.
Sunrise over a cube (the PrognosticSkin case), the test that caught a mirrored solar azimuth in the provider:
Also in this PR, both opt-in, both bit-identical when off
erf.if_snap_partial_cells: the immersed forcing on a height-map building grows a vertical checkerboard on the sliver cells of the reader's one-cell ramp over about two hours (with the balance off, a no-slip ground, the forcing outside the substeps and a coarser small-cell threshold as well). The switch reads the blanking snapped to whole cells and uses the point-implicit drag: a momentum face between a solid and a fluid cell gets the drag toward zero (no penetration), a face between two solid cells the roof or wall law of its row or the interior drag; buildings only, terrain by immersed forcing keeps its fractions. RegtestExec/RegTests/ImmersedForcingTest/PartialCells.erf.pbl_ib_aware(MRF and YSUNew only): the schemes measure height from the terrain surface and the ground surface layer evaluates u*, theta* and L on cells inside a building, so over an immersed building MRF fills the domain with NaN and YSUNew drives the density negative at its second step. The switch makes each column's surface the first fluid cell above the solid; the stored PBL height stays the absolute height in both schemes, the PBL height smoothing runs on that height and keeps every column at or above its own floor. RegtestsExec/RegTests/ImmersedForcingTest/PBL_IBAware, CTestsPBL_IBAware_MRF_Tiling,PBL_IBAware_YSUNew_Tiling,PBL_IBAware_MRF_Smoothing.Limitations and the effort to remove them
"Session" is one working day with this workflow (implement, regtest, document).
two_streamprovider behind the existingerf.ibseb.radiationinput from the now-merged two-stream model (#3950, #4004); the balance keeps its own solar geometry (Spencer's declination and equation of time, the azimuth the CESM orbital form of the radiation models does not give) and only asks for the sun vector, the beam, the horizontal diffuse and the sky longwaveParallelForover device vectors) but have not run on oneQ_ext, ignition from the skin temperature) is not hereReview rounds
utc_offset_hoursinput removed, range checks on latitude, longitude and day of year, theibseb_*plotfile variables documented.PBL_IBAware_MRF_Smoothinggold test; the report and CSV step numbers as completed steps; a restart with anothern_slab_layersrejected against the checkpoint's field width.erf.if_z0changes the flow); the checkpoint field on column blocks around the buildings with a transfer layer so a restart works on any rank count (FaceStorage restarts on one rank from a four-rank checkpoint); the ray cast's column map over the built bounding box; the bulk Richardson shear from the wind vector; the slab response guarded like its sibling.Findings to report separately
🤖 Generated with Claude Code