Surface energy balance on immersed-boundary building faces (erf.ibseb), with immersed-boundary fixes for the forcing and the MRF/YSUNew schemes - #3960
Open
hgopalan wants to merge 21 commits into
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>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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. The work was done in eight phases, each with a regtest 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; the living plan with every design decision and finding isSource/ImmersedBoundarySEB/IBSEB_DEVELOPMENT.md.Q_exthook for external radiation; bounds as inputsDemo
The isolated building over a solstice day (phase 7): 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 (phase 8): 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 (phase 6 regtest), 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. 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. RegtestExec/RegTests/ImmersedForcingTest/PBL_IBAware.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; the balance 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 hereFindings to report separately
🤖 Generated with Claude Code