diff --git a/TESTING.md b/TESTING.md index 396b7ea..9655416 100644 --- a/TESTING.md +++ b/TESTING.md @@ -214,6 +214,15 @@ absent. `run.bash` runs it after the examples. python3 tests/python/test_validate.py ``` +`tests/python/test_farsite_reference.py` compares the Farsite model with the +FARSITE reference code (`mbedward/farsite`) over the 13 fuel models. It needs +a C++ compiler and is not run by CI. Gaps on slope and at flank/back are +reported; `--strict` makes the slope gap fail. + +```bash +./.venv/bin/python tests/python/test_farsite_reference.py +``` + `tests/python/` also holds `idealizedwind.py`, `farsite_flat.py` and `percolation.py`. They are examples rather than tests — they produce plots and assert nothing — and are not run anywhere. They are the closest thing to diff --git a/tests/python/README.md b/tests/python/README.md index a3635d2..045431a 100644 --- a/tests/python/README.md +++ b/tests/python/README.md @@ -55,6 +55,20 @@ python tests/python/test_moisture_invariants.py # -v for every probe `--model NAME` and `--test NAME`, both repeatable, narrow it down. +### `test_farsite_reference.py` + +Compares the `Farsite` model with FARSITE's own code for the 13 standard fuel +models, four midflame winds, three slopes and three moisture sets. The +reference is `farsite_reference.cpp`, FARSITE's spread rate and ellipse copied +from mbedward/farsite and compiled by the test with the system C++ compiler. +Fails when the head rate of spread differs from the reference; flank, back and +the gap on a slope are printed, not asserted (`--strict` asserts the slope +too). Needs a C++ compiler besides `pyforefire`. Manually: + +```bash +python tests/python/test_farsite_reference.py # -v for every case +``` + ### `test_runner.py` Checks `pyforefire.ForeFireRunner`: eight simulations, each in its own worker diff --git a/tests/python/farsite_reference.cpp b/tests/python/farsite_reference.cpp new file mode 100644 index 0000000..f2ec682 --- /dev/null +++ b/tests/python/farsite_reference.cpp @@ -0,0 +1,500 @@ +/** + * @file farsite_reference.cpp + * @brief Surface spread of the FARSITE reference, for comparison with ForeFire. + * @license Public domain, as the FARSITE source it is copied from (17 USC 105), + * except contributions from outside the RMRS lab, under the BSD + * license below. + * + * The arithmetic below is copied from mbedward/farsite, commit + * 4537d60ab013fa91d0b49d8efb385ffb1b3ddd13 (https://github.com/mbedward/farsite). + * Every block names the file and lines it comes from. Only the plumbing around + * it changed: the fuel comes from a local table instead of Farsite5, the slope + * arrives as a tangent instead of whole degrees, and the wind always blows + * straight upslope, so the wind/slope vector is a plain sum. + * + * Reads one case per line on stdin: + * fuel ones tens hundreds liveh livew wind_mph slope_tan + * (moistures as fractions of dry weight, wind at midflame) and writes one line + * per case on stdout: + * rateo phiw phis phiew head lb flank back + * rateo is the no-wind no-slope spread, head/flank/back the spread rates at + * the head, flank and back of the ellipse, all in m/min. + * + * Build: c++ -O2 -o farsite_reference farsite_reference.cpp + * + * Copyright (c) 2016 The FARSITE Authors. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above + * copyright notice, this list of conditions and the following disclaimer + * in the documentation and/or other materials provided with the + * distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include +#include +#include + +using namespace std; + +static const double PI = acos(-1.0); // Farsite5.cpp:48 + +inline double pow2(double in) { return in * in; } // Farsite5.h:40-43 + +// Farsite5.h:183-202 +typedef struct +{ + long number; + char code[8]; + double h1; + double h10; + double h100; + double lh; + double lw; + long dynamic; + long sav1; + long savlh; + long savlw; + double depth; + double xmext; + double heatd; + double heatl; + char desc[256]; +}NewFuel; + +// Farsite5.cpp:6907-7128, Farsite5::InitializeNewFuel, the 13 standard models +static const NewFuel newfuels[14] = { + {0, "", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ""}, + {1, "FM1", 0.74, 0, 0, 0, 0, 0, 3500, 1800, 1500, 1.0, 0.12, 8000, 8000, "Short Grass"}, + {2, "FM2", 2.0, 1.0, 0.5, 0, 0.5, 0, 3000, 1800, 1500, 1.0, 0.15, 8000, 8000, "Timber Grass/Understory"}, + {3, "FM3", 3.01, 0, 0, 0, 0, 0, 1500, 1800, 1500, 2.5, 0.25, 8000, 8000, "Tall Grass"}, + {4, "FM4", 5.01, 4.01, 2.0, 0, 5.01, 0, 2000, 1800, 1500, 6.0, 0.20, 8000, 8000, "Chaparral"}, + {5, "FM5", 1.0, 0.5, 0, 0, 2.0, 0, 2000, 1800, 1500, 2.0, 0.20, 8000, 8000, "Short Brush"}, + {6, "FM6", 1.5, 2.5, 2.0, 0, 0, 0, 1750, 1800, 1500, 2.5, 0.25, 8000, 8000, "Dormant Brush"}, + {7, "FM7", 1.13, 1.87, 1.5, 0, 0.37, 0, 1550, 1800, 1500, 2.5, 0.40, 8000, 8000, "Southern Rough"}, + {8, "FM8", 1.5, 1.0, 2.5, 0, 0, 0, 2000, 1800, 1500, 0.2, 0.30, 8000, 8000, "Closed Timber Litter"}, + {9, "FM9", 2.92, 0.41, 0.15, 0, 0, 0, 2500, 1800, 1500, 0.2, 0.25, 8000, 8000, "Hardwood Litter"}, + {10, "FM10", 3.01, 2.0, 5.01, 0, 2.0, 0, 2000, 1800, 1500, 1.0, 0.25, 8000, 8000, "Timber Litter/Understory"}, + {11, "FM11", 1.5, 4.51, 5.51, 0, 0, 0, 1500, 1800, 1500, 1.0, 0.15, 8000, 8000, "Light Slash"}, + {12, "FM12", 4.01, 14.03, 16.53, 0, 0, 0, 1500, 1800, 1500, 2.3, 0.20, 8000, 8000, "Medium Slash"}, + {13, "FM13", 7.01, 23.04, 28.05, 0, 0, 0, 1500, 1800, 1500, 3.0, 0.25, 8000, 8000, "Heavy Slash"}, +}; + +static bool GetNewFuel(long number, NewFuel* newfuel) +{ + if (number < 1 || number > 13) + return false; + memcpy(newfuel, &newfuels[number], sizeof(NewFuel)); + return true; +} + +// The members of Mechanix (fsx4.hpp:483-522) that the functions below use. +struct Mechanix +{ + double lb_ratio, hb_ratio, rateo; + double b, part1; + double m_ones, m_tens, m_hundreds, m_livew, m_liveh; + double phiw, phis, phiew, LocalWindDir; + double slopespd; + double vecspeed, ivecspeed, m_windspd; + double fros, vecros; + double react, savx; + double head, back, flank; + + double spreadrate(double slope_tan, double windspd, int fuel); + void VecSurf(double slope_tan, double aspectf); + double vectorspd(double* ivecspeed, double aspectf, long CrownWinds); + double headback(); + void ellipse(double iros, double wspeed); +}; + +// fsxwmech.cpp:784-1045, Mechanix::spreadrate (the one compiled in; the two +// versions after it in that file are commented out). Changed: the slope +// argument is already a tangent, and GetNewFuel is the local table. +double Mechanix::spreadrate(double slope_tan, double windspd, int fuel) +{// Rothermel spread equation based on BEHAVE source code, + // support for dynamic fuel models added 10/13/2004 + long i, j, ndead=0, nlive=0; + double depth; + double seff[4][2]={{0.01,0.01},{0.01,0.01},{0.01,0},{0.01,0.0}}; //mineral content + double wtfact, fined=0.0, finel=0.0, wmfd=0.0, fdmois=0.0, w=0.0, wo=0.0, beta; + double rm, sigma=0.0, rhob, sum3=0.0, betaop, rat, aa, gammax=0.0, gamma=0, wind=0; + double xir, rbqig=0, xi=0, c, e, slopex=0; + double ewind=0, wlim, sum1=0, sum2=0; + + double mois[4][2]= // fraction of oven-dry weight + { {m_ones, m_liveh}, // use Mechanix Class copies of FE data + {m_tens, m_livew}, + {m_hundreds, 0.0}, + {m_ones, 0.0} + }; + double gx[5]={0.0, 0.0, 0.0, 0.0, 0.0}; + double wn[4][2]={{0,0},{0,0},{0,0},{0,0}}; + double qig[4][2]={{0,0},{0,0},{0,0},{0,0}}; + double a[4][2]={{0,0},{0,0},{0,0},{0,0}}; + double f[4][2]={{0,0},{0,0},{0,0},{0,0}}; + double g[4][2]={{0,0},{0,0},{0,0},{0,0}}; + double ai[2]={0,0}; + double fi[2]={0,0}; + double hi[2]={0,0}; + double se[2]={0,0}; + double xmf[2]={0,0}; + double si[2]={0,0}; + double wni[2]={0,0}; + double etam[2]={0,0}; + double etas[2]={0,0}; + double rir[2]={0,0}; + + NewFuel newfuel; + + memset(&newfuel, 0x0, sizeof(NewFuel)); + if(!GetNewFuel(fuel, &newfuel)) // will get all models, 13, 40, and custom + { rateo=0.0; + + return rateo; + } + + // count number of fuels + if(newfuel.h1) ndead++; + if(newfuel.h10) ndead++; + if(newfuel.h100) ndead++; + if(newfuel.lh) nlive++; + if(newfuel.lw) nlive++; + + if(nlive==0) + newfuel.dynamic=0; + + if(nlive>0) + nlive=2; // boost to max number + if(ndead>0) + ndead=4; + + double nclas[2]={(double) ndead,(double) nlive}; // # of dead & live fuel classes + double xmext[2]={newfuel.xmext, 0}; + + double load[4][2]= // tons per acre, later converted to lb/ft2 + { {newfuel.h1, newfuel.lh}, + {newfuel.h10, newfuel.lw}, + {newfuel.h100, 0.0}, + {0.0, 0.0}, + }; + + depth=newfuel.depth; + + //------------------------------------------------------------------------- + // do the dynamic load transfer + //------------------------------------------------------------------------- + if(newfuel.dynamic) + { if(mois[0][1]<0.30) // if live herbaceous is less than 30.0 + { load[3][0]=load[0][1]; + load[0][1]=0.0; + } + else if(mois[0][1]<1.20) + { load[3][0]=load[0][1]*(1.20-mois[0][1])/0.9; + load[0][1]-=load[3][0]; + } + } + //------------------------------------------------------------------------- + + double sav[4][2]= // 1/ft + { {(double) newfuel.sav1, (double) newfuel.savlh}, + {109.0, (double) newfuel.savlw}, + {30.0, 0.0}, + {(double) newfuel.savlh, 0.0}, + }; + + double heat[4][2]= + { {newfuel.heatd, newfuel.heatl}, + {newfuel.heatd, newfuel.heatl}, + {newfuel.heatd, 0.0}, + {newfuel.heatd, 0.0}, + }; + + wind=windspd*88.0; // ft/minute + slopex=slope_tan; // was tan((double) slope/180.0*PI), slope in whole degrees + + // fuel weighting factors + for(i=0; i<2; i++) + { for(j=0; j0.0) + f[j][i]=a[j][i]/ai[i]; + else + f[j][i]=0.0; + } + memset(gx, 0x0, 5*sizeof(double)); + for(j=0; j=1200.0) + gx[0]+=f[j][i]; + else if(sav[j][i]>=192.0) + gx[1]+=f[j][i]; + else if(sav[j][i]>=96.0) + gx[2]+=f[j][i]; + else if(sav[j][i]>=48.0) + gx[3]+=f[j][i]; + else if(sav[j][i]>=16.0) + gx[4]+=f[j][i]; + } + for(j=0; j=1200.0) + g[j][i]=gx[0]; + else if(sav[j][i]>=192.0) + g[j][i]=gx[1]; + else if(sav[j][i]>=96.0) + g[j][i]=gx[2]; + else if(sav[j][i]>=48.0) + g[j][i]=gx[3]; + else if(sav[j][i]>=16.0) + g[j][i]=gx[4]; + else + g[j][i]=0.0; + } + + } + } + fi[0]=ai[0]/(ai[0]+ai[1]); + fi[1]=1.0-fi[0]; + + /* no need for this, because extinction moistures are assigned */ + /* as on last page of Burgan and Rothermel 1984 */ + /* rhob=(wo/depth); + beta=rhob/32; + xmext[0]=.12+4.*beta; + */ + + //moisture of extinction + if(nclas[1]!=0) + { for(j=0; j= xmext[i]) + etam[i]=0; + etas[i]=0.174/(pow(se[i],0.19)); + if(etas[i]>1.0) + etas[i]=1.0; + sigma=sigma+fi[i]*si[i]; + rir[i]=wni[i]*hi[i]*etas[i]*etam[i]; + } + } + + /* final calculations */ + rhob=sum1/depth; + beta=sum2/depth; + betaop=3.348/pow(sigma,0.8189); + rat=beta/betaop; + aa=133.0/pow(sigma,0.7913); + //gammax=pow(sigma,1.5)/(495.0+0.0594*pow(sigma,1.5)); + gammax=(sigma*sqrt(sigma))/(495.0+0.0594*sigma*sqrt(sigma)); + gamma=gammax*pow(rat,aa)*exp(aa*(1.0-rat)); + xir=gamma*(rir[0]+rir[1]); + rbqig=rhob*sum3; + xi=exp((0.792+0.681*sqrt(sigma))*(beta+0.1))/(192.0+0.2595*sigma); +/* flux=xi*xir;*/ + rateo=xir*xi/rbqig; /* this is in English units */ + + phis=5.275*pow(beta,-0.3)*pow2(slopex); + c=7.47*exp(-0.133*pow(sigma,0.55)); + b=0.02526*pow(sigma,0.54); + e=0.715*exp(-0.000359*sigma); + part1=c*pow(rat,-e); + phiw=pow(wind,b)*part1; + + wlim=0.9*xir; + + if(phis>0.0) + { if(phis>wlim) // can't have inifinite windspeed + phis=wlim; + slopespd=pow((phis/part1),1.0/b)/88.0; // converts phis to windspd in mph + } + else + slopespd=0.0; + +/* maximum windspeed effect on ros*/ + phiew=phiw+phis; + ewind=pow(((phiew*pow(rat,e))/c),1.0/b); + + if(ewind>wlim) + { ewind=wlim; + phiew=c*pow(wlim,b)*pow(rat,-e); + } + savx=sigma; // SAVX, REACT IN MECHANIX.PUBLIC: + react=xir*0.189275; // convert btu/f2/min to kW/m2 + rateo=rateo*0.30480060960; // convert from f/min to m/min + + return rateo; +} + +// fsxwmech.cpp:443-466, MechCalls::VecSurf. Changed: ld.slope > 0 tests the +// tangent, and TransformWindDir is skipped because the wind blows upslope +// (fsxwmech.cpp:96, a wind direction of -1 means "aligned with the aspect"). +void Mechanix::VecSurf(double slope_tan, double aspectf) +{ + // calculates wind/slope vector for surface fire + if (slope_tan > 0) + { + LocalWindDir = aspectf; + vecspeed = vectorspd(&ivecspeed, aspectf, 0); + if (vecspeed > phiew) + vecspeed = phiew; // vecspeed is still dimensionless + vecros = (fros * (1.0 + vecspeed)); // wind/slope vectored spread rate + } + else + { + ivecspeed = m_windspd; + if (phiw > phiew) + phiw = phiew; + vecspeed = phiw; + vecros = (fros * (1 + phiw)); // winddriven spread rate in wind directn + } +} + +// fsxwmech.cpp:1726-1760, Mechanix::vectorspd, surface branch only +double Mechanix::vectorspd(double* ivecspeed, double aspectf, long CrownWinds) +{ + // determines the direction and strength of windXslope vector + double angle, vecspeed; + + angle = fabs(LocalWindDir - aspectf); + if (angle != PI) // dimensionless coefs + { + vecspeed = sqrt(pow2(phiw) + + pow2(phis) + + (2.0 * phiw * phis * cos(angle))); + *ivecspeed = sqrt(pow2(m_windspd) + + pow2(slopespd) + + (2.0 * m_windspd * slopespd * cos(angle))); + } // actual windspeed vector, not dimensionless coefs + else + { + vecspeed = fabs(phis - phiw); + *ivecspeed = fabs(slopespd - m_windspd); + } + (void) CrownWinds; + + return vecspeed; +} + +// fsxwmech.cpp:1630-1639, Mechanix::headback +double Mechanix::headback() +{ + // calculates Alexanders head/backing ratio assuming fire origin at focus equation [16] + double hb_ratio, part; + + part = sqrt(pow2(lb_ratio) - 1); + hb_ratio = (lb_ratio + part) / (lb_ratio - part); + + return hb_ratio; +} + +// fsxwmech.cpp:1642-1666, Mechanix::ellipse. ConstantBackingSpreadRate is +// false by default (Farsite5.cpp:291). +void Mechanix::ellipse(double iros, double wspeed) +{ + // ANDERSON 1983 FOR MIDFLAME WINDSPEED + lb_ratio = .936 * exp(.1147 * wspeed) + + .461 * exp(-.0692 * wspeed) - + .397; + + if (lb_ratio > 8.0) // maximum eccentricity + lb_ratio = 8.0; + hb_ratio = headback(); // subt 0.397 to make circle at 0 windspd (my own addition) + head = iros; + back = head / hb_ratio; // backspread from hb_ratio from Alexander 1985 + flank = ((head + back) / lb_ratio) / 2.0; + head = (head + back) / 2.0; + back = head - back; // modifies rates of spread for elliptical dimensions +} + +int main() +{ + int fuel; + double ones, tens, hundreds, liveh, livew, wind_mph, slope_tan; + while (scanf("%d %lf %lf %lf %lf %lf %lf %lf", &fuel, &ones, &tens, + &hundreds, &liveh, &livew, &wind_mph, &slope_tan) == 8) + { + Mechanix m; + memset(&m, 0, sizeof(m)); + m.m_ones = ones; + m.m_tens = tens; + m.m_hundreds = hundreds; + m.m_liveh = liveh; + m.m_livew = livew; + m.m_windspd = wind_mph; + + // fsxwburn4.cpp:279-285, BurnThread::SurfaceFire + m.fros = m.spreadrate(slope_tan, m.m_windspd, fuel); + double phiw = m.phiw, phis = m.phis, phiew = m.phiew; + double head = 0, lb = 1, flank = 0, back = 0; + if (m.fros > 0.0) + { + m.VecSurf(slope_tan, 0.0); + head = m.vecros; + m.ellipse(m.vecros, m.ivecspeed); + lb = m.lb_ratio; + // fsxwmech.cpp:509-512, SpreadCorrect at Normal = pi/2 and pi + flank = m.flank; + back = m.head - m.back; + } + printf("%.12g %.12g %.12g %.12g %.12g %.12g %.12g %.12g\n", + m.fros, phiw, phis, phiew, head, lb, flank, back); + } + return 0; +} diff --git a/tests/python/test_farsite_reference.py b/tests/python/test_farsite_reference.py new file mode 100644 index 0000000..6589729 --- /dev/null +++ b/tests/python/test_farsite_reference.py @@ -0,0 +1,285 @@ +#!/usr/bin/env python3 +"""Compare ForeFire's Farsite model with the FARSITE reference code. + +The reference is ``farsite_reference.cpp`` next to this file: FARSITE's own +spread rate, wind/slope vector and ellipse, copied from mbedward/farsite with +only the plumbing changed. It is compiled with the system C++ compiler when +the test starts. + +The ForeFire side runs real simulations. Each case lays a circle of 72 nodes +on uniform fuel, uniform wind blowing north and a plane rising to the north, +and reads back what ``Farsite::getSpeed`` was called with and returned, from +the CSV written by ``FarsiteLoggerCSVPath``. The node facing north is the +head, the one facing south the back, the ones facing east or west the flanks. + +Two comparisons, both against the reference, never against ForeFire itself: + +model The reference is fed the wind and slope ForeFire handed to the + model at the head node. Any gap is in FarsitePropagationModel.cpp + or in the fuel table and moistures ForeFire feeds it. +pipeline The reference is fed the case as FARSITE sees it: the midflame wind + and the tangent of the terrain slope. Any gap beyond the model one + is in how ForeFire derives the normal wind and slope it hands over. + +The test fails when the head rate of spread differs beyond REL_TOL in the +model comparison, or in the pipeline comparison on flat ground. On a slope the +pipeline gap is known and only reported, unless ``--strict``. Flank, back and +length-to-breadth are only reported: ForeFire has no ellipse, so a node not +facing the wind spreads at the no-wind rate. + +Every case runs in its own interpreter, because a second ``ForeFire()`` in a +process inherits the first one's parameters (see test_moisture_invariants.py). + +Usage +----- + python3 tests/python/test_farsite_reference.py # summary + python3 tests/python/test_farsite_reference.py -v # every case + python3 tests/python/test_farsite_reference.py --strict # sloped pipeline too +""" + +import argparse +import json +import math +import os +import shutil +import subprocess +import sys +import tempfile +from concurrent.futures import ThreadPoolExecutor + +HERE = os.path.dirname(os.path.abspath(__file__)) +REFERENCE_SRC = os.path.join(HERE, "farsite_reference.cpp") + +FUELS = range(1, 14) +WINDS = (0.0, 2.0, 5.0, 10.0) # m/s at midflame +SLOPES = (0.0, 10.0, 30.0) # degrees +MOISTURES = { # fractions of dry weight: ones, tens, hundreds, liveh, livew + "dry": (0.03, 0.04, 0.05, 0.60, 0.90), + "moderate": (0.06, 0.07, 0.08, 0.90, 1.20), + "moist": (0.09, 0.10, 0.11, 1.20, 1.50), +} + +MPH = 0.44704 # m/s +# The ForeFire log prints 6 significant digits, so the comparison cannot be +# tighter than about 1e-5; anything a porting bug produces is far larger. +REL_TOL = 1e-4 + +DOMAIN = 400.0 # m, square +GRID = 80 # cells per side +RADIUS = 40.0 # m, the initial circle +NODES = 72 # node 0 faces north, NODES/4 east, NODES/2 south + +_SENTINEL = "@@PROBE@@" + + +# -------------------------------------------------------------------------- +# Reference +# -------------------------------------------------------------------------- + +def build_reference(workdir): + exe = os.path.join(workdir, "farsite_reference") + cxx = os.environ.get("CXX") or shutil.which("c++") or shutil.which("g++") + if cxx is None: + raise RuntimeError("no C++ compiler found to build farsite_reference.cpp") + subprocess.run([cxx, "-O2", "-o", exe, REFERENCE_SRC], check=True) + return exe + + +def run_reference(exe, rows): + """rows: (fuel, moistures, wind m/s, slope tangent). Returns dicts, m/s.""" + lines = [] + for fuel, moist, wind, slope in rows: + args = (fuel,) + tuple(moist) + (wind / MPH, slope) + lines.append(" ".join(repr(float(a)) if i else str(a) + for i, a in enumerate(args))) + out = subprocess.run([exe], input="\n".join(lines) + "\n", + capture_output=True, text=True, check=True).stdout + keys = ("rateo", "phiw", "phis", "phiew", "head", "lb", "flank", "back") + results = [] + for line in out.splitlines(): + r = dict(zip(keys, map(float, line.split()))) + for k in ("rateo", "head", "flank", "back"): + r[k] /= 60.0 # m/min -> m/s + results.append(r) + if len(results) != len(rows): + raise RuntimeError("reference returned %d rows for %d cases" + % (len(results), len(rows))) + return results + + +# -------------------------------------------------------------------------- +# ForeFire +# -------------------------------------------------------------------------- + +def probe(case): + proc = subprocess.run( + [sys.executable, os.path.abspath(__file__), "--probe", json.dumps(case)], + capture_output=True, text=True, timeout=300, + ) + for line in proc.stdout.splitlines(): + if line.startswith(_SENTINEL): + return json.loads(line[len(_SENTINEL):]) + raise RuntimeError( + f"probe {case} produced no result (exit {proc.returncode})\n" + f"stdout: {proc.stdout[-2000:]}\nstderr: {proc.stderr[-2000:]}" + ) + + +def _probe_in_process(case): + """One case in a fresh interpreter: every getSpeed call at ignition.""" + import numpy as np + import pyforefire as forefire + + ones, tens, hundreds, liveh, livew = MOISTURES[case["moisture"]] + wind, slope_deg = case["wind"], case["slope"] + + with tempfile.TemporaryDirectory() as tmp: + log = os.path.join(tmp, "farsite.csv") + ff = forefire.ForeFire() + ff["propagationModel"] = "Farsite" + ff["fuelsTable"] = "STDfarsiteFuelsTable" + ff["FarsiteLoggerCSVPath"] = log + ff["moistures.ones"] = ones + ff["moistures.tens"] = tens + ff["moistures.hundreds"] = hundreds + ff["moistures.liveh"] = liveh + ff["moistures.livew"] = livew + ff["minSpeed"] = 0.0 + ff["relax"] = 1.0 + ff["spatialIncrement"] = 1.0 + ff["perimeterResolution"] = 10.0 + ff["minimalPropagativeFrontDepth"] = 10.0 + ff["bmapLayer"] = 1 + ff["SWx"] = 0.0 + ff["SWy"] = 0.0 + ff["Lx"] = DOMAIN + ff["Ly"] = DOMAIN + + ff.execute(f"FireDomain[sw=(0.,0.,0.);ne=({DOMAIN},{DOMAIN},0);t=0]") + ff.addLayer("propagation", "Farsite", "propagationModel") + ff.addIndexLayer("table", "fuel", 0.0, 0.0, 0, DOMAIN, DOMAIN, 0, + np.full((1, 1, GRID, GRID), case["fuel"], dtype=np.int32)) + u = np.zeros((1, 2, GRID, GRID)) + v = np.zeros((1, 2, GRID, GRID)) + u[0, 0] = 1.0 + v[0, 1] = 1.0 + ff.addScalarLayer("windScalDir", "windU", 0.0, 0.0, 0, DOMAIN, DOMAIN, 0, u) + ff.addScalarLayer("windScalDir", "windV", 0.0, 0.0, 0, DOMAIN, DOMAIN, 0, v) + # A plane rising to the north, sampled at cell centres. + y = (np.arange(GRID) + 0.5) * DOMAIN / GRID + altitude = np.zeros((1, 1, GRID, GRID)) + altitude[0, 0] = (y * math.tan(math.radians(slope_deg)))[:, None] + ff.addScalarLayer("data", "altitude", 0.0, 0.0, 0, DOMAIN, DOMAIN, 0, + altitude) + ff.execute(f"trigger[wind;loc=(0.,0.,0.);vel=(0.,{wind},0.);t=0]") + + ff.execute("FireFront[id=2;domain=0;t=0]") + for k in range(NODES): + a = math.pi / 2 - 2 * math.pi * k / NODES # clockwise from north + x = DOMAIN / 2 + RADIUS * math.cos(a) + yy = DOMAIN / 2 + RADIUS * math.sin(a) + ff.execute( + f"FireNode[domain=0;id={4 + 2 * k};fdepth=5;kappa=0;" + f"loc=({x},{yy},0);vel=({0.01 * math.cos(a)}," + f"{0.01 * math.sin(a)},0);t=0;state=init;frontId=2]") + ff.execute("goTo[t=0.01]") + del ff # log lines end in endl, so the file is already flushed + + with open(log) as f: + header = f.readline().strip().split(";") + rows = [dict(zip(header, map(float, line.split(";")))) + for line in f if line.strip()] + keep = ("ROS", "rateo", "phiew", "slope", "normalWind") + return [{k: r[k] for k in keep} for r in rows] + + +# -------------------------------------------------------------------------- +# Comparison +# -------------------------------------------------------------------------- + +def pick(rows, which): + """The logged call at the head, back or flank of the circle.""" + def facing(r): # +1 north (head), -1 south (back), 0 east/west (flank) + return r["normalWind"] + r["slope"] + if which == "head": + return max(rows, key=facing) + if which == "back": + return min(rows, key=facing) + return min(rows, key=lambda r: abs(r["normalWind"]) + abs(r["slope"])) + + +def rel(a, b): + if a == b: + return 0.0 + return abs(a - b) / max(abs(a), abs(b)) + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__.split("\n")[0]) + parser.add_argument("-v", "--verbose", action="store_true", + help="print every case, not only the mismatches") + parser.add_argument("--strict", action="store_true", + help="also fail on the pipeline gap on a slope") + parser.add_argument("--probe", help=argparse.SUPPRESS) + args = parser.parse_args(argv) + + if args.probe: + print(_SENTINEL + json.dumps(_probe_in_process(json.loads(args.probe)))) + return 0 + + cases = [{"fuel": f, "moisture": m, "wind": w, "slope": s} + for f in FUELS for m in MOISTURES for w in WINDS for s in SLOPES] + with ThreadPoolExecutor(max_workers=min(8, os.cpu_count() or 2)) as pool: + logs = list(pool.map(probe, cases)) + + heads = [pick(rows, "head") for rows in logs] + with tempfile.TemporaryDirectory() as tmp: + exe = build_reference(tmp) + model = run_reference(exe, [ + (c["fuel"], MOISTURES[c["moisture"]], h["normalWind"], h["slope"]) + for c, h in zip(cases, heads)]) + pipeline = run_reference(exe, [ + (c["fuel"], MOISTURES[c["moisture"]], c["wind"], + math.tan(math.radians(c["slope"]))) + for c in cases]) + + failures = [] + columns = ("fuel moist U slope | R0 ff/ref phiew ff/ref " + "| head ff model pipe | LB flank ff/ref back ff/ref") + print(columns) + for c, rows, h, mo, pi in zip(cases, logs, heads, model, pipeline): + flank, back = pick(rows, "flank"), pick(rows, "back") + errs = { + "R0": rel(h["rateo"], mo["rateo"]), + "phiew": rel(h["phiew"], mo["phiew"]), + "head(model)": rel(h["ROS"], mo["head"]), + "head(pipeline)": rel(h["ROS"], pi["head"]), + } + bad = [k for k in ("R0", "phiew", "head(model)") if errs[k] > REL_TOL] + if errs["head(pipeline)"] > REL_TOL and (args.strict or c["slope"] == 0): + bad.append("head(pipeline)") + if bad: + failures.append((c, bad)) + if args.verbose or bad: + print(f"FM{c['fuel']:<3} {c['moisture']:<8} {c['wind']:4.0f} " + f"{c['slope']:4.0f} | {h['rateo']:.4g}/{mo['rateo']:.4g} " + f"{h['phiew']:.4g}/{mo['phiew']:.4g} " + f"| {h['ROS']:.4g} {mo['head']:.4g} {pi['head']:.4g} " + f"| {pi['lb']:.2f} {flank['ROS']:.4g}/{pi['flank']:.4g} " + f"{back['ROS']:.4g}/{pi['back']:.4g}" + + (" <- " + ", ".join(bad) if bad else "")) + + sloped = [(c, h, pi) for c, h, pi in zip(cases, heads, pipeline) + if c["slope"] > 0 and pi["head"] > 0] + print(f"\n{len(cases)} cases, {len(failures)} failing (rel. tol. {REL_TOL})") + for deg in sorted({c["slope"] for c, _, _ in sloped}): + ratios = [h["ROS"] / pi["head"] for c, h, pi in sloped + if c["slope"] == deg] + print(f"slope {deg:.0f} deg, ForeFire head / FARSITE head: " + f"{min(ratios):.3f} to {max(ratios):.3f}" + + ("" if args.strict else " (reported, not asserted)")) + return 1 if failures else 0 + + +if __name__ == "__main__": + sys.exit(main())