Subversion Repositories public iLand

Rev

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