Commit 1f52e511 authored by jouke's avatar jouke
Browse files

Merging three classes for leaf area 

- removing code duplication
- fixing bug that causes a run with stress but no carbon model to have
potential leaf area, rather than stress reduced leaf area.
parent fbe96ca5

OpenSimRoot/src/modules/ShootModule/LeafArea.cpp

+110 −127
Original line number Diff line number Diff line
@@ -29,31 +29,73 @@ NOTE: The GPL.v3 license requires that all derivative work is distributed under 
#include <algorithm>

#endif



PotentialLeafArea::PotentialLeafArea(SimulaDynamic* const pSV):

	DerivativeBase(pSV),LASimulator(nullptr)

DerivativeBase * newInstantiationPotentialLeafArea(SimulaDynamic* const pSV){

	msg::warning("Input files specify deprecated potentialLeafGrowthRate. Using 'leafArea' instead", 2);

   return new LeafArea(pSV);

}



DerivativeBase * newInstantiationStressAdjustedPotentialLeafArea(SimulaDynamic* const pSV){

	msg::warning("Input files specify deprecated stressAdjustedPotentialLeafGrowthRate. Using 'leafArea' instead", 2);

   return new LeafArea(pSV);

}





LeafArea::LeafArea(SimulaDynamic* const pSD):

	DerivativeBase(pSD), SLASimulator(nullptr), c2lSimulator(nullptr), potentialLA(nullptr), rgrCoefficient(nullptr), stress(nullptr), CinDryWeight(nullptr), LASimulator(nullptr), SLAisLAR_(true)

{

	//check if unit given in input file agrees with this function

	pSD->checkUnit("cm2");

	//planting time

	plantingTime=pSD->getStartTime();

	if(pSD->getUnit()!="cm2") msg::error("LeafArea: Expecting unit cm2 but found "+pSD->getUnit().name);

	//plant type

	std::string plantType;

	PLANTTYPE(plantType,pSD)

	// shoot params

	SimulaBase* shootParam=GETSHOOTPARAMETERS(plantType);

	SimulaBase* shootLA=shootParam->existingChild("leafAreaExpansionRate", "cm2/day");

	if(!shootLA){

		//backup because of spelling error in older files t should be s

		shootLA=shootParam->getChild("leafAreaExpantionRate", "cm2/day");

	//see if carbon model is run

	std::string name=pSD->getName();

	//planting time

	pSD->getParent(3)->getChild("plantingTime")->get(plantingTime);

	//plant parameters

	SimulaBase *parameters(GETSHOOTPARAMETERS(plantType));



	//modus

	bool addStress=false;

	if(name.find("potential")!=std::string::npos){

		potentialLA=nullptr;

	}else if(name.find("stress")!=std::string::npos){

		potentialLA=pSD->existingSibling("potentialLeafArea");

		if(!potentialLA) msg::error("LeafArea: stressAdjustedPotentialLeafArea can not be computed if the potentialLeafArea is not defined");

		addStress=true;

	}else{

		c2lSimulator=pSD->existingSibling("carbonAllocation2Leafs","g");

		if(c2lSimulator){

			SLASimulator=parameters->existingChild("specificLeafArea");

			if(!SLASimulator) SLASimulator=parameters->getChild("leafAreaRatio");

			if(SLASimulator->getUnit()=="cm2/g"){

				SLAisLAR_=false;

			}else{

				if(SLASimulator->getUnit()!="g/cm2") msg::error("LeafArea: Unkown unit for "+SLASimulator->getPath()+". Use g/cm2 or cm2/g");

			}

			msg::warning("LeafArea: using carbon allocation to compute "+name, 2);

		}else{

			potentialLA=pSD->existingSibling("stressAdjustedPotentialLeafArea");

			if(!potentialLA){

				potentialLA=pSD->existingSibling("potentialLeafArea");

				addStress=true;

			}

		}

	}

	SimulaBase* tillerParam=GETSHOOTPARAMETERS(plantType)->existingChild("tillers");





	if(!potentialLA){

			// shoot params

			SimulaBase* shootLA=parameters->getOneOfThesePaths("leafAreaExpansionRate","leafAreaExpantionRate");

			shootLA->checkUnit("cm2/day");

			msg::warning("LeafArea: using leafAreaExpansionRate from parameter section to compute "+name,2);

			//plant parameters

			std::string parent=pSD->getParent()->getName();

			if ( parent == "shoot"){

				SimulaBase *p(ORIGIN->existingChild("tillerTemplate"));

				if(p && p->existingChild("potentialLeafArea")){

					//add up tillers;

			msg::warning("PotentialLeafArea::leafArea is sum of tillers");

					msg::warning("LeafArea::leafArea is sum of tillers", 2);

					LASimulator=nullptr;

				}else{

					LASimulator=shootLA;
@@ -61,126 +103,45 @@ PotentialLeafArea::PotentialLeafArea(SimulaDynamic* const pSV): 
			}else if ( parent == "mainShoot"){

				LASimulator = shootLA;

			}else {

				SimulaBase* tillerParam=parameters->existingChild("tillers");

				if (tillerParam){

					LASimulator=tillerParam->existingChild("leafAreaExpansionRatePerTiller", "cm2/day");

					if (!LASimulator) {

						//backup because of typos in older files

						LASimulator=tillerParam->existingChild("leafAreaExpantionRatePerTiller", "cm2/day");

						msg::warning("LeafArea: Using deprecated 'leafAreaExpantionRatePerTiller', consider renaming it into 'leafAreaExpansionRatePerTiller'.",2);

					}

				}

				if (!LASimulator) {

					LASimulator=shootLA;

			msg::warning("PotentialLeafArea: leafAreaExpantionRatePerTiller not found, using for tillers the leafAreaExpantionRate");

		}

	}

}

void PotentialLeafArea::calculate(const Time &t, double& r){

	//get rate from parameter set

	if(LASimulator){

		LASimulator->get(t-plantingTime,r);

	}else{

		SimulaBase::List l;

		pSD->getSibling("tillers")->getAllChildren(l,t);

		pSD->getSibling("mainShoot")->getChild("potentialLeafArea")->getRate(t,r);

		for(auto &i:l){

			double la;

			i->getChild("potentialLeafArea")->getRate(t,la);

			r+=la;

		}

					msg::warning("PotentialLeafArea: leafAreaExpansionRatePerTiller not found, using for tillers the leafAreaExpansionRate",2);

				}

			}

std::string PotentialLeafArea::getName()const{

	return "potentialLeafGrowthRate";

}



DerivativeBase * newInstantiationPotentialLeafArea(SimulaDynamic* const pSV){

   return new PotentialLeafArea(pSV);

	}



//if we are simulating stress - switch to relative leafAreaExpantionRate as soon as stress factor has reduced growth

StressAdjustedPotentialLeafArea::StressAdjustedPotentialLeafArea(SimulaDynamic* const pSV):

	DerivativeBase(pSV)

{

	//simulators

	if(addStress){

		rgrCoefficient=pSD->existingSibling("leafAreaReductionCoefficient","cm2/cm2");

	potential=pSD->getSibling("potentialLeafArea",pSD->getUnit());

	//stress factor

		stress=pSD->getParent(3)->existingChild("stressFactor:impactOn:leafAreaExpansionRate");

		if (!stress) stress=pSD->getParent(3)->existingChild("stressFactor:impactOn:leafAreaExpantionRate");

		if(!stress) {

			stress=pSD->getParent(3)->getChild("stressFactor");

		msg::warning("StressAdjustedPotentialLeafArea:: \"stressFactor:impactOn:leafAreaExpantionRate\" or \"stressFactor:impactOn:leafAreaExpansionRate\" not found, using raw stress factor values");

	}

}

void StressAdjustedPotentialLeafArea::calculate(const Time &t, double& r){

	//get rate from parameter set

	potential->getRate(t,r);

	//relative growth rate for potential growth

	double c;

	if(rgrCoefficient){

		rgrCoefficient->get(t,c);

		r*=c;

	}

	//stress

	stress->get(t,c);

	if(c<1){

		r*=c;

	}

	if(r<0) {

		msg::warning("StressAdjustedPotentialLeafArea: fixing negative growth");

		r=0;

	}

			msg::warning("LeafArea: \"stressFactor:impactOn:leafAreaExpansionRate\" not found, using raw stress factor values");

		}

std::string StressAdjustedPotentialLeafArea::getName()const{

	return "stressAdjustedPotentialLeafGrowthRate";

}





DerivativeBase * newInstantiationStressAdjustedPotentialLeafArea(SimulaDynamic* const pSV){

   return new StressAdjustedPotentialLeafArea(pSV);

	}





LeafArea::LeafArea(SimulaDynamic* const pSD):

	DerivativeBase(pSD), SLASimulator(nullptr), c2lSimulator(nullptr), CinDryWeight(nullptr), SLAisLAR_(true)

{

	//check if unit given in input file agrees with this function

	if(pSD->getUnit()!="cm2") msg::error("LeafArea: Expecting unit cm2 but found "+pSD->getUnit().name);

	//plant type

	std::string plantType;

	PLANTTYPE(plantType,pSD)

	//see if carbon model is run

	c2lSimulator=pSD->existingSibling("carbonAllocation2Leafs","g");

	//planting time

	pSD->getParent(3)->getChild("plantingTime")->get(plantingTime);

	//plant parameters

	SimulaBase *parameters(GETSHOOTPARAMETERS(plantType));

	if(c2lSimulator){

		SLASimulator=parameters->existingChild("specificLeafArea");

		if(!SLASimulator) SLASimulator=parameters->getChild("leafAreaRatio");

		if(SLASimulator->getUnit()=="cm2/g"){

			SLAisLAR_=false;

		}else{

			if(SLASimulator->getUnit()!="g/cm2") msg::error("LeafArea: Unkown unit for "+SLASimulator->getPath()+". Use g/cm2 or cm2/g");

		}

	}else{

		//plant parameters

		SLASimulator=parameters->existingChild("leafAreaExpansionRate");

		if (!SLASimulator) SLASimulator=parameters->getChild("leafAreaExpantionRate");

	}

	//carbon conversion factor.

	CinDryWeight=pSD->getParent(3)->getChild("carbonToDryWeightRatio","100%");

}

void LeafArea::calculate(const Time &t, double&r){

	//note that the original version starts with an exponential growth curve instead of using photosynthesis

	//local time

	Time localTime(t - plantingTime);

	if(c2lSimulator){

		//use carbon allocation to compute leaf growth

		//get portion of shoot allocated carbon that is used for leafs.

		c2lSimulator->getRate(t,r);

		//get specific leaf area SLA

		double sla;

		SLASimulator->get(localTime,sla);

		SLASimulator->get(t - plantingTime,sla);

		if(sla<1e-10) msg::error("LeafArea: SLA is too small.");

		if(SLAisLAR_){

			r/=sla;
@@ -192,15 +153,37 @@ void LeafArea::calculate(const Time &t, double&r){ 
		CinDryWeight->get(t,cdw);

		//multiply and return the result - assume that carbon to drymatter ratio is 0.54

		r/=cdw;

		if(r<0) msg::error("LeafArea:Negative leaf growth");

		if(std::isnan(r)) {

			msg::error("LeafArea: area is NaN");

		}

	}else if(potentialLA){

		//use potential leaf growth and add stress if available.

		potentialLA->getRate(t,r);

	}else{

		//no carbon model, use potential growth

		//code duplication of the potential leaf area growth function.

		//this is only here for backward compatibility reasons, as no we use the replace option in the inputfiles to switch from potential to actual (carbon driven).

		SLASimulator->get(localTime,r);

		//compute potential growth from the input parameters, or the tillers.

		if(LASimulator){

			LASimulator->get(t-plantingTime,r);

		}else{

			SimulaBase::List l;

			pSD->getSibling("tillers")->getAllChildren(l,t);

			pSD->getSibling("mainShoot")->getChild("potentialLeafArea")->getRate(t,r);

			for(auto &i:l){

				double la;

				i->getChild("potentialLeafArea")->getRate(t,la);

				r+=la;

			}

		}

	}

	//relative growth rate for potential growth

	double c;

	if(rgrCoefficient){

		rgrCoefficient->get(t,c);

		r*=c;

	}

	if(stress){

		stress->get(t,c);

		if(c<1)	r*=c;

	}

	if(r<0) {

		msg::warning("LeafArea: fixing negative growth");

		r=0;

	}

	if(std::isnan(r)) {

		msg::error("LeafArea: area is NaN");

OpenSimRoot/src/modules/ShootModule/LeafArea.hpp

+1 −18
Original line number Diff line number Diff line
@@ -29,7 +29,7 @@ public: 
	std::string getName()const;

protected:

	void calculate(const Time &t, double&);

	SimulaBase *SLASimulator, *c2lSimulator, *CinDryWeight;

	SimulaBase *SLASimulator, *c2lSimulator, *potentialLA, *rgrCoefficient, *stress, *CinDryWeight, *LASimulator;

	Time plantingTime;

	bool SLAisLAR_;

};
@@ -42,24 +42,7 @@ protected: 
	SimulaBase *pSunlitLeafAreaIndex;

	double areaPerPlant, cachedTime, cachedSLAI;

};

class PotentialLeafArea:public DerivativeBase {

public:

	PotentialLeafArea(SimulaDynamic* const pSV);

	std::string getName()const;

protected:

	void calculate(const Time &t, double&);

	SimulaBase *LASimulator;

	Time plantingTime;

};



class StressAdjustedPotentialLeafArea:public DerivativeBase {

public:

	StressAdjustedPotentialLeafArea(SimulaDynamic* const pSV);

	std::string getName()const;

protected:

	void calculate(const Time &t, double&);

	SimulaBase *rgrCoefficient, *potential, *stress;

};

class LeafAreaIndex:public DerivativeBase {

public:

	LeafAreaIndex(SimulaDynamic* const pSD);