From 110262e225c28ce21140f8374ca8afb223ecbaa7 Mon Sep 17 00:00:00 2001 From: Raymond Harding Date: Fri, 11 Sep 2026 14:13:30 -0700 Subject: [PATCH] Add gated WRF-style Rothermel wind coupling; prescribed wind and wind-limit bypass for testing Wind coupling (opt-in, backward compatible): The default directional projection formula evaluates Rothermel's model with the wind projected onto the front normal before raising it to the B power, R(theta) ~ (U cos theta)^B. This is non-convex once phi_w(B-1) > 1, so a finite ignition line's endpoints grow a Wulff wedge whose tip falls well below the head rate over time -- verified against a Munoz-Esparza-style finite-line test case, where the head rate degraded from ~100% to ~40% of the theoretical Rf over 2100s. erf.fire.directional_wind_coupling = "wrf" (default: "projection", unchanged) instead exponentiates Rothermel's wind/slope factor from the raw, unprojected wind speed and scales the result by the cosine to the spread direction afterward, matching WRF-Fire's fire_ros (module_fr_fire_phys.F). This keeps R(theta) linear in cos(theta) -- the convex support function of a circle -- so the wedge does not form; the same test case now tracks the theoretical head rate to within 0.5% for the full run. Default behavior is verified bit-for-bit identical to the prior code (max|delta phi| = 0 over the whole domain on the same test case). Test infrastructure: - erf.fire.prescribed_wind[_x/_y]: bypass atmospheric wind interpolation with an exact constant wind, for controlled idealized-wind test cases. - erf.fire.use_wind_limit: allow disabling the MEWS wind-speed cap (default on, unchanged), needed when a test's wind speed exceeds the cap and the cap would otherwise mask the ROS model's true behavior. Co-Authored-By: Claude Sonnet 5 --- Source/Fire/ERF_DirectionalRos.H | 91 ++++++++++++++++++++++++++------ Source/Fire/ERF_FireLayer.cpp | 30 ++++++----- Source/Fire/ERF_FireParams.H | 31 +++++++++++ Source/Fire/ERF_Rothermel.H | 16 +++++- 4 files changed, 139 insertions(+), 29 deletions(-) diff --git a/Source/Fire/ERF_DirectionalRos.H b/Source/Fire/ERF_DirectionalRos.H index e0975c3d25..85dfb779ea 100644 --- a/Source/Fire/ERF_DirectionalRos.H +++ b/Source/Fire/ERF_DirectionalRos.H @@ -135,18 +135,46 @@ inline constexpr int DIRECTIONAL_ROS_FBP = 4; /// isotropic wind and slope magnitudes are used instead. inline constexpr amrex::Real DIRECTIONAL_ROS_GRAD_MIN = 1.0e-12; +/// Angular coupling of the wind/slope factor into the projected rate, +/// erf.fire.directional_wind_coupling. Rothermel only; every other model +/// ignores this and directional_ros_cell's U_mag/s_mag arguments. +inline constexpr int DIRECTIONAL_WIND_COUPLING_PROJECTION = 0; ///< default: project wind, then exponentiate (Rothermel's U^B) +inline constexpr int DIRECTIONAL_WIND_COUPLING_WRF = 1; ///< opt-in: exponentiate on the raw wind, then scale by cos(theta) (WRF-Fire's fire_ros) + /** * @brief Rate of spread along one direction * - * @param[in] st Precomputed coefficients for the selected model - * @param[in] U_n Wind component along the spread direction [m/s] - * @param[in] s_n Terrain slope along the spread direction [-] + * @param[in] st Precomputed coefficients for the selected model + * @param[in] U_n Wind component along the spread direction [m/s] + * @param[in] s_n Terrain slope along the spread direction [-] + * @param[in] wind_coupling DIRECTIONAL_WIND_COUPLING_PROJECTION (default) or + * _WRF; Rothermel only, see below. + * @param[in] U_mag Raw (unprojected) wind speed [m/s]; read only + * when wind_coupling == _WRF. + * @param[in] s_mag Raw (unprojected) slope magnitude [-]; read only + * when wind_coupling == _WRF. * @return ROS [m/s] + * + * DIRECTIONAL_WIND_COUPLING_PROJECTION (default, backward compatible): the + * Rothermel branch projects the wind onto the spread direction and then + * raises it to Rothermel's B power, R(theta) ~ (U cos theta)^B -- the + * original behaviour, byte-for-byte, since U_mag/s_mag are not read on this + * path. This is non-convex once phi_w(B-1) > 1, which is what lets a finite + * line's ends grow a Wulff wedge that falls behind the head rate over time + * (see the file comment). DIRECTIONAL_WIND_COUPLING_WRF (opt-in) instead + * exponentiates the wind/slope factor from the raw magnitude U_mag/s_mag and + * only afterward scales by the cosine to the spread direction (U_n/U_mag), + * matching WRF-Fire's fire_ros (module_fr_fire_phys.F). That keeps R(theta) + * linear in cos(theta) -- the convex support function of a circle -- so the + * wedge does not form. Selected via erf.fire.directional_wind_coupling. */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real directional_ros_cell(const DirectionalRosState& st, amrex::Real U_n, - amrex::Real s_n) noexcept + amrex::Real s_n, + int wind_coupling = DIRECTIONAL_WIND_COUPLING_PROJECTION, + amrex::Real U_mag = amrex::Real(0.0), + amrex::Real s_mag = amrex::Real(0.0)) noexcept { // Only a wind pushing into the unburned fuel drives the front, and only an // upslope component enhances it; the models take magnitudes, so the @@ -164,6 +192,17 @@ amrex::Real directional_ros_cell(const DirectionalRosState& st, case DIRECTIONAL_ROS_FBP: return fbp_ros(st.fbp, U, s); default: + if (wind_coupling == DIRECTIONAL_WIND_COUPLING_WRF) { + // The raw magnitude is at least the projected component (it is + // the same vector, just not dotted with the normal); guard + // against a caller passing 0 by mistake so cos_w falls back + // to 1 rather than dividing by (near) zero. + const amrex::Real Um = amrex::max(U_mag, U); + const amrex::Real sm = amrex::max(s_mag, s); + const amrex::Real cos_w = (Um > amrex::Real(1.0e-12)) ? U / Um : amrex::Real(1.0); + const amrex::Real cos_s = (sm > amrex::Real(1.0e-12)) ? s / sm : amrex::Real(1.0); + return rothermel_ros_cell(Um, amrex::Real(0.0), sm, amrex::Real(0.0), st.rc, cos_w, cos_s); + } return rothermel_ros_cell(U, amrex::Real(0.0), s, amrex::Real(0.0), st.rc); } } @@ -308,6 +347,9 @@ amrex::Real directional_ellipse_speed(const DirectionalEllipse& d, * or DIRECTIONAL_SHAPE_ELLIPSE (the model's spread ellipse) * @param[in] ellipse_lw Flank rule of the ellipse, DIRECTIONAL_ELLIPSE_LW_* * @param[in] ellipse_lw_max Cap on Anderson's length-to-width ratio + * @param[in] wind_coupling DIRECTIONAL_WIND_COUPLING_PROJECTION (default) or + * _WRF; Rothermel-with-DIRECTIONAL_SHAPE_PROJECTION only, + * see directional_ros_cell() */ inline void fill_directional_ros(amrex::MultiFab& R_mf, const amrex::MultiFab& phi, @@ -319,7 +361,9 @@ inline void fill_directional_ros(amrex::MultiFab& R_mf, bool wall_extrap = false, int shape = DIRECTIONAL_SHAPE_PROJECTION, int ellipse_lw = DIRECTIONAL_ELLIPSE_LW_MODEL, - amrex::Real ellipse_lw_max = amrex::Real(8.0)) + amrex::Real ellipse_lw_max = amrex::Real(8.0), + const fire_levelset::LevelSetGradient& grad_scheme = fire_levelset::LevelSetGradient{}, + int wind_coupling = DIRECTIONAL_WIND_COUPLING_PROJECTION) { const amrex::Real dx = geom.CellSize(0); const amrex::Real dy = geom.CellSize(1); @@ -345,12 +389,17 @@ inline void fill_directional_ros(amrex::MultiFab& R_mf, const amrex::Real dzdx = slopes_arr(i, j, k, 0); const amrex::Real dzdy = slopes_arr(i, j, k, 1); - // Front normal from central differences; it points toward increasing - // phi, that is into the unburned fuel. - const amrex::Real dpdx = (wall.at(phi_arr, i, j, k, 1, 0) - - wall.at(phi_arr, i, j, k, -1, 0)) / (2.0 * dx); - const amrex::Real dpdy = (wall.at(phi_arr, i, j, k, 0, 1) - - wall.at(phi_arr, i, j, k, 0, -1)) / (2.0 * dy); + // Front normal from the same one-sided derivatives (first order or + // HJ-WENO5-Z, per grad_scheme) as the magnitude in + // compute_levelset_rhs, averaged into a centred estimate; this + // reduces exactly to the prior central difference when + // grad_scheme.scheme == LEVELSET_GRAD_UPWIND1, and points toward + // increasing phi, that is into the unburned fuel. + amrex::Real dxm, dxp, dym, dyp; + fire_levelset::one_sided_derivatives(phi_arr, dx, i, j, k, 1, 0, grad_scheme, wall, dxm, dxp); + fire_levelset::one_sided_derivatives(phi_arr, dy, i, j, k, 0, 1, grad_scheme, wall, dym, dyp); + const amrex::Real dpdx = amrex::Real(0.5) * (dxm + dxp); + const amrex::Real dpdy = amrex::Real(0.5) * (dym + dyp); const amrex::Real grad = std::sqrt(dpdx*dpdx + dpdy*dpdy); if (shape == DIRECTIONAL_SHAPE_ELLIPSE) { @@ -364,6 +413,13 @@ inline void fill_directional_ros(amrex::MultiFab& R_mf, return; } + // Raw (unprojected) magnitudes: passed through so the Rothermel + // branch of directional_ros_cell can exponentiate from these, + // WRF-style, then scale by the cosine to the front normal below, + // rather than projecting first. Ignored by every other model. + const amrex::Real U_mag = std::sqrt(u*u + v*v); + const amrex::Real s_mag = std::sqrt(dzdx*dzdx + dzdy*dzdy); + amrex::Real U_n, s_n; if (grad > grad_min) { const amrex::Real nx = dpdx / grad; @@ -371,11 +427,11 @@ inline void fill_directional_ros(amrex::MultiFab& R_mf, U_n = u*nx + v*ny; s_n = dzdx*nx + dzdy*ny; } else { - U_n = std::sqrt(u*u + v*v); - s_n = std::sqrt(dzdx*dzdx + dzdy*dzdy); + U_n = U_mag; + s_n = s_mag; } - R_arr(i, j, k) = directional_ros_cell(st, U_n, s_n); + R_arr(i, j, k) = directional_ros_cell(st, U_n, s_n, wind_coupling, U_mag, s_mag); }); } } @@ -398,6 +454,8 @@ inline void fill_directional_ros(amrex::MultiFab& R_mf, * @param[in] ellipse_lw Flank rule of the ellipse, DIRECTIONAL_ELLIPSE_LW_* * @param[in] ellipse_lw_max Cap on Anderson's length-to-width ratio * @param[in] ros_scale Optional per-cell factor on the rebuilt ROS (front acceleration clock) + * @param[in] wind_coupling DIRECTIONAL_WIND_COUPLING_PROJECTION (default) or _WRF; + * Rothermel-with-DIRECTIONAL_SHAPE_PROJECTION only, see directional_ros_cell() */ inline void advect_levelset_directional_rk3(amrex::MultiFab& phi, const amrex::MultiFab& wind, @@ -412,12 +470,13 @@ inline void advect_levelset_directional_rk3(amrex::MultiFab& phi, int shape = DIRECTIONAL_SHAPE_PROJECTION, int ellipse_lw = DIRECTIONAL_ELLIPSE_LW_MODEL, amrex::Real ellipse_lw_max = amrex::Real(8.0), - const amrex::MultiFab* ros_scale = nullptr) + const amrex::MultiFab* ros_scale = nullptr, + int wind_coupling = DIRECTIONAL_WIND_COUPLING_PROJECTION) { fire_levelset::advect_levelset_rk3_with_fill(phi, slopes, geom, dt, eps_visc, [&](amrex::MultiFab& R, const amrex::MultiFab& phi_stage) { fill_directional_ros(R, phi_stage, wind, slopes, geom, st, nonburnable, wall_extrap, shape, - ellipse_lw, ellipse_lw_max); + ellipse_lw, ellipse_lw_max, grad_scheme, wind_coupling); }, nonburnable, wall_extrap, grad_scheme, ros_scale); } diff --git a/Source/Fire/ERF_FireLayer.cpp b/Source/Fire/ERF_FireLayer.cpp index e81392ad08..4d08e1dd0f 100644 --- a/Source/Fire/ERF_FireLayer.cpp +++ b/Source/Fire/ERF_FireLayer.cpp @@ -622,17 +622,22 @@ void FireLayer::advance(Real time, Real dt, SurfaceLayer& surface_layer, if (m_params.fire_debug) amrex::Print() << "[FIRE DEBUG] Starting fire advance step with dt=" << dt << std::endl; - const bool wind_open = m_params.structures.wind_open_columns && m_open_frac_atm && m_roof_h_atm; - fill_fire_wind_from_interpolation(*fire_wind_ref, *fire_wind_extract_z, xvel, yvel, z_phys_cc, - *fire_surface_z, *fire_col_ground, - m_fg, m_params.wind_ref_ht, m_nz, - m_use_per_fuel_wind_ht ? fire_fuel_model.get() : nullptr, - m_use_per_fuel_wind_ht ? m_d_fcwh.data() : nullptr, - m_use_per_fuel_wind_ht ? FUEL_SLOT_COUNT - 1 : 0, - m_params.wind_interp, - wind_open ? m_open_frac_atm.get() : nullptr, - wind_open ? m_roof_h_atm.get() : nullptr, - m_params.wind_sample_ht, m_params.wind_sample_z0); + if (m_params.prescribed_wind) { + fire_wind_ref->setVal(m_params.prescribed_wind_x, 0, 1, 0); + fire_wind_ref->setVal(m_params.prescribed_wind_y, 1, 1, 0); + } else { + const bool wind_open = m_params.structures.wind_open_columns && m_open_frac_atm && m_roof_h_atm; + fill_fire_wind_from_interpolation(*fire_wind_ref, *fire_wind_extract_z, xvel, yvel, z_phys_cc, + *fire_surface_z, *fire_col_ground, + m_fg, m_params.wind_ref_ht, m_nz, + m_use_per_fuel_wind_ht ? fire_fuel_model.get() : nullptr, + m_use_per_fuel_wind_ht ? m_d_fcwh.data() : nullptr, + m_use_per_fuel_wind_ht ? FUEL_SLOT_COUNT - 1 : 0, + m_params.wind_interp, + wind_open ? m_open_frac_atm.get() : nullptr, + wind_open ? m_roof_h_atm.get() : nullptr, + m_params.wind_sample_ht, m_params.wind_sample_z0); + } if (m_params.fire_debug) { if (m_params.wind_sample_ht > 0.0) { amrex::Print() << "[FIRE DEBUG] Wind sampled at " << m_params.wind_sample_ht @@ -1002,7 +1007,8 @@ void FireLayer::advance(Real time, Real dt, SurfaceLayer& surface_layer, ls_grad, m_params.directional_shape, m_params.directional_ellipse_lw, m_params.directional_ellipse_lw_max, - accel_factor.get()); + accel_factor.get(), + m_params.directional_wind_coupling); } else if (m_params.levelset_ellipse) { // Huygens ellipse: the model's rate is the head rate and the // normal speed follows the ellipse set by the midflame wind. diff --git a/Source/Fire/ERF_FireParams.H b/Source/Fire/ERF_FireParams.H index dd6419e04d..64e974cc8e 100644 --- a/Source/Fire/ERF_FireParams.H +++ b/Source/Fire/ERF_FireParams.H @@ -32,6 +32,9 @@ struct FireParams amrex::Real ignition_y = 0.0; ///< Ignition center y [m] amrex::Real ignition_r = 20.0; ///< Ignition radius [m] amrex::Real wind_ref_ht = 6.1; ///< MOST wind reference height [m] + bool prescribed_wind = false; ///< Bypass atmospheric wind interpolation + amrex::Real prescribed_wind_x = 0.0; ///< Constant reference wind, x [m/s] + amrex::Real prescribed_wind_y = 0.0; ///< Constant reference wind, y [m/s] // Phase 13A: Per-fuel wind height bool use_per_fuel_wind_ht = false; ///< Use per-fuel fcwh table for wind height amrex::Real waf_fcz0_scale = 1.0; ///< Scale factor on fcz0 roughness for WAF [dimensionless] @@ -385,6 +388,21 @@ struct FireParams // macarthur, cheney_gould and fbp. int directional_shape = 0; + // Angular coupling of the wind/slope factor into that projection, + // erf.fire.directional_wind_coupling. "projection" (0, default, backward + // compatible) projects the wind onto the spread direction before + // Rothermel's power-law wind factor, R(theta) ~ (U cos theta)^B -- this is + // the original behaviour, unchanged. "wrf" (1, opt-in) instead computes + // the wind/slope factor from the raw, unprojected speed and scales the + // result by the cosine afterward, matching WRF-Fire's fire_ros + // (module_fr_fire_phys.F): R(theta) is then linear in cos(theta), the + // convex support function of a circle, rather than cos^B(theta), which + // stops being convex once phi_w(B-1) > 1 and lets a finite line's ends + // grow a Wulff wedge that falls behind the head rate over time (see + // ERF_DirectionalRos.H). Rothermel with directional_shape = "projection" + // only; ignored otherwise (the ellipse path is already convex). + int directional_wind_coupling = 0; + // Flank rule of that ellipse, erf.fire.directional_ellipse_lw. "model" (0, // default) keeps the model's no-wind, no-slope rate at the flanks, so the // length-to-width ratio (R_h/R0 + 1)/2 is not capped. "anderson" (1) sets the @@ -664,6 +682,9 @@ struct FireParams pp.query("ignition_y", ignition_y); pp.query("ignition_r", ignition_r); pp.query("wind_ref_ht", wind_ref_ht); + pp.query("prescribed_wind", prescribed_wind); + pp.query("prescribed_wind_x", prescribed_wind_x); + pp.query("prescribed_wind_y", prescribed_wind_y); pp.query("wind_sample_ht", wind_sample_ht); pp.query("wind_sample_z0", wind_sample_z0); if (wind_sample_ht > 0.0) { @@ -929,6 +950,16 @@ struct FireParams amrex::Abort("erf.fire.directional_shape must be \"projection\" or \"ellipse\""); } + std::string coupling_str = (directional_wind_coupling == 1) ? "wrf" : "projection"; + pp.query("directional_wind_coupling", coupling_str); + if (coupling_str == "projection") { + directional_wind_coupling = 0; + } else if (coupling_str == "wrf") { + directional_wind_coupling = 1; + } else { + amrex::Abort("erf.fire.directional_wind_coupling must be \"projection\" or \"wrf\""); + } + std::string lw_str = (directional_ellipse_lw == 1) ? "anderson" : "model"; pp.query("directional_ellipse_lw", lw_str); pp.query("directional_ellipse_lw_max", directional_ellipse_lw_max); diff --git a/Source/Fire/ERF_Rothermel.H b/Source/Fire/ERF_Rothermel.H index 7e9fbe4cbc..7315767344 100644 --- a/Source/Fire/ERF_Rothermel.H +++ b/Source/Fire/ERF_Rothermel.H @@ -296,7 +296,9 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real rothermel_ros_cell( amrex::Real ux_eff, amrex::Real uy_eff, amrex::Real sx, amrex::Real sy, - const RothermelComputed& rc) noexcept + const RothermelComputed& rc, + amrex::Real cos_theta_wind = amrex::Real(1.0), + amrex::Real cos_theta_slope = amrex::Real(1.0)) noexcept { // Effective wind speed [m/s] amrex::Real U_eff = std::sqrt(ux_eff*ux_eff + uy_eff*uy_eff); @@ -321,6 +323,18 @@ amrex::Real rothermel_ros_cell( // Wind factor amrex::Real phi_w = rc.C * std::pow(U_eff_ftmin, rc.B) * rc.beta_ratio_E; + // WRF-Fire (fire_ros, module_fr_fire_phys.F) couples direction in by + // scaling the already-exponentiated wind/slope excess by the cosine to + // the spread direction, rather than projecting (ux_eff,uy_eff) before + // raising it to the B power. cos_theta_* = 1 (the default) reproduces the + // omnidirectional rate exactly, unchanged from before this parameter + // existed; a caller doing directional spread may pass the raw wind speed + // as (ux_eff,uy_eff) and the true cosine here instead, which keeps R + // linear in cos(theta) and so a convex (circle) support function instead + // of cos^B(theta) -- see ERF_DirectionalRos.H's file comment. + phi_w *= amrex::max(cos_theta_wind, amrex::Real(0.0)); + phi_s *= amrex::max(cos_theta_slope, amrex::Real(0.0)); + // Rate of spread: R = R0 * (1 + φ_w + φ_s) amrex::Real ROS = rc.R0 * (1.0 + phi_w + phi_s);