From d5d761032967bb89493e9f49b839966b6d04d567 Mon Sep 17 00:00:00 2001 From: FILIPPI Jean-Baptiste Date: Thu, 12 Sep 2024 16:39:09 +0200 Subject: [PATCH] modified learn from run --- src/RothermelAndrews2018.cpp | 93 +++++++++++++++++------------------- 1 file changed, 44 insertions(+), 49 deletions(-) diff --git a/src/RothermelAndrews2018.cpp b/src/RothermelAndrews2018.cpp index 110f4e8..c4ff250 100644 --- a/src/RothermelAndrews2018.cpp +++ b/src/RothermelAndrews2018.cpp @@ -45,17 +45,17 @@ class RothermelAndrews2018: public PropagationModel { static int isInitialized; double windReductionFactor; /*! properties needed by the model */ - 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 st_; - size_t se_; - size_t me_; + 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 st; + size_t se; + size_t me; /*! result of the model */ double getSpeed(double*); @@ -89,18 +89,18 @@ RothermelAndrews2018::RothermelAndrews2018(const int & mindex, DataBroker* db) /* defining the properties needed for the model */ windReductionFactor = params->getDouble("windReductionFactor"); - 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"); - pp_ = registerProperty("fuel.fuelDens_lbft3"); - h_ = registerProperty("fuel.H_BTUlb"); - me_ = registerProperty("fuel.Dme_pc"); - st_ = registerProperty("fuel.totMineral_r"); - se_ = registerProperty("fuel.effectMineral_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"); + pp = registerProperty("fuel.fuelDens_lbft3"); + h = registerProperty("fuel.H_BTUlb"); + me = registerProperty("fuel.Dme_pc"); + 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]; @@ -124,7 +124,7 @@ string RothermelAndrews2018::getName(){ /* Model for the propagation velovity of the front */ /* *********************************************** */ -double RothermelAndrews2018::getSpeed(double* valueOf){ +double RothermelAndrews2018::getSpeed(double* Z){ // Constants double msToftmin = 196.85039; @@ -133,40 +133,35 @@ double RothermelAndrews2018::getSpeed(double* valueOf){ 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_] * tacTokgm2 / lbft2Tokgm2; // convert US short tons per acre to pounds per square foot - double h = valueOf[h_]; - double st = valueOf[st_]; - double se = valueOf[se_]; - - if (wv < 0) wv = 0; - - if(wo > 0){ - double Beta_op = 3.348 * pow(fpsa, -0.8189); // Optimum packing ratio - double ODBD = wo / fd; // Ovendry bulk density + double wvC = msToftmin * Z[wv]; // convert m/s to feat/min + double meC = Z[me] * pcTor; // convert percentage to ratio + double woC = Z[wo] * tacTokgm2 / lbft2Tokgm2; // convert US short tons per acre to pounds per square foot + + if (wvC < 0) wvC = 0; + + if(woC > 0){ + double Beta_op = 3.348 * pow(Z[fpsa], -0.8189); // Optimum packing ratio + double ODBD = woC / Z[fd]; // Ovendry bulk density double Beta = ODBD / pp; // Packing ratio - double WN = wo / (1 + st); // Net fuel loading + double WN = woC / (1 + Z[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 / 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 NM = 1. - 2.59 * (Z[mf] / meC) + 5.11 * pow(mf / meC, 2.) - 3.52 * pow(mf / meC,3.); // Moisture damping coeff. + double NS = 0.174 * pow(Z[se], -0.19); // Mineral damping coefficient + double RI = T * WN * Z[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 Coefficien t double B = 0.02526 * pow(fpsa, 0.54); double C = 7.47 * exp(-0.1333 * pow(fpsa, 0.55)); 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); + if (wvC > (0.9 * RI)) wvC = 0.9 * RI; + double WC = (C * pow(wvC, B)) * pow((Beta / Beta_op), (-E)); + + double SC = 0; + if (Z[slope] >0){ // Rothermel only for upslope, assuming no slope if negative slope is encountered + 5.275*(pow(Beta, -0.3))*pow(Z[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)