Skip to content

Fire: opt-in WRF-style Rothermel wind coupling for the directional projection formula #410

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
RaymndH wants to merge 1 commit into
base: ERF-Hazard
Choose a base branch
from
Open

Fire: opt-in WRF-style Rothermel wind coupling for the directional projection formula #410

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 75 additions & 16 deletions Source/Fire/ERF_DirectionalRos.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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);
}
}
Expand Down Expand Up @@ -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,
Expand All @@ -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);
Expand All @@ -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) {
Expand All @@ -364,18 +413,25 @@ 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;
const amrex::Real ny = dpdy / grad;
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);
});
}
}
Expand All @@ -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,
Expand All @@ -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);
}

Expand Down
30 changes: 18 additions & 12 deletions Source/Fire/ERF_FireLayer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down
31 changes: 31 additions & 0 deletions Source/Fire/ERF_FireParams.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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);
Expand Down
16 changes: 15 additions & 1 deletion Source/Fire/ERF_Rothermel.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);

Expand Down