merge with ff farsite

.gitignore

@@ -1,25 +1,70 @@
-#Ignoring files and folders by github
+slurm*
+output/
+#dir
+src/__pycache__/
+output/*
+old/
+*.png
+*.npy
+#gvim 
+*.bak
+*.swp
+*.swo
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
 
-#build directories
-bin
-build
-lib
-sconsbuild
-forefire2
+# C extensions
+*.so
 
-#python
-__pycache__
+# Distribution / packaging
+.Python
+env/
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
 
-# compilation files
-*.dblite
-*.o
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
 
-# editors
-.vscode
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
 
-#mac files
-**/.DS_Store
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*,cover
 
-# jupyter notebooks checkpoints
-.ipynb_checkpoints
-CMakeLists.txt
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/

src/Andrews.cpp

@@ -0,0 +1,214 @@
+/*
+
+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 <math.h>
+using namespace std;
+namespace libforefire {
+
+class Andrews: 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 Md;
+	size_t normalWind;
+	size_t Sigmad;
+	size_t sd;
+	size_t Rhod;
+	size_t e;
+	size_t Mchi;
+	size_t DeltaH;
+
+	/*! coefficients needed by the model */
+
+	/*! local variables */
+
+	/*! result of the model */
+	double getSpeed(double*);
+
+public:
+	Andrews(const int& = 0, DataBroker* db=0);
+	virtual ~Andrews();
+
+	string getName();
+
+};
+
+PropagationModel* getAndrewsModel(const int& = 0, DataBroker* db=0);
+
+
+/* name of the model */
+const string Andrews::name = "Andrews";
+
+/* instantiation */
+PropagationModel* getAndrewsModel(const int & mindex, DataBroker* db) {
+	return new Andrews(mindex, db);
+}
+
+/* registration */
+int Andrews::isInitialized =
+		FireDomain::registerPropagationModelInstantiator(name, getAndrewsModel );
+
+/* constructor */
+Andrews::Andrews(const int & mindex, DataBroker* db)
+: PropagationModel(mindex, db) {
+	/* defining the properties needed for the model */
+	windReductionFactor = params->getDouble("windReductionFactor");
+
+	slope = registerProperty("slope");
+	Md = registerProperty("moisture");
+
+    normalWind = registerProperty("normalWind");
+	Sigmad = registerProperty("fuel.Sigmad");
+	sd = registerProperty("fuel.sd");
+	Rhod = registerProperty("fuel.Rhod");
+	e = registerProperty("fuel.e");
+	Mchi = registerProperty("fuel.Mchi");
+    DeltaH = registerProperty("fuel.DeltaH");
+
+	/* 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) */
+Andrews::~Andrews() {
+}
+
+/* accessor to the name of the model */
+string Andrews::getName(){
+	return name;
+}
+
+/* *********************************************** */
+/* Model for the propagation velovity of the front */
+/* *********************************************** */
+
+double Andrews::getSpeed(double* valueOf){
+
+	double lRhod = valueOf[Rhod] * 0.06; // conversion kg/m^3 -> lb/ft^3
+	double lMd  = valueOf[Md];
+	double lMchi  = valueOf[Mchi];
+	double lsd  = valueOf[sd] / 3.2808399; // conversion 1/m -> 1/ft
+	double le   = valueOf[e] * 3.2808399; // conversion m -> ft
+	if (le==0) return 0;
+	double lSigmad = valueOf[Sigmad] * 0.2048; // conversion kg/m^2 -> lb/ft^2
+	double lDeltaH = valueOf[DeltaH] / 2326.0;// conversion J/kg -> BTU/lb
+	double normal_wind  = valueOf[normalWind] * 196.850394 ; //conversion m/s -> ft/min
+	double localngle =  valueOf[slope];
+	//	if (normal_wind  > 0) cout<<"wind is"<<valueOf[normalWind]<<endl;
+
+/*	lRhod = 30.0;
+	 lMd =0.1;
+	  lsd =1523.99999768352;
+	  le =1.0;
+	  lSigmad =1044.27171; */
+
+	normal_wind *= windReductionFactor; // factor in the data seen in 2013
+
+	if (normal_wind < 0) normal_wind = 0;
+
+
+	double tanangle = localngle;
+	if (tanangle<0) tanangle=0;
+
+
+	//double Mchi = 0.3; // Moisture of extinction
+
+	double Etas = 1; // no mineral damping
+
+	double Wn = lSigmad;
+
+	double Mratio = lMd / lMchi;
+
+	double Etam = 1  + Mratio * (-2.59 + Mratio * (5.11 - 3.52 * Mratio));
+
+	double A = 1 / (4.774 * pow(lsd, 0.1) - 7.27);
+
+	double lRhobulk = Wn / le; // Dead bulk density = Dead fuel load / Fuel height
+
+	double Beta = lRhobulk / lRhod;  // Packing ratio = Bulk density / Particle density
+
+	double Betaop = 3.348 * pow(lsd, -0.8189);
+
+	double RprimeMax = pow(lsd, 1.5) * (1 / (495 + 0.0594 * pow(lsd, 1.5)));
+
+	double Rprime = RprimeMax * pow((Beta/Betaop), A) *  exp(A*(1-(Beta/Betaop))) ;
+
+	double chi = pow(192 + 0.259*lsd, -1) * exp((0.792 + 0.681*pow(lsd, 0.5)) * (Beta+0.1));
+
+	double epsilon  = exp(-138 / lsd);
+
+	double Qig = 250 + 1116 * lMd; // "1,116 in Andrews" != 1.116
+
+	double C = 7.47 * exp(-0.133*pow(lsd, 0.55));
+
+	double B = 0.02526*pow(lsd, 0.54);
+
+	double E = 0.715* exp(-3.59*(1.E-4 * lsd));
+
+	double Ir = Rprime*Wn*lDeltaH*Etam*Etas;
+	double Uf = 0.9*Ir;
+
+	// wind limit 2013 10.1071/WF12122 andrews/cruz/rothermel
+	if(windReductionFactor < 1.0){
+		Uf = 96.81*pow(Ir, 1./3);
+	}
+
+
+	if (normal_wind>Uf) {
+		normal_wind = Uf;
+	}
+
+	double phiV = C * pow((Beta/Betaop), -E) * pow(normal_wind,B) ;
+
+	double phiP = 5.275 * pow(Beta, -0.3) * pow(tanangle, 2);
+
+	double R0 = (Ir * chi) / (lRhobulk * epsilon * Qig);
+
+	double R = R0 * (1 + phiV + phiP);
+
+	if(R < R0)  R = R0;
+
+	if(R > 0.0) {
+		//cout << "Andrews =" << R * 0.00508 << endl; 
+		return R * 0.005588 ; // ft/min -> m/s
+	}else{
+		//cout << " Rhod "<< lRhod << " lMd "<< lMd << " lsd "<< lsd << " le "<< le << " lSigmad "<< lSigmad <<endl;
+		//cout << " R "<< R << " R0 "<< R0 << " phiv " << phiV <<" phiP" << phiP <<endl;
+	}
+	return 0;
+}
+
+} /* namespace libforefire */

src/CLibForeFire.cpp

@@ -237,17 +237,15 @@ void FFPutDoubleArray(const char* mname, double* x,
 	string tmpname(mname);
 	// searching for concerned layer
 
-	//cout<<session->fd->getDomainID()<<" is putting "<<tmpname<<endl;
+    //cout<<session->fd->getDomainID()<<" is putting "<<tmpname<<endl;
 	DataLayer<double>* myLayer = session->fd->getDataLayer(tmpname);
-
 	if ( myLayer ){
 		FFArray<double>* myMatrix;
 		// getting the pointer
 		myLayer->getMatrix(&myMatrix, executor.getTime());
 
 		myMatrix->copyDataToFortran(x);
 	} else {
-		cout<<"Error trying to put data from unknown layer "<<tmpname<<endl;
 	}
 }
 

src/DataBroker.cpp

@@ -278,8 +278,6 @@ void DataBroker::registerLayer(string name, DataLayer<double>* layer) {
 
 
 	/* inserting the layer into the map of layers */
-
-
 	ilayer = layersMap.find(name);
 	if (ilayer != layersMap.end()) {
 
@@ -292,6 +290,7 @@ void DataBroker::registerLayer(string name, DataLayer<double>* layer) {
 		}
 	}
 
+
 	layersMap.insert(make_pair(name, layer));
 	layers.push_back(layer);
 
@@ -318,7 +317,6 @@ void DataBroker::registerLayer(string name, DataLayer<double>* layer) {
 		temperatureLayer = layer;
 	}
 	if (name.find("windU") != string::npos){
-
 			windULayer = layer;
 	}
 	if (name.find("windV") != string::npos){
@@ -912,7 +910,6 @@ void DataBroker::loadFromNCFile(string filename) {
 		NcVarAtt layerTypeNC;
 		string layerType;
 		for (auto it = allVariables.begin(); it != allVariables.end(); ++it){
-
 				try {
                     layerTypeNC = it->second.getAtt("type");
                     if (layerTypeNC.isNull()) {
@@ -936,8 +933,6 @@ void DataBroker::loadFromNCFile(string filename) {
 
 				string varName = it->first;
 
-
-
 				if (varName == "wind") {
 					XYZTDataLayer<double> * wul = constructXYZTLayerFromFile(&dataFile, "wind",0);
 					registerLayer("windU", wul);
@@ -1073,7 +1068,7 @@ FFPoint DataBroker::getNetCDFSWCorner(NcVar* var) {
 			params->getDouble("SHIFT_ALL_DATA_ORDINATES_BY"));
 	return shiftPos + relSWC;
 }
-double DataBroker::getNetCDFTimeOrigin(NcVar* var) {
+double DataBroker::getNetCDFTimeSpan(NcVar* var) {
 	double Lt =0;
  	NcVarAtt att = var->getAtt("Lt");
 	if(!att.isNull()) att.getValues(&Lt);
@@ -1091,7 +1086,7 @@ FFPoint DataBroker::getNetCDFSpatialSpan(NcVar* var) {
 	if(!att.isNull()) att.getValues(&Lz);
 	return FFPoint(Lx, Ly, Lz);
 }
-double DataBroker::getNetCDFTimeSpan(NcVar* var) {
+double DataBroker::getNetCDFTimeOrigin(NcVar* var) {
 	double timeOrigin =0;
  	NcVarAtt att = var->getAtt("t0");
 	if(!att.isNull()) att.getValues(&timeOrigin);
@@ -1143,9 +1138,9 @@ XYZTDataLayer<double>* DataBroker::constructXYZTLayerFromFile(	NcFile* NcdataFil
 
 	if (isRelevantData(SWCorner, spatialExtent)) {
 //
-		//cout << "Variable " << property <<" " <<data[0]<< " have "<<timeOrigin<<" "<< values.getDimCount()<< " " << Lt<< " " << nx<< " " << ny<< " " << nz<< " " << nt<< " " << " SW" <<SWCorner.print()<< " ext:"<<spatialExtent.print()<< endl;
+		cout << "Variable " << property <<" " <<data[0]<< " have "<<values.getDimCount()<<" "<< timeOrigin<< " " << Lt<< " " << nx<< " " << ny<< " " << nz<< " " << nt<< " " << " SW" <<SWCorner.print()<< " ext:"<<spatialExtent.print()<< endl;
 
-		XYZTDataLayer<double>* newlayer = new XYZTDataLayer<double>(property,	SWCorner, timeOrigin, spatialExtent, Lt, nx, ny, nz, nt, data);
+		XYZTDataLayer<double>* newlayer = new XYZTDataLayer<double>(property, SWCorner, timeOrigin, spatialExtent, Lt, nx, ny, nz, nt, data);
 
 		delete[] data;
 		return newlayer;