Commit 62a99ece authored by jouke's avatar jouke
Browse files

Fixes for the fertilization code. Mostly surface area was not correctly

calculated. But also reduced scope of the variables, reducing memory
usage.
parent 202aca84
Loading
Loading
Loading
Loading
+24 −37
Original line number Diff line number Diff line
@@ -20,13 +20,9 @@
#include "Fertilization.hpp"


Fertilization::Fertilization(const Mesh_base &mesh, const std::string &nutrient)
Fertilization::Fertilization(const Mesh_base &mesh, const std::string &nutrient):
  rate(nullptr), depthMultipliers(mesh.getNumNP())
{
	// Z is the vertical coordinate in this mesh!
	auto &z(mesh.getCordZ());
	numNodes = z.size();
	tVector meshVols = mesh.getvol();

	SimulaBase* p;
	p = ORIGIN->getPath("environment/soil/" + nutrient + "/fertilization");
	rate = p->getChild("fertilizerRate");
@@ -38,14 +34,6 @@ Fertilization::Fertilization(const Mesh_base &mesh, const std::string &nutrient)
	else if (rate->getUnit() == "g/cm2/day"){
		SimulaBase *probe;
		double area;
		probe = ORIGIN->existingPath("environment/plantingScheme/rowSpacing");
		if (probe){
			double wid, len;
			probe->get(wid);
			probe = ORIGIN->getPath("environment/plantingScheme/inTheRowSpacing");
			probe->get(len);
			area = wid*len;
		} else{
		Coordinate minCorner, maxCorner;
		probe = ORIGIN->getPath("environment/dimensions/minCorner");
		probe->get(minCorner);
@@ -53,31 +41,30 @@ Fertilization::Fertilization(const Mesh_base &mesh, const std::string &nutrient)
		probe->get(maxCorner);
		// Here y is the vertical coordinate!
		area = fabs(minCorner.x - maxCorner.x)*fabs(minCorner.z - maxCorner.z);
		}
		unitCorrection = area;
	} else{
		msg::error("Unknown unit");
	}

	depthDistribution = p->getChild("distribution");
	SimulaBase* depthDistribution = p->existingChild("distribution");
	const size_t numNP=mesh.getNumNP();
	tVector nodeFactors(0.,numNP);
	if(depthDistribution) {
		for (unsigned int i = 0; i != numNP; ++i) {
			double temp;
	depthFactor = 0;
	tVector nodeFactors;
	nodeFactors.resize(numNodes);
	for (unsigned int i = 0; i != numNodes; ++i) {
		depthDistribution->get(z[i], temp);
		nodeFactors[i] = temp*meshVols[i];
		depthFactor += temp*meshVols[i];
			depthDistribution->get(mesh.getCordZ()[i], temp);
			nodeFactors[i] = temp*mesh.getvol()[i];
		}
	}else{
		for (auto i : mesh.getTopBoundaryNodes()) {
			nodeFactors[i] = mesh.getvol()[i];
		}
	depthMultipliers.resize(numNodes);
	for (unsigned int i = 0; i < numNodes; ++i){
		depthMultipliers[i] = unitCorrection*nodeFactors[i]/depthFactor;
	}
	double depthFactor = nodeFactors.sum();
	depthMultipliers = unitCorrection*nodeFactors/depthFactor;
}

void Fertilization::fertilization(const double & t0, tVector &fertil){
	// output:  fertil
	if(fertil.size()!=numNodes) fertil.resize(numNodes);
	double currentRate;
	rate->get(t0, currentRate);
	fertil = depthMultipliers*currentRate;
+1 −4
Original line number Diff line number Diff line
@@ -40,11 +40,8 @@ public:
	*/
	void fertilization(const double & t0, tVector &fert);
private:
	SimulaBase *rate, *depthDistribution;
	double depthFactor;
	unsigned int numNodes;
	SimulaBase *rate;
	tVector depthMultipliers;

};

#endif