From 1ea362e6c011066f3205c49742d768d3af76abbf Mon Sep 17 00:00:00 2001 From: Thoreau Romain Date: Thu, 5 Sep 2024 16:10:47 +0200 Subject: [PATCH 1/6] Add Rothermel-Andrews-2018 ROS model --- src/RothermelAndrews2018.cpp | 171 +++++++++++++++++++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100644 src/RothermelAndrews2018.cpp diff --git a/src/RothermelAndrews2018.cpp b/src/RothermelAndrews2018.cpp new file mode 100644 index 0000000..618573a --- /dev/null +++ b/src/RothermelAndrews2018.cpp @@ -0,0 +1,171 @@ +/* + +Copyright (C) 2012 ForeFire Team, SPE, Universit� de Corse. + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU Lesser General Public +License as published by the Free Software Foundation; either +version 2.1 of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public +License along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 US + + */ + + + +#include "PropagationModel.h" +#include "FireDomain.h" +#include +using namespace std; +namespace libforefire { + +class RothermelAndrews2018: public PropagationModel { + + /*! name the model */ + static const string name; + + /*! boolean for initialization */ + static int isInitialized; + double windReductionFactor; + /*! properties needed by the model */ + size_t slope; + size_t normalWind; + size_t Rhod; + size_t Md; + size_t sd; + size_t e; + size_t Sigmad; + size_t DeltaH; + + /*! result of the model */ + double getSpeed(double*); + +public: + RothermelAndrews2018(const int& = 0, DataBroker* db=0); + virtual ~RothermelAndrews2018(); + + string getName(); + +}; + +PropagationModel* getRothermelAndrews2018Model(const int& = 0, DataBroker* db=0); + + +/* name of the model */ +const string RothermelAndrews2018::name = "RothermelAndrews2018"; + +/* instantiation */ +PropagationModel* getRothermelAndrews2018Model(const int & mindex, DataBroker* db) { + return new RothermelAndrews2018(mindex, db); +} + +/* registration */ +int RothermelAndrews2018::isInitialized = + FireDomain::registerPropagationModelInstantiator(name, getRothermelAndrews2018Model ); + +/* constructor */ +RothermelAndrews2018::RothermelAndrews2018(const int & mindex, DataBroker* db) +: PropagationModel(mindex, db) { + /* defining the properties needed for the model */ + windReductionFactor = params->getDouble("windReductionFactor"); + + wv = registerProperty("normalWind_ms"); + slope = registerProperty("slope_deg"); + + wo = registerProperty("fl1h_tac"); + fd = registerProperty("fd_ft"); + fpsa = registerProperty("SAVcar_ftinv"); + mf = registerProperty("mdOnDry1h_r"); + pp = registerProperty("fuelDens_lbft3"); + h = registerProperty("H_BTUlb"); + + /* allocating the vector for the values of these properties */ + if ( numProperties > 0 ) properties = new double[numProperties]; + + /* registering the model in the data broker */ + dataBroker->registerPropagationModel(this); + + /* Definition of the coefficients */ +} + +/* destructor (shoudn't be modified) */ +RothermelAndrews2018::~RothermelAndrews2018() { +} + +/* accessor to the name of the model */ +string RothermelAndrews2018::getName(){ + return name; +} + +/* *********************************************** */ +/* Model for the propagation velovity of the front */ +/* *********************************************** */ + +double RothermelAndrews2018::getSpeed(double* valueOf){ + + double mois_ext = 0.3; + double st = 0.0555; + double se = 0.01; + double Pi = 3.14159265358; + double ms2ftmin = 196.85039; + double ftmin2ms = 0.00508; + + double slope = Pi * valueOf[slope] / 180; // convert degrees to radians + double wv = ms2ftmin * valueOf[wv]; // convert m/s to feat/min + double pp = valueOf[pp]; + double mf = valueOf[mf]; + double fpsa = valueOf[fpsa]; + double fd = valueOf[fd]; + double wo = valueOf[wo]; + double h = valueOf[h]; + + if (wv < 0) wv = 0; + + if(wv < 0){ + double tan_slope = tan(slope); + double Beta_op = 3.348 * pow(fpsa, -0.8189); // Optimum packing ratio + double ODBD = wo / fd; // Ovendry bulk density + double Beta = ODBD / pp; // Packing ratio + double WN = wo / (1 + st); // Net fuel loading + double A = 133.0 / pow(fpsa, 0.7913); // updated A + double T_max = pow(fpsa,1.5) * pow(495.0 + 0.0594 * pow(fpsa, 1.5),-1.0); // Maximum reaction velocity + double T = T_max * pow((Beta / Beta_op), A) * exp(A * (1 - Beta / Beta_op)); // Optimum reaction velocity + double NM = 1. - 2.59 * (mf / mois_ext) + 5.11 * pow(mf / mois_ext, 2.) - 3.52 * pow(mf / mois_ext,3.); // Moisture damping coeff. + double NS = 0.174 * pow(se, -0.19); // Mineral damping coefficient + double RI = T * WN * h * NM * NS; + double PFR = pow(192.0 + 0.2595 * fpsa, -1) * exp((0.792 + 0.681 * pow(fpsa, 0.5)) * (Beta + 0.1)); // Propogating flux ratio + // Wind Coefficient + double B = 0.02526 * pow(fpsa, 0.54); + double C = 7.47 * exp(-0.1333 * pow(fpsa, 0.55)); + double = 0.715 * exp(-3.59 * pow(10, -4) * fpsa); + if (wv > (0.9 * RI)) wv = 0.9 * RI; + double WC = (C * pow(wv, B)) * pow((Beta / Beta_op), (-E)); + double SC = 5.275*(pow(Beta, -0.3))*pow(tan_slope, 2); + // Heat sink + double EHN = exp(-138. / fpsa); // Effective Heating Number = f(surface are volume ratio) + double QIG = 250. + 1116. * mf; // Heat of preignition= f(moisture content) + // rate of spread (ft per minute) // RI = BTU/ft^2 + double numerator = (RI * PFR * (1 + WC + SC)); + double denominator = (ODBD * EHN * QIG); + double R = numerator / denominator; // WC and SC will be zero at slope = wind = 0 + double RT = 384.0/fpsa; + double HA = RI*RT; + // fireline intensity as described by Albini via USDA Forest Service RMRS-GTR-371. 2018 + double FI = (384.0/fpsa)*RI*(R); // Uses Reaction Intensity in BTU / ft/ min + + if(R <= 0.0) { + return 0; + }else{ + return R * ftmin2ms; + } + }else{ + return 0 + } +} /* namespace libforefire */ From dfbbfd702f9149ca77668f865dd792125c21bf0f Mon Sep 17 00:00:00 2001 From: Thoreau Romain Date: Thu, 5 Sep 2024 16:57:18 +0200 Subject: [PATCH 2/6] update RothermelAndrews2018 --- src/RothermelAndrews2018.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/RothermelAndrews2018.cpp b/src/RothermelAndrews2018.cpp index 618573a..4e73706 100644 --- a/src/RothermelAndrews2018.cpp +++ b/src/RothermelAndrews2018.cpp @@ -85,6 +85,9 @@ RothermelAndrews2018::RothermelAndrews2018(const int & mindex, DataBroker* db) mf = registerProperty("mdOnDry1h_r"); pp = registerProperty("fuelDens_lbft3"); h = registerProperty("H_BTUlb"); + mois_ext = registerProperty("mois_ext_r"); + st = registerProperty("st_r"); + se = registerProperty("se_r"); /* allocating the vector for the values of these properties */ if ( numProperties > 0 ) properties = new double[numProperties]; From 07328e8d1dd8795fc5a3847b19ed1f064ea6e935 Mon Sep 17 00:00:00 2001 From: Thoreau Romain Date: Thu, 5 Sep 2024 17:21:08 +0200 Subject: [PATCH 3/6] update RothermelAndrews2018.cpp --- src/RothermelAndrews2018.cpp | 62 ++++++++++++++++++------------------ 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/src/RothermelAndrews2018.cpp b/src/RothermelAndrews2018.cpp index 4e73706..b589dfc 100644 --- a/src/RothermelAndrews2018.cpp +++ b/src/RothermelAndrews2018.cpp @@ -35,14 +35,17 @@ class RothermelAndrews2018: public PropagationModel { static int isInitialized; double windReductionFactor; /*! properties needed by the model */ - size_t slope; - size_t normalWind; - size_t Rhod; - size_t Md; - size_t sd; - size_t e; - size_t Sigmad; - size_t DeltaH; + size_t wv_; + size_t slope_; + size_t wo_; + size_t fd_; + size_t fpsa_; + size_t mf_; + size_t pp_; + size_t h_; + size_t mois_ext; + size_t st; + size_t se; /*! result of the model */ double getSpeed(double*); @@ -76,15 +79,15 @@ RothermelAndrews2018::RothermelAndrews2018(const int & mindex, DataBroker* db) /* defining the properties needed for the model */ windReductionFactor = params->getDouble("windReductionFactor"); - wv = registerProperty("normalWind_ms"); - slope = registerProperty("slope_deg"); + wv_ = registerProperty("normalWind_ms"); + slope_ = registerProperty("slope_deg"); - wo = registerProperty("fl1h_tac"); - fd = registerProperty("fd_ft"); - fpsa = registerProperty("SAVcar_ftinv"); - mf = registerProperty("mdOnDry1h_r"); - pp = registerProperty("fuelDens_lbft3"); - h = registerProperty("H_BTUlb"); + wo_ = registerProperty("fl1h_tac"); + fd_ = registerProperty("fd_ft"); + fpsa_ = registerProperty("SAVcar_ftinv"); + mf_ = registerProperty("mdOnDry1h_r"); + pp_ = registerProperty("fuelDens_lbft3"); + h_ = registerProperty("H_BTUlb"); mois_ext = registerProperty("mois_ext_r"); st = registerProperty("st_r"); se = registerProperty("se_r"); @@ -120,14 +123,14 @@ double RothermelAndrews2018::getSpeed(double* valueOf){ double ms2ftmin = 196.85039; double ftmin2ms = 0.00508; - double slope = Pi * valueOf[slope] / 180; // convert degrees to radians - double wv = ms2ftmin * valueOf[wv]; // convert m/s to feat/min - double pp = valueOf[pp]; - double mf = valueOf[mf]; - double fpsa = valueOf[fpsa]; - double fd = valueOf[fd]; - double wo = valueOf[wo]; - double h = valueOf[h]; + double slope = Pi * valueOf[slope_] / 180; // convert degrees to radians + double wv = ms2ftmin * valueOf[wv_]; // convert m/s to feat/min + double pp = valueOf[pp_]; + double mf = valueOf[mf_]; + double fpsa = valueOf[fpsa_]; + double fd = valueOf[fd_]; + double wo = valueOf[wo_]; + double h = valueOf[h_]; if (wv < 0) wv = 0; @@ -144,10 +147,10 @@ double RothermelAndrews2018::getSpeed(double* valueOf){ double NS = 0.174 * pow(se, -0.19); // Mineral damping coefficient double RI = T * WN * h * NM * NS; double PFR = pow(192.0 + 0.2595 * fpsa, -1) * exp((0.792 + 0.681 * pow(fpsa, 0.5)) * (Beta + 0.1)); // Propogating flux ratio - // Wind Coefficient + // Wind Coefficien t double B = 0.02526 * pow(fpsa, 0.54); double C = 7.47 * exp(-0.1333 * pow(fpsa, 0.55)); - double = 0.715 * exp(-3.59 * pow(10, -4) * fpsa); + double E = 0.715 * exp(-3.59 * pow(10, -4) * fpsa); if (wv > (0.9 * RI)) wv = 0.9 * RI; double WC = (C * pow(wv, B)) * pow((Beta / Beta_op), (-E)); double SC = 5.275*(pow(Beta, -0.3))*pow(tan_slope, 2); @@ -158,10 +161,6 @@ double RothermelAndrews2018::getSpeed(double* valueOf){ double numerator = (RI * PFR * (1 + WC + SC)); double denominator = (ODBD * EHN * QIG); double R = numerator / denominator; // WC and SC will be zero at slope = wind = 0 - double RT = 384.0/fpsa; - double HA = RI*RT; - // fireline intensity as described by Albini via USDA Forest Service RMRS-GTR-371. 2018 - double FI = (384.0/fpsa)*RI*(R); // Uses Reaction Intensity in BTU / ft/ min if(R <= 0.0) { return 0; @@ -169,6 +168,7 @@ double RothermelAndrews2018::getSpeed(double* valueOf){ return R * ftmin2ms; } }else{ - return 0 + return 0; } } /* namespace libforefire */ +} \ No newline at end of file From 68e2eddd2fe1d95aaa231feae5b5b5952cc2d86a Mon Sep 17 00:00:00 2001 From: Thoreau Romain Date: Thu, 5 Sep 2024 20:25:09 +0200 Subject: [PATCH 4/6] correct angle --- src/RothermelAndrews2018.cpp | 68 ++++++++++++++++++++---------------- 1 file changed, 38 insertions(+), 30 deletions(-) diff --git a/src/RothermelAndrews2018.cpp b/src/RothermelAndrews2018.cpp index b589dfc..a7cea8b 100644 --- a/src/RothermelAndrews2018.cpp +++ b/src/RothermelAndrews2018.cpp @@ -43,9 +43,9 @@ class RothermelAndrews2018: public PropagationModel { size_t mf_; size_t pp_; size_t h_; - size_t mois_ext; size_t st; size_t se; + size_t me_; /*! result of the model */ double getSpeed(double*); @@ -79,18 +79,18 @@ RothermelAndrews2018::RothermelAndrews2018(const int & mindex, DataBroker* db) /* defining the properties needed for the model */ windReductionFactor = params->getDouble("windReductionFactor"); - wv_ = registerProperty("normalWind_ms"); - slope_ = registerProperty("slope_deg"); - - wo_ = registerProperty("fl1h_tac"); - fd_ = registerProperty("fd_ft"); - fpsa_ = registerProperty("SAVcar_ftinv"); - mf_ = registerProperty("mdOnDry1h_r"); - pp_ = registerProperty("fuelDens_lbft3"); - h_ = registerProperty("H_BTUlb"); - mois_ext = registerProperty("mois_ext_r"); - st = registerProperty("st_r"); - se = registerProperty("se_r"); + wv_ = registerProperty("normalWind"); + slope_ = registerProperty("slope"); + + wo_ = registerProperty("fuel.fl1h_tac"); + fd_ = registerProperty("fuel.fd_ft"); + fpsa_ = registerProperty("fuel.SAVcar_ftinv"); + // mf_ = registerProperty("fuel.mdOnDry1h_r"); // is an environmental property that do not depends on fuel? = 0.06 + // pp_ = registerProperty("fuel.fuelDens_lbft3"); // is an environmental property that do not depends on fuel? = 32 + // h_ = registerProperty("fuel.H_BTUlb"); // is an environmental property that do not depends on fuel? = 8000 + me_ = registerProperty("fuel.Dme_pc"); + // st = registerProperty("fuel.st_r"); + // se = registerProperty("fuel.se_r"); /* allocating the vector for the values of these properties */ if ( numProperties > 0 ) properties = new double[numProperties]; @@ -116,26 +116,34 @@ string RothermelAndrews2018::getName(){ double RothermelAndrews2018::getSpeed(double* valueOf){ - double mois_ext = 0.3; - double st = 0.0555; - double se = 0.01; - double Pi = 3.14159265358; - double ms2ftmin = 196.85039; - double ftmin2ms = 0.00508; - - double slope = Pi * valueOf[slope_] / 180; // convert degrees to radians - double wv = ms2ftmin * valueOf[wv_]; // convert m/s to feat/min - double pp = valueOf[pp_]; - double mf = valueOf[mf_]; + // Constants + double msToftmin = 196.85039; + double ftminToms = 0.00508; + double lbft2Tokgm2 = 4.88243; + double tacTokgm2 = 0.224; + double pcTor = 0.01; + + // Fuel properties + double tan_slope = valueOf[slope_]; + double wv = msToftmin * valueOf[wv_]; // convert m/s to feat/min + double me = valueOf[me_] * pcTor; // convert percentage to ratio + // double pp = valueOf[pp_]; + // double mf = valueOf[mf_]; double fpsa = valueOf[fpsa_]; double fd = valueOf[fd_]; - double wo = valueOf[wo_]; - double h = valueOf[h_]; + double wo = valueOf[wo_] * tacTokgm2 / lbft2Tokgm2; // convert US short tons per acre to pounds per square foot + // double h = valueOf[h_]; + + // Fuel properties that are not in the fuel table + double st = 0.0555; + double se = 0.01; + double mf = 0.06; + double pp = 32; + double h = 8000; if (wv < 0) wv = 0; - if(wv < 0){ - double tan_slope = tan(slope); + if(wo > 0){ double Beta_op = 3.348 * pow(fpsa, -0.8189); // Optimum packing ratio double ODBD = wo / fd; // Ovendry bulk density double Beta = ODBD / pp; // Packing ratio @@ -143,7 +151,7 @@ double RothermelAndrews2018::getSpeed(double* valueOf){ double A = 133.0 / pow(fpsa, 0.7913); // updated A double T_max = pow(fpsa,1.5) * pow(495.0 + 0.0594 * pow(fpsa, 1.5),-1.0); // Maximum reaction velocity double T = T_max * pow((Beta / Beta_op), A) * exp(A * (1 - Beta / Beta_op)); // Optimum reaction velocity - double NM = 1. - 2.59 * (mf / mois_ext) + 5.11 * pow(mf / mois_ext, 2.) - 3.52 * pow(mf / mois_ext,3.); // Moisture damping coeff. + double NM = 1. - 2.59 * (mf / me) + 5.11 * pow(mf / me, 2.) - 3.52 * pow(mf / me,3.); // Moisture damping coeff. double NS = 0.174 * pow(se, -0.19); // Mineral damping coefficient double RI = T * WN * h * NM * NS; double PFR = pow(192.0 + 0.2595 * fpsa, -1) * exp((0.792 + 0.681 * pow(fpsa, 0.5)) * (Beta + 0.1)); // Propogating flux ratio @@ -165,7 +173,7 @@ double RothermelAndrews2018::getSpeed(double* valueOf){ if(R <= 0.0) { return 0; }else{ - return R * ftmin2ms; + return R * ftminToms; } }else{ return 0; From cf7bdc07a54022b8cd0f34f99bc404b980ad4160 Mon Sep 17 00:00:00 2001 From: Thoreau Romain Date: Fri, 6 Sep 2024 11:37:28 +0200 Subject: [PATCH 5/6] update RothermelAndrews2018.cpp --- src/RothermelAndrews2018.cpp | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/src/RothermelAndrews2018.cpp b/src/RothermelAndrews2018.cpp index a7cea8b..b1f1738 100644 --- a/src/RothermelAndrews2018.cpp +++ b/src/RothermelAndrews2018.cpp @@ -43,8 +43,8 @@ class RothermelAndrews2018: public PropagationModel { size_t mf_; size_t pp_; size_t h_; - size_t st; - size_t se; + size_t st_; + size_t se_; size_t me_; /*! result of the model */ @@ -85,12 +85,12 @@ RothermelAndrews2018::RothermelAndrews2018(const int & mindex, DataBroker* db) wo_ = registerProperty("fuel.fl1h_tac"); fd_ = registerProperty("fuel.fd_ft"); fpsa_ = registerProperty("fuel.SAVcar_ftinv"); - // mf_ = registerProperty("fuel.mdOnDry1h_r"); // is an environmental property that do not depends on fuel? = 0.06 - // pp_ = registerProperty("fuel.fuelDens_lbft3"); // is an environmental property that do not depends on fuel? = 32 - // h_ = registerProperty("fuel.H_BTUlb"); // is an environmental property that do not depends on fuel? = 8000 + mf_ = registerProperty("fuel.mdOnDry1h_r"); + pp_ = registerProperty("fuel.fuelDens_lbft3"); + h_ = registerProperty("fuel.H_BTUlb"); me_ = registerProperty("fuel.Dme_pc"); - // st = registerProperty("fuel.st_r"); - // se = registerProperty("fuel.se_r"); + st_ = registerProperty("fuel.totMineral_r"); + se_ = registerProperty("fuel.effectMineral_r"); /* allocating the vector for the values of these properties */ if ( numProperties > 0 ) properties = new double[numProperties]; @@ -127,19 +127,14 @@ double RothermelAndrews2018::getSpeed(double* valueOf){ double tan_slope = valueOf[slope_]; double wv = msToftmin * valueOf[wv_]; // convert m/s to feat/min double me = valueOf[me_] * pcTor; // convert percentage to ratio - // double pp = valueOf[pp_]; - // double mf = valueOf[mf_]; + double pp = valueOf[pp_]; + double mf = valueOf[mf_]; double fpsa = valueOf[fpsa_]; double fd = valueOf[fd_]; double wo = valueOf[wo_] * tacTokgm2 / lbft2Tokgm2; // convert US short tons per acre to pounds per square foot - // double h = valueOf[h_]; - - // Fuel properties that are not in the fuel table - double st = 0.0555; - double se = 0.01; - double mf = 0.06; - double pp = 32; - double h = 8000; + double h = valueOf[h_]; + double st = valueOf[st_]; + double se = valueOf[se_]; if (wv < 0) wv = 0; From 744e7c76d70d2d991f93cb5badeea6182025277a Mon Sep 17 00:00:00 2001 From: Thoreau Romain Date: Sat, 7 Sep 2024 11:27:26 +0200 Subject: [PATCH 6/6] add Andrews2018 ref --- src/RothermelAndrews2018.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/RothermelAndrews2018.cpp b/src/RothermelAndrews2018.cpp index b1f1738..110f4e8 100644 --- a/src/RothermelAndrews2018.cpp +++ b/src/RothermelAndrews2018.cpp @@ -27,6 +27,16 @@ using namespace std; namespace libforefire { class RothermelAndrews2018: public PropagationModel { + + // ROS model as defined in Table 5 and Table 6 in Andrews, 2018. + // @book{Andrews_2018, + // title={The Rothermel surface fire spread model and associated developments: A comprehensive explanation}, + // url={http://dx.doi.org/10.2737/RMRS-GTR-371}, + // DOI={10.2737/rmrs-gtr-371}, + // institution={U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station}, + // author={Andrews, Patricia L.}, + // year={2018} + // } /*! name the model */ static const string name;