Subversion Repositories public iLand

Rev

Rev 882 | Rev 890 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1
 
671 werner 2
/********************************************************************************************
3
**    iLand - an individual based forest landscape and disturbance model
4
**    http://iland.boku.ac.at
5
**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl
6
**
7
**    This program is free software: you can redistribute it and/or modify
8
**    it under the terms of the GNU General Public License as published by
9
**    the Free Software Foundation, either version 3 of the License, or
10
**    (at your option) any later version.
11
**
12
**    This program is distributed in the hope that it will be useful,
13
**    but WITHOUT ANY WARRANTY; without even the implied warranty of
14
**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15
**    GNU General Public License for more details.
16
**
17
**    You should have received a copy of the GNU General Public License
18
**    along with this program.  If not, see <http://www.gnu.org/licenses/>.
19
********************************************************************************************/
20
 
92 Werner 21
/** @class Model
247 werner 22
  Main object of the iLand model composited of various sub models / sub components.
697 werner 23
  @ingroup core
24
  The class Model is the top level container of iLand. The Model holds a collection of ResourceUnits, links to SpeciesSet and Climate.
25
  ResourceUnit are grid cells with (currently) a size of 1 ha (100x100m). Many stand level processes (NPP produciton, WaterCycle) operate on this
26
  level.
27
  The Model also contain the landscape-wide 2m LIF-grid (http://iland.boku.ac.at/competition+for+light).
28
 
92 Werner 29
  */
30
#include "global.h"
31
#include "model.h"
286 werner 32
#include "sqlhelper.h"
92 Werner 33
 
105 Werner 34
#include "xmlhelper.h"
808 werner 35
#include "debugtimer.h"
281 werner 36
#include "environment.h"
340 werner 37
#include "timeevents.h"
106 Werner 38
#include "helper.h"
189 iland 39
#include "resourceunit.h"
208 werner 40
#include "climate.h"
92 Werner 41
#include "speciesset.h"
106 Werner 42
#include "standloader.h"
43
#include "tree.h"
185 werner 44
#include "management.h"
200 werner 45
#include "modelsettings.h"
615 werner 46
#include "standstatistics.h"
543 werner 47
#include "mapgrid.h"
632 werner 48
#include "modelcontroller.h"
641 werner 49
#include "modules.h"
654 werner 50
#include "dem.h"
92 Werner 51
 
202 werner 52
#include "outputmanager.h"
176 werner 53
 
105 Werner 54
#include <QtCore>
55
#include <QtXml>
56
 
107 Werner 57
/** iterate over all trees of the model. return NULL if all trees processed.
58
  Usage:
59
  @code
60
  AllTreeIterator trees(model);
61
  while (Tree *tree = trees.next()) { // returns NULL when finished.
62
     tree->something(); // do something
63
  }
281 werner 64
  @endcode  */
107 Werner 65
Tree *AllTreeIterator::next()
66
{
143 Werner 67
 
107 Werner 68
    if (!mTreeEnd) {
69
        // initialize to first ressource unit
70
        mRUIterator = mModel->ruList().constBegin();
314 werner 71
        // fast forward to the first RU with trees
72
        while (mRUIterator!=mModel->ruList().constEnd()) {
73
            if ((*mRUIterator)->trees().count()>0)
74
                break;
753 werner 75
            ++mRUIterator;
314 werner 76
        }
77
            // finished if all RU processed
78
        if (mRUIterator == mModel->ruList().constEnd())
79
            return NULL;
143 Werner 80
        mTreeEnd = &((*mRUIterator)->trees().back()) + 1; // let end point to "1 after end" (STL-style)
107 Werner 81
        mCurrent = &((*mRUIterator)->trees().front());
82
    }
83
    if (mCurrent==mTreeEnd) {
753 werner 84
        ++mRUIterator; // switch to next RU (loop until RU with trees is found)
314 werner 85
        while (mRUIterator!=mModel->ruList().constEnd()) {
86
            if ((*mRUIterator)->trees().count()>0) {
87
                break;
88
            }
753 werner 89
            ++mRUIterator;
314 werner 90
        }
107 Werner 91
        if (mRUIterator == mModel->ruList().constEnd()) {
143 Werner 92
            mCurrent = NULL;
93
            return NULL; // finished!!
107 Werner 94
        }else {
143 Werner 95
            mTreeEnd = &((*mRUIterator)->trees().back()) + 1;
107 Werner 96
            mCurrent = &((*mRUIterator)->trees().front());
97
        }
98
    }
143 Werner 99
 
100
    return mCurrent++;
107 Werner 101
}
157 werner 102
Tree *AllTreeIterator::nextLiving()
103
{
104
    while (Tree *t = next())
158 werner 105
        if (!t->isDead()) return t;
157 werner 106
    return NULL;
107
}
108
Tree *AllTreeIterator::current() const
109
{
110
    return mCurrent?mCurrent-1:NULL;
111
}
107 Werner 112
 
113
 
200 werner 114
ModelSettings Model::mSettings;
92 Werner 115
Model::Model()
116
{
117
    initialize();
137 Werner 118
    GlobalSettings::instance()->setModel(this);
767 werner 119
    GlobalSettings::instance()->resetScriptEngine(); // clear the script
130 Werner 120
    QString dbg="running in release mode.";
442 werner 121
    DBGMODE( dbg="running in debug mode."; );
130 Werner 122
    qDebug() << dbg;
92 Werner 123
}
124
 
125
Model::~Model()
126
{
127
    clear();
137 Werner 128
    GlobalSettings::instance()->setModel(NULL);
92 Werner 129
}
130
 
131
/** Initial setup of the Model.
132
  */
133
void Model::initialize()
134
{
151 iland 135
   mSetup = false;
162 werner 136
   GlobalSettings::instance()->setCurrentYear(0);
151 iland 137
   mGrid = 0;
138
   mHeightGrid = 0;
185 werner 139
   mManagement = 0;
281 werner 140
   mEnvironment = 0;
340 werner 141
   mTimeEvents = 0;
549 werner 142
   mStandGrid = 0;
641 werner 143
   mModules = 0;
654 werner 144
   mDEM = 0;
92 Werner 145
}
146
 
261 werner 147
/** sets up the simulation space.
148
*/
103 Werner 149
void Model::setupSpace()
150
{
194 werner 151
    XmlHelper xml(GlobalSettings::instance()->settings().node("model.world"));
192 werner 152
    double cellSize = xml.value("cellSize", "2").toDouble();
153
    double width = xml.value("width", "100").toDouble();
154
    double height = xml.value("height", "100").toDouble();
549 werner 155
    double buffer = xml.value("buffer", "60").toDouble();
156
    mModelRect = QRectF(0., 0., width, height);
157
 
103 Werner 158
    qDebug() << QString("setup of the world: %1x%2m with cell-size=%3m and %4m buffer").arg(width).arg(height).arg(cellSize).arg(buffer);
159
 
160
    QRectF total_grid(QPointF(-buffer, -buffer), QPointF(width+buffer, height+buffer));
161
    qDebug() << "setup grid rectangle:" << total_grid;
162
 
151 iland 163
    if (mGrid)
164
        delete mGrid;
103 Werner 165
    mGrid = new FloatGrid(total_grid, cellSize);
156 werner 166
    mGrid->initialize(1.f);
151 iland 167
    if (mHeightGrid)
168
        delete mHeightGrid;
169
    mHeightGrid = new HeightGrid(total_grid, cellSize*5);
285 werner 170
    mHeightGrid->wipe(); // set all to zero
156 werner 171
    Tree::setGrid(mGrid, mHeightGrid);
105 Werner 172
 
569 werner 173
    // setup the spatial location of the project area
174
    if (xml.hasNode("location")) {
175
        // setup of spatial location
176
        double loc_x = xml.valueDouble("location.x");
177
        double loc_y = xml.valueDouble("location.y");
178
        double loc_z = xml.valueDouble("location.z");
179
        double loc_rot = xml.valueDouble("location.rotation");
180
        setupGISTransformation(loc_x, loc_y, loc_z, loc_rot);
181
        qDebug() << "setup of spatial location: x/y/z" << loc_x << loc_y << loc_z << "rotation:" << loc_rot;
182
    } else {
183
        setupGISTransformation(0., 0., 0., 0.);
184
    }
185
 
567 werner 186
    // load environment (multiple climates, speciesSets, ...
187
    if (mEnvironment)
188
        delete mEnvironment;
189
    mEnvironment = new Environment();
281 werner 190
 
567 werner 191
    if (xml.valueBool("environmentEnabled", false)) {
192
        QString env_file = GlobalSettings::instance()->path(xml.value("environmentFile"));
193
        bool grid_mode = (xml.value("environmentMode")=="grid");
194
        QString grid_file = GlobalSettings::instance()->path(xml.value("environmentGrid"));
195
        if (grid_mode && QFile::exists(grid_file))
196
            mEnvironment->setGridMode(grid_file);
197
 
198
        if (!mEnvironment->loadFromFile(env_file))
199
            return;
200
        // retrieve species sets and climates:
201
        mSpeciesSets << mEnvironment->speciesSetList();
202
        mClimates << mEnvironment->climateList();
203
    } else {
204
        // load and prepare default values
205
        // (2) SpeciesSets: currently only one a global species set.
206
        SpeciesSet *speciesSet = new SpeciesSet();
207
        mSpeciesSets.push_back(speciesSet);
208
        speciesSet->setup();
209
        // Climate...
210
        Climate *c = new Climate();
211
        mClimates.push_back(c);
212
        mEnvironment->setDefaultValues(c, speciesSet);
213
    } // environment?
214
 
215
    // time series data
216
    if (xml.valueBool(".timeEventsEnabled", false)) {
217
        mTimeEvents = new TimeEvents();
584 werner 218
        mTimeEvents->loadFromFile(GlobalSettings::instance()->path(xml.value("timeEventsFile"), "script"));
567 werner 219
    }
220
 
646 werner 221
 
105 Werner 222
    // simple case: create ressource units in a regular grid.
194 werner 223
    if (xml.valueBool("resourceUnitsAsGrid")) {
881 werner 224
 
281 werner 225
        mRUmap.setup(QRectF(0., 0., width, height),100.); // Grid, that holds positions of resource units
881 werner 226
        mRUmap.wipe();
646 werner 227
 
539 werner 228
        bool mask_is_setup = false;
229
        if (xml.valueBool("standGrid.enabled")) {
230
            QString fileName = GlobalSettings::instance()->path(xml.value("standGrid.fileName"));
882 werner 231
            mStandGrid = new MapGrid(fileName,false); // create stand grid index later
543 werner 232
 
549 werner 233
            if (mStandGrid->isValid()) {
714 werner 234
                for (int i=0;i<mStandGrid->grid().count();i++) {
235
                    const int &grid_value = mStandGrid->grid().constValueAtIndex(i);
236
                    mHeightGrid->valueAtIndex(i).setValid( grid_value > -1 );
881 werner 237
                    if (grid_value>-1)
238
                        mRUmap.valueAt(mStandGrid->grid().cellCenterPoint(i)) = (ResourceUnit*)1;
714 werner 239
                    if (grid_value < -1)
718 werner 240
                        mHeightGrid->valueAtIndex(i).setForestOutside(true);
714 werner 241
                }
539 werner 242
            }
243
            mask_is_setup = true;
802 werner 244
        } else {
245
            if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
246
                // in the case we have no stand grid but only a large rectangle (without the torus option)
247
                // we assume a forest outside
248
                for (int i=0;i<mHeightGrid->count();++i) {
249
                    const QPointF &p = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i));
250
                    if (p.x() < 0. || p.x()>width || p.y()<0. || p.y()>height) {
251
                        mHeightGrid->valueAtIndex(i).setForestOutside(true);
252
                        mHeightGrid->valueAtIndex(i).setValid(false);
253
                    }
254
                }
255
 
256
            }
539 werner 257
        }
258
 
881 werner 259
        ResourceUnit **p; // ptr to ptr!
260
        ResourceUnit *new_ru;
261
 
262
        int ru_index = 0;
263
        for (p=mRUmap.begin(); p!=mRUmap.end(); ++p) {
264
            QRectF r = mRUmap.cellRect(mRUmap.indexOf(p));
265
            mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used.
889 werner 266
            if (!mStandGrid || !mStandGrid->isValid() || *p>0) {
881 werner 267
                // create resource units for valid positions only
268
                new_ru = new ResourceUnit(ru_index++); // create resource unit
269
                new_ru->setClimate( mEnvironment->climate() );
270
                new_ru->setSpeciesSet( mEnvironment->speciesSet() );
271
                new_ru->setup();
272
                new_ru->setID( mEnvironment->currentID() ); // set id of resource unit in grid mode
273
                new_ru->setBoundingBox(r);
274
                mRU.append(new_ru);
889 werner 275
                *p = new_ru; // save in the RUmap grid
881 werner 276
            }
277
        }
278
 
889 werner 279
        if (mStandGrid && mStandGrid->isValid())
882 werner 280
            mStandGrid->createIndex();
881 werner 281
        // now store the pointers in the grid.
282
        // Important: This has to be done after the mRU-QList is complete - otherwise pointers would
283
        // point to invalid memory when QList's memory is reorganized (expanding)
284
//        ru_index = 0;
285
//        for (p=mRUmap.begin();p!=mRUmap.end(); ++p) {
286
//            *p = mRU.value(ru_index++);
287
//        }
288
        qDebug() << "created a grid of ResourceUnits: count=" << mRU.count() << "number of RU-map-cells:" << mRUmap.count();
289
 
290
 
574 werner 291
        calculateStockableArea();
292
 
285 werner 293
        // setup of the project area mask
539 werner 294
        if (!mask_is_setup && xml.valueBool("areaMask.enabled", false) && xml.hasNode("areaMask.imageFile")) {
285 werner 295
            // to be extended!!! e.g. to load ESRI-style text files....
296
            // setup a grid with the same size as the height grid...
297
            FloatGrid tempgrid((int)mHeightGrid->cellsize(), mHeightGrid->sizeX(), mHeightGrid->sizeY());
298
            QString fileName = GlobalSettings::instance()->path(xml.value("areaMask.imageFile"));
299
            loadGridFromImage(fileName, tempgrid); // fetch from image
300
            for (int i=0;i<tempgrid.count(); i++)
301
                mHeightGrid->valueAtIndex(i).setValid( tempgrid.valueAtIndex(i)>0.99 );
302
            qDebug() << "loaded project area mask from" << fileName;
303
        }
304
 
590 werner 305
        // list of "valid" resource units
306
        QList<ResourceUnit*> valid_rus;
307
        foreach(ResourceUnit* ru, mRU)
308
            if (ru->id()!=-1)
309
                valid_rus.append(ru);
310
 
654 werner 311
        // setup of the digital elevation map (if present)
312
        QString dem_file = xml.value("DEM");
313
        if (!dem_file.isEmpty()) {
656 werner 314
            mDEM = new DEM(GlobalSettings::instance()->path(dem_file));
315
            // add them to the visuals...
316
            GlobalSettings::instance()->controller()->addGrid(mDEM, "DEM height", GridViewRainbow, 0, 1000);
317
            GlobalSettings::instance()->controller()->addGrid(mDEM->slopeGrid(), "DEM slope", GridViewRainbow, 0, 3);
318
            GlobalSettings::instance()->controller()->addGrid(mDEM->aspectGrid(), "DEM aspect", GridViewRainbow, 0, 360);
319
            GlobalSettings::instance()->controller()->addGrid(mDEM->viewGrid(), "DEM view", GridViewGray, 0, 1);
320
 
654 werner 321
        }
322
 
646 werner 323
        // setup of external modules
324
        mModules->setup();
325
        if (mModules->hasSetupResourceUnits()) {
326
            for (ResourceUnit **p=mRUmap.begin(); p!=mRUmap.end(); ++p) {
327
                QRectF r = mRUmap.cellRect(mRUmap.indexOf(p));
328
                mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used.
329
                mModules->setupResourceUnit( *p );
330
            }
331
        }
332
 
767 werner 333
        // setup of scripting environment
334
        ScriptGlobal::setupGlobalScripting();
335
 
123 Werner 336
        // setup the helper that does the multithreading
590 werner 337
        threadRunner.setup(valid_rus);
194 werner 338
        threadRunner.setMultithreading(GlobalSettings::instance()->settings().valueBool("system.settings.multithreading"));
123 Werner 339
        threadRunner.print();
105 Werner 340
 
194 werner 341
    } else  {
342
        throw IException("resourceUnitsAsGrid MUST be set to true - at least currently :)");
105 Werner 343
    }
120 Werner 344
    mSetup = true;
103 Werner 345
}
346
 
347
 
92 Werner 348
/** clear() frees all ressources allocated with the run of a simulation.
349
 
350
  */
351
void Model::clear()
352
{
143 Werner 353
    mSetup = false;
151 iland 354
    qDebug() << "Model clear: attempting to clear" << mRU.count() << "RU, " << mSpeciesSets.count() << "SpeciesSets.";
92 Werner 355
    // clear ressource units
103 Werner 356
    qDeleteAll(mRU); // delete ressource units (and trees)
92 Werner 357
    mRU.clear();
103 Werner 358
 
359
    qDeleteAll(mSpeciesSets); // delete species sets
92 Werner 360
    mSpeciesSets.clear();
103 Werner 361
 
208 werner 362
    // delete climate data
363
    qDeleteAll(mClimates);
364
 
151 iland 365
    // delete the grids
366
    if (mGrid)
367
        delete mGrid;
368
    if (mHeightGrid)
369
        delete mHeightGrid;
185 werner 370
    if (mManagement)
371
        delete mManagement;
281 werner 372
    if (mEnvironment)
373
        delete mEnvironment;
340 werner 374
    if (mTimeEvents)
375
        delete mTimeEvents;
549 werner 376
    if (mStandGrid)
377
        delete mStandGrid;
641 werner 378
    if (mModules)
379
        delete mModules;
654 werner 380
    if (mDEM)
381
        delete mDEM;
103 Werner 382
 
185 werner 383
    mGrid = 0;
384
    mHeightGrid = 0;
385
    mManagement = 0;
281 werner 386
    mEnvironment = 0;
340 werner 387
    mTimeEvents = 0;
549 werner 388
    mStandGrid  = 0;
641 werner 389
    mModules = 0;
654 werner 390
    mDEM = 0;
185 werner 391
 
583 werner 392
    GlobalSettings::instance()->outputManager()->close();
393
 
92 Werner 394
    qDebug() << "Model ressources freed.";
395
}
396
 
397
/** Setup of the Simulation.
398
  This really creates the simulation environment and does the setup of various aspects.
399
  */
102 Werner 400
void Model::loadProject()
92 Werner 401
{
109 Werner 402
    DebugTimer dt("load project");
93 Werner 403
    GlobalSettings *g = GlobalSettings::instance();
679 werner 404
    g->printDirecories();
102 Werner 405
    const XmlHelper &xml = g->settings();
189 iland 406
 
93 Werner 407
    g->clearDatabaseConnections();
92 Werner 408
    // database connections: reset
94 Werner 409
    GlobalSettings::instance()->clearDatabaseConnections();
286 werner 410
    // input and climate connection
411
    // see initOutputDatabase() for output database
191 werner 412
    QString dbPath = g->path( xml.value("system.database.in"), "database");
194 werner 413
    GlobalSettings::instance()->setupDatabaseConnection("in", dbPath, true);
193 werner 414
    dbPath = g->path( xml.value("system.database.climate"), "database");
194 werner 415
    GlobalSettings::instance()->setupDatabaseConnection("climate", dbPath, true);
92 Werner 416
 
200 werner 417
    mSettings.loadModelSettings();
418
    mSettings.print();
416 werner 419
    // random seed: if stored value is <> 0, use this as the random seed (and produce hence always an equal sequence of random numbers)
420
    uint seed = xml.value("system.settings.randomSeed","0").toUInt();
708 werner 421
    RandomGenerator::setup(RandomGenerator::ergMersenneTwister, seed); // use the MersenneTwister as default
428 werner 422
    // linearization of expressions: if true *and* linearize() is explicitely called, then
423
    // function results will be cached over a defined range of values.
424
    bool do_linearization = xml.valueBool("system.settings.expressionLinearizationEnabled", false);
425
    Expression::setLinearizationEnabled(do_linearization);
431 werner 426
    if (do_linearization)
427
        qDebug() << "The linearization of certains expressions is enabled (performance optimization).";
428
 
429
    // log level
430
    QString log_level = xml.value("system.settings.logLevel", "debug").toLower();
431
    if (log_level=="debug") setLogLevel(0);
432
    if (log_level=="info") setLogLevel(1);
433
    if (log_level=="warning") setLogLevel(2);
434
    if (log_level=="error") setLogLevel(3);
435
 
475 werner 436
    // snag dynamics / soil model enabled? (info used during setup of world)
437
    changeSettings().carbonCycleEnabled = xml.valueBool("model.settings.carbonCycleEnabled", false);
528 werner 438
    // class size of snag classes
439
    Snag::setupThresholds(xml.valueDouble("model.settings.soil.swdDBHClass12"),
440
                          xml.valueDouble("model.settings.soil.swdDBHClass23"));
475 werner 441
 
646 werner 442
    // setup of modules
443
    if (mModules)
444
        delete mModules;
445
    mModules = new Modules();
446
 
103 Werner 447
    setupSpace();
281 werner 448
    if (mRU.isEmpty())
449
        throw IException("Setup of Model: no resource units present!");
185 werner 450
 
451
    // (3) additional issues
837 werner 452
    // (3.1) load javascript code into the engine
453
    ScriptGlobal::loadScript(g->path(xml.value("system.javascript.fileName"),"script"));
454
 
455
    // (3.2) setup of regeneration
391 werner 456
    changeSettings().regenerationEnabled = xml.valueBool("model.settings.regenerationEnabled", false);
457
    if (settings().regenerationEnabled) {
458
        foreach(SpeciesSet *ss, mSpeciesSets)
459
            ss->setupRegeneration();
387 werner 460
    }
493 werner 461
    Sapling::setRecruitmentVariation(xml.valueDouble("model.settings.seedDispersal.recruitmentDimensionVariation",0.1));
468 werner 462
 
837 werner 463
    // (3.3) management
191 werner 464
    QString mgmtFile = xml.value("model.management.file");
465
    if (!mgmtFile.isEmpty() && xml.valueBool("model.management.enabled")) {
185 werner 466
        mManagement = new Management();
467
        QString path = GlobalSettings::instance()->path(mgmtFile, "script");
468
        mManagement->loadScript(path);
469
        qDebug() << "setup management using script" << path;
470
    }
468 werner 471
 
472
 
92 Werner 473
}
105 Werner 474
 
475
 
543 werner 476
ResourceUnit *Model::ru(QPointF coord)
105 Werner 477
{
106 Werner 478
    if (!mRUmap.isEmpty() && mRUmap.coordValid(coord))
479
        return mRUmap.valueAt(coord);
281 werner 480
    return ru(); // default RU if there is only one
105 Werner 481
}
482
 
286 werner 483
void Model::initOutputDatabase()
484
{
485
    GlobalSettings *g = GlobalSettings::instance();
486
    QString dbPath = g->path(g->settings().value("system.database.out"), "output");
487
    // create run-metadata
488
    int maxid = SqlHelper::queryValue("select max(id) from runs", g->dbin()).toInt();
489
 
490
    maxid++;
491
    QString timestamp = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss");
287 werner 492
    SqlHelper::executeSql(QString("insert into runs (id, timestamp) values (%1, '%2')").arg(maxid).arg(timestamp), g->dbin());
286 werner 493
    // replace path information
494
    dbPath.replace("$id$", QString::number(maxid));
495
    dbPath.replace("$date$", timestamp);
496
    // setup final path
497
   g->setupDatabaseConnection("out", dbPath, false);
498
 
499
}
500
 
824 werner 501
/// multithreaded running function for the resource unit level establishment
502
ResourceUnit *nc_sapling_growth_establishment(ResourceUnit *unit)
440 werner 503
{
632 werner 504
    try {
824 werner 505
        // define a height map for the current resource unit on the stack
632 werner 506
        float sapling_map[cPxPerRU*cPxPerRU];
824 werner 507
        // set the map and initialize it:
632 werner 508
        unit->setSaplingHeightMap(sapling_map);
444 werner 509
 
450 werner 510
 
451 werner 511
        // (1) calculate the growth of (already established) saplings (populate sapling map)
824 werner 512
        QList<ResourceUnitSpecies*>::const_iterator rus;
513
        for (rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)
514
            (*rus)->calclulateSaplingGrowth();
615 werner 515
 
824 werner 516
        // (2) calculate the establishment probabilities of new saplings
517
        for (rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)
864 werner 518
            (*rus)->calculateEstablishment();
615 werner 519
 
520
 
440 werner 521
    } catch (const IException& e) {
632 werner 522
        GlobalSettings::instance()->controller()->throwError(e.message());
440 werner 523
    }
632 werner 524
 
824 werner 525
    unit->setSaplingHeightMap(0); // invalidate again
440 werner 526
    return unit;
824 werner 527
 
440 werner 528
}
529
 
824 werner 530
 
475 werner 531
/// multithreaded execution of the carbon cycle routine
532
ResourceUnit *nc_carbonCycle(ResourceUnit *unit)
533
{
611 werner 534
    try {
632 werner 535
        // (1) do calculations on snag dynamics for the resource unit
536
        unit->calculateCarbonCycle();
537
        // (2) do the soil carbon and nitrogen dynamics calculations (ICBM/2N)
538
    } catch (const IException& e) {
539
        GlobalSettings::instance()->controller()->throwError(e.message());
611 werner 540
    }
541
 
475 werner 542
    return unit;
543
}
544
 
440 werner 545
/// beforeRun performs several steps before the models starts running.
546
/// inter alia: * setup of the stands
547
///             * setup of the climates
105 Werner 548
void Model::beforeRun()
549
{
395 werner 550
    // setup outputs
551
    // setup output database
539 werner 552
    if (GlobalSettings::instance()->dbout().isOpen())
553
        GlobalSettings::instance()->dbout().close();
395 werner 554
    initOutputDatabase();
555
    GlobalSettings::instance()->outputManager()->setup();
556
    GlobalSettings::instance()->clearDebugLists();
557
 
105 Werner 558
    // initialize stands
251 werner 559
    {
106 Werner 560
    DebugTimer loadtrees("load trees");
561
    StandLoader loader(this);
562
    loader.processInit();
251 werner 563
    }
106 Werner 564
 
214 werner 565
    // load climate
251 werner 566
    {
734 werner 567
    if (logLevelDebug()) qDebug() << "attempting to load climate..." ;
251 werner 568
    DebugTimer loadclim("load climate");
214 werner 569
    foreach(Climate *c, mClimates)
280 werner 570
        if (!c->isSetup())
571
            c->setup();
214 werner 572
 
251 werner 573
    }
574
 
575
    { DebugTimer loadinit("load standstatistics");
734 werner 576
    if (logLevelDebug()) qDebug() << "attempting to calculate initial stand statistics (incl. apply and read pattern)..." ;
135 Werner 577
    Tree::setGrid(mGrid, mHeightGrid);
737 werner 578
    // debugCheckAllTrees(); // introduced for debugging session (2012-04-06)
135 Werner 579
    applyPattern();
580
    readPattern();
106 Werner 581
 
376 werner 582
    // force the compilation of initial stand statistics
241 werner 583
    createStandStatistics();
251 werner 584
    }
286 werner 585
 
264 werner 586
    // outputs to create with inital state (without any growth) are called here:
257 werner 587
    GlobalSettings::instance()->outputManager()->execute("stand"); // year=0
837 werner 588
    GlobalSettings::instance()->outputManager()->execute("landscape"); // year=0
504 werner 589
    GlobalSettings::instance()->outputManager()->execute("sapling"); // year=0
264 werner 590
    GlobalSettings::instance()->outputManager()->execute("tree"); // year=0
395 werner 591
    GlobalSettings::instance()->outputManager()->execute("dynamicstand"); // year=0
264 werner 592
 
257 werner 593
    GlobalSettings::instance()->setCurrentYear(1); // set to first year
105 Werner 594
}
595
 
331 werner 596
/** Main model runner.
597
  The sequence of actions is as follows:
598
  (1) Load the climate of the new year
599
  (2) Reset statistics for resource unit as well as for dead/managed trees
714 werner 600
  (3) Invoke Management.
331 werner 601
  (4) *after* that, calculate Light patterns
602
  (5) 3PG on stand level, tree growth. Clear stand-statistcs before they are filled by single-tree-growth. calculate water cycle (with LAIs before management)
440 werner 603
  (6) execute Regeneration
714 werner 604
  (7) invoke disturbance modules
605
  (8) calculate statistics for the year
606
  (9) write database outputs
331 werner 607
  */
105 Werner 608
void Model::runYear()
609
{
223 werner 610
    DebugTimer t("Model::runYear()");
615 werner 611
    GlobalSettings::instance()->systemStatistics()->reset();
707 werner 612
    RandomGenerator::checkGenerator(); // see if we need to generate new numbers...
649 werner 613
    // initalization at start of year for external modules
614
    mModules->yearBegin();
341 werner 615
    // execute scheduled events for the current year
616
    if (mTimeEvents)
617
        mTimeEvents->run();
618
 
619
    // load the next year of the climate database
214 werner 620
    foreach(Climate *c, mClimates)
621
        c->nextYear();
622
 
278 werner 623
    // reset statistics
624
    foreach(ResourceUnit *ru, mRU)
625
        ru->newYear();
626
 
391 werner 627
    foreach(SpeciesSet *set, mSpeciesSets)
628
        set->newYear();
278 werner 629
    // management
615 werner 630
    if (mManagement) {
631
        DebugTimer t("management");
278 werner 632
        mManagement->run();
615 werner 633
        GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed();
634
    }
278 werner 635
 
214 werner 636
    // process a cycle of individual growth
277 werner 637
    applyPattern(); // create Light Influence Patterns
638
    readPattern(); // readout light state of individual trees
639
    grow(); // let the trees grow (growth on stand-level, tree-level, mortality)
106 Werner 640
 
391 werner 641
    // regeneration
642
    if (settings().regenerationEnabled) {
440 werner 643
        // seed dispersal
483 werner 644
        DebugTimer tseed("Regeneration and Establishment");
391 werner 645
        foreach(SpeciesSet *set, mSpeciesSets)
475 werner 646
            set->regeneration(); // parallel execution for each species set
440 werner 647
 
615 werner 648
        GlobalSettings::instance()->systemStatistics()->tSeedDistribution+=tseed.elapsed();
440 werner 649
        // establishment
824 werner 650
        { DebugTimer t("saplingGrowthEstablishment");
651
        executePerResourceUnit( nc_sapling_growth_establishment, false /* true: force single thraeded operation */);
864 werner 652
        GlobalSettings::instance()->systemStatistics()->tSaplingAndEstablishment+=t.elapsed();
615 werner 653
        }
824 werner 654
 
615 werner 655
    }
440 werner 656
 
468 werner 657
    // calculate soil / snag dynamics
475 werner 658
    if (settings().carbonCycleEnabled) {
659
        DebugTimer ccycle("carbon cylce");
660
        executePerResourceUnit( nc_carbonCycle, false /* true: force single thraeded operation */);
615 werner 661
        GlobalSettings::instance()->systemStatistics()->tCarbonCycle+=ccycle.elapsed();
662
 
475 werner 663
    }
649 werner 664
 
665
    // external modules/disturbances
666
    mModules->run();
667
 
664 werner 668
    foreach(ResourceUnit *ru, mRU) {
720 werner 669
        // remove trees that died during disturbances.
670
        if (ru->hasDiedTrees()) {
671
            ru->cleanTreeList(); // clean up the died trees
672
            ru->recreateStandStatistics(); // re-evaluate the stand statistics for this resource unit
673
        }
664 werner 674
    }
675
 
615 werner 676
    DebugTimer toutput("outputs");
278 werner 677
    // calculate statistics
678
    foreach(ResourceUnit *ru, mRU)
679
        ru->yearEnd();
124 Werner 680
 
277 werner 681
    // create outputs
184 werner 682
    OutputManager *om = GlobalSettings::instance()->outputManager();
285 werner 683
    om->execute("tree"); // single tree output
278 werner 684
    om->execute("stand"); //resource unit level x species
837 werner 685
    om->execute("landscape"); //landscape x species
504 werner 686
    om->execute("sapling"); // sapling layer per RU x species
278 werner 687
    om->execute("production_month"); // 3pg responses growth per species x RU x month
285 werner 688
    om->execute("dynamicstand"); // output with user-defined columns (based on species x RU)
278 werner 689
    om->execute("standdead"); // resource unit level x species
690
    om->execute("management"); // resource unit level x species
587 werner 691
    om->execute("carbon"); // resource unit level, carbon pools above and belowground
609 werner 692
    om->execute("carbonflow"); // resource unit level, GPP, NPP and total carbon flows (atmosphere, harvest, ...)
185 werner 693
 
615 werner 694
    GlobalSettings::instance()->systemStatistics()->tWriteOutput+=toutput.elapsed();
695
    GlobalSettings::instance()->systemStatistics()->tTotalYear+=t.elapsed();
696
    GlobalSettings::instance()->systemStatistics()->writeOutput();
697
 
162 werner 698
    GlobalSettings::instance()->setCurrentYear(GlobalSettings::instance()->currentYear()+1);
105 Werner 699
}
700
 
278 werner 701
 
702
 
105 Werner 703
void Model::afterStop()
704
{
705
    // do some cleanup
706
}
106 Werner 707
 
369 werner 708
/// multithreaded running function for LIP printing
187 iland 709
ResourceUnit* nc_applyPattern(ResourceUnit *unit)
106 Werner 710
{
711
 
712
    QVector<Tree>::iterator tit;
118 Werner 713
    QVector<Tree>::iterator tend = unit->trees().end();
107 Werner 714
 
632 werner 715
    try {
107 Werner 716
 
632 werner 717
        // light concurrence influence
718
        if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
719
            // height dominance grid
720
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
721
                (*tit).heightGrid(); // just do it ;)
155 werner 722
 
632 werner 723
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
724
                (*tit).applyLIP(); // just do it ;)
155 werner 725
 
632 werner 726
        } else {
727
            // height dominance grid
728
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
729
                (*tit).heightGrid_torus(); // just do it ;)
155 werner 730
 
632 werner 731
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
732
                (*tit).applyLIP_torus(); // do it the wraparound way
733
        }
734
        return unit;
735
    } catch (const IException &e) {
736
        GlobalSettings::instance()->controller()->throwError(e.message());
118 Werner 737
    }
738
    return unit;
739
}
106 Werner 740
 
369 werner 741
/// multithreaded running function for LIP value extraction
187 iland 742
ResourceUnit *nc_readPattern(ResourceUnit *unit)
118 Werner 743
{
744
    QVector<Tree>::iterator tit;
745
    QVector<Tree>::iterator  tend = unit->trees().end();
632 werner 746
    try {
747
        if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
748
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
749
                (*tit).readLIF(); // multipliactive approach
750
        } else {
751
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
752
                (*tit).readLIF_torus(); // do it the wraparound way
753
        }
754
    } catch (const IException &e) {
755
        GlobalSettings::instance()->controller()->throwError(e.message());
118 Werner 756
    }
757
    return unit;
106 Werner 758
}
759
 
369 werner 760
/// multithreaded running function for the growth of individual trees
187 iland 761
ResourceUnit *nc_grow(ResourceUnit *unit)
118 Werner 762
{
763
    QVector<Tree>::iterator tit;
764
    QVector<Tree>::iterator  tend = unit->trees().end();
632 werner 765
    try {
766
        unit->beforeGrow(); // reset statistics
767
        // calculate light responses
768
        // responses are based on *modified* values for LightResourceIndex
769
        for (tit=unit->trees().begin(); tit!=tend; ++tit) {
770
            (*tit).calcLightResponse();
771
        }
118 Werner 772
 
632 werner 773
        unit->calculateInterceptedArea();
251 werner 774
 
632 werner 775
        for (tit=unit->trees().begin(); tit!=tend; ++tit) {
776
            (*tit).grow(); // actual growth of individual trees
777
        }
778
    } catch (const IException &e) {
779
        GlobalSettings::instance()->controller()->throwError(e.message());
118 Werner 780
    }
615 werner 781
 
782
    GlobalSettings::instance()->systemStatistics()->treeCount+=unit->trees().count();
118 Werner 783
    return unit;
784
}
785
 
369 werner 786
/// multithreaded running function for the resource level production
787
ResourceUnit *nc_production(ResourceUnit *unit)
788
{
632 werner 789
    try {
790
        unit->production();
791
    } catch (const IException &e) {
792
        GlobalSettings::instance()->controller()->throwError(e.message());
793
    }
369 werner 794
    return unit;
795
}
796
 
440 werner 797
 
124 Werner 798
void Model::test()
799
{
800
    // Test-funktion: braucht 1/3 time von readGrid()
801
    DebugTimer t("test");
802
    FloatGrid averaged = mGrid->averaged(10);
803
    int count = 0;
804
    float *end = averaged.end();
805
    for (float *p=averaged.begin(); p!=end; ++p)
806
        if (*p > 0.9)
807
            count++;
808
    qDebug() << count << "LIF>0.9 of " << averaged.count();
809
}
810
 
734 werner 811
void Model::debugCheckAllTrees()
812
{
813
    AllTreeIterator at(this);
814
    bool has_errors = false; double dummy=0.;
815
    while (Tree *t = at.next()) {
816
        // plausibility
817
        if (t->dbh()<0 || t->dbh()>10000. || t->biomassFoliage()<0. || t->height()>1000. || t->height() < 0.
818
                || t->biomassFoliage() <0.)
819
            has_errors = true;
820
        // check for objects....
821
        dummy = t->stamp()->offset() + t->ru()->ruSpecies()[1]->statistics().count();
822
    }
823
    if (has_errors)
824
        qDebug() << "model: debugCheckAllTrees found problems" << dummy;
825
}
826
 
118 Werner 827
void Model::applyPattern()
828
{
829
 
830
    DebugTimer t("applyPattern()");
831
    // intialize grids...
720 werner 832
    initializeGrid();
551 werner 833
 
406 werner 834
    // initialize height grid with a value of 4m. This is the height of the regeneration layer
551 werner 835
    for (HeightGridValue *h=mHeightGrid->begin();h!=mHeightGrid->end();++h) {
836
        h->resetCount(); // set count = 0, but do not touch the flags
837
        h->height = 4.f;
838
    }
118 Werner 839
 
123 Werner 840
    threadRunner.run(nc_applyPattern);
615 werner 841
    GlobalSettings::instance()->systemStatistics()->tApplyPattern+=t.elapsed();
118 Werner 842
}
843
 
106 Werner 844
void Model::readPattern()
845
{
846
    DebugTimer t("readPattern()");
123 Werner 847
    threadRunner.run(nc_readPattern);
615 werner 848
    GlobalSettings::instance()->systemStatistics()->tReadPattern+=t.elapsed();
849
 
106 Werner 850
}
851
 
331 werner 852
/** Main function for the growth of stands and trees.
853
   This includes several steps.
854
   (1) calculate the stocked area (i.e. count pixels in height grid)
855
   (2) 3PG production (including response calculation, water cycle)
856
   (3) single tree growth (including mortality)
857
   (4) cleanup of tree lists (remove dead trees)
858
  */
106 Werner 859
void Model::grow()
860
{
151 iland 861
 
369 werner 862
    if (!settings().growthEnabled)
863
        return;
864
    { DebugTimer t("growRU()");
865
    calculateStockedArea();
113 Werner 866
 
374 werner 867
    // multithreaded: mutex for the message handler in mainwindow solved the crashes.
868
    threadRunner.run(nc_production);
370 werner 869
    }
870
 
871
    DebugTimer t("growTrees()");
251 werner 872
    threadRunner.run(nc_grow); // actual growth of individual trees
159 werner 873
 
187 iland 874
    foreach(ResourceUnit *ru, mRU) {
159 werner 875
       ru->cleanTreeList();
376 werner 876
       ru->afterGrow();
168 werner 877
       //qDebug() << (b-n) << "trees died (of" << b << ").";
159 werner 878
   }
615 werner 879
   GlobalSettings::instance()->systemStatistics()->tTreeGrowth+=t.elapsed();
123 Werner 880
}
151 iland 881
 
240 werner 882
/** calculate for each resource unit the fraction of area which is stocked.
883
  This is done by checking the pixels of the global height grid.
151 iland 884
  */
885
void Model::calculateStockedArea()
886
{
887
    // iterate over the whole heightgrid and count pixels for each ressource unit
888
    HeightGridValue *end = mHeightGrid->end();
889
    QPointF cp;
187 iland 890
    ResourceUnit *ru;
151 iland 891
    for (HeightGridValue *i=mHeightGrid->begin(); i!=end; ++i) {
892
        cp = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i));
893
        if (mRUmap.coordValid(cp)) {
894
            ru = mRUmap.valueAt(cp);
895
            if (ru) {
285 werner 896
                ru->countStockedPixel( (*i).count()>0 );
151 iland 897
            }
898
        }
899
 
900
    }
901
}
240 werner 902
 
664 werner 903
/** calculate for each resource unit the stockable area.
574 werner 904
  "stockability" is determined by the isValid flag of resource units which in turn
905
  is derived from stand grid values.
906
  */
907
void Model::calculateStockableArea()
908
{
909
 
910
    foreach(ResourceUnit *ru, mRU) {
720 werner 911
        // //
912
        //        if (ru->id()==-1) {
913
        //            ru->setStockableArea(0.);
914
        //            continue;
915
        //        }
574 werner 916
        GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox());
917
        int valid=0, total=0;
918
        while (runner.next()) {
919
            if ( runner.current()->isValid() )
920
                valid++;
921
            total++;
922
        }
575 werner 923
        if (total) {
574 werner 924
            ru->setStockableArea( cHeightPixelArea * valid);
734 werner 925
            if (valid==0 && ru->id()>-1) {
926
                // invalidate this resource unit
927
                ru->setID(-1);
928
            }
575 werner 929
            if (valid>0 && ru->id()==-1) {
930
                qDebug() << "Warning: a resource unit has id=-1 but stockable area (id was set to 0)!!! ru: " << ru->boundingBox() << "with index" << ru->index();
931
                ru->setID(0);
664 werner 932
                // test-code
933
                //GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox());
934
                //while (runner.next()) {
935
                //    qDebug() << mHeightGrid->cellCenterPoint(mHeightGrid->indexOf( runner.current() )) << ": " << runner.current()->isValid();
936
                //}
585 werner 937
 
575 werner 938
            }
939
        } else
574 werner 940
            throw IException("calculateStockableArea: resource unit without pixels!");
941
 
942
    }
720 werner 943
    // mark those pixels that are at the edge of a "forest-out-of-area"
944
    GridRunner<HeightGridValue> runner(*mHeightGrid, mHeightGrid->metricRect());
945
    HeightGridValue* neighbors[8];
946
    while (runner.next()) {
947
        if (runner.current()->isForestOutside()) {
948
            // if the current pixel is a "radiating" border pixel,
949
            // then check the neighbors and set a flag if the pixel is a neighbor of a in-project-area pixel.
950
            runner.neighbors8(neighbors);
951
            for (int i=0;i<8;++i)
952
                if (neighbors[i] &&  neighbors[i]->isValid())
953
                    runner.current()->setIsRadiating();
954
 
955
        }
956
    }
957
 
574 werner 958
}
959
 
720 werner 960
void Model::initializeGrid()
961
{
962
    // fill the whole grid with a value of "1."
963
    mGrid->initialize(1.f);
574 werner 964
 
720 werner 965
    // apply special values for grid cells border regions where out-of-area cells
966
    // radiate into the main LIF grid.
967
    QPoint p;
968
    int ix_min, ix_max, iy_min, iy_max, ix_center, iy_center;
969
    const int px_offset = cPxPerHeight / 2; // for 5 px per height grid cell, the offset is 2
970
    const int max_radiate_distance = 7;
971
    const float step_width = 1.f / (float)max_radiate_distance;
972
    int c_rad = 0;
973
    for (HeightGridValue *hgv=mHeightGrid->begin(); hgv!=mHeightGrid->end(); ++hgv) {
974
        if (hgv->isRadiating()) {
975
            p=mHeightGrid->indexOf(hgv);
976
            ix_min = p.x() * cPxPerHeight - max_radiate_distance + px_offset;
977
            ix_max = ix_min + 2*max_radiate_distance + 1;
978
            ix_center = ix_min + max_radiate_distance;
979
            iy_min = p.y() * cPxPerHeight - max_radiate_distance + px_offset;
980
            iy_max = iy_min + 2*max_radiate_distance + 1;
981
            iy_center = iy_min + max_radiate_distance;
982
            for (int y=iy_min; y<=iy_max; ++y) {
983
                for (int x=ix_min; x<=ix_max; ++x) {
802 werner 984
                    if (!mGrid->isIndexValid(x,y) ||  !(*mHeightGrid)(x/cPxPerHeight, y/cPxPerHeight).isValid())
720 werner 985
                        continue;
986
                    float value = qMax(qAbs(x-ix_center), qAbs(y-iy_center)) * step_width;
987
                    float &v = mGrid->valueAtIndex(x, y);
988
                    if (value>=0.f && v>value)
989
                        v = value;
990
                }
991
            }
992
            c_rad++;
993
        }
994
    }
721 werner 995
    if (logLevelDebug())
996
        qDebug() << "initialize grid:" << c_rad << "radiating pixels...";
720 werner 997
 
998
}
999
 
1000
 
241 werner 1001
/// Force the creation of stand statistics.
1002
/// - stocked area (for resourceunit-areas)
1003
/// - ru - statistics
240 werner 1004
void Model::createStandStatistics()
1005
{
241 werner 1006
    calculateStockedArea();
482 werner 1007
    foreach(ResourceUnit *ru, mRU) {
1008
        ru->addTreeAgingForAllTrees();
240 werner 1009
        ru->createStandStatistics();
482 werner 1010
    }
240 werner 1011
}
584 werner 1012
 
1013