diff --git a/doc/Configuration.rst b/doc/Configuration.rst index 519d515..0dc1186 100644 --- a/doc/Configuration.rst +++ b/doc/Configuration.rst @@ -292,6 +292,36 @@ output fields retain their existing single-precision NetCDF representation. Rate of Spread (ROS) parameterization option. 0: Rothermel model (only option currently implemented) + For burnable cells, effective ground moisture ``fmc_g`` must be finite and + nonnegative, and fuel-table extinction moisture ``fuelmce`` must be finite + and positive. Both use kg water per kg dry fuel. Invalid inputs stop with + the cell indices or fuel category in the diagnostic. Fuel-table values are + checked during initialization, and prescribed ``fuelmc_g`` is checked in + the namelist (zero moisture is allowed). The moisture component validates + fixed averaging weights during initialization and checks computed classes + before averaging and updating spread parameters. The restart reader checks + restored class and effective moisture separately: ``fmc_g`` is read directly, + not recomputed from ``fmc_gc``, and prescribed-moisture restarts have no class + state. Spatial moisture checks run outside OpenMP tile regions and exclude + nonburnable cells. + Existing restart checks also require available class state and forcing history. + + Aggregate moisture damping retains the existing polynomial below extinction + and is zero at or above extinction. Reaction intensity, ``iboros`` (reaction + intensity times residence time), and the Rothermel base spread rate are + therefore zero at extinction. The stored physical moisture is not capped. + The independent empirical chaparral spread law for Anderson category 4 is + unchanged; its reaction-based flame length is zero at or above extinction. + + Restart files retain physical moisture. The reader reconstructs spread + parameters from that state using the current damping calculation, so no + new restart field or namelist option is required. Existing PR #39 files + remain readable when their inputs satisfy the same validation. A legacy + checkpoint can contain flame diagnostics calculated with the old damping; + these are replaced at the next fire advance. Readability does not imply + numerical equivalence across the extinction correction. Continuity must be + assessed with continuous and restarted runs of the same model revision. + ``fmc_opt``: *integer* (Default: ``-1``) :term:`FMC` model -1 = Constant fuel moisture (only option currently implemented) diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index 1befbc0..e50265f 100644 --- a/io/namelist_mod.F90 +++ b/io/namelist_mod.F90 @@ -396,6 +396,9 @@ subroutine Check_nml (this, require_atm_interval) if (this%ideal_opt /= 0 .and. this%fmoist_run) & call Stop_simulation ('ideal runs do not support a FMC model') + ! Prescribed moisture may describe dry fuel, but must be a physical mass ratio. + if (.not. ieee_is_finite (this%fuelmc_g)) call Stop_simulation ('fuelmc_g must be finite') + if (this%fuelmc_g < 0.0) call Stop_simulation ('fuelmc_g must be nonnegative') if (this%fire_is_real_perim .and. this%fire_num_ignitions /= 1) & call Stop_simulation ('perimeter mode requires one ignition record for its ignition time') if (.not. this%fire_is_real_perim .and. this%fire_num_ignitions <= 0) & diff --git a/physics/fmc_wrffire_mod.F90 b/physics/fmc_wrffire_mod.F90 index b5bd5c9..163421f 100644 --- a/physics/fmc_wrffire_mod.F90 +++ b/physics/fmc_wrffire_mod.F90 @@ -55,6 +55,7 @@ module fmc_wrffire_mod procedure, public :: Advance_fmc_model => Advance_fmc_model procedure, public :: Advance_moisture_classes => Advance_moisture_classes procedure, public :: Average_moisture_classes => Average_moisture_classes + procedure, public :: Check_moisture_classes => Check_moisture_classes procedure, public :: Add_restart_vars => Add_restart_vars procedure, public :: Init => Init_fmc_wrffire end type fmc_wrffire_t @@ -170,11 +171,20 @@ subroutine Init_fmc_wrffire (this, fuels, fuelmc_g, fuelmc_g_live, ifms, ifme, j this%fmc_gc_initial_value(4) = fuelmc_g this%fmc_gc_initial_value(5) = fuelmc_g_live - ! Calc averaging weights + ! Fixed weights must form a nonnegative, normalized mean for each loaded + ! fuel category. Validate once here; subsequent updates check class state. do i = 1, mfuelcats fgi_t(i) = 0.0 do k = 1, MOISTURE_CLASSES + if (.not. ieee_is_finite (fgi_c(i, k))) then + write (msg, '(a,2(1x,i0))') 'Moisture averaging fuel load must be finite for category,class:', i, k + call Stop_simulation (msg) + end if if (fgi_c(i, k) >= 0.0) then + if (fgi_c(i, k) > huge (fgi_t(i)) - fgi_t(i)) then + write (msg, '(a,1x,i0)') 'Moisture averaging total fuel load exceeds representable range for category:', i + call Stop_simulation (msg) + end if fgi_t(i) = fgi_t(i) + fgi_c(i, k) else write (msg, *) 'fuel load in category', i, ' fuel class ', k, ' is ', fgi_c(i, k), ',must be nonegative.' @@ -189,6 +199,17 @@ subroutine Init_fmc_wrffire (this, fuels, fuelmc_g, fuelmc_g_live, ifms, ifme, j this%fmc_gw(i, k) = 0.0 end if end do + if (.not. all (ieee_is_finite (this%fmc_gw(i,:))) .or. any (this%fmc_gw(i,:) < 0.0)) then + write (msg, '(a,1x,i0)') 'Moisture averaging weights must be finite and nonnegative for category:', i + call Stop_simulation (msg) + end if + if (fgi_t(i) > 0.0) then + ! Allow rounding from summing MOISTURE_CLASSES default-real fractions. + if (abs (sum (this%fmc_gw(i,:)) - 1.0) > MOISTURE_CLASSES * epsilon (1.0)) then + write (msg, '(a,1x,i0)') 'Moisture averaging weights must sum to one for category:', i + call Stop_simulation (msg) + end if + end if end do ! moisture model derived scalars @@ -260,6 +281,13 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms, !$OMP END PARALLEL DO this%initialized = .true. + ! Validate the completed update on the initial thread, before averaging. + ! Fatal input diagnostics must not race across OpenMP tiles in MPI_Abort. + do ij = 1, num_tiles + call this%Check_moisture_classes (ifms, ifme, jfms, jfme, i_start(ij), i_end(ij), j_start(ij), j_end(ij), & + nfuel_cat, fuels%no_fuel_cat, 'Computed') + end do + !$OMP PARALLEL DO & !$OMP PRIVATE (ij) do ij = 1, num_tiles @@ -278,6 +306,35 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms, end subroutine Advance_fmc_model + subroutine Check_moisture_classes (this, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, nfuel_cat, no_fuel_cat, source) + + implicit none + + class (fmc_wrffire_t), intent (in) :: this + integer, intent (in) :: ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, no_fuel_cat + real, dimension(ifms:ifme, jfms:jfme), intent (in) :: nfuel_cat + character (len = *), intent (in) :: source + integer :: i, j, k + character (len = 160) :: msg + + + do j = jfts, jfte + do k = 1, MOISTURE_CLASSES + do i = ifts, ifte + if (int (nfuel_cat(i, j)) == no_fuel_cat) cycle + if (.not. ieee_is_finite (this%fmc_gc(i, k, j))) then + write (msg, '(a,a,3(1x,i0))') source, ' fmc_gc must be finite at i,class,j:', i, k, j + call Stop_simulation (msg) + else if (this%fmc_gc(i, k, j) < 0.0) then + write (msg, '(a,a,3(1x,i0))') source, ' fmc_gc must be nonnegative at i,class,j:', i, k, j + call Stop_simulation (msg) + end if + end do + end do + end do + + end subroutine Check_moisture_classes + subroutine Average_moisture_classes (this, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, nfuel_cat, fmc_g) implicit none diff --git a/physics/fuel_anderson_mod.F90 b/physics/fuel_anderson_mod.F90 index e9ec224..ba11b49 100644 --- a/physics/fuel_anderson_mod.F90 +++ b/physics/fuel_anderson_mod.F90 @@ -117,6 +117,8 @@ subroutine Init_anderson_fuel_model (this, fuelmc_c) end if end do + call this%Check_fuels () + ! Initialize wind adjustment factor call this%Calc_wind_adjustment_factor() diff --git a/physics/ros_wrffire_mod.F90 b/physics/ros_wrffire_mod.F90 index b8cc25f..2b0b2bc 100644 --- a/physics/ros_wrffire_mod.F90 +++ b/physics/ros_wrffire_mod.F90 @@ -176,9 +176,16 @@ subroutine Set_ros_parameters_wrffire (this, ifms, ifme, jfms, jfme, ifts, ifte, ! net fuel loading, lb/ft^2 wn = fuelload/(1 + fuels%st(k)) - rtemp1 = fmc_g(i, j) / fuels%fuelmce(k) ! moist damp coef - etam = 1.0 - 2.59 * rtemp1 + 5.11 * rtemp1 ** 2 - 3.52 * rtemp1 ** 3 + ! Zero reaction at extinction without capping physical moisture. + ! Polynomial cancellation can leave a positive residual at equality. + ! Retain the existing polynomial and arithmetic below extinction. + if (fmc_g(i, j) >= fuels%fuelmce(k)) then + etam = 0.0 + else + rtemp1 = fmc_g(i, j) / fuels%fuelmce(k) + etam = 1.0 - 2.59 * rtemp1 + 5.11 * rtemp1 ** 2 - 3.52 * rtemp1 ** 3 + end if ! mineral damping coef etas = 0.174 * fuels%se(k) ** (-0.19) !rxn intensity,btu/ft^2 min diff --git a/share/fmc_mod.F90 b/share/fmc_mod.F90 index 165ee5f..e26b966 100644 --- a/share/fmc_mod.F90 +++ b/share/fmc_mod.F90 @@ -1,7 +1,9 @@ module fmc_mod + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite use fuel_mod, only : fuel_t use ros_mod, only : ros_t + use stderrout_mod, only : Stop_simulation implicit none @@ -16,6 +18,7 @@ module fmc_mod procedure (Advance_fmc_model), deferred :: Advance_fmc_model procedure (Add_restart_vars), deferred :: Add_restart_vars procedure (Init), deferred :: Init + procedure, public :: Check_fmc => Check_fmc end type fmc_t abstract interface @@ -57,4 +60,37 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms, end subroutine Advance_fmc_model end interface + contains + + subroutine Check_fmc (this, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, nfuel_cat, no_fuel_cat, fmc_g, source) + + implicit none + + class (fmc_t), intent (in) :: this + integer, intent (in) :: ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, no_fuel_cat + real, dimension(ifms:ifme, jfms:jfme), intent (in) :: nfuel_cat, fmc_g + character (len = *), intent (in) :: source + integer :: i, j + character (len = 160) :: msg + + + ! Check computed or restored moisture before it reaches spread physics. + ! Prescribed moisture has an FMC object but no initialized class arrays; + ! this common check therefore uses only the supplied effective moisture. + ! Only owned burnable cells require moisture; halo values may be unavailable. + do j = jfts, jfte + do i = ifts, ifte + if (int (nfuel_cat(i, j)) == no_fuel_cat) cycle + if (.not. ieee_is_finite (fmc_g(i, j))) then + write (msg, '(a,a,2(1x,i0))') source, ' fmc_g must be finite at i,j:', i, j + call Stop_simulation (msg) + else if (fmc_g(i, j) < 0.0) then + write (msg, '(a,a,2(1x,i0))') source, ' fmc_g must be nonnegative at i,j:', i, j + call Stop_simulation (msg) + end if + end do + end do + + end subroutine Check_fmc + end module fmc_mod diff --git a/share/fuel_mod.F90 b/share/fuel_mod.F90 index 529da46..6398dc9 100644 --- a/share/fuel_mod.F90 +++ b/share/fuel_mod.F90 @@ -1,5 +1,8 @@ module fuel_mod + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite + use stderrout_mod, only : Stop_simulation + implicit none private @@ -40,6 +43,7 @@ module fuel_mod real, dimension(:), allocatable :: waf contains procedure (Initialization), deferred :: Initialization + procedure, public :: Check_fuels => Check_fuels procedure, public :: Calc_wind_adjustment_factor => Calc_wind_adjustment_factor end type fuel_t @@ -53,6 +57,29 @@ end subroutine Initialization contains + subroutine Check_fuels (this) + + implicit none + + class (fuel_t), intent (in) :: this + integer :: k + character (len = 128) :: msg + + + ! Extinction moisture is a fixed table property, checked once at initialization. + do k = 1, this%n_fuel_cat + if (k == this%no_fuel_cat) cycle + if (.not. ieee_is_finite (this%fuelmce(k))) then + write (msg, '(a,1x,i0)') 'fuelmce must be finite for fuel category:', k + call Stop_simulation (msg) + else if (this%fuelmce(k) <= 0.0) then + write (msg, '(a,1x,i0)') 'fuelmce must be positive for fuel category:', k + call Stop_simulation (msg) + end if + end do + + end subroutine Check_fuels + pure function Crosswalk_from_scottburgan_to_anderson (fuel_cat) result (return_value) implicit none diff --git a/state/state_restart_submod.F90 b/state/state_restart_submod.F90 index c1bf51d..3aef76b 100644 --- a/state/state_restart_submod.F90 +++ b/state/state_restart_submod.F90 @@ -520,6 +520,11 @@ module subroutine Read_restart (this, config_flags) if (allocated (this%ros_param) .and. allocated (this%fuels)) then do ij = 1, this%num_tiles + ! fmc_g is read independently from the checkpoint, not recomputed from + ! fmc_gc. Valid classes cannot establish its validity, and prescribed + ! moisture restarts have no class state; check it before rebuilding ROS. + call this%fmc_param%Check_fmc (this%ifms, this%ifme, this%jfms, this%jfme, this%i_start(ij), this%i_end(ij), & + this%j_start(ij), this%j_end(ij), this%nfuel_cat, this%fuels%no_fuel_cat, this%fmc_g, 'Restart') call this%ros_param%Set_params (this%ifms, this%ifme, this%jfms, this%jfme, this%i_start(ij), this%i_end(ij), & this%j_start(ij), this%j_end(ij), this%fuels, this%nfuel_cat, this%fmc_g) end do @@ -667,6 +672,9 @@ subroutine Read_restart_fmc (file_name) this%ifps, this%ifpe, this%jfps, this%jfpe, NAME_VAR_FMC_GC, & fmc_param%fmc_gc(this%ifps:this%ifpe, 1:n_moisture_classes, this%jfps:this%jfpe)) + call fmc_param%Check_moisture_classes (this%ifms, this%ifme, this%jfms, this%jfme, & + this%ifps, this%ifpe, this%jfps, this%jfpe, this%nfuel_cat, this%fuels%no_fuel_cat, 'Restart') + fmoist_lasttime = 0.0 fmoist_nexttime = 0.0 initialization_time = 0.0 diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index d228ab7..c3cb835 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -22,6 +22,91 @@ foreach(test IN ITEMS tiles_unit tile_namelist_unit observed_perimeter_unit target_link_libraries(${test} PRIVATE firelib) endforeach() +add_executable(moisture_extinction_unit moisture_extinction_unit.F90) +target_link_libraries(moisture_extinction_unit PRIVATE firelib) +foreach(scenario IN ITEMS below edge at above) + add_cfbm_unit_test(NAME moisture_extinction_${scenario}_unit + COMMAND moisture_extinction_unit ${scenario} LABELS fuel moisture ros + EXPECT "moisture extinction checks passed: ${scenario}") + if(OPENMP) + set_tests_properties(moisture_extinction_${scenario}_unit PROPERTIES + ENVIRONMENT "OMP_NUM_THREADS=2;OMP_DYNAMIC=FALSE") + endif() +endforeach() +add_cfbm_unit_test(NAME moisture_extinction_negative_moisture_unit + COMMAND moisture_extinction_unit negative_moisture LABELS fuel moisture ros + EXPECT "STOP:Restart fmc_g must be nonnegative") +foreach(scenario IN ITEMS nan_moisture inf_moisture neginf_moisture) + add_cfbm_unit_test(NAME moisture_extinction_${scenario}_unit + COMMAND moisture_extinction_unit ${scenario} LABELS fuel moisture ros + EXPECT "STOP:Restart fmc_g must be finite") +endforeach() +foreach(scenario IN ITEMS zero_extinction negative_extinction) + add_cfbm_unit_test(NAME moisture_extinction_${scenario}_unit + COMMAND moisture_extinction_unit ${scenario} LABELS fuel moisture ros + EXPECT "STOP:fuelmce must be positive") +endforeach() +foreach(scenario IN ITEMS nan_extinction inf_extinction neginf_extinction) + add_cfbm_unit_test(NAME moisture_extinction_${scenario}_unit + COMMAND moisture_extinction_unit ${scenario} LABELS fuel moisture ros + EXPECT "STOP:fuelmce must be finite") +endforeach() + +add_cfbm_unit_test(NAME namelist_fuelmc_zero_unit + COMMAND namelist_validation_unit fuelmc_zero LABELS namelist moisture + EXPECT "namelist validation accepted: fuelmc_zero") +add_cfbm_unit_test(NAME namelist_fuelmc_negative_unit + COMMAND namelist_validation_unit fuelmc_negative LABELS namelist moisture + EXPECT "STOP:fuelmc_g must be nonnegative") +foreach(scenario IN ITEMS nan inf neginf) + add_cfbm_unit_test(NAME namelist_fuelmc_${scenario}_unit + COMMAND namelist_validation_unit fuelmc_${scenario} LABELS namelist moisture + EXPECT "STOP:fuelmc_g must be finite") +endforeach() + +add_executable(fmc_validation_unit fmc_validation_unit.F90) +target_link_libraries(fmc_validation_unit PRIVATE firelib) +foreach(scenario IN ITEMS valid zero_load restored_zero restored_wet) + add_cfbm_unit_test(NAME fmc_validation_${scenario}_unit + COMMAND fmc_validation_unit ${scenario} LABELS fuel moisture + EXPECT "moisture validation checks passed: ${scenario}") +endforeach() +foreach(quantity IN ITEMS class restored) + foreach(scenario IN ITEMS negative nan inf neginf) + if(quantity STREQUAL "class") + set(field fmc_gc) + set(source Computed) + else() + set(field fmc_g) + set(source Restart) + endif() + if(scenario STREQUAL "negative") + set(requirement nonnegative) + else() + set(requirement finite) + endif() + add_cfbm_unit_test(NAME fmc_validation_${quantity}_${scenario}_unit + COMMAND fmc_validation_unit ${quantity}_${scenario} LABELS fuel moisture + EXPECT "STOP:${source} ${field} must be ${requirement}") + endforeach() +endforeach() +add_cfbm_unit_test(NAME fmc_validation_load_negative_unit + COMMAND fmc_validation_unit load_negative LABELS fuel moisture + EXPECT "STOP:.*fuel load in category.*must be nonegative") +foreach(scenario IN ITEMS nan inf neginf) + add_cfbm_unit_test(NAME fmc_validation_load_${scenario}_unit + COMMAND fmc_validation_unit load_${scenario} LABELS fuel moisture + EXPECT "STOP:Moisture averaging fuel load must be finite for category,class:") +endforeach() +add_cfbm_unit_test(NAME fmc_validation_load_overflow_unit + COMMAND fmc_validation_unit load_overflow LABELS fuel moisture + EXPECT "STOP:Moisture averaging total fuel load exceeds representable range for category:") +if(OPENMP) + get_property(fmc_tests DIRECTORY PROPERTY TESTS) + list(FILTER fmc_tests INCLUDE REGEX "^fmc_validation_") + set_tests_properties(${fmc_tests} PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=2;OMP_DYNAMIC=FALSE") +endif() + # Check the production wind routine's array bounds even in Release builds. # Keep its test module separate from firelib's normal compilation products. target_sources(wind_order_unit PRIVATE ${PROJECT_SOURCE_DIR}/io/coupling_mod.F90) diff --git a/tests/unit/fmc_validation_unit.F90 b/tests/unit/fmc_validation_unit.F90 new file mode 100644 index 0000000..fdb1c87 --- /dev/null +++ b/tests/unit/fmc_validation_unit.F90 @@ -0,0 +1,127 @@ +program fmc_validation_unit + + ! Exercise validation at the moisture-component boundary with two owned tiles. + ! Invalid fuel loads must be rejected when fixed averaging weights are built. + ! Invalid class initial values must be diagnosed after the parallel update and + ! before averaging, including a live class with zero weight in Anderson fuel 1. + ! Independently restored fmc_g must be checked even when every class is valid. + ! Restored zero and wetter-than-extinction states remain valid and unchanged; + ! unavailable halos and nonburnable cells must not enter the validity checks. + + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite, ieee_value, ieee_quiet_nan, & + ieee_positive_inf, ieee_negative_inf + use fmc_wrffire_mod, only : fmc_wrffire_t + use fuel_anderson_mod, only : fuel_anderson_t + use ros_wrffire_mod, only : ros_wrffire_t +#ifdef DM_PARALLEL + use mpi +#endif + + implicit none + + type (fmc_wrffire_t) :: moisture + type (fuel_anderson_t) :: fuels + type (ros_wrffire_t) :: ros + real :: rain(0:5,0:3), t(0:5,0:3), q(0:5,0:3), p(0:5,0:3) + real :: rain_old(0:5,0:3), t_old(0:5,0:3), q_old(0:5,0:3), p_old(0:5,0:3) + real :: rh(0:5,0:3), fmc(0:5,0:3), category(0:5,0:3), nan, invalid, restored + integer :: i, j, ierr + integer :: i_start(2) = [1,1], i_end(2) = [4,4], j_start(2) = [1,2], j_end(2) = [1,2] + character (len = 40) :: scenario + +#ifdef DM_PARALLEL + call MPI_Init (ierr) +#endif + call get_command_argument (1, scenario) + nan = ieee_value (0.0, ieee_quiet_nan) + select case (trim (scenario)) + case ('valid', 'zero_load', 'load_overflow', 'restored_zero', 'restored_wet') + invalid = 0.0 + case ('class_negative', 'load_negative', 'restored_negative') + invalid = -1.0 + case ('class_nan', 'load_nan', 'restored_nan') + invalid = nan + case ('class_inf', 'load_inf', 'restored_inf') + invalid = ieee_value (0.0, ieee_positive_inf) + case ('class_neginf', 'load_neginf', 'restored_neginf') + invalid = ieee_value (0.0, ieee_negative_inf) + case default + error stop 'unknown moisture validation scenario' + end select + call fuels%Initialization (1.0) + ! Perturb the input to initialization rather than corrupting fixed weights + ! after they have been validated. No moisture advance should be needed to fail. + if (index (scenario, 'load_') == 1) fuels%fgi_1h(1) = invalid + if (scenario == 'load_overflow') then + fuels%fgi_1h(1) = huge (1.0) + fuels%fgi_10h(1) = huge (1.0) + else if (scenario == 'zero_load') then + fuels%fgi_1h(1) = 0.0 + fuels%fgi_10h(1) = 0.0 + fuels%fgi_100h(1) = 0.0 + fuels%fgi_1000h(1) = 0.0 + fuels%fgi_live(1) = 0.0 + end if + call moisture%Init (fuels, 0.08, 0.30, 0, 5, 0, 3, 0, 4.0) + if (index (scenario, 'load_') == 1) error stop 'invalid averaging fuel load was accepted' + call ros%Init (0, 5, 0, 3) + category = nan + category(1:4,1:2) = 1.0 + category(4,2) = real (fuels%no_fuel_cat) + fmc = nan + moisture%fmc_gc = nan + + if (index (scenario, 'restored_') == 1) then + restored = 0.0 + if (scenario == 'restored_wet') restored = 2.0 + fmc(1:4,1:2) = restored + moisture%fmc_gc(1:4,:,1:2) = restored + fmc(4,2) = nan + moisture%fmc_gc(4,:,2) = nan + call moisture%Check_moisture_classes (0,5,0,3,1,4,1,2,category,fuels%no_fuel_cat,'Restart') + ! The restart reader loads fmc_g independently of the valid class array and + ! validates it through the common binding inherited from the FMC base type. + if (scenario /= 'restored_zero' .and. scenario /= 'restored_wet') fmc(2,1) = invalid + call moisture%Check_fmc (0,5,0,3,1,4,1,2,category,fuels%no_fuel_cat,fmc,'Restart') + if (scenario /= 'restored_zero' .and. scenario /= 'restored_wet') error stop 'invalid restored moisture was accepted' + if (any (fmc(1:3,1:2) /= restored)) error stop 'restored moisture was modified' + if (any (moisture%fmc_gc(1:3,:,1:2) /= restored)) error stop 'restored class moisture was modified' + else + rain = nan + t = nan + q = nan + p = nan + rain_old = nan + t_old = nan + q_old = nan + p_old = nan + rain(1:4,1:2) = 0.0 + t(1:4,1:2) = 300.0 + q(1:4,1:2) = 0.008 + p(1:4,1:2) = 90000.0 + + if (index (scenario, 'class_') == 1) moisture%fmc_gc_initial_value(5) = invalid + + call moisture%Advance_fmc_model (0,4.0,1,4.0,0,5,0,3,i_start,i_end,j_start,j_end,2, & + rain,t,q,p,rain_old,t_old,q_old,p_old,rh,0.08,fmc,category,fuels,ros) + if (scenario /= 'valid' .and. scenario /= 'zero_load') error stop 'invalid moisture update was accepted' + if (scenario == 'zero_load') then + ! An empty category retains the existing all-zero weights and mean. + if (any (moisture%fmc_gw(1,:) /= 0.0)) error stop 'nonzero weights for empty fuel category' + if (any (fmc(1:3,1:2) /= 0.0)) error stop 'nonzero mean for empty fuel category' + end if + do j = 1, 2 + do i = 1, 4 + if (int (category(i,j)) == fuels%no_fuel_cat) cycle + if (.not. ieee_is_finite (fmc(i,j)) .or. fmc(i,j) < 0.0) error stop 'invalid effective moisture' + if (.not. ieee_is_finite (ros%iboros(i,j)) .or. ros%iboros(i,j) < 0.0) error stop 'invalid reaction' + end do + end do + end if + + print *, 'moisture validation checks passed: ', trim (scenario) +#ifdef DM_PARALLEL + call MPI_Finalize (ierr) +#endif + +end program fmc_validation_unit diff --git a/tests/unit/moisture_extinction_unit.F90 b/tests/unit/moisture_extinction_unit.F90 new file mode 100644 index 0000000..7665be5 --- /dev/null +++ b/tests/unit/moisture_extinction_unit.F90 @@ -0,0 +1,217 @@ +program moisture_extinction_unit + + ! Check the Rothermel extinction correction for all 13 Anderson fuel categories. + ! Moisture is water mass per dry fuel mass and must remain unchanged, including + ! valid values above extinction. Reaction intensity, its residence-time product + ! (iboros), and base spread must be exactly zero at and above extinction, while + ! valid values below extinction retain the existing damping polynomial. + ! Downstream checks cover spread and flame length, including the independent + ! chaparral spread law and nonburnable category 14. + + use, intrinsic :: iso_fortran_env, only : REAL64 + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite, ieee_value, ieee_quiet_nan, & + ieee_positive_inf, ieee_negative_inf + use constants_mod, only : CMBCNST, CONVERT_J_PER_KG_TO_BTU_PER_POUND + use fuel_anderson_mod, only : fuel_anderson_t + use fmc_wrffire_mod, only : fmc_wrffire_t + use ros_wrffire_mod, only : ros_wrffire_t + use fire_physics_mod, only : Calc_flame_length +#ifdef DM_PARALLEL + use mpi +#endif + + implicit none + + type(fuel_anderson_t) :: fuels + type(fmc_wrffire_t) :: moisture_params + type(ros_wrffire_t) :: params + real :: category(0:15,0:3), moisture(0:15,0:3), original(1:14,1:2) + real :: spread(0:15,0:3), flame(0:15,0:3), front(0:15,0:3), area(0:15,0:3) + real :: ir, expected_ir, expected_r0, expected_iboros, factor, nan, wind, slope + integer :: i, j, tile, sample, direction, ierr + character(len=40) :: scenario + +#ifdef DM_PARALLEL + call MPI_Init(ierr) +#endif + call get_command_argument(1, scenario) + ! Prescribed moisture allocates an FMC object without initializing class state. + ! Its inherited Check_fmc binding must validate the supplied array in this case. + call fuels%Initialization(1.0) + call params%Init(0,15,0,3) + nan = ieee_value(0.0, ieee_quiet_nan) + ! Unavailable halos and the nonburnable cell must not enter fuel physics. + category = nan + moisture = nan + do i = 1, 14 + category(i,1:2) = real(i) + end do + area = 0.5 + + ! Sample dry-to-moist fuel below extinction and the nearest representable values + ! on each side of the threshold. These boundary cases distinguish an exact + ! extinction cutoff from polynomial cancellation or an unintended moisture cap. + do sample = 1, 4 + select case (trim(scenario)) + case ('below') + factor = 0.25 * (sample - 1) + case ('edge') + factor = 1.0 + case ('at') + factor = 1.0 + case ('above') + factor = real(sample + 1) + case default + factor = 0.5 + end select + do i = 1, 13 + moisture(i,1:2) = factor * fuels%fuelmce(i) + if (scenario == 'edge') moisture(i,1:2) = nearest(fuels%fuelmce(i), -1.0) + if (scenario == 'above' .and. sample == 1) moisture(i,1:2) = nearest(fuels%fuelmce(i), 1.0) + end do + ! Invalid inputs are rejected by the owning fuel/moisture layer, before ROS. + ! CTest requires the specific diagnostic; a generic failure is insufficient. + select case (trim(scenario)) + case ('negative_moisture') + moisture(1,1) = -0.01 + case ('nan_moisture') + moisture(1,1) = nan + case ('inf_moisture') + moisture(1,1) = ieee_value(0.0, ieee_positive_inf) + case ('neginf_moisture') + moisture(1,1) = ieee_value(0.0, ieee_negative_inf) + case ('zero_extinction') + fuels%fuelmce(1) = 0.0 + case ('negative_extinction') + fuels%fuelmce(1) = -0.1 + case ('nan_extinction') + fuels%fuelmce(1) = nan + case ('inf_extinction') + fuels%fuelmce(1) = ieee_value(0.0, ieee_positive_inf) + case ('neginf_extinction') + fuels%fuelmce(1) = ieee_value(0.0, ieee_negative_inf) + case ('below', 'edge', 'at', 'above') + case default + error stop 'unknown extinction scenario' + end select + original = moisture(1:14,1:2) + if (index(scenario, 'moisture') > 0 .or. index(scenario, 'extinction') > 0) then + ! Match initialization/restart validation on the initial thread. Fuel-table + ! mutations emulate an invalid definition before Check_fuels is called. + call fuels%Check_fuels () + call moisture_params%Check_fmc (0,15,0,3,1,14,1,2,category,fuels%no_fuel_cat,moisture,'Restart') + error stop 'invalid input was accepted' + end if + call moisture_params%Check_fmc (0,15,0,3,1,14,1,2,category,fuels%no_fuel_cat,moisture,'Restart') + ! Exercise disjoint owned tiles, concurrently in OpenMP builds. Set_params + ! also reconstructs derived spread parameters after restart; checkpoint I/O + ! and continuous-versus-restarted evolution are covered by regression cases. + !$omp parallel do private(tile) + do tile = 1, 2 + call params%Set_params(0,15,0,3,1,14,tile,tile,fuels,category,moisture) + end do + !$omp end parallel do + if (any(moisture(1:13,1:2) /= original(1:13,1:2))) error stop 'physical moisture was modified' + + do i = 1, 13 + call Reference_reaction(i, moisture(i,1), expected_ir, expected_r0, expected_iboros) + do j = 1, 2 + ! Recover reaction intensity from its stored residence-time product. + ir = params%iboros(i,j) / (1055.0 / (0.3048**2 * 60.0) * 1.e-3 * (60.0 * 384.0 / fuels%savr(i))) + call Check_value(ir, expected_ir, 'reaction intensity') + call Check_value(params%iboros(i,j), expected_iboros, 'iboros') + call Check_value(params%r_0(i,j), expected_r0, 'base spread') + if (scenario == 'at' .or. scenario == 'above') then + if (ir /= 0.0 .or. params%iboros(i,j) /= 0.0 .or. params%r_0(i,j) /= 0.0) & + error stop 'nonzero reaction at or above extinction' + end if + end do + ! Hexadecimal values permit an exact before/after comparison below extinction. + if (scenario == 'below') write(*,'(a,2i3,2(1x,z8.8))') 'below signature ', sample, i, & + transfer(params%iboros(i,1),0), transfer(params%r_0(i,1),0) + end do + if (any(params%iboros(14,1:2) /= 0.0) .or. any(params%r_0(14,1:2) /= 0.0)) & + error stop 'nonburnable cell has a reaction' + + do direction = 1, 3 + wind = 0.0 + slope = 0.0 + if (direction > 1) wind = 3.0 + if (direction == 3) slope = 0.2 + do j = 1, 2 + do i = 1, 14 + spread(i,j) = params%Calc_ros(0,15,0,3,i,j,1.0,0.0,wind,0.0,slope,0.0) + end do + end do + call Calc_flame_length(1,14,1,2,0,15,0,3,spread,params%iboros,flame,front,area) + if (any(.not. ieee_is_finite(spread(1:14,1:2))) .or. any(spread(1:14,1:2) < 0.0)) & + error stop 'invalid spread rate' + if (any(.not. ieee_is_finite(flame(1:14,1:2))) .or. any(flame(1:14,1:2) < 0.0)) & + error stop 'invalid flame length' + if (any(front(1:14,1:2) /= spread(1:14,1:2))) error stop 'front spread differs' + do i = 1, 13 + if (moisture(i,1) >= fuels%fuelmce(i)) then + if (any(flame(i,1:2) /= 0.0)) error stop 'nonzero extinction flame length' + ! Anderson 4 uses the independent empirical chaparral spread law. + if (fuels%ichap(i) == 0 .and. any(spread(i,1:2) /= 0.0)) error stop 'Rothermel spread at extinction' + else + if (any(flame(i,1:2) <= 0.0 .and. spread(i,1:2) > 0.0)) error stop 'missing below-extinction flame' + end if + if (fuels%ichap(i) /= 0) call Check_value(spread(i,1), & + min(max(1.2974 * wind**1.41, 0.03333) * wind / (wind + tiny(wind)), 6.0), 'chaparral spread') + end do + if (any(flame(14,1:2) /= 0.0) .or. any(spread(14,1:2) /= 0.0)) error stop 'nonburnable propagation' + end do + end do + print *, 'moisture extinction checks passed: ', trim(scenario) +#ifdef DM_PARALLEL + call MPI_Finalize(ierr) +#endif + +contains + + subroutine Reference_reaction(k, m, intensity, base_ros, integrated) + integer, intent(in) :: k + real, intent(in) :: m + real, intent(out) :: intensity, base_ros, integrated + real(REAL64) :: load, depth, packing, optimum, sigma, a, rate, damping, flux, heat + real :: ratio + ! Evaluate the dimensional Rothermel expressions in double precision. + ! Retain single-precision damping to test the existing polynomial even + ! one representable moisture value below extinction, where it cancels. + ratio = m / fuels%fuelmce(k) + damping = 0.0_REAL64 + if (m < fuels%fuelmce(k)) damping = 1.0 - 2.59 * ratio + 5.11 * ratio**2 - 3.52 * ratio**3 + sigma = fuels%savr(k) + load = real(fuels%fgi(k), REAL64) / (1.0_REAL64 + m) * 0.3048_REAL64**2 * 2.205_REAL64 + depth = real(fuels%fueldepthm(k), REAL64) / 0.3048_REAL64 + packing = load / (depth * fuels%fueldens(k)) + optimum = 3.348_REAL64 * sigma**(-0.8189_REAL64) + a = 1.0_REAL64 / (4.774_REAL64 * sigma**0.1_REAL64 - 7.27_REAL64) + rate = sigma**1.5_REAL64 / (495.0_REAL64 + 0.0594_REAL64 * sigma**1.5_REAL64) * & + (packing / optimum)**a * exp(a * (1.0_REAL64 - packing / optimum)) + heat = CMBCNST * CONVERT_J_PER_KG_TO_BTU_PER_POUND + intensity = rate * load / (1.0_REAL64 + fuels%st(k)) * heat * damping * & + 0.174_REAL64 * real(fuels%se(k), REAL64)**(-0.19_REAL64) + integrated = intensity * 1055.0_REAL64 / (0.3048_REAL64**2 * 60.0_REAL64) * & + 1.e-3_REAL64 * (60.0_REAL64 * 384.0_REAL64 / sigma) + flux = exp((0.792_REAL64 + 0.681_REAL64 * sqrt(sigma)) * (packing + 0.1_REAL64)) / & + (192.0_REAL64 + 0.2595_REAL64 * sigma) + base_ros = intensity * flux / (load / depth * exp(-138.0_REAL64 / sigma) * & + (250.0_REAL64 + 1116.0_REAL64 * m)) + end subroutine Reference_reaction + + subroutine Check_value(value, expected, quantity) + real, intent(in) :: value, expected + character(len=*), intent(in) :: quantity + if (.not. ieee_is_finite(value)) error stop 'nonfinite reaction quantity' + if (value < 0.0) error stop 'negative reaction quantity' + ! Bound accumulated single-precision rounding against the independent + ! double-precision calculation, without a nonzero absolute-error floor. + if (abs(value - expected) > 64.0 * epsilon(value) * abs(expected)) then + print *, quantity, value, expected + error stop 'reaction reference mismatch' + end if + end subroutine Check_value + +end program moisture_extinction_unit diff --git a/tests/unit/namelist_validation_unit.F90 b/tests/unit/namelist_validation_unit.F90 index 56e716b..19e234b 100644 --- a/tests/unit/namelist_validation_unit.F90 +++ b/tests/unit/namelist_validation_unit.F90 @@ -1,10 +1,11 @@ program namelist_validation_unit ! Exercise configuration checks for atmospheric intervals, ignition counts, - ! and perimeter activation times. CTest selects valid or invalid scenarios + ! perimeter activation times, and prescribed moisture. CTest selects scenarios ! and checks the expected rejection diagnostics. Atmospheric-only scenarios ! validate their timing without requiring the fire component's settings. + use, intrinsic :: ieee_arithmetic, only : ieee_value, ieee_quiet_nan, ieee_positive_inf, ieee_negative_inf use namelist_mod, only : namelist_t, FIRE_MAX_IGNITIONS_IN_NAMELIST #ifdef DM_PARALLEL use mpi @@ -34,6 +35,16 @@ program namelist_validation_unit if (scenario == 'atm_missing') config%interval_atm = -1 if (scenario == 'atm_zero_dt') config%dt = 0.0 case ('aligned') + case ('fuelmc_zero') + config%fuelmc_g = 0.0 + case ('fuelmc_negative') + config%fuelmc_g = -0.01 + case ('fuelmc_nan') + config%fuelmc_g = ieee_value (0.0, ieee_quiet_nan) + case ('fuelmc_inf') + config%fuelmc_g = ieee_value (0.0, ieee_positive_inf) + case ('fuelmc_neginf') + config%fuelmc_g = ieee_value (0.0, ieee_negative_inf) case ('reversed') config%dt = 8.0 config%interval_atm = 4