WebSVN - public iLand - Diff - Rev 278 and 280


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':
/** @class ResourceUnit
  ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components
  (Climate, Soil, Water, ...).

  */
#include <QtCore>
#include "global.h"

#include "resourceunit.h"
#include "resourceunitspecies.h"
#include "speciesset.h"
#include "species.h"
#include "production3pg.h"
#include "model.h"
#include "climate.h"
#include "watercycle.h"
#include "helper.h"

ResourceUnit::~ResourceUnit()
{
    if (mWater)
        delete mWater;
}

ResourceUnit::ResourceUnit(const int index)
{
    mSpeciesSet = 0;
    mClimate = 0;
    mIndex = index;
    mWater = new WaterCycle();

    mTrees.reserve(100); // start with space for 100 trees.
}

void ResourceUnit::setup()
{
    mWater->setup(this);
}

/// set species and setup the species-per-RU-data
void ResourceUnit::setSpeciesSet(SpeciesSet *set)
{
    mSpeciesSet = set;
    mRUSpecies.clear();
    mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated
    for (int i=0;i<set->count();i++) {
        Species *s = const_cast<Species*>(mSpeciesSet->species(i));
        if (!s)
            throw IException("ResourceUnit::setSpeciesSet: invalid index!");

        /* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.
           If the container memory is relocated (QVector), the pointer gets invalid!!!
           Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */
        mRUSpecies[i].setup(s,this); // setup this element

    }
}

ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)
{
    return mRUSpecies[species->index()];
}

Tree &ResourceUnit::newTree()
{
    // start simple: just append to the vector...
    mTrees.append(Tree());
    return mTrees.back();
}

/// remove dead trees from tree list
/// reduce size of vector if lots of space is free
/// tests showed that this way of cleanup is very fast,
/// because no memory allocations are performed (simple memmove())
/// when trees are moved.
void ResourceUnit::cleanTreeList()
{
    QVector<Tree>::iterator last=mTrees.end()-1;
    QVector<Tree>::iterator current = mTrees.begin();
    while (last>=current && (*last).isDead())
        --last;

    while (current<last) {
        if ((*current).isDead()) {
            *current = *last; // copy data!
            --last; //
            while (last>=current && (*last).isDead())
                --last;
        }
        ++current;
    }
    ++last; // last points now to the first dead tree

    // free ressources
    if (last!=mTrees.end()) {
        mTrees.erase(last, mTrees.end());
        if (mTrees.capacity()>100) {
            if (mTrees.count() / double(mTrees.capacity()) < 0.2) {
                int target_size = mTrees.count()*2;
                qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;
                mTrees.reserve(qMax(target_size, 100));
            }
        }
    }
}

void ResourceUnit::newYear()
{
    mAggregatedWLA = 0.;
    mAggregatedLA = 0.;
    mAggregatedLR = 0.;
    mEffectiveArea = 0.;
    mPixelCount = mStockedPixelCount = 0;
    // clear statistics global and per species...
    ResourceUnitSpecies *i;
    QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end();
    mStatistics.clear();
    for (i=mRUSpecies.begin(); i!=iend; ++i) {
        i->statistics().clear();
        i->statisticsDead().clear();
        i->statisticsMgmt().clear();
    }

}

/** production() is the "stand-level" part of the biomass production (3PG).
    - The amount of radiation intercepted by the stand is calculated
    - The 3PG production for each species and ressource unit is invoked  */
void ResourceUnit::production()
{

    if (mAggregatedWLA==0 || mPixelCount==0) {
        // nothing to do...
        return;
    }

    // the pixel counters are filled during the height-grid-calculations
    mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)

    // calculate the leaf area index (LAI)
    double LAI = mAggregatedLA / mStockedArea;
    // calculate the intercepted radiation fraction using the law of Beer Lambert
    const double k = Model::settings().lightExtinctionCoefficient;
    double interception_fraction = 1. - exp(-k * LAI);
    mEffectiveArea = mStockedArea * interception_fraction; // m2

    // calculate the total weighted leaf area on this RU:
    mLRI_modification = interception_fraction *  mStockedArea / mAggregatedWLA;
    if (mLRI_modification == 0.)
        qDebug() << "lri modifaction==0!";


    DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")
            .arg(LAI)
            .arg(interception_fraction)
            .arg(mLRI_modification)
            .arg(mStockedArea);
    );
    // soil water model - this determines soil water contents needed for response calculations
    {
    DebugTimer tw("water:run");
    mWater->run();
    }

    // invoke species specific calculation (3PG)
    ResourceUnitSpecies *i;
    QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end();

    for (i=mRUSpecies.begin(); i!=iend; ++i) {
        i->calculate();
 
        qDebug() << "species" << (*i).species()->id() << "raw_gpp_m2" << i->prod3PG().GPPperArea() << "area:" << productiveArea() << "gpp:" << productiveArea()*i->prod3PG().GPPperArea();
    }
}

void ResourceUnit::calculateInterceptedArea()
{
    if (mAggregatedLR==0) {
        mEffectiveArea_perWLA = 0.;
        return;
    }
    Q_ASSERT(mAggregatedLR>0.);
    mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;
    qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR  << "eff.area./wla:" << mEffectiveArea_perWLA;
}

void ResourceUnit::yearEnd()
{
    // calculate statistics for all tree species of the ressource unit
    int c = mRUSpecies.count();
    for (int i=0;i<c; i++) {
        mRUSpecies[i].statisticsDead().calculate(); // calculate the dead trees
        mRUSpecies[i].statisticsMgmt().calculate(); // stats of removed trees
        mRUSpecies[i].updateGWL(); // get sum of dead trees (died + removed)
        mRUSpecies[i].statistics().calculate(); // calculate the living (and add removed volume to gwl)
        mStatistics.add(mRUSpecies[i].statistics());
    }
    mStatistics.calculate(); // aggreagte on stand level
}

/// refresh of tree based statistics.
void ResourceUnit::createStandStatistics()
{
    // clear statistics (ru-level and ru-species level)
    mStatistics.clear();
    for (int i=0;i<mRUSpecies.count();i++) {
        mRUSpecies[i].statistics().clear();
        mRUSpecies[i].statisticsDead().clear();
        mRUSpecies[i].statisticsMgmt().clear();
    }

    // add all trees to the statistics objects of the species
    foreach(const Tree &t, mTrees) {
        if (!t.isDead())
            resourceUnitSpecies(t.species()).statistics().add(&t, 0);
    }
    // summarize statistics for the whole resource unit
    for (int i=0;i<mRUSpecies.count();i++) {
        mRUSpecies[i].statistics().calculate();
        mStatistics.add(mRUSpecies[i].statistics());
    }
}
 

Rev 278	Rev 280
1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':	1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':
2	/** @class ResourceUnit	2	/** @class ResourceUnit
3	ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components	3	ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components
4	(Climate, Soil, Water, ...).	4	(Climate, Soil, Water, ...).
5		5
6	*/	6	*/
7	#include <QtCore>	7	#include <QtCore>
8	#include "global.h"	8	#include "global.h"
9		9
10	#include "resourceunit.h"	10	#include "resourceunit.h"
11	#include "resourceunitspecies.h"	11	#include "resourceunitspecies.h"
12	#include "speciesset.h"	12	#include "speciesset.h"
13	#include "species.h"	13	#include "species.h"
14	#include "production3pg.h"	14	#include "production3pg.h"
15	#include "model.h"	15	#include "model.h"
16	#include "climate.h"	16	#include "climate.h"
17	#include "watercycle.h"	17	#include "watercycle.h"
18	#include "helper.h"	18	#include "helper.h"
19		19
20	ResourceUnit::~ResourceUnit()	20	ResourceUnit::~ResourceUnit()
21	{	21	{
22	if (mWater)	22	if (mWater)
23	delete mWater;	23	delete mWater;
24	}	24	}
25		25
26	ResourceUnit::ResourceUnit(const int index)	26	ResourceUnit::ResourceUnit(const int index)
27	{	27	{
28	mSpeciesSet = 0;	28	mSpeciesSet = 0;
29	mClimate = 0;	29	mClimate = 0;
30	mIndex = index;	30	mIndex = index;
31	mWater = new WaterCycle();	31	mWater = new WaterCycle();
32		32
33	mTrees.reserve(100); // start with space for 100 trees.	33	mTrees.reserve(100); // start with space for 100 trees.
34	}	34	}
35		35
36	void ResourceUnit::setup()	36	void ResourceUnit::setup()
37	{	37	{
38	mWater->setup(this);	38	mWater->setup(this);
39	}	39	}
40		40
41	/// set species and setup the species-per-RU-data	41	/// set species and setup the species-per-RU-data
42	void ResourceUnit::setSpeciesSet(SpeciesSet *set)	42	void ResourceUnit::setSpeciesSet(SpeciesSet *set)
43	{	43	{
44	mSpeciesSet = set;	44	mSpeciesSet = set;
45	mRUSpecies.clear();	45	mRUSpecies.clear();
46	mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated	46	mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated
47	for (int i=0;i<set->count();i++) {	47	for (int i=0;i<set->count();i++) {
48	Species s = const_cast<Species>(mSpeciesSet->species(i));	48	Species s = const_cast<Species>(mSpeciesSet->species(i));
49	if (!s)	49	if (!s)
50	throw IException("ResourceUnit::setSpeciesSet: invalid index!");	50	throw IException("ResourceUnit::setSpeciesSet: invalid index!");
51		51
52	/* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.	52	/* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.
53	If the container memory is relocated (QVector), the pointer gets invalid!!!	53	If the container memory is relocated (QVector), the pointer gets invalid!!!
54	Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */	54	Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */
55	mRUSpecies[i].setup(s,this); // setup this element	55	mRUSpecies[i].setup(s,this); // setup this element
56		56
57	}	57	}
58	}	58	}
59		59
60	ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)	60	ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)
61	{	61	{
62	return mRUSpecies[species->index()];	62	return mRUSpecies[species->index()];
63	}	63	}
64		64
65	Tree &ResourceUnit::newTree()	65	Tree &ResourceUnit::newTree()
66	{	66	{
67	// start simple: just append to the vector...	67	// start simple: just append to the vector...
68	mTrees.append(Tree());	68	mTrees.append(Tree());
69	return mTrees.back();	69	return mTrees.back();
70	}	70	}
71		71
72	/// remove dead trees from tree list	72	/// remove dead trees from tree list
73	/// reduce size of vector if lots of space is free	73	/// reduce size of vector if lots of space is free
74	/// tests showed that this way of cleanup is very fast,	74	/// tests showed that this way of cleanup is very fast,
75	/// because no memory allocations are performed (simple memmove())	75	/// because no memory allocations are performed (simple memmove())
76	/// when trees are moved.	76	/// when trees are moved.
77	void ResourceUnit::cleanTreeList()	77	void ResourceUnit::cleanTreeList()
78	{	78	{
79	QVector<Tree>::iterator last=mTrees.end()-1;	79	QVector<Tree>::iterator last=mTrees.end()-1;
80	QVector<Tree>::iterator current = mTrees.begin();	80	QVector<Tree>::iterator current = mTrees.begin();
81	while (last>=current && (*last).isDead())	81	while (last>=current && (*last).isDead())
82	--last;	82	--last;
83		83
84	while (current<last) {	84	while (current<last) {
85	if ((*current).isDead()) {	85	if ((*current).isDead()) {
86	current = last; // copy data!	86	current = last; // copy data!
87	--last; //	87	--last; //
88	while (last>=current && (*last).isDead())	88	while (last>=current && (*last).isDead())
89	--last;	89	--last;
90	}	90	}
91	++current;	91	++current;
92	}	92	}
93	++last; // last points now to the first dead tree	93	++last; // last points now to the first dead tree
94		94
95	// free ressources	95	// free ressources
96	if (last!=mTrees.end()) {	96	if (last!=mTrees.end()) {
97	mTrees.erase(last, mTrees.end());	97	mTrees.erase(last, mTrees.end());
98	if (mTrees.capacity()>100) {	98	if (mTrees.capacity()>100) {
99	if (mTrees.count() / double(mTrees.capacity()) < 0.2) {	99	if (mTrees.count() / double(mTrees.capacity()) < 0.2) {
100	int target_size = mTrees.count()*2;	100	int target_size = mTrees.count()*2;
101	qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;	101	qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;
102	mTrees.reserve(qMax(target_size, 100));	102	mTrees.reserve(qMax(target_size, 100));
103	}	103	}
104	}	104	}
105	}	105	}
106	}	106	}
107		107
108	void ResourceUnit::newYear()	108	void ResourceUnit::newYear()
109	{	109	{
110	mAggregatedWLA = 0.;	110	mAggregatedWLA = 0.;
111	mAggregatedLA = 0.;	111	mAggregatedLA = 0.;
112	mAggregatedLR = 0.;	112	mAggregatedLR = 0.;
113	mEffectiveArea = 0.;	113	mEffectiveArea = 0.;
114	mPixelCount = mStockedPixelCount = 0;	114	mPixelCount = mStockedPixelCount = 0;
115	// clear statistics global and per species...	115	// clear statistics global and per species...
116	ResourceUnitSpecies *i;	116	ResourceUnitSpecies *i;
117	QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end();	117	QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end();
118	mStatistics.clear();	118	mStatistics.clear();
119	for (i=mRUSpecies.begin(); i!=iend; ++i) {	119	for (i=mRUSpecies.begin(); i!=iend; ++i) {
120	i->statistics().clear();	120	i->statistics().clear();
121	i->statisticsDead().clear();	121	i->statisticsDead().clear();
122	i->statisticsMgmt().clear();	122	i->statisticsMgmt().clear();
123	}	123	}
124		124
125	}	125	}
126		126
127	/** production() is the "stand-level" part of the biomass production (3PG).	127	/** production() is the "stand-level" part of the biomass production (3PG).
128	- The amount of radiation intercepted by the stand is calculated	128	- The amount of radiation intercepted by the stand is calculated
129	- The 3PG production for each species and ressource unit is invoked */	129	- The 3PG production for each species and ressource unit is invoked */
130	void ResourceUnit::production()	130	void ResourceUnit::production()
131	{	131	{
132		132
133	if (mAggregatedWLA==0 \|\| mPixelCount==0) {	133	if (mAggregatedWLA==0 \|\| mPixelCount==0) {
134	// nothing to do...	134	// nothing to do...
135	return;	135	return;
136	}	136	}
137		137
138	// the pixel counters are filled during the height-grid-calculations	138	// the pixel counters are filled during the height-grid-calculations
139	mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)	139	mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)
140		140
141	// calculate the leaf area index (LAI)	141	// calculate the leaf area index (LAI)
142	double LAI = mAggregatedLA / mStockedArea;	142	double LAI = mAggregatedLA / mStockedArea;
143	// calculate the intercepted radiation fraction using the law of Beer Lambert	143	// calculate the intercepted radiation fraction using the law of Beer Lambert
144	const double k = Model::settings().lightExtinctionCoefficient;	144	const double k = Model::settings().lightExtinctionCoefficient;
145	double interception_fraction = 1. - exp(-k * LAI);	145	double interception_fraction = 1. - exp(-k * LAI);
146	mEffectiveArea = mStockedArea * interception_fraction; // m2	146	mEffectiveArea = mStockedArea * interception_fraction; // m2
147		147
148	// calculate the total weighted leaf area on this RU:	148	// calculate the total weighted leaf area on this RU:
149	mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA;	149	mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA;
150	if (mLRI_modification == 0.)	150	if (mLRI_modification == 0.)
151	qDebug() << "lri modifaction==0!";	151	qDebug() << "lri modifaction==0!";
152		152
153		153
154	DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")	154	DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")
155	.arg(LAI)	155	.arg(LAI)
156	.arg(interception_fraction)	156	.arg(interception_fraction)
157	.arg(mLRI_modification)	157	.arg(mLRI_modification)
158	.arg(mStockedArea);	158	.arg(mStockedArea);
159	);	159	);
160	// soil water model - this determines soil water contents needed for response calculations	160	// soil water model - this determines soil water contents needed for response calculations
161	{	161	{
162	DebugTimer tw("water:run");	162	DebugTimer tw("water:run");
163	mWater->run();	163	mWater->run();
164	}	164	}
165		165
166	// invoke species specific calculation (3PG)	166	// invoke species specific calculation (3PG)
167	ResourceUnitSpecies *i;	167	ResourceUnitSpecies *i;
168	QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end();	168	QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end();
169		169
170	for (i=mRUSpecies.begin(); i!=iend; ++i) {	170	for (i=mRUSpecies.begin(); i!=iend; ++i) {
171	i->calculate();	171	i->calculate();
172	qDebug() << "species" << (*i).species()->id() << "raw_gpp_m2" << i->prod3PG().GPPperArea();	-
-		172	qDebug() << "species" << (i).species()->id() << "raw_gpp_m2" << i->prod3PG().GPPperArea() << "area:" << productiveArea() << "gpp:" << productiveArea()i->prod3PG().GPPperArea();
173	}	173	}
174	}	174	}
175		175
176	void ResourceUnit::calculateInterceptedArea()	176	void ResourceUnit::calculateInterceptedArea()
177	{	177	{
178	if (mAggregatedLR==0) {	178	if (mAggregatedLR==0) {
179	mEffectiveArea_perWLA = 0.;	179	mEffectiveArea_perWLA = 0.;
180	return;	180	return;
181	}	181	}
182	Q_ASSERT(mAggregatedLR>0.);	182	Q_ASSERT(mAggregatedLR>0.);
183	mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;	183	mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;
184	qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA;	184	qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA;
185	}	185	}
186		186
187	void ResourceUnit::yearEnd()	187	void ResourceUnit::yearEnd()
188	{	188	{
189	// calculate statistics for all tree species of the ressource unit	189	// calculate statistics for all tree species of the ressource unit
190	int c = mRUSpecies.count();	190	int c = mRUSpecies.count();
191	for (int i=0;i<c; i++) {	191	for (int i=0;i<c; i++) {
192	mRUSpecies[i].statisticsDead().calculate(); // calculate the dead trees	192	mRUSpecies[i].statisticsDead().calculate(); // calculate the dead trees
193	mRUSpecies[i].statisticsMgmt().calculate(); // stats of removed trees	193	mRUSpecies[i].statisticsMgmt().calculate(); // stats of removed trees
194	mRUSpecies[i].updateGWL(); // get sum of dead trees (died + removed)	194	mRUSpecies[i].updateGWL(); // get sum of dead trees (died + removed)
195	mRUSpecies[i].statistics().calculate(); // calculate the living (and add removed volume to gwl)	195	mRUSpecies[i].statistics().calculate(); // calculate the living (and add removed volume to gwl)
196	mStatistics.add(mRUSpecies[i].statistics());	196	mStatistics.add(mRUSpecies[i].statistics());
197	}	197	}
198	mStatistics.calculate(); // aggreagte on stand level	198	mStatistics.calculate(); // aggreagte on stand level
199	}	199	}
200		200
201	/// refresh of tree based statistics.	201	/// refresh of tree based statistics.
202	void ResourceUnit::createStandStatistics()	202	void ResourceUnit::createStandStatistics()
203	{	203	{
204	// clear statistics (ru-level and ru-species level)	204	// clear statistics (ru-level and ru-species level)
205	mStatistics.clear();	205	mStatistics.clear();
206	for (int i=0;i<mRUSpecies.count();i++) {	206	for (int i=0;i<mRUSpecies.count();i++) {
207	mRUSpecies[i].statistics().clear();	207	mRUSpecies[i].statistics().clear();
208	mRUSpecies[i].statisticsDead().clear();	208	mRUSpecies[i].statisticsDead().clear();
209	mRUSpecies[i].statisticsMgmt().clear();	209	mRUSpecies[i].statisticsMgmt().clear();
210	}	210	}
211		211
212	// add all trees to the statistics objects of the species	212	// add all trees to the statistics objects of the species
213	foreach(const Tree &t, mTrees) {	213	foreach(const Tree &t, mTrees) {
214	if (!t.isDead())	214	if (!t.isDead())
215	resourceUnitSpecies(t.species()).statistics().add(&t, 0);	215	resourceUnitSpecies(t.species()).statistics().add(&t, 0);
216	}	216	}
217	// summarize statistics for the whole resource unit	217	// summarize statistics for the whole resource unit
218	for (int i=0;i<mRUSpecies.count();i++) {	218	for (int i=0;i<mRUSpecies.count();i++) {
219	mRUSpecies[i].statistics().calculate();	219	mRUSpecies[i].statistics().calculate();
220	mStatistics.add(mRUSpecies[i].statistics());	220	mStatistics.add(mRUSpecies[i].statistics());
221	}	221	}
222	}	222	}
223		223

Subversion Repositories public iLand

(root)/src/core/resourceunit.cpp @ 1222 - Rev 278 → 280