Subversion Repositories public iLand

Rev

Rev 893 | Rev 902 | 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;
890 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));
272
            mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used.
889 werner 273
            if (!mStandGrid || !mStandGrid->isValid() || *p>0) {
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;
890 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;
890 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
477
        mABEManagement = new AMIE::ForestManagementEngine();
478
        // setup of ABE after loading of trees.
479
 
480
    } else {
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
        }
185 werner 489
    }
468 werner 490
 
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
 
824 werner 520
/// multithreaded running function for the resource unit level establishment
521
ResourceUnit *nc_sapling_growth_establishment(ResourceUnit *unit)
440 werner 522
{
632 werner 523
    try {
824 werner 524
        // define a height map for the current resource unit on the stack
632 werner 525
        float sapling_map[cPxPerRU*cPxPerRU];
824 werner 526
        // set the map and initialize it:
632 werner 527
        unit->setSaplingHeightMap(sapling_map);
444 werner 528
 
450 werner 529
 
451 werner 530
        // (1) calculate the growth of (already established) saplings (populate sapling map)
824 werner 531
        QList<ResourceUnitSpecies*>::const_iterator rus;
532
        for (rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)
533
            (*rus)->calclulateSaplingGrowth();
615 werner 534
 
824 werner 535
        // (2) calculate the establishment probabilities of new saplings
536
        for (rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)
864 werner 537
            (*rus)->calculateEstablishment();
615 werner 538
 
539
 
440 werner 540
    } catch (const IException& e) {
632 werner 541
        GlobalSettings::instance()->controller()->throwError(e.message());
440 werner 542
    }
632 werner 543
 
824 werner 544
    unit->setSaplingHeightMap(0); // invalidate again
440 werner 545
    return unit;
824 werner 546
 
440 werner 547
}
548
 
824 werner 549
 
475 werner 550
/// multithreaded execution of the carbon cycle routine
551
ResourceUnit *nc_carbonCycle(ResourceUnit *unit)
552
{
611 werner 553
    try {
632 werner 554
        // (1) do calculations on snag dynamics for the resource unit
555
        unit->calculateCarbonCycle();
556
        // (2) do the soil carbon and nitrogen dynamics calculations (ICBM/2N)
557
    } catch (const IException& e) {
558
        GlobalSettings::instance()->controller()->throwError(e.message());
611 werner 559
    }
560
 
475 werner 561
    return unit;
562
}
563
 
440 werner 564
/// beforeRun performs several steps before the models starts running.
565
/// inter alia: * setup of the stands
566
///             * setup of the climates
105 Werner 567
void Model::beforeRun()
568
{
395 werner 569
    // setup outputs
570
    // setup output database
539 werner 571
    if (GlobalSettings::instance()->dbout().isOpen())
572
        GlobalSettings::instance()->dbout().close();
395 werner 573
    initOutputDatabase();
574
    GlobalSettings::instance()->outputManager()->setup();
575
    GlobalSettings::instance()->clearDebugLists();
576
 
105 Werner 577
    // initialize stands
251 werner 578
    {
106 Werner 579
    DebugTimer loadtrees("load trees");
580
    StandLoader loader(this);
581
    loader.processInit();
251 werner 582
    }
106 Werner 583
 
214 werner 584
    // load climate
251 werner 585
    {
734 werner 586
    if (logLevelDebug()) qDebug() << "attempting to load climate..." ;
251 werner 587
    DebugTimer loadclim("load climate");
214 werner 588
    foreach(Climate *c, mClimates)
280 werner 589
        if (!c->isSetup())
590
            c->setup();
214 werner 591
 
251 werner 592
    }
593
 
594
    { DebugTimer loadinit("load standstatistics");
734 werner 595
    if (logLevelDebug()) qDebug() << "attempting to calculate initial stand statistics (incl. apply and read pattern)..." ;
135 Werner 596
    Tree::setGrid(mGrid, mHeightGrid);
737 werner 597
    // debugCheckAllTrees(); // introduced for debugging session (2012-04-06)
135 Werner 598
    applyPattern();
599
    readPattern();
106 Werner 600
 
376 werner 601
    // force the compilation of initial stand statistics
241 werner 602
    createStandStatistics();
251 werner 603
    }
286 werner 604
 
890 werner 605
    // initalization of ABE
606
    if (mABEManagement) {
607
        mABEManagement->setup();
608
    }
609
 
610
 
264 werner 611
    // outputs to create with inital state (without any growth) are called here:
257 werner 612
    GlobalSettings::instance()->outputManager()->execute("stand"); // year=0
837 werner 613
    GlobalSettings::instance()->outputManager()->execute("landscape"); // year=0
504 werner 614
    GlobalSettings::instance()->outputManager()->execute("sapling"); // year=0
264 werner 615
    GlobalSettings::instance()->outputManager()->execute("tree"); // year=0
395 werner 616
    GlobalSettings::instance()->outputManager()->execute("dynamicstand"); // year=0
264 werner 617
 
257 werner 618
    GlobalSettings::instance()->setCurrentYear(1); // set to first year
105 Werner 619
}
620
 
331 werner 621
/** Main model runner.
622
  The sequence of actions is as follows:
623
  (1) Load the climate of the new year
624
  (2) Reset statistics for resource unit as well as for dead/managed trees
714 werner 625
  (3) Invoke Management.
331 werner 626
  (4) *after* that, calculate Light patterns
627
  (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 628
  (6) execute Regeneration
714 werner 629
  (7) invoke disturbance modules
630
  (8) calculate statistics for the year
631
  (9) write database outputs
331 werner 632
  */
105 Werner 633
void Model::runYear()
634
{
223 werner 635
    DebugTimer t("Model::runYear()");
615 werner 636
    GlobalSettings::instance()->systemStatistics()->reset();
707 werner 637
    RandomGenerator::checkGenerator(); // see if we need to generate new numbers...
649 werner 638
    // initalization at start of year for external modules
639
    mModules->yearBegin();
341 werner 640
    // execute scheduled events for the current year
641
    if (mTimeEvents)
642
        mTimeEvents->run();
643
 
644
    // load the next year of the climate database
214 werner 645
    foreach(Climate *c, mClimates)
646
        c->nextYear();
647
 
278 werner 648
    // reset statistics
649
    foreach(ResourceUnit *ru, mRU)
650
        ru->newYear();
651
 
391 werner 652
    foreach(SpeciesSet *set, mSpeciesSets)
653
        set->newYear();
890 werner 654
    // management classic
615 werner 655
    if (mManagement) {
656
        DebugTimer t("management");
278 werner 657
        mManagement->run();
615 werner 658
        GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed();
659
    }
890 werner 660
    // ... or ABE (the agent based variant)
661
    if (mABEManagement) {
901 werner 662
        DebugTimer t("ABE:run");
890 werner 663
        mABEManagement->run();
664
        GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed();
665
    }
278 werner 666
 
214 werner 667
    // process a cycle of individual growth
277 werner 668
    applyPattern(); // create Light Influence Patterns
669
    readPattern(); // readout light state of individual trees
670
    grow(); // let the trees grow (growth on stand-level, tree-level, mortality)
106 Werner 671
 
391 werner 672
    // regeneration
673
    if (settings().regenerationEnabled) {
440 werner 674
        // seed dispersal
483 werner 675
        DebugTimer tseed("Regeneration and Establishment");
391 werner 676
        foreach(SpeciesSet *set, mSpeciesSets)
475 werner 677
            set->regeneration(); // parallel execution for each species set
440 werner 678
 
615 werner 679
        GlobalSettings::instance()->systemStatistics()->tSeedDistribution+=tseed.elapsed();
440 werner 680
        // establishment
824 werner 681
        { DebugTimer t("saplingGrowthEstablishment");
682
        executePerResourceUnit( nc_sapling_growth_establishment, false /* true: force single thraeded operation */);
864 werner 683
        GlobalSettings::instance()->systemStatistics()->tSaplingAndEstablishment+=t.elapsed();
615 werner 684
        }
824 werner 685
 
615 werner 686
    }
440 werner 687
 
468 werner 688
    // calculate soil / snag dynamics
475 werner 689
    if (settings().carbonCycleEnabled) {
690
        DebugTimer ccycle("carbon cylce");
691
        executePerResourceUnit( nc_carbonCycle, false /* true: force single thraeded operation */);
615 werner 692
        GlobalSettings::instance()->systemStatistics()->tCarbonCycle+=ccycle.elapsed();
693
 
475 werner 694
    }
649 werner 695
 
696
    // external modules/disturbances
697
    mModules->run();
698
 
664 werner 699
    foreach(ResourceUnit *ru, mRU) {
720 werner 700
        // remove trees that died during disturbances.
701
        if (ru->hasDiedTrees()) {
702
            ru->cleanTreeList(); // clean up the died trees
703
            ru->recreateStandStatistics(); // re-evaluate the stand statistics for this resource unit
704
        }
664 werner 705
    }
706
 
615 werner 707
    DebugTimer toutput("outputs");
278 werner 708
    // calculate statistics
709
    foreach(ResourceUnit *ru, mRU)
710
        ru->yearEnd();
124 Werner 711
 
277 werner 712
    // create outputs
184 werner 713
    OutputManager *om = GlobalSettings::instance()->outputManager();
285 werner 714
    om->execute("tree"); // single tree output
278 werner 715
    om->execute("stand"); //resource unit level x species
837 werner 716
    om->execute("landscape"); //landscape x species
504 werner 717
    om->execute("sapling"); // sapling layer per RU x species
278 werner 718
    om->execute("production_month"); // 3pg responses growth per species x RU x month
285 werner 719
    om->execute("dynamicstand"); // output with user-defined columns (based on species x RU)
278 werner 720
    om->execute("standdead"); // resource unit level x species
721
    om->execute("management"); // resource unit level x species
587 werner 722
    om->execute("carbon"); // resource unit level, carbon pools above and belowground
609 werner 723
    om->execute("carbonflow"); // resource unit level, GPP, NPP and total carbon flows (atmosphere, harvest, ...)
185 werner 724
 
615 werner 725
    GlobalSettings::instance()->systemStatistics()->tWriteOutput+=toutput.elapsed();
726
    GlobalSettings::instance()->systemStatistics()->tTotalYear+=t.elapsed();
727
    GlobalSettings::instance()->systemStatistics()->writeOutput();
728
 
162 werner 729
    GlobalSettings::instance()->setCurrentYear(GlobalSettings::instance()->currentYear()+1);
105 Werner 730
}
731
 
278 werner 732
 
733
 
105 Werner 734
void Model::afterStop()
735
{
736
    // do some cleanup
737
}
106 Werner 738
 
369 werner 739
/// multithreaded running function for LIP printing
187 iland 740
ResourceUnit* nc_applyPattern(ResourceUnit *unit)
106 Werner 741
{
742
 
743
    QVector<Tree>::iterator tit;
118 Werner 744
    QVector<Tree>::iterator tend = unit->trees().end();
107 Werner 745
 
632 werner 746
    try {
107 Werner 747
 
632 werner 748
        // light concurrence influence
749
        if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
750
            // height dominance grid
751
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
752
                (*tit).heightGrid(); // just do it ;)
155 werner 753
 
632 werner 754
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
755
                (*tit).applyLIP(); // just do it ;)
155 werner 756
 
632 werner 757
        } else {
758
            // height dominance grid
759
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
760
                (*tit).heightGrid_torus(); // just do it ;)
155 werner 761
 
632 werner 762
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
763
                (*tit).applyLIP_torus(); // do it the wraparound way
764
        }
765
        return unit;
766
    } catch (const IException &e) {
767
        GlobalSettings::instance()->controller()->throwError(e.message());
118 Werner 768
    }
769
    return unit;
770
}
106 Werner 771
 
369 werner 772
/// multithreaded running function for LIP value extraction
187 iland 773
ResourceUnit *nc_readPattern(ResourceUnit *unit)
118 Werner 774
{
775
    QVector<Tree>::iterator tit;
776
    QVector<Tree>::iterator  tend = unit->trees().end();
632 werner 777
    try {
778
        if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
779
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
780
                (*tit).readLIF(); // multipliactive approach
781
        } else {
782
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
783
                (*tit).readLIF_torus(); // do it the wraparound way
784
        }
785
    } catch (const IException &e) {
786
        GlobalSettings::instance()->controller()->throwError(e.message());
118 Werner 787
    }
788
    return unit;
106 Werner 789
}
790
 
369 werner 791
/// multithreaded running function for the growth of individual trees
187 iland 792
ResourceUnit *nc_grow(ResourceUnit *unit)
118 Werner 793
{
794
    QVector<Tree>::iterator tit;
795
    QVector<Tree>::iterator  tend = unit->trees().end();
632 werner 796
    try {
797
        unit->beforeGrow(); // reset statistics
798
        // calculate light responses
799
        // responses are based on *modified* values for LightResourceIndex
800
        for (tit=unit->trees().begin(); tit!=tend; ++tit) {
801
            (*tit).calcLightResponse();
802
        }
118 Werner 803
 
632 werner 804
        unit->calculateInterceptedArea();
251 werner 805
 
632 werner 806
        for (tit=unit->trees().begin(); tit!=tend; ++tit) {
807
            (*tit).grow(); // actual growth of individual trees
808
        }
809
    } catch (const IException &e) {
810
        GlobalSettings::instance()->controller()->throwError(e.message());
118 Werner 811
    }
615 werner 812
 
813
    GlobalSettings::instance()->systemStatistics()->treeCount+=unit->trees().count();
118 Werner 814
    return unit;
815
}
816
 
369 werner 817
/// multithreaded running function for the resource level production
818
ResourceUnit *nc_production(ResourceUnit *unit)
819
{
632 werner 820
    try {
821
        unit->production();
822
    } catch (const IException &e) {
823
        GlobalSettings::instance()->controller()->throwError(e.message());
824
    }
369 werner 825
    return unit;
826
}
827
 
440 werner 828
 
124 Werner 829
void Model::test()
830
{
831
    // Test-funktion: braucht 1/3 time von readGrid()
832
    DebugTimer t("test");
833
    FloatGrid averaged = mGrid->averaged(10);
834
    int count = 0;
835
    float *end = averaged.end();
836
    for (float *p=averaged.begin(); p!=end; ++p)
837
        if (*p > 0.9)
838
            count++;
839
    qDebug() << count << "LIF>0.9 of " << averaged.count();
840
}
841
 
734 werner 842
void Model::debugCheckAllTrees()
843
{
844
    AllTreeIterator at(this);
845
    bool has_errors = false; double dummy=0.;
846
    while (Tree *t = at.next()) {
847
        // plausibility
848
        if (t->dbh()<0 || t->dbh()>10000. || t->biomassFoliage()<0. || t->height()>1000. || t->height() < 0.
849
                || t->biomassFoliage() <0.)
850
            has_errors = true;
851
        // check for objects....
852
        dummy = t->stamp()->offset() + t->ru()->ruSpecies()[1]->statistics().count();
853
    }
854
    if (has_errors)
855
        qDebug() << "model: debugCheckAllTrees found problems" << dummy;
856
}
857
 
118 Werner 858
void Model::applyPattern()
859
{
860
 
861
    DebugTimer t("applyPattern()");
862
    // intialize grids...
720 werner 863
    initializeGrid();
551 werner 864
 
406 werner 865
    // initialize height grid with a value of 4m. This is the height of the regeneration layer
551 werner 866
    for (HeightGridValue *h=mHeightGrid->begin();h!=mHeightGrid->end();++h) {
867
        h->resetCount(); // set count = 0, but do not touch the flags
868
        h->height = 4.f;
869
    }
118 Werner 870
 
123 Werner 871
    threadRunner.run(nc_applyPattern);
615 werner 872
    GlobalSettings::instance()->systemStatistics()->tApplyPattern+=t.elapsed();
118 Werner 873
}
874
 
106 Werner 875
void Model::readPattern()
876
{
877
    DebugTimer t("readPattern()");
123 Werner 878
    threadRunner.run(nc_readPattern);
615 werner 879
    GlobalSettings::instance()->systemStatistics()->tReadPattern+=t.elapsed();
880
 
106 Werner 881
}
882
 
331 werner 883
/** Main function for the growth of stands and trees.
884
   This includes several steps.
885
   (1) calculate the stocked area (i.e. count pixels in height grid)
886
   (2) 3PG production (including response calculation, water cycle)
887
   (3) single tree growth (including mortality)
888
   (4) cleanup of tree lists (remove dead trees)
889
  */
106 Werner 890
void Model::grow()
891
{
151 iland 892
 
369 werner 893
    if (!settings().growthEnabled)
894
        return;
895
    { DebugTimer t("growRU()");
896
    calculateStockedArea();
113 Werner 897
 
374 werner 898
    // multithreaded: mutex for the message handler in mainwindow solved the crashes.
899
    threadRunner.run(nc_production);
370 werner 900
    }
901
 
902
    DebugTimer t("growTrees()");
251 werner 903
    threadRunner.run(nc_grow); // actual growth of individual trees
159 werner 904
 
187 iland 905
    foreach(ResourceUnit *ru, mRU) {
159 werner 906
       ru->cleanTreeList();
376 werner 907
       ru->afterGrow();
168 werner 908
       //qDebug() << (b-n) << "trees died (of" << b << ").";
159 werner 909
   }
615 werner 910
   GlobalSettings::instance()->systemStatistics()->tTreeGrowth+=t.elapsed();
123 Werner 911
}
151 iland 912
 
240 werner 913
/** calculate for each resource unit the fraction of area which is stocked.
914
  This is done by checking the pixels of the global height grid.
151 iland 915
  */
916
void Model::calculateStockedArea()
917
{
918
    // iterate over the whole heightgrid and count pixels for each ressource unit
919
    HeightGridValue *end = mHeightGrid->end();
920
    QPointF cp;
187 iland 921
    ResourceUnit *ru;
151 iland 922
    for (HeightGridValue *i=mHeightGrid->begin(); i!=end; ++i) {
923
        cp = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i));
924
        if (mRUmap.coordValid(cp)) {
925
            ru = mRUmap.valueAt(cp);
926
            if (ru) {
285 werner 927
                ru->countStockedPixel( (*i).count()>0 );
151 iland 928
            }
929
        }
930
 
931
    }
932
}
240 werner 933
 
664 werner 934
/** calculate for each resource unit the stockable area.
574 werner 935
  "stockability" is determined by the isValid flag of resource units which in turn
936
  is derived from stand grid values.
937
  */
938
void Model::calculateStockableArea()
939
{
940
 
941
    foreach(ResourceUnit *ru, mRU) {
720 werner 942
        // //
943
        //        if (ru->id()==-1) {
944
        //            ru->setStockableArea(0.);
945
        //            continue;
946
        //        }
574 werner 947
        GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox());
948
        int valid=0, total=0;
949
        while (runner.next()) {
950
            if ( runner.current()->isValid() )
951
                valid++;
952
            total++;
953
        }
575 werner 954
        if (total) {
574 werner 955
            ru->setStockableArea( cHeightPixelArea * valid);
734 werner 956
            if (valid==0 && ru->id()>-1) {
957
                // invalidate this resource unit
958
                ru->setID(-1);
959
            }
575 werner 960
            if (valid>0 && ru->id()==-1) {
961
                qDebug() << "Warning: a resource unit has id=-1 but stockable area (id was set to 0)!!! ru: " << ru->boundingBox() << "with index" << ru->index();
962
                ru->setID(0);
664 werner 963
                // test-code
964
                //GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox());
965
                //while (runner.next()) {
966
                //    qDebug() << mHeightGrid->cellCenterPoint(mHeightGrid->indexOf( runner.current() )) << ": " << runner.current()->isValid();
967
                //}
585 werner 968
 
575 werner 969
            }
970
        } else
574 werner 971
            throw IException("calculateStockableArea: resource unit without pixels!");
972
 
973
    }
720 werner 974
    // mark those pixels that are at the edge of a "forest-out-of-area"
975
    GridRunner<HeightGridValue> runner(*mHeightGrid, mHeightGrid->metricRect());
976
    HeightGridValue* neighbors[8];
977
    while (runner.next()) {
978
        if (runner.current()->isForestOutside()) {
979
            // if the current pixel is a "radiating" border pixel,
980
            // then check the neighbors and set a flag if the pixel is a neighbor of a in-project-area pixel.
981
            runner.neighbors8(neighbors);
982
            for (int i=0;i<8;++i)
983
                if (neighbors[i] &&  neighbors[i]->isValid())
984
                    runner.current()->setIsRadiating();
985
 
986
        }
987
    }
988
 
574 werner 989
}
990
 
720 werner 991
void Model::initializeGrid()
992
{
993
    // fill the whole grid with a value of "1."
994
    mGrid->initialize(1.f);
574 werner 995
 
720 werner 996
    // apply special values for grid cells border regions where out-of-area cells
997
    // radiate into the main LIF grid.
998
    QPoint p;
999
    int ix_min, ix_max, iy_min, iy_max, ix_center, iy_center;
1000
    const int px_offset = cPxPerHeight / 2; // for 5 px per height grid cell, the offset is 2
1001
    const int max_radiate_distance = 7;
1002
    const float step_width = 1.f / (float)max_radiate_distance;
1003
    int c_rad = 0;
1004
    for (HeightGridValue *hgv=mHeightGrid->begin(); hgv!=mHeightGrid->end(); ++hgv) {
1005
        if (hgv->isRadiating()) {
1006
            p=mHeightGrid->indexOf(hgv);
1007
            ix_min = p.x() * cPxPerHeight - max_radiate_distance + px_offset;
1008
            ix_max = ix_min + 2*max_radiate_distance + 1;
1009
            ix_center = ix_min + max_radiate_distance;
1010
            iy_min = p.y() * cPxPerHeight - max_radiate_distance + px_offset;
1011
            iy_max = iy_min + 2*max_radiate_distance + 1;
1012
            iy_center = iy_min + max_radiate_distance;
1013
            for (int y=iy_min; y<=iy_max; ++y) {
1014
                for (int x=ix_min; x<=ix_max; ++x) {
802 werner 1015
                    if (!mGrid->isIndexValid(x,y) ||  !(*mHeightGrid)(x/cPxPerHeight, y/cPxPerHeight).isValid())
720 werner 1016
                        continue;
1017
                    float value = qMax(qAbs(x-ix_center), qAbs(y-iy_center)) * step_width;
1018
                    float &v = mGrid->valueAtIndex(x, y);
1019
                    if (value>=0.f && v>value)
1020
                        v = value;
1021
                }
1022
            }
1023
            c_rad++;
1024
        }
1025
    }
721 werner 1026
    if (logLevelDebug())
1027
        qDebug() << "initialize grid:" << c_rad << "radiating pixels...";
720 werner 1028
 
1029
}
1030
 
1031
 
241 werner 1032
/// Force the creation of stand statistics.
1033
/// - stocked area (for resourceunit-areas)
1034
/// - ru - statistics
240 werner 1035
void Model::createStandStatistics()
1036
{
241 werner 1037
    calculateStockedArea();
482 werner 1038
    foreach(ResourceUnit *ru, mRU) {
1039
        ru->addTreeAgingForAllTrees();
240 werner 1040
        ru->createStandStatistics();
482 werner 1041
    }
240 werner 1042
}
584 werner 1043
 
1044