Subversion Repositories public iLand

Rev

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