ERF_IBFaceSet.H

Go to the documentation of this file.

#ifndef ERF_IB_FACE_SET_H
#define ERF_IB_FACE_SET_H
 
#include <AMReX_MultiFab.H>
#include <AMReX_BoxArray.H>
#include <AMReX_DistributionMapping.H>
#include <AMReX_Geometry.H>
#include <AMReX_GpuContainers.H>
#include <cmath>
#include <vector>
#include "ERF_IBSEBParams.H"
#include "ERF_IBSEBMaterials.H"
#include <ERF_TwoStreamCanopyForcing.H>
 
/**
 * \file ERF_IBFaceSet.H
 * \brief Compact list of the wall faces of resolved buildings on one AMR level.
 *
 * **What a face is.** The immersed forcing marks every cell of a level as
 * fluid or solid through the cell-centred blanking ``terrain_blanking``
 * (one minus the embedded-boundary volume fraction, so 0 in the air and 1
 * inside a building). A wall face is the cell face between a fluid cell
 * (blanking below 0.5) and a solid cell (blanking at or above 0.5). Walls
 * are x or y faces, roofs are z faces with the fluid above the solid.
 *
 * **Face convention.** A face is identified by its *fluid* cell ``(i, j, k)``,
 * its direction ``dir`` (0 = x, 1 = y, 2 = z) and its ``side`` (-1 when the
 * solid neighbour is at ``i-1`` / ``j-1`` / ``k-1``, +1 when it is at
 * ``i+1`` / ``j+1`` / ``k+1``). The outward unit normal of the wall, pointing
 * into the fluid, is therefore ``-side`` along ``dir``. The fluid cell is where
 * every atmospheric input of the balance (wind, temperature, humidity,
 * radiation) is read and where the wall's heat flux is deposited.
 *
 * **Ownership and layout.** Each face is stored exactly once, on the rank
 * that owns its fluid cell, so the per-step update of the balance is a kernel
 * over this rank's faces with no communication; only the reports reduce
 * across ranks. The faces are a struct of device arrays (one
 * ``amrex::Gpu::DeviceVector`` per quantity, all indexed by the face number)
 * and are stored contiguously per local fab in MFIter order:
 * ``fab_start()[mfi.LocalIndex()]`` to ``fab_start()[mfi.LocalIndex()+1]``
 * are the faces whose fluid cell is in that fab, which is what lets a kernel
 * scatter into that fab's arrays without a search.
 *
 * **Building ids.** The solid columns of the blanking (any solid cell in the
 * column) are labelled by 4-connectivity in the horizontal and numbered in
 * scan order over ``(i, j)``, identically on every rank from a reduced column
 * mask; a face carries the id of the column of its solid neighbour. Ids are
 * 1-based, 0 is never assigned.
 *
 * **Output.** The list is never written directly. scatter_diagnostics() puts
 * the number of faces and their mean skin temperature into cell-centred
 * fields for the plotfile; save_state() / load_state() move the skin
 * temperature, the sensible flux and the slab temperatures through a
 * cell-centred field with one slot per face of a cell (see state_ncomp())
 * for the checkpoint. On restart the list is rebuilt from the blanking and
 * refilled from that field, so a restart does not depend on the rank count
 * or on the order the faces were found in.
 *
 * **Per-step order.** build() detects and stores the faces and
 * compute_view_fractions() fills ``d_f_sky`` / ``d_f_ground`` / ``d_f_bldg``
 * once at initialisation; every step then runs compute_shortwave(),
 * compute_longwave(), compute_sensible() and, with ``erf.ibseb.prognostic``,
 * solve_balance(), which finds the skin temperature closing the balance,
 * advances the slab with it and rewrites the fluxes at the solution; without
 * it compute_ground() advances the slab under the fixed skin. ``d_Q_ext`` is
 * the hook through which an external incident flux (a fire's radiation)
 * enters without the balance knowing its source: whoever owns it writes the
 * array before solve_balance() runs.
 */
class IBFaceSet
{
public:
    /**
     * Construct an empty face set for one level; build() fills it.
     * @param params Inputs of the balance (copied; the set never reads ParmParse).
     * @param lev    AMR level the set belongs to, for reports only.
     */
    IBFaceSet (const IBSEBParams& params, int lev) : m_params(params), m_lev(lev) {}
 
    /**
     * Detect this rank's faces from the cell-centred blanking and allocate
     * every per-face array.
     *
     * Runs on the host once at initialisation: the blanking is copied fab by
     * fab (a pinned copy on GPU builds), the valid cells are scanned, and the
     * arrays are uploaded to the device at the end. Must be called after the
     * blanking has been built and its ghost cells filled, since the solid
     * neighbour of a fluid cell on a box edge lives in a ghost cell. Stops
     * when a solid cell outside the level's grids neighbours a cell of the
     * level, fluid or solid: a building across the edge of a refined level.
     *
     * @param blanking Cell-centred blanking of the level (needs at least one
     *                 ghost cell, filled).
     * @param geom     Geometry of the level: cell sizes for the face areas,
     *                 domain and periodicity for the neighbour tests.
     */
    void build (const amrex::MultiFab& blanking, const amrex::Geometry& geom);
 
    /**
     * Add the buildings this refined level does not hold to the column map of
     * its ray cast, from the map of the level below.
     *
     * build() takes the map from the level's own cells, so on a refined level
     * that covers part of the domain the rays of the shortwave and the view
     * fractions would pass through every building outside it. Every column of
     * the coarser map whose top cell this level's grids do not cover
     * (ERF::ibseb_check_refined_levels() has made sure a building is wholly
     * inside or wholly outside) is copied onto the ratio x ratio columns of
     * this level beneath it, at the same height, and the map grows to the
     * bounding box of both. Only the ray cast reads the result: the building
     * ids, the face heights of the wall function and the reports keep the
     * level's own buildings. Call after build() and before
     * compute_view_fractions(); the coarser set must have had its own call
     * already, so buildings outside several levels pass down the hierarchy.
     *
     * @param coarser   Face set of the level below.
     * @param ratio     Refinement ratio from the level below to this one.
     * @param grids_crse This level's grids coarsened by ``ratio``.
     */
    void add_outside_occluders (const IBFaceSet& coarser, const amrex::IntVect& ratio,
                                const amrex::BoxArray& grids_crse);
 
    /**
     * Give the sun of the next compute_shortwave() with
     * ``erf.ibseb.sun_mode = two_stream``: the two-stream radiation's
     * declination and hour angle at that time, its site latitude and its
     * top-of-atmosphere irradiance (the date's distance factor included).
     * ERF::ibseb_set_two_stream_sun() calls it before every call; the call
     * uses it up, so a shortwave without a fresh sun stops rather than reusing
     * an old one.
     */
    void set_two_stream_sun (amrex::Real declination, amrex::Real hour_angle,
                             amrex::Real latitude_deg, amrex::Real toa_irradiance)
    {
        m_ts_sun = {true, declination, hour_angle, latitude_deg, toa_irradiance};
    }
 
    /**
     * Number this refined level's buildings as level 0 numbers them.
     *
     * Each level labels its buildings from its own cells, so the ids of a
     * refined level differ from level 0's wherever the levels see different
     * buildings: a building outside the refined level, or two that level 0
     * merges into one. Inputs that go by building (erf.ibseb.material_by_building)
     * and the report's ``building_level0`` column use level 0's numbering, so
     * every building of this level is looked up in the labels of the level
     * below, and through that level's own map in level 0's: each of its
     * columns votes with the id under it and the majority wins. A building
     * only this level resolves maps to 0 (the default material), with a
     * warning.
     * Call after build() of both levels; the level below keeps its labels until
     * release_labels().
     *
     * @param coarser Face set of the level below, already mapped.
     * @param ratio   Refinement ratio from the level below to this one.
     */
    void map_buildings_to_level0 (const IBFaceSet& coarser, const amrex::IntVect& ratio);
 
    /** Level 0's number for building ``b`` of this level (1-based; 0 for one only this level resolves). */
    int building_level0 (int b) const { return m_to_lev0.at(b); }
 
    /** Free the column labels once the level above has been mapped (map_buildings_to_level0()). */
    void release_labels () { std::vector<int>().swap(m_lab); }
 
    /** Number of faces owned by this rank (not the global count; see report()). */
    int n_faces () const { return m_nface; }
    /** Number of buildings of the level, the same on every rank. */
    int n_buildings () const { return m_nbld; }
    /** Conduction layers per face, from ``erf.ibseb.n_slab_layers``. */
    int n_layers () const { return m_params.n_slab_layers; }
 
    /**
     * Number of components of the checkpoint state field.
     *
     * The field carries n_slots() slots of ``2 + n_layers()`` values per cell:
     * component ``slot * (2 + n_layers())`` is the skin temperature of that
     * slot, the next its sensible flux (read as the previous step's by the
     * convective velocity scale), and the following ``n_layers()`` components
     * are its slab layers.
     *
     * A cell carries one slot per face it owns, numbered in the ``(dir, side)``
     * order in which build() finds them, so the width is the largest face
     * count on any one cell of the level rather than the six a cell could hold
     * in principle. Around ordinary buildings that maximum is two or three (an
     * outside edge or corner); six needs a one-cell slot. build() rebuilds the
     * numbering from the same blanking on restart, so the slots line up
     * whatever the rank count.
     */
    int state_ncomp () const { return m_nslots * (2 + n_layers()); }
 
    /** Slots per cell of the checkpoint field: the largest face count on a cell of the level. */
    int n_slots () const { return m_nslots; }
 
    /**
     * Boxes of the checkpoint field: 4 x 4 column blocks clipped to the
     * k-range that actually owns faces, so the field follows the shell of
     * face-owning cells rather than the built volume. A block over the
     * interior of a wide footprint keeps only the layer above the roof, while
     * a block on a wall still spans its height. Built identically on every
     * rank by build() from a reduced block k-range map; the same buildings
     * give the same boxes on restart, whatever the rank count. Empty when the
     * level has no faces.
     */
    const amrex::BoxArray& state_boxarray () const { return m_state_ba; }
    /** Whether the level has a checkpoint field to write (any faces at all). */
    bool has_state () const { return !m_state_ba.empty(); }
    /** A field on state_boxarray() with state_ncomp() components, for the checkpoint. */
    amrex::MultiFab make_state () const;
 
    /**
     * View fractions: sample a cosine-weighted hemisphere around
     * every face's outward normal with ``view_n_az x view_n_el`` rays through
     * the column walk and count where they end: sky, ground or a building.
     * The three fractions sum to one and are view factors. Called once at
     * initialisation; the geometry is static.
     */
    void compute_view_fractions ();
 
    /**
     * Longwave of the current step. Incoming on a face:
     * ``f_sky LW_sky + f_ground LW_ground + f_bldg sigma T_skin^4``,
     * the sky term fixed or gray (``sky_emissivity sigma T_air^4`` with the
     * air temperature of the face's fluid cell) with the ground term
     * ``eps_g sigma T_g^4``, or with ``lw_mode = two_stream`` the longwave
     * down and up at the face's height in the face's column
     * (``canopy``); the building term is the isothermal-surroundings
     * approximation. Net: ``eps (LW_in - sigma T_skin^4)``, positive into the face.
     * @param cons   Conserved state of the level, for the air temperature.
     * @param canopy With ``lw_mode = two_stream``, this step's two-stream view of the
     *               level (ERF_TwoStreamCanopyForcing.H), or an empty one before the
     *               first sweep, when the incoming and net longwave are zero.
     */
    void compute_longwave (const amrex::MultiFab& cons, const TwoStreamCanopyView& canopy = TwoStreamCanopyView{});
 
    /**
     * Assign every face its material properties: from the material
     * library by building id, or the uniform inputs. Called once after
     * build(); fills d_mat, d_albedo, d_emis, d_k, d_rhocp, d_thick and
     * initialises the slab layers to the skin temperature.
     */
    void assign_materials ();
 
    /**
     * Ground heat flux of the current step: advance every face's
     * slab by dt with the skin temperature at the top and the interior
     * temperature at the bottom (implicit, ERF_IBSEBSlab.H) and store the
     * conduction into the slab, G, positive into the wall.
     */
    void compute_ground (amrex::Real dt);
 
    /**
     * Prognostic balance of the current step. Per face: the
     * slab's linear response to the skin temperature for this step, the
     * Newton solve of ERF_IBSEBBalance.H with the shortwave, the external
     * flux, the sky and ground longwave, the frozen wall-function
     * coefficient and the latent flux, the slab advanced with the solution,
     * and the fluxes (LW_in, LW_net, H, G) rewritten at the solution. Fills
     * d_resid (|balance| at the end, non-zero only at a bound or the
     * iteration cap) and d_niter. Needs compute_shortwave(),
     * compute_longwave() and compute_sensible() of the same step first.
     * @param dt Atmospheric time step [s], the slab's step.
     */
    void solve_balance (amrex::Real dt);
 
    /**
     * Sensible heat of the current step through a wall function on
     * every face: 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 gives theta* with z0h_wall, and
     * H = rho c_p u* theta*, positive out of the face. Neutral, or with the
     * surface layer's stability functions on roofs when asked. The latent
     * flux is not modelled and stays zero.
     *
     * Beyond neutral: with ``convective_velocity = deardorff`` the
     * wind is ``sqrt(U_tan^2 + (beta w*)^2)`` with w* from the previous
     * step's H (zero when H is into the face) and a depth that is the mixed
     * layer above a roof or the building height for a wall; with
     * ``stability_correction`` the roofs iterate the similarity functions on
     * the face's own Obukhov length, seeded from the ground surface layer's
     * 2D field at the face's column when given.
     *
     * Fills d_H, d_ustar, d_U_tan, d_rho, d_theta_air, d_H_coeff (the
     * coefficient of the skin-to-air potential-temperature difference that
     * the balance keeps frozen), d_w_star, d_olen and d_z_i; T_air was read
     * by compute_longwave(), which must run first in a step.
     * @param cons  Conserved state of the level (density, rho theta).
     * @param xvel, yvel, zvel  Face-centred velocities of the level.
     * @param c_p   Specific heat of air [J/kg/K].
     * @param olen_ground  Ground surface layer's Obukhov length (2D field on
     *                     the level's box array collapsed to k = 0), or null.
     * @param pblh_ground  Its boundary-layer height field, or null.
     * @param z_i_bulk     Mixed-layer depth of the level from the bulk
     *                     Richardson diagnostic or the fixed input [m].
     */
    void compute_sensible (const amrex::MultiFab& cons, const amrex::MultiFab& xvel,
                           const amrex::MultiFab& yvel, const amrex::MultiFab& zvel, amrex::Real c_p,
                           const amrex::MultiFab* olen_ground = nullptr,
                           const amrex::MultiFab* pblh_ground = nullptr,
                           amrex::Real z_i_bulk = 0.0);
 
    /**
     * Add the face heat flux to the temperature equation: every face
     * deposits ``H A / (c_p V Pi)`` into the rho-theta source of its fluid
     * cell (Pi the Exner function, so the flux heats the cell by H A / (c_p V)
     * in temperature). Called after the sources are rebuilt at every slow
     * stage; adds, never overwrites, and uses atomic adds because a corner
     * cell receives several faces.
     */
    void add_heat_flux_to_source (amrex::MultiFab& source, const amrex::MultiFab& cons,
                                  const amrex::Geometry& geom, amrex::Real c_p, amrex::Real rdOcp) const;
 
    /**
     * Shortwave of the current step: the sun's position (sun_mode), the
     * irradiances of the prescribed provider or, with ``radiation =
     * two_stream``, those of the face's own column at its height (``canopy``),
     * then per face
     * the direct beam on the face (zero when the ray toward the sun hits a
     * building), the diffuse light through the view fractions, and the
     * absorbed sum.
     *
     * Fills d_shadow, d_SW_direct_in, d_SW_diffuse_in and d_SW_abs; stores
     * the sun vector and irradiances for the reports (with two_stream the
     * means over all the level's faces of the direct-normal and diffuse
     * irradiance each takes, written to every row of the level).
     * @param time   Simulation time [s]; with sun_mode = solar it is added to
     *               time_zero_utc_s (with sun_mode = two_stream the sun comes
     *               from set_two_stream_sun()).
     * @param canopy With ``radiation = two_stream``, this step's two-stream view of
     *               the level (see compute_longwave()), or an empty one before the
     *               first sweep, when every face gets no shortwave.
     */
    void compute_shortwave (amrex::Real time, const TwoStreamCanopyView& canopy = TwoStreamCanopyView{});
 
    /**
     * The highest interface of the level's vertical grid that any of its faces
     * samples with radiation = two_stream (a roof the bottom of its fluid cell, a
     * ceiling and a wall the top of its cell too), counted from the domain bottom, over all
     * ranks; -1 when the level has no faces. The two-stream beam the faces read
     * holds the interfaces up to there.
     */
    int top_sample_interface () const { return m_top_sample; }
 
    /** Whether the level has faces on any rank. */
    bool has_faces () const { return m_top_sample >= 0; }
 
    /**
     * Stop unless the fields of ``canopy`` (this level's two-stream view) lie on
     * the columns of ``cons``, with the same distribution: the two-stream kernels
     * read the column of a face from the fabs with the local index of the face's
     * own fab.
     */
    void check_canopy_layout (const amrex::MultiFab& cons, const TwoStreamCanopyView& canopy) const;
 
    // The two-stream branches of compute_shortwave() (once the sun is placed) and
    // compute_longwave(). Public only because a CUDA extended lambda cannot sit in
    // a private member function; call the two functions above instead.
    void compute_shortwave_two_stream (const TwoStreamCanopyView& canopy);
    void compute_longwave_two_stream (const amrex::MultiFab& cons, const TwoStreamCanopyView& canopy);
 
    /**
     * Scatter one per-face array into a cell-centred field as the mean over
     * the faces touching each cell (zero where there are none).
     */
    void scatter_field (const amrex::Gpu::DeviceVector<amrex::Real>& v, amrex::MultiFab& out) const;
 
    /**
     * Write every face of this rank to ``<prefix>.rank<N>.csv`` with its
     * geometry, view fractions, shadow flag, shortwave terms and skin
     * temperature, for the regression tests and for users who want the raw list.
     */
    void dump_faces (const std::string& prefix) const;
 
    /**
     * Scatter the faces into two cell-centred diagnostics for the plotfile.
     *
     * @param[out] nfaces Number of wall faces touching each fluid cell (0 away
     *                    from buildings; up to 3 at an outside corner, 6 in a
     *                    one-cell slot).
     * @param[out] tskin  Mean skin temperature of those faces, 0 where there
     *                    are none.
     *
     * Both fields must be on the level's box array with no ghost cells. Faces
     * of one fab are accumulated with atomic adds because a corner cell
     * receives several faces.
     */
    void scatter_diagnostics (amrex::MultiFab& nfaces, amrex::MultiFab& tskin) const;
 
    /**
     * Write the face state (skin temperature, sensible flux and slab
     * temperatures) into a field from make_state(), for the checkpoint.
     *
     * The faces of this rank are written into a transfer field on this
     * rank's grids cut by the state boxes, which a ParallelCopy moves onto
     * the state's own distribution. Each face writes its own slot (see
     * state_ncomp()), so no two faces touch the same component of the same
     * cell and no atomics are needed; cells without a face in a slot keep
     * zero there.
     */
    void save_state (amrex::MultiFab& state) const;
 
    /**
     * Read the face state (skin, flux, slab) back from a field written by
     * save_state(), after build() has recreated the list from the blanking:
     * the reverse ParallelCopy onto the transfer field, then each face reads
     * its slot. Only the slots that correspond to a face of this list are read.
     */
    void load_state (const amrex::MultiFab& state);
 
    /**
     * Print the ``[IBSEB]`` summary line and, when asked, append one row per
     * building to ``erf.ibseb.csv_file``.
     *
     * The summary carries the global face counts per direction, the number of
     * buildings, the total face area and the range of the skin temperature;
     * the CSV rows carry per building the face count, the area, the
     * area-weighted mean fluxes and skin temperature with its range, the
     * largest balance residual, and the sun the step used, from its start
     * (zenith, azimuth, direct-normal and diffuse irradiance), on a row
     * stamped with the time at the end of the step. All numbers are reduced over the
     * ranks; the file is written by the I/O rank only, with a header only
     * when the file does not exist yet, so a restarted run appends. Every
     * level appends to the same file, its rows marked by the level column.
     * With ``erf.ibseb.dump_faces_file`` the rows come with a face dump,
     * dump_faces() to ``<prefix>[.levN][.stepNNNNNN].rank<N>.csv``: a refined
     * level N adds ``.levN`` so it does not overwrite the level below, and
     * ``erf.ibseb.dump_faces_tag_step`` the step.
     *
     * @param time      Simulation time [s] written to the rows.
     * @param step      Step number written to the rows.
     * @param write_csv Whether to append rows, or only print the summary.
     */
    void report (amrex::Real time, int step, bool write_csv) const;
 
    // ---- Per-face device arrays, indexed by face number ----------------------
    // Static geometry, set by build():
    amrex::Gpu::DeviceVector<int>         d_i, d_j, d_k;   ///< Fluid cell of the face
    amrex::Gpu::DeviceVector<int>         d_dir;           ///< Face direction: 0 x, 1 y, 2 z
    amrex::Gpu::DeviceVector<int>         d_side;          ///< Solid neighbour at -1 (low) or +1 (high) side
    amrex::Gpu::DeviceVector<int>         d_slot;          ///< Slot of the face within its cell's checkpoint record
    amrex::Gpu::DeviceVector<int>         d_bid;           ///< Building id (1-based)
    amrex::Gpu::DeviceVector<int>         d_mat;           ///< Material id (0 = uniform inputs)
    amrex::Gpu::DeviceVector<amrex::Real> d_albedo, d_emis;///< Optical properties of the face
    amrex::Gpu::DeviceVector<amrex::Real> d_kth, d_rhocp, d_thick; ///< Slab conductivity [W/m/K], heat capacity [J/m3/K], thickness [m]
    amrex::Gpu::DeviceVector<amrex::Real> d_area;          ///< Face area [m2]
    amrex::Gpu::DeviceVector<amrex::Real> d_xf, d_yf, d_zf;///< Face centre [m]
    // State, advanced by the balance and checkpointed:
    amrex::Gpu::DeviceVector<amrex::Real> d_T_skin;        ///< Skin temperature [K]
    amrex::Gpu::DeviceVector<amrex::Real> d_T_slab;        ///< Slab layers [K], face-major: [face * n_layers + layer]
    // Static radiative geometry:
    amrex::Gpu::DeviceVector<amrex::Real> d_f_sky, d_f_ground, d_f_bldg;  ///< View fractions, sum to 1
    // Fluxes of the current step [W/m2]: SW_abs and LW_net positive into the
    // face, H, LE and G positive out of it:
    amrex::Gpu::DeviceVector<amrex::Real> d_SW_abs, d_LW_net, d_H, d_LE, d_G;
    amrex::Gpu::DeviceVector<amrex::Real> d_Q_ext;         ///< External incident flux [W/m2], e.g. fire radiation
    // Shortwave detail, all [W/m2] except the flag:
    amrex::Gpu::DeviceVector<amrex::Real> d_shadow;        ///< 1 when the ray to the sun hits a building, else 0
    amrex::Gpu::DeviceVector<amrex::Real> d_SW_direct_in;  ///< Direct beam incident on the face
    amrex::Gpu::DeviceVector<amrex::Real> d_SW_diffuse_in; ///< Sky and ground-reflected diffuse incident
    // Longwave detail [W/m2] and the air temperature read [K]:
    amrex::Gpu::DeviceVector<amrex::Real> d_LW_down_in;   ///< Incoming longwave on the face (sky, ground and walls)
    amrex::Gpu::DeviceVector<amrex::Real> d_T_air;        ///< Air temperature of the fluid cell
    // Sensible heat detail:
    amrex::Gpu::DeviceVector<amrex::Real> d_theta_air;    ///< Potential temperature of the fluid cell [K]
    amrex::Gpu::DeviceVector<amrex::Real> d_rho;          ///< Density of the fluid cell [kg/m3]
    amrex::Gpu::DeviceVector<amrex::Real> d_U_tan;        ///< Tangential wind speed at the fluid cell [m/s]
    amrex::Gpu::DeviceVector<amrex::Real> d_ustar;        ///< Friction velocity of the wall function [m/s]
    amrex::Gpu::DeviceVector<amrex::Real> d_H_coeff;      ///< H = H_coeff (theta_skin - theta_air) [W/m2/K]
    // Wall function beyond neutral:
    amrex::Gpu::DeviceVector<amrex::Real> d_hbld;         ///< Height of the face's building column [m]
    amrex::Gpu::DeviceVector<amrex::Real> d_w_star;       ///< Convective velocity scale of the step [m/s]
    amrex::Gpu::DeviceVector<amrex::Real> d_olen;         ///< Obukhov length used [m] (1e30 when neutral)
    amrex::Gpu::DeviceVector<amrex::Real> d_z_i;          ///< Depth in w* [m]: mixed layer above a roof, building height for a wall
    // Prognostic balance detail:
    amrex::Gpu::DeviceVector<amrex::Real> d_LW_ext;       ///< Incoming longwave from sky and ground only [W/m2]
    amrex::Gpu::DeviceVector<amrex::Real> d_resid;        ///< |balance| after the solve [W/m2]
    amrex::Gpu::DeviceVector<int>         d_niter;        ///< Newton iterations of the last solve
    // Column tops over the bounding box of the built columns (origin
    // m_col_i0, m_col_j0 in 0-based domain columns, m_col_nx x m_col_ny
    // values, index (ci - i0) * m_col_ny + (cj - j0)), replicated on every
    // rank for the ray cast: 4 bytes per built column per rank. Each entry is
    // the domain-relative k of the highest solid cell of its column, -1 where
    // the column is fluid; ibseb::column_top() turns it into a height, and
    // columns outside the box are open ground.
    amrex::Gpu::DeviceVector<int> d_col_top;
    int m_col_i0 = 0, m_col_j0 = 0, m_col_nx = 0, m_col_ny = 0;
 
    // Checkpoint layout (see state_boxarray()) and the transfer layer of
    // save_state() / load_state(): this rank's grids cut by the state boxes,
    // gathered into one global BoxArray owned box by box by the rank whose
    // grid it came from, with the local grid index each box came from.
    amrex::BoxArray m_state_ba;
    amrex::BoxArray m_xfer_ba;
    amrex::DistributionMapping m_xfer_dm;
    std::vector<int> m_xfer_src;
    amrex::Long m_domain_cells = 0;            ///< Cells of the level, for the debug summary of the field
 
    /// Sun and irradiances of the last compute_shortwave() call, for reports.
    struct SunState {
        amrex::Real zenith = 0.0, azimuth = 0.0;        ///< [rad]
        amrex::Real sx = 0.0, sy = 0.0, sz = 1.0;       ///< Unit vector toward the sun
        amrex::Real dni = 0.0;                          ///< Direct-normal irradiance [W/m2]
        amrex::Real diffuse_h = 0.0;                    ///< Diffuse on a horizontal surface [W/m2]
    };
    const SunState& sun () const { return m_sun; }
 
    /**
     * Print the ``[IBSEB DEBUG]`` description of the set: the inputs, every
     * rank's face count and per-fab ranges, the buildings with their footprint
     * columns and bounding boxes, and the device memory of the arrays.
     * Called once from build() when ``erf.ibseb.debug`` is set.
     */
    void print_debug_summary () const;
 
    /** Add the wall-clock time of one per-step update of this rank, for the cost line of report(). */
    void add_cost (double seconds) { m_cost_s += seconds; ++m_cost_n; }
    /** Record the wall-clock time of the initialisation (build and view fractions). */
    void set_init_cost (double seconds) { m_init_cost_s = seconds; }
 
    /**
     * First face of each local fab in MFIter order; size number of fabs + 1.
     * Faces ``fab_start()[n]`` to ``fab_start()[n+1]-1`` belong to the fab with
     * local index ``n``, which is what the scatter kernels rely on.
     */
    const std::vector<int>& fab_start () const { return m_fab_start; }
 
private:
    IBSEBParams m_params;          ///< Copy of the inputs
    int m_lev   = 0;               ///< AMR level, for reports
    int m_nface = 0;               ///< Faces owned by this rank
    int m_nslots = 0;              ///< Slots per cell of the checkpoint field (largest face count on a cell)
    int m_top_sample = -1;         ///< Highest interface a face samples, over all ranks (-1: no faces)
    int m_nbld  = 0;               ///< Buildings of the level (global)
    std::vector<int> m_fab_start;  ///< See fab_start()
 
    // Static totals over all ranks, computed once in build(); building index 0 unused.
    std::vector<amrex::Long> m_bld_nface;      ///< Faces per building
    std::vector<amrex::Real> m_bld_area;       ///< Face area per building [m2]
    amrex::Long m_nface_dir[3] = {0, 0, 0};    ///< Faces per direction
    amrex::Real m_area_total = 0.0;            ///< Total face area [m2]
    // Footprint of each building (global): columns, and bounding box in cell indices.
    std::vector<amrex::Long> m_bld_ncol;
    std::vector<int>  m_bld_ilo, m_bld_ihi, m_bld_jlo, m_bld_jhi;
    // Level-0 numbering (map_buildings_to_level0()): the level-0 id of each
    // building of this level (identity on level 0), and the column labels over
    // the level's own built bounding box, kept until the level above has been
    // mapped.
    std::vector<int>  m_to_lev0;
    std::vector<int>  m_lab;
    int m_lab_i0 = 0, m_lab_j0 = 0, m_lab_nx = 0, m_lab_ny = 0;
 
    // Geometry kept for the ray cast and the reports.
    int  m_nx = 0, m_ny = 0;                   ///< Columns of the domain
    amrex::IntVect m_dom_lo{AMREX_D_DECL(0, 0, 0)}; ///< Low corner of the level's domain
    amrex::Real m_x_lo = 0.0, m_y_lo = 0.0;    ///< Domain origin [m]
    amrex::Real m_dx[3] = {0.0, 0.0, 0.0};     ///< Cell sizes [m]
    bool m_per_x = false, m_per_y = false;     ///< Periodicity in x and y
    amrex::Real m_col_top_max = 0.0;           ///< Tallest column [m]
    amrex::Real m_z_ground = 0.0;              ///< Ground height [m] (flat, the domain bottom)
    amrex::Real m_max_path = 0.0;              ///< Ray path cap [m]
    SunState m_sun;
    /// The two-stream sun of the next compute_shortwave() (sun_mode = two_stream).
    struct TwoStreamSun { bool set = false; amrex::Real decl = 0.0, hour_angle = 0.0, lat_deg = 0.0, S0 = 0.0; };
    TwoStreamSun m_ts_sun;
    double m_cost_s = 0.0;        ///< Wall-clock seconds spent in the per-step update on this rank
    amrex::Long m_cost_n = 0;     ///< Number of per-step updates timed
    double m_init_cost_s = 0.0;   ///< Wall-clock seconds of build() and compute_view_fractions()
};
 
#endif