ERF_IBSEBParams.H

Go to the documentation of this file.

#ifndef ERF_IBSEB_PARAMS_H
#define ERF_IBSEB_PARAMS_H
 
#include <AMReX_ParmParse.H>
#include <AMReX_REAL.H>
#include <AMReX.H>
#include <string>
#include <vector>
 
/**
 * \file ERF_IBSEBParams.H
 * \brief Inputs of the immersed-boundary surface energy balance (``erf.ibseb.*``).
 *
 * The balance runs on the faces of buildings represented by immersed forcing
 * (``erf.buildings_type = ImmersedForcing``). Every input lives under the
 * ``erf.ibseb`` prefix and everything is off unless ``erf.ibseb.enable`` is set,
 * in which case ERF::init_ibseb() builds one IBFaceSet per level.
 *
 * The parameters are parsed once, in ERF::ReadParameters(), and copied into
 * each IBFaceSet, so a face set never reads ParmParse itself. Every input of
 * the balance lives here with its default, so that the whole input surface
 * is in one place; the inputs are documented in Docs/sphinx_doc/Inputs.rst.
 */
struct IBSEBParams
{
    bool        enable        = false;                 ///< Master switch
    int         n_slab_layers = 4;                     ///< Conduction layers per face
    amrex::Real T_skin_init   = 300.0;                 ///< Initial skin temperature [K]
    amrex::Real T_interior    = 293.0;                 ///< Interior (deep) temperature [K]
    std::string csv_file      = "ibseb_buildings.csv"; ///< Per-building report
    int         csv_int       = 100;                   ///< Steps between report rows; <= 0 disables
    bool        debug         = false;                 ///< [IBSEB DEBUG] prints: set description at build, per-step summary, per-building rows
    std::string dump_faces_file;                       ///< Per-rank face dump at report time (empty = off); a refined level N adds ".levN"
    bool        dump_faces_tag_step = false;           ///< Name each dump by its step (keeps every one) instead of overwriting
 
    // ---- Radiation provider ------------------------------------------------
    // Where the downwelling shortwave and longwave come from.
    // "prescribed" is built in: clear-sky formulas with a fixed or computed
    // sun. "two_stream": the column fluxes of the TwoStream radiation model
    // (erf.radiation_model = TwoStream) of each face's own column at the face's
    // height (ERF_TwoStreamCanopyForcing.H): the direct beam (as a direct-normal
    // irradiance, the beam over the cosine of the zenith angle), the diffuse
    // light from the sky and the light the ground reflects; and the longwave
    // with lw_mode = two_stream, its default there. The sun's position still
    // comes from sun_mode.
    std::string radiation     = "prescribed";
    // "fixed": the sun stays at sun_zenith_deg / sun_azimuth_deg; with the
    // prescribed provider, with the direct-normal and diffuse irradiances given
    // below (analytic tests).
    // "solar": sun from the site and time; irradiances from the clear-sky
    // formulas.
    // "two_stream": the sun of the two-stream radiation (erf.radiation_model =
    // TwoStream following start_datetime): its declination, hour angle (no
    // equation of time), site (erf.rad_cons_lat / lon) and top-of-atmosphere
    // irradiance; the clear-sky formulas of "solar" then give the irradiances.
    std::string sun_mode      = "fixed";
    amrex::Real sun_zenith_deg    = 45.0;   ///< fixed: zenith angle [deg]
    amrex::Real sun_azimuth_deg   = 180.0;  ///< fixed: azimuth clockwise from north [deg]; 180 = sun in the south
    amrex::Real sw_direct_normal  = 800.0;  ///< fixed: direct-normal irradiance [W/m2]
    amrex::Real sw_diffuse        = 100.0;  ///< fixed: diffuse irradiance on a horizontal surface [W/m2]
    amrex::Real latitude_deg      = 40.0;   ///< solar: site latitude, north positive
    amrex::Real longitude_deg     = -105.0; ///< solar: site longitude, east positive
    amrex::Real day_of_year       = 172.0;  ///< solar: day of year
    amrex::Real time_zero_utc_s   = 12.0 * 3600.0; ///< solar: seconds after 00:00 UTC at simulation time zero (the sun follows UTC and the longitude)
    amrex::Real solar_constant    = 1361.0; ///< solar: [W/m2]
    amrex::Real sw_transmission   = 0.7;    ///< solar, two_stream: bulk clear-sky transmission (Bird)
    amrex::Real sw_diffuse_coeff  = 0.5;    ///< solar, two_stream: share of the attenuated beam that arrives as diffuse light
    // The lowest sun elevation [deg] at which the shadow test is exact everywhere: the
    // ray map each rank keeps (IBFaceSet::build_ray_map()) reaches as far as a ray that
    // low can travel before it clears the tallest building, (tallest height) / tan(elevation),
    // beyond this rank's columns. With the sun lower, a building farther than that from a
    // face casts no shadow on it. The clear-sky direct beam at 2 deg is 0.05 W/m2 with the
    // default transmission (2.3 W/m2 at sw_transmission = 0.8; a fixed beam keeps its
    // value at any elevation), and a small domain is covered whole whatever the value.
    amrex::Real sun_min_elev_deg = 2.0;
    // Optical properties, uniform unless the material library below overrides them.
    amrex::Real albedo            = 0.3;    ///< Face shortwave albedo
    amrex::Real albedo_ground     = 0.2;    ///< Ground albedo for the reflected diffuse term
    amrex::Real emissivity        = 0.9;    ///< Face longwave emissivity
    amrex::Real emissivity_ground = 0.95;   ///< Ground emissivity for the ground term
 
    // ---- View fractions and longwave ----------------------------------------
    // Rays per face of the cosine-weighted hemisphere sampled once at
    // initialisation: n_az azimuths x n_el elevations.
    int         view_n_az         = 16;
    int         view_n_el         = 8;
    // Downwelling longwave from the sky: "fixed" uses lw_down as given,
    // "gray" uses sky_emissivity x sigma x T_air^4 with the air temperature
    // of the face's fluid cell; with both the ground term is
    // emissivity_ground x sigma x T_ground^4. "two_stream" (radiation =
    // two_stream only, and its default) takes the sky term from the column's
    // longwave down at the face's height and the ground term from its
    // longwave up at the same height (the ground's emission and its reflection
    // of the sky, seen through the air in between).
    std::string lw_mode           = "gray";
    amrex::Real lw_down           = 300.0;  ///< fixed: sky longwave [W/m2]
    amrex::Real sky_emissivity    = 0.83;   ///< gray: effective clear-sky emissivity
    amrex::Real T_ground          = 300.0;  ///< Ground temperature for the ground term [K]
 
    // ---- Sensible and latent heat -------------------------------------------
    // Wall function between the face skin temperature and the adjacent fluid
    // cell: u* from the tangential wind at half a cell from the wall with the
    // roughness length z0_wall, theta* from the skin-to-air difference with
    // z0h_wall; H = rho c_p u* theta* (positive out of the face). Neutral on
    // every face unless stability_correction is set, which applies the
    // surface layer's similarity functions on roofs. couple_heat = false
    // diagnoses H without adding it to the temperature equation.
    amrex::Real z0_wall           = 0.01;   ///< Momentum roughness of the faces [m]
    amrex::Real z0h_wall          = 0.001;  ///< Heat roughness of the faces [m]
    bool        stability_correction = false;
    bool        couple_heat       = true;   ///< Add the face flux to the temperature equation
 
    // ---- Ground heat flux and materials -------------------------------------
    // Every face carries a slab of n_slab_layers layers between the skin and
    // the building interior at T_interior; the conduction into it is G. The
    // slab's conductivity, heat capacity and thickness, and the face's albedo
    // and emissivity, come from the material library (SLUCM schema) when
    // material_file is given: material_default for every building, or one
    // id per building in material_by_building (indexed by level 0's building
    // id from 1; every level maps its buildings to level 0's). Without a file
    // the uniform values below apply to every face.
    std::string material_file;
    int         material_default   = 1;
    std::vector<int> material_by_building;
    amrex::Real k_therm            = 1.0;    ///< Uniform conductivity [W/m/K]
    amrex::Real rho_cp             = 1.6e6;  ///< Uniform volumetric heat capacity [J/m3/K]
    amrex::Real thickness          = 0.3;    ///< Uniform slab thickness [m]
 
    // ---- Prognostic balance -------------------------------------------------
    // With prognostic on, every step solves the massless-skin balance
    // SW_abs + eps Q_ext + LW_net - H - LE - G = 0 for the skin temperature
    // of every face by Newton's method (ERF_IBSEBBalance.H), consistent with
    // the implicit slab step, and stores the fluxes at the solution. Off, the
    // skin stays at T_skin_init (or the restart value) and the terms are
    // diagnosed around it, which is what the term-by-term regression tests
    // use to check each term on its own.
    bool        prognostic         = true;
    amrex::Real T_skin_min         = 260.0;  ///< Lower bound of the skin temperature [K]
    amrex::Real T_skin_max         = 380.0;  ///< Upper bound [K]; raise it for faces under a fire
    int         newton_max_iter    = 20;     ///< Newton iteration cap per face and step
    amrex::Real newton_tol_K       = 1.0e-3; ///< Newton stops when a step changes T_skin by less than this [K]
    amrex::Real newton_max_step_K  = 20.0;   ///< Largest Newton step [K]
    amrex::Real Q_ext_uniform      = 0.0;    ///< Uniform external incident flux on every face [W/m2] (test hook)
 
    // ---- Wall function beyond neutral ----------------------------------------
    // stability_correction (above) applies the surface layer's similarity
    // functions on roofs with the face's own Obukhov length, seeded from the
    // ground surface layer's 2D Obukhov field at the face's column when
    // obukhov_seed = "ground" (or from neutral). convective_velocity =
    // "deardorff" adds a convective velocity scale to the wind the wall
    // function sees, U_eff^2 = U_tan^2 + (beta w*)^2 with
    // w* = (g/theta H/(rho c_p) depth)^(1/3) from the previous step's H when it
    // is out of the face, so u* does not collapse when the wind dies. The
    // depth is the mixed layer above a roof (z_i - z_face, floored at the
    // building height; z_i from the bulk Richardson diagnostic on the
    // horizontal-mean profile, from ERF's surface-layer pblh, or fixed) and
    // the building height for a wall. All off by default.
    std::string obukhov_seed        = "ground";   ///< "ground" or "neutral"
    std::string convective_velocity = "none";     ///< "none" or "deardorff"
    amrex::Real beta_conv           = 1.2;        ///< Gustiness factor on w*
    std::string z_i_mode            = "bulk_ri";  ///< "bulk_ri", "pblh" or "fixed"
    amrex::Real z_i                 = 1000.0;     ///< Fixed mixed-layer depth [m] (z_i_mode = fixed, and the fallback)
    amrex::Real ri_crit             = 0.25;       ///< Critical bulk Richardson number of the diagnostic
    amrex::Real obukhov_relax       = 0.5;        ///< Under-relaxation of the roof's u*, theta* iteration (1 = none), as erf-model #3486
 
    /**
     * Read the ``erf.ibseb.*`` inputs and validate them.
     *
     * Only queried, never required: a deck without any ``erf.ibseb`` key gets
     * the defaults above with the balance off. Validation is skipped when the
     * balance is off so an unrelated deck cannot abort on an unused key.
     */
    void init_params ()
    {
        amrex::ParmParse pp("erf.ibseb");
        pp.query("enable",        enable);
        pp.query("n_slab_layers", n_slab_layers);
        pp.query("T_skin_init",   T_skin_init);
        pp.query("T_interior",    T_interior);
        pp.query("csv_file",      csv_file);
        pp.query("csv_int",       csv_int);
        pp.query("debug",         debug);
        pp.query("dump_faces_file", dump_faces_file);
        pp.query("dump_faces_tag_step", dump_faces_tag_step);
        pp.query("radiation",        radiation);
        pp.query("sun_mode",         sun_mode);
        pp.query("sun_zenith_deg",   sun_zenith_deg);
        pp.query("sun_azimuth_deg",  sun_azimuth_deg);
        pp.query("sw_direct_normal", sw_direct_normal);
        pp.query("sw_diffuse",       sw_diffuse);
        pp.query("latitude_deg",     latitude_deg);
        pp.query("longitude_deg",    longitude_deg);
        pp.query("day_of_year",      day_of_year);
        pp.query("time_zero_utc_s",  time_zero_utc_s);
        pp.query("solar_constant",   solar_constant);
        pp.query("sw_transmission",  sw_transmission);
        pp.query("sw_diffuse_coeff", sw_diffuse_coeff);
        pp.query("sun_min_elev_deg", sun_min_elev_deg);
        pp.query("albedo",           albedo);
        pp.query("albedo_ground",    albedo_ground);
        pp.query("emissivity",       emissivity);
        pp.query("emissivity_ground", emissivity_ground);
        pp.query("view_n_az",        view_n_az);
        pp.query("view_n_el",        view_n_el);
        pp.query("lw_mode",          lw_mode);
        pp.query("lw_down",          lw_down);
        pp.query("sky_emissivity",   sky_emissivity);
        pp.query("T_ground",         T_ground);
        pp.query("z0_wall",          z0_wall);
        pp.query("z0h_wall",         z0h_wall);
        pp.query("stability_correction", stability_correction);
        pp.query("couple_heat",      couple_heat);
        pp.query("material_file",    material_file);
        pp.query("material_default", material_default);
        pp.queryarr("material_by_building", material_by_building);
        pp.query("k_therm",          k_therm);
        pp.query("rho_cp",           rho_cp);
        pp.query("thickness",        thickness);
        pp.query("prognostic",       prognostic);
        pp.query("T_skin_min",       T_skin_min);
        pp.query("T_skin_max",       T_skin_max);
        pp.query("newton_max_iter",  newton_max_iter);
        pp.query("newton_tol_K",     newton_tol_K);
        pp.query("newton_max_step_K", newton_max_step_K);
        pp.query("Q_ext_uniform",    Q_ext_uniform);
        pp.query("obukhov_seed",        obukhov_seed);
        pp.query("convective_velocity", convective_velocity);
        pp.query("beta_conv",           beta_conv);
        pp.query("z_i_mode",            z_i_mode);
        pp.query("z_i",                 z_i);
        pp.query("ri_crit",             ri_crit);
        pp.query("obukhov_relax",       obukhov_relax);
        if (enable) {
            if (radiation != "prescribed" && radiation != "two_stream") {
                amrex::Abort("erf.ibseb.radiation must be \"prescribed\" or \"two_stream\"");
            }
            // The column's longwave goes with its shortwave unless the deck says otherwise.
            if (radiation == "two_stream" && !pp.contains("lw_mode")) { lw_mode = "two_stream"; }
            // The columns bring the irradiances and the ground's reflection; the
            // clear-sky formulas and the uniform ground albedo are not used.
            if (radiation == "two_stream") {
                for (const char* k : {"sw_direct_normal", "sw_diffuse", "sw_transmission", "sw_diffuse_coeff",
                                      "solar_constant", "albedo_ground"}) {
                    if (pp.contains(k)) {
                        amrex::Abort(std::string("erf.ibseb.") + k + " is not used with erf.ibseb.radiation = two_stream:"
                                     " the faces take the two-stream columns' irradiances and the ground's reflection;"
                                     " remove it");
                    }
                }
            }
            if (sun_mode != "fixed" && sun_mode != "solar" && sun_mode != "two_stream") {
                amrex::Abort("erf.ibseb.sun_mode must be \"fixed\", \"solar\" or \"two_stream\"");
            }
            // The two-stream sun brings its own position, site, date and irradiance.
            if (sun_mode == "two_stream") {
                for (const char* k : {"latitude_deg", "longitude_deg", "day_of_year", "time_zero_utc_s", "solar_constant",
                                      "sun_zenith_deg", "sun_azimuth_deg", "sw_direct_normal", "sw_diffuse"}) {
                    if (pp.contains(k)) {
                        amrex::Abort(std::string("erf.ibseb.") + k + " is not used with erf.ibseb.sun_mode = two_stream:"
                                     " the faces take the two-stream sun (start_datetime, erf.rad_cons_lat / lon,"
                                     " erf.fixed_total_solar_irradiance); remove it");
                    }
                }
            }
            if (sun_zenith_deg < 0.0 || sun_zenith_deg > 180.0) { amrex::Abort("erf.ibseb.sun_zenith_deg must be in [0, 180]"); }
            if (latitude_deg < -90.0 || latitude_deg > 90.0) { amrex::Abort("erf.ibseb.latitude_deg must be in [-90, 90]"); }
            if (longitude_deg < -180.0 || longitude_deg > 180.0) { amrex::Abort("erf.ibseb.longitude_deg must be in [-180, 180]"); }
            if (day_of_year < 1.0 || day_of_year > 366.0) { amrex::Abort("erf.ibseb.day_of_year must be in [1, 366]"); }
            // The sun is computed from UTC (time_zero_utc_s) and the longitude; a
            // time-zone offset has no place in that and used to cancel silently.
            if (pp.contains("utc_offset_hours")) {
                amrex::Abort("erf.ibseb.utc_offset_hours is not an input: the sun follows time_zero_utc_s (UTC) and longitude_deg; remove it");
            }
            if (sw_direct_normal < 0.0 || sw_diffuse < 0.0) { amrex::Abort("erf.ibseb.sw_direct_normal and sw_diffuse must be >= 0"); }
            if (sw_transmission <= 0.0 || sw_transmission > 1.0) { amrex::Abort("erf.ibseb.sw_transmission must be in (0, 1]"); }
            if (sw_diffuse_coeff < 0.0 || sw_diffuse_coeff > 1.0) { amrex::Abort("erf.ibseb.sw_diffuse_coeff must be in [0, 1]"); }
            if (solar_constant <= 0.0) { amrex::Abort("erf.ibseb.solar_constant must be > 0"); }
            if (sun_min_elev_deg <= 0.0 || sun_min_elev_deg >= 90.0) {
                amrex::Abort("erf.ibseb.sun_min_elev_deg must be in (0, 90) [deg]");
            }
            if (albedo < 0.0 || albedo > 1.0 || albedo_ground < 0.0 || albedo_ground > 1.0) {
                amrex::Abort("erf.ibseb.albedo and albedo_ground must be in [0, 1]");
            }
            if (emissivity <= 0.0 || emissivity > 1.0 || emissivity_ground <= 0.0 || emissivity_ground > 1.0 ||
                sky_emissivity <= 0.0 || sky_emissivity > 1.0) {
                amrex::Abort("erf.ibseb emissivities must be in (0, 1]");
            }
            if (view_n_az < 1 || view_n_el < 1) { amrex::Abort("erf.ibseb.view_n_az and view_n_el must be >= 1"); }
            if (lw_mode != "fixed" && lw_mode != "gray" && lw_mode != "two_stream") {
                amrex::Abort("erf.ibseb.lw_mode must be \"fixed\", \"gray\" or \"two_stream\"");
            }
            if (lw_mode == "two_stream" && radiation != "two_stream") {
                amrex::Abort("erf.ibseb.lw_mode = two_stream needs erf.ibseb.radiation = two_stream");
            }
            if (lw_mode == "two_stream") {
                for (const char* k : {"lw_down", "sky_emissivity", "T_ground", "emissivity_ground"}) {
                    if (pp.contains(k)) {
                        amrex::Abort(std::string("erf.ibseb.") + k + " is not used with erf.ibseb.lw_mode = two_stream:"
                                     " the faces take the two-stream columns' longwave down and the ground's longwave up;"
                                     " remove it");
                    }
                }
            }
            if (lw_down < 0.0 || T_ground <= 0.0) { amrex::Abort("erf.ibseb.lw_down must be >= 0 and T_ground > 0"); }
            if (z0_wall <= 0.0 || z0h_wall <= 0.0) { amrex::Abort("erf.ibseb.z0_wall and z0h_wall must be positive [m]"); }
            if (k_therm <= 0.0 || rho_cp <= 0.0 || thickness <= 0.0) { amrex::Abort("erf.ibseb.k_therm, rho_cp and thickness must be positive"); }
            if (n_slab_layers > 32) { amrex::Abort("erf.ibseb.n_slab_layers must be <= 32 (the slab kernels keep their layers on the stack)"); }
            if (n_slab_layers < 1) { amrex::Abort("erf.ibseb.n_slab_layers must be >= 1"); }
            if (T_skin_init <= 0.0 || T_interior <= 0.0) {
                amrex::Abort("erf.ibseb.T_skin_init and T_interior must be positive [K]");
            }
            if (T_skin_min <= 0.0 || T_skin_max <= T_skin_min) {
                amrex::Abort("erf.ibseb.T_skin_min must be positive and below T_skin_max");
            }
            if (prognostic && (T_skin_init < T_skin_min || T_skin_init > T_skin_max)) {
                amrex::Abort("erf.ibseb.T_skin_init must lie within [T_skin_min, T_skin_max]");
            }
            if (newton_max_iter < 1 || newton_tol_K <= 0.0 || newton_max_step_K <= 0.0) {
                amrex::Abort("erf.ibseb.newton_max_iter must be >= 1, newton_tol_K and newton_max_step_K positive");
            }
            if (Q_ext_uniform < 0.0) { amrex::Abort("erf.ibseb.Q_ext_uniform must be >= 0"); }
            if (obukhov_seed != "ground" && obukhov_seed != "neutral") { amrex::Abort("erf.ibseb.obukhov_seed must be \"ground\" or \"neutral\""); }
            if (convective_velocity != "none" && convective_velocity != "deardorff") {
                amrex::Abort("erf.ibseb.convective_velocity must be \"none\" or \"deardorff\"");
            }
            if (z_i_mode != "bulk_ri" && z_i_mode != "pblh" && z_i_mode != "fixed") {
                amrex::Abort("erf.ibseb.z_i_mode must be \"bulk_ri\", \"pblh\" or \"fixed\"");
            }
            if (beta_conv <= 0.0 || z_i <= 0.0 || ri_crit <= 0.0) { amrex::Abort("erf.ibseb.beta_conv, z_i and ri_crit must be positive"); }
            if (obukhov_relax <= 0.0 || obukhov_relax > 1.0) { amrex::Abort("erf.ibseb.obukhov_relax must be in (0, 1]"); }
        }
    }
};
 
#endif