From 3b2c4c2123e14fd319555bc388d06182a89a46a6 Mon Sep 17 00:00:00 2001 From: ftomei Date: Mon, 13 Jan 2025 17:47:39 +0100 Subject: [PATCH 01/12] update 3D --- soilFluxes3D/header/solver.h | 2 ++ soilFluxes3D/solver.cpp | 36 ++++++++++++++++++++++++++++++++++++ soilFluxes3D/water.cpp | 1 - 3 files changed, 38 insertions(+), 1 deletion(-) diff --git a/soilFluxes3D/header/solver.h b/soilFluxes3D/header/solver.h index 52195724..9c0129ea 100644 --- a/soilFluxes3D/header/solver.h +++ b/soilFluxes3D/header/solver.h @@ -13,5 +13,7 @@ bool GaussSeidelRelaxation (int myApproximation, double myResidualTolerance, int myProcess); + bool JacobiRelaxation (int approximation, double toleranceThreshold, int process); + #endif // SOLVER_H diff --git a/soilFluxes3D/solver.cpp b/soilFluxes3D/solver.cpp index dbcb9736..cc17334f 100644 --- a/soilFluxes3D/solver.cpp +++ b/soilFluxes3D/solver.cpp @@ -233,3 +233,39 @@ bool GaussSeidelRelaxation (int approximation, double toleranceThreshold, int pr return true; } + + +bool JacobiRelaxation (int approximation, double toleranceThreshold, int process) +{ + double currentNorm = 1.0; + double bestNorm = currentNorm; + + int maxIterationsNr = getMaxIterationsNr(approximation); + int iteration = 0; + + while ((currentNorm > toleranceThreshold) && (iteration < maxIterationsNr)) + { + if (process == PROCESS_HEAT) + { + currentNorm = GaussSeidelIterationHeat(); + } + if (process == PROCESS_WATER) + { + // TODO Jacobi method + + if (currentNorm < bestNorm) + { + bestNorm = currentNorm; + } + else if (currentNorm > (bestNorm * 10.0)) + { + // non-convergent system + return false; + } + } + + iteration++; + } + + return true; +} diff --git a/soilFluxes3D/water.cpp b/soilFluxes3D/water.cpp index 13bb92d2..709c590b 100644 --- a/soilFluxes3D/water.cpp +++ b/soilFluxes3D/water.cpp @@ -201,7 +201,6 @@ double redistribution(long i, TlinkedNode *link, int linkType) } - bool computeFlux(long i, int matrixIndex, TlinkedNode *link, double deltaT, unsigned long myApprox, int linkType) { if ((*link).index == NOLINK) return false; From c86a2c69029de02bd2a72eff959037461bdf280f Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Tue, 14 Jan 2025 17:03:31 +0100 Subject: [PATCH 02/12] bug fixed on saving derived variables --- pragaProject/pragaProject.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 8e68ec61..cc6e86a5 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2947,6 +2947,10 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (! loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem())) return false; } + // save also derived variables + foreach (myVar, derivedVariables) + varToSave.push_back(myVar); + // save also time aggregated variables foreach (myVar, aggrVariables) varToSave.push_back(myVar); From b1914c2ba1f820ed4c51b9ef30e78be47b12832f Mon Sep 17 00:00:00 2001 From: ftomei Date: Wed, 15 Jan 2025 18:53:22 +0100 Subject: [PATCH 03/12] update --- dbMeteoGrid/dbMeteoGrid.cpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/dbMeteoGrid/dbMeteoGrid.cpp b/dbMeteoGrid/dbMeteoGrid.cpp index 06b2d3d8..6dca62ea 100644 --- a/dbMeteoGrid/dbMeteoGrid.cpp +++ b/dbMeteoGrid/dbMeteoGrid.cpp @@ -3427,7 +3427,9 @@ bool Crit3DMeteoGridDbHandler::saveGridData(QString *errorStr, QDateTime firstTi if (lastTime.time().hour() == 0) lastDate = lastDate.addDays(-1); for (int row = 0; row < gridStructure().header().nrRows; row++) + { for (int col = 0; col < gridStructure().header().nrCols; col++) + { if (meteoGrid()->getMeteoPointActiveId(row, col, &id)) { if (! gridStructure().isFixedFields()) @@ -3441,6 +3443,8 @@ bool Crit3DMeteoGridDbHandler::saveGridData(QString *errorStr, QDateTime firstTi if (isDaily) saveCellGridDailyDataFF(errorStr, QString::fromStdString(id), row, col, firstTime.date(), lastDate, meteoSettings); } } + } + } return true; } @@ -3506,7 +3510,7 @@ bool Crit3DMeteoGridDbHandler::saveCellGridHourlyData(QString *errorStr, QString QString statement = QString("CREATE TABLE IF NOT EXISTS `%1` " "(`%2` datetime, VariableCode tinyint(3) UNSIGNED, Value float(6,1), PRIMARY KEY(`%2`,VariableCode))").arg(tableH, _tableHourly.fieldTime); - if( !qry.exec(statement) ) + if(! qry.exec(statement) ) { *errorStr = qry.lastError().text(); return false; @@ -3516,6 +3520,7 @@ bool Crit3DMeteoGridDbHandler::saveCellGridHourlyData(QString *errorStr, QString statement = QString(("REPLACE INTO `%1` VALUES")).arg(tableH); foreach (meteoVariable meteoVar, meteoVariableList) + { if (getVarFrequency(meteoVar) == hourly) { for (QDateTime myTime = firstTime; myTime <= lastTime; myTime = myTime.addSecs(3600)) @@ -3528,10 +3533,11 @@ bool Crit3DMeteoGridDbHandler::saveCellGridHourlyData(QString *errorStr, QString statement += QString(" ('%1','%2',%3),").arg(myTime.toString("yyyy-MM-dd hh:mm")).arg(varCode).arg(valueS); } } + } statement = statement.left(statement.length() - 1); - if( !qry.exec(statement) ) + if(! qry.exec(statement)) { *errorStr = qry.lastError().text(); return false; @@ -3562,6 +3568,7 @@ bool Crit3DMeteoGridDbHandler::saveCellGridHourlyDataEnsemble(QString *errorStr, statement = QString(("REPLACE INTO `%1` (%2, VariableCode, Value, MemberNr) VALUES ")).arg(tableH, _tableHourly.fieldTime); foreach (meteoVariable meteoVar, meteoVariableList) + { if (getVarFrequency(meteoVar) == hourly) { for (QDateTime myTime = firstTime; myTime <= lastTime; myTime = myTime.addSecs(3600)) @@ -3574,10 +3581,11 @@ bool Crit3DMeteoGridDbHandler::saveCellGridHourlyDataEnsemble(QString *errorStr, statement += QString(" ('%1','%2',%3,'%4'),").arg(myTime.toString("yyyy-MM-dd hh:mm")).arg(varCode).arg(valueS).arg(memberNr); } } + } statement = statement.left(statement.length() - 1); - if( !qry.exec(statement) ) + if(! qry.exec(statement)) { *errorStr = qry.lastError().text(); return false; @@ -3587,13 +3595,13 @@ bool Crit3DMeteoGridDbHandler::saveCellGridHourlyDataEnsemble(QString *errorStr, return true; } + bool Crit3DMeteoGridDbHandler::saveCellGridHourlyDataFF(QString *errorStr, QString meteoPointID, int row, int col, QDateTime firstTime, QDateTime lastTime) { QSqlQuery qry(_db); QString tableH = _tableHourly.prefix + meteoPointID + _tableHourly.postFix; QString tableFields; - for (unsigned int i=0; i < _tableHourly.varcode.size(); i++) { QString var = _tableHourly.varcode[i].varPragaName; From 917225c015d177d476a0c2df7a6a6ea668efa24a Mon Sep 17 00:00:00 2001 From: Atena Niazi Nasihati Date: Thu, 16 Jan 2025 12:04:51 +0100 Subject: [PATCH 04/12] aggiunte statistical summary --- commonDialogs/formSelectionSource.cpp | 8 +++++++- commonDialogs/formSelectionSource.h | 1 + 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/commonDialogs/formSelectionSource.cpp b/commonDialogs/formSelectionSource.cpp index d25b4256..73130e78 100644 --- a/commonDialogs/formSelectionSource.cpp +++ b/commonDialogs/formSelectionSource.cpp @@ -13,10 +13,12 @@ FormSelectionSource::FormSelectionSource() gridButton = new QRadioButton(tr("Meteo Grid")); pointButton =new QRadioButton(tr("Meteo Points")); + interpolationButton =new QRadioButton(tr("Interpolation Raster")); QHBoxLayout *sourceLayout = new QHBoxLayout; sourceLayout->addWidget(gridButton); sourceLayout->addWidget(pointButton); + sourceLayout->addWidget(interpolationButton); QGroupBox *sourceGroupBox = new QGroupBox("Source"); sourceGroupBox->setLayout(sourceLayout); @@ -39,7 +41,7 @@ void FormSelectionSource::done(int res) { if (res == QDialog::Accepted) // ok { - if (!pointButton->isChecked() && !gridButton->isChecked()) + if (!pointButton->isChecked() && !gridButton->isChecked() && !interpolationButton->isChecked()) { QMessageBox::information(nullptr, "Missing source selection.", "Please choose a data source."); return; @@ -65,6 +67,10 @@ int FormSelectionSource::getSourceSelectionId() { return 2; } + else if (interpolationButton->isChecked()) + { + return 3; + } else { return NODATA; diff --git a/commonDialogs/formSelectionSource.h b/commonDialogs/formSelectionSource.h index abdd0de6..efef9e4c 100644 --- a/commonDialogs/formSelectionSource.h +++ b/commonDialogs/formSelectionSource.h @@ -15,6 +15,7 @@ private: QRadioButton* pointButton; QRadioButton* gridButton; + QRadioButton* interpolationButton; void done(int res); }; From 1d2dab81c14dc79d0993b8f3ebfa4713626fcbc2 Mon Sep 17 00:00:00 2001 From: ftomei Date: Fri, 17 Jan 2025 11:03:38 +0100 Subject: [PATCH 05/12] v2.0.1 --- pragaProject/pragaShell.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pragaProject/pragaShell.cpp b/pragaProject/pragaShell.cpp index e8fa0cb1..8f641dda 100644 --- a/pragaProject/pragaShell.cpp +++ b/pragaProject/pragaShell.cpp @@ -900,7 +900,7 @@ int pragaBatch(PragaProject* myProject, QString scriptFileName) attachOutputToConsole(); #endif - myProject->logInfo("\nPRAGA v2.0.0"); + myProject->logInfo("\nPRAGA v2.0.1"); myProject->logInfo("Execute script: " + scriptFileName); if (scriptFileName == "") From 6d37d8a3213557d69a6426ce0060cbdfc80f5203 Mon Sep 17 00:00:00 2001 From: Caterina Toscano Date: Fri, 17 Jan 2025 14:58:46 +0100 Subject: [PATCH 06/12] created hydrall library --- hydrall/hydrall.cpp | 40 +++++++++++++++++++++++++++++++++ hydrall/hydrall.h | 31 +++++++++++++++++++++++++ hydrall/hydrall.pro | 35 +++++++++++++++++++++++++++++ mathFunctions/commonConstants.h | 2 ++ 4 files changed, 108 insertions(+) create mode 100644 hydrall/hydrall.cpp create mode 100644 hydrall/hydrall.h create mode 100644 hydrall/hydrall.pro diff --git a/hydrall/hydrall.cpp b/hydrall/hydrall.cpp new file mode 100644 index 00000000..bb6341d7 --- /dev/null +++ b/hydrall/hydrall.cpp @@ -0,0 +1,40 @@ +/*! + \name hydrall.cpp + \brief + \authors Antonio Volta, Caterina Toscano + +*/ + + +#include +#include +#include "crit3dDate.h" +#include "commonConstants.h" +#include "hydrall.h" + + +bool computeHydrall(Crit3DDate myDate, double myTemperature, double myElevation) +{ + getCO2(myDate, myTemperature, myElevation); + return true; +} + +double getCO2(Crit3DDate myDate, double myTemperature, double myElevation) +{ + double atmCO2 ; // fitting from data of Mauna Loa,Hawaii + if (myDate.year < 1990) + { + atmCO2= 280 * exp(0.0014876*(myDate.year -1840));//exponential change in CO2 concentration (ppm) + } + else + { + atmCO2= 350 * exp(0.00630*(myDate.year - 1990)); + } + atmCO2 += 3*cos(2*PI*getDoyFromDate(myDate)/365.0); // to consider the seasonal effects + return atmCO2*getPressureFromElevation(myTemperature, myElevation)/1000000 ; // [Pa] in +- ppm/10 +} + +double getPressureFromElevation(double myTemperature, double myElevation) +{ + return SEA_LEVEL_PRESSURE * exp((- GRAVITY * M_AIR * myElevation) / (R_GAS * myTemperature)); +} diff --git a/hydrall/hydrall.h b/hydrall/hydrall.h new file mode 100644 index 00000000..2882e175 --- /dev/null +++ b/hydrall/hydrall.h @@ -0,0 +1,31 @@ +#ifndef HYDRALL_H +#define HYDRALL_H + + #ifndef COMMONCONSTANTS_H + #include "commonConstants.h" + #endif + #ifndef CRIT3DDATE_H + #include "crit3dDate.h" + #endif + + #define UPSCALINGFUNC(z,LAI) ((1.0 - exp(-(z)*(LAI))) / (z)) + + // Tree-plant properties + #define FORM 0.5 // stem form factor + #define RHOF 0.1 // [KgDM m-3] foliage density + #define RHOS 750 // [KgDM m-3] wood-stem density + + // Hydraulic properties + #define H50 0.4 // height for 50% maturation of xylem cells (m) [not relevant] + #define KR 4.0E-7 // root specific conductance (m3 MPa-1 s-1 kg-1) [not relevant] + #define KSMAX 2.5E-3 // max. sapwood specific conductivity (m2 MPa-1 s-1) [not relevant] + #define PSITHR -2.5 // water potential threshold for cavitation (MPa) [not relevant] + + + #define NOT_INITIALIZED_VINE -1 + + double getCO2(Crit3DDate myDate, double myTemperature, double myElevation); + double getPressureFromElevation(double myTemperature, double myElevation); + + +#endif // HYDRALL_H diff --git a/hydrall/hydrall.pro b/hydrall/hydrall.pro new file mode 100644 index 00000000..8d0724a7 --- /dev/null +++ b/hydrall/hydrall.pro @@ -0,0 +1,35 @@ +#--------------------------------------------------- +# +# hydrall library +# This project is part of CRITERIA3D distribution +# +#--------------------------------------------------- + +QT -= core gui + +TEMPLATE = lib +CONFIG += staticlib + +DEFINES += _CRT_SECURE_NO_WARNINGS + +CONFIG += debug_and_release + +unix:{ + CONFIG(debug, debug|release) { + TARGET = debug/hydrall + } else { + TARGET = release/hydrall + } +} +win32:{ + TARGET = hydrall +} + +INCLUDEPATH += ../crit3dDate ../mathFunctions ../soil ../crop + +SOURCES += hydrall.cpp + + +HEADERS += hydrall.h + + diff --git a/mathFunctions/commonConstants.h b/mathFunctions/commonConstants.h index f37fd5ec..cce50ef9 100644 --- a/mathFunctions/commonConstants.h +++ b/mathFunctions/commonConstants.h @@ -159,6 +159,8 @@ #define MO2 0.032 // [kg mol-1] mass of molecular nitrogen (N2) #define MN2 0.028 + // [kg mol-1] mass of air + #define M_AIR 0.029 // [K] zero Celsius #define ZEROCELSIUS 273.15 // [] ratio molecular weight of water vapour/dry air From 3c36b0e92817ce4d282ca813b8991e9fa43c7d53 Mon Sep 17 00:00:00 2001 From: avolta Date: Fri, 17 Jan 2025 16:05:12 +0100 Subject: [PATCH 07/12] aggiunta libreria hydrall in criteria3d --- hydrall/hydrall.pro | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/hydrall/hydrall.pro b/hydrall/hydrall.pro index 8d0724a7..86998cee 100644 --- a/hydrall/hydrall.pro +++ b/hydrall/hydrall.pro @@ -5,14 +5,16 @@ # #--------------------------------------------------- -QT -= core gui +QT -= core gui TEMPLATE = lib CONFIG += staticlib -DEFINES += _CRT_SECURE_NO_WARNINGS - CONFIG += debug_and_release +CONFIG += c++11 c++14 c++17 + +#DEFINES += _CRT_SECURE_NO_WARNINGS + unix:{ CONFIG(debug, debug|release) { @@ -32,4 +34,3 @@ SOURCES += hydrall.cpp HEADERS += hydrall.h - From cedb5d4f046cd70e5f1fe46655bca493bd656bb7 Mon Sep 17 00:00:00 2001 From: ftomei Date: Tue, 21 Jan 2025 13:07:18 +0100 Subject: [PATCH 08/12] version --- pragaProject/pragaProject.cpp | 6 ++++++ pragaProject/pragaProject.h | 2 ++ pragaProject/pragaShell.cpp | 18 ++++++++++++++++-- pragaProject/pragaShell.h | 1 + 4 files changed, 25 insertions(+), 2 deletions(-) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index cc6e86a5..4f38073d 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -71,6 +71,12 @@ void PragaProject::clearPragaProject() } } + +QString PragaProject::getVersion() +{ + return "PRAGA V2.0.1"; +} + void PragaProject::createPragaProject(QString path_, QString name_, QString description_) { createProject(path_, name_, description_); diff --git a/pragaProject/pragaProject.h b/pragaProject/pragaProject.h index b0276ef5..790dd3df 100644 --- a/pragaProject/pragaProject.h +++ b/pragaProject/pragaProject.h @@ -90,6 +90,8 @@ void initializePragaProject(); void clearPragaProject(); + QString getVersion(); + void createPragaProject(QString path_, QString name_, QString description_); void savePragaProject(); void savePragaParameters(); diff --git a/pragaProject/pragaShell.cpp b/pragaProject/pragaShell.cpp index 8f641dda..3c797a63 100644 --- a/pragaProject/pragaShell.cpp +++ b/pragaProject/pragaShell.cpp @@ -13,6 +13,7 @@ QList getPragaCommandList() // praga commands cmdList.append("List | ListCommands"); + cmdList.append("Version | PragaVersion"); cmdList.append("Proj | OpenProject"); cmdList.append("Download | Download"); cmdList.append("AggrOnZones | GridAggregationOnZones"); @@ -50,6 +51,14 @@ int cmdList(PragaProject* myProject) } +int pragaVersion(PragaProject* myProject) +{ + myProject->logInfo(myProject->getVersion()); + + return PRAGA_OK; +} + + int PragaProject::executePragaCommand(QList argumentList, bool* isCommandFound) { *isCommandFound = false; @@ -57,11 +66,16 @@ int PragaProject::executePragaCommand(QList argumentList, bool* isComma QString command = argumentList[0].toUpper(); - if (command == "?" || command == "LIST" || command == "LISTCOMMANDS") + if (command == "?" || command == "LS" || command == "LIST" || command == "LISTCOMMANDS") { *isCommandFound = true; return cmdList(this); } + if (command == "VERSION" || command == "PRAGAVERSION") + { + *isCommandFound = true; + return pragaVersion(this); + } else if (command == "PROJ" || command == "OPENPROJECT") { *isCommandFound = true; @@ -900,7 +914,7 @@ int pragaBatch(PragaProject* myProject, QString scriptFileName) attachOutputToConsole(); #endif - myProject->logInfo("\nPRAGA v2.0.1"); + myProject->logInfo(myProject->getVersion()); myProject->logInfo("Execute script: " + scriptFileName); if (scriptFileName == "") diff --git a/pragaProject/pragaShell.h b/pragaProject/pragaShell.h index 8a2981f9..f2b2267a 100644 --- a/pragaProject/pragaShell.h +++ b/pragaProject/pragaShell.h @@ -7,6 +7,7 @@ QList getPragaCommandList(); int cmdList(PragaProject* myProject); + int pragaVersion(PragaProject* myProject); int executeCommand(QList argumentList, PragaProject* myProject); int pragaShell(PragaProject* myProject); From a550c9271ace571719eb105910a9b4e2af8e4c0f Mon Sep 17 00:00:00 2001 From: ftomei Date: Tue, 21 Jan 2025 14:26:13 +0100 Subject: [PATCH 09/12] help --- pragaProject/pragaProject.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 4f38073d..fdfca60a 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -74,7 +74,7 @@ void PragaProject::clearPragaProject() QString PragaProject::getVersion() { - return "PRAGA V2.0.1"; + return "PRAGA V2.0.1 (2025)"; } void PragaProject::createPragaProject(QString path_, QString name_, QString description_) From c5c57c9046cd1c6ba8ac93b1c1433c98287cc63a Mon Sep 17 00:00:00 2001 From: avolta Date: Tue, 21 Jan 2025 14:32:41 +0100 Subject: [PATCH 10/12] aggiunta mappa temperature medie del mese precedente per hydrall --- grapevine/grapevine.cpp | 2 +- hydrall/hydrall.cpp | 39 +++++++++++++++++++++++++++++++++++---- hydrall/hydrall.h | 4 +++- 3 files changed, 39 insertions(+), 6 deletions(-) diff --git a/grapevine/grapevine.cpp b/grapevine/grapevine.cpp index 83e04c5d..dd5ebd7f 100644 --- a/grapevine/grapevine.cpp +++ b/grapevine/grapevine.cpp @@ -27,7 +27,7 @@ bool Vine3D_Grapevine::compute(bool computeDaily, int secondsPerStep, Crit3DMode { simulationStepInSeconds = double(secondsPerStep); isAmphystomatic = true; - myLeafWidth = 0.2; // [m] + myLeafWidth = 0.2; // [cm] // Stomatal conductance Adjust stom conductance-photosynth ratio for soil water (Pa) alphaLeuning = modelCase->cultivar->parameterWangLeuning.alpha; getFixSimulationParameters(); diff --git a/hydrall/hydrall.cpp b/hydrall/hydrall.cpp index bb6341d7..5770584d 100644 --- a/hydrall/hydrall.cpp +++ b/hydrall/hydrall.cpp @@ -6,29 +6,45 @@ */ -#include +//#include #include #include "crit3dDate.h" #include "commonConstants.h" #include "hydrall.h" -bool computeHydrall(Crit3DDate myDate, double myTemperature, double myElevation) +bool computeHydrall(Crit3DDate myDate, double myTemperature, double myElevation, int secondPerStep) { getCO2(myDate, myTemperature, myElevation); + double actualLAI = getLAI(); + /* necessaria per ogni specie: + * il contenuto di clorofilla (g cm-2) il default è 500 + * lo spessore della foglia 0.2 cm default + * un booleano che indichi se la specie è anfistomatica oppure no + * parametro alpha del modello di Leuning + * + */ + // la temperatura del mese precedente arriva da fuori + + + return true; } double getCO2(Crit3DDate myDate, double myTemperature, double myElevation) { - double atmCO2 ; // fitting from data of Mauna Loa,Hawaii + double atmCO2 ; //https://www.eea.europa.eu/data-and-maps/daviz/atmospheric-concentration-of-carbon-dioxide-5/download.table + double year[24] = {1750,1800,1850,1900,1910,1920,1930,1940,1950,1960,1970,1980,1990,2000,2010,2020,2030,2040,2050,2060,2070,2080,2090,2100}; + double valueCO2[24] = {278,283,285,296,300,303,307,310,311,317,325,339,354,369,389,413,443,473,503,530,550,565,570,575}; + + // exponential fitting Mauna Loa if (myDate.year < 1990) { atmCO2= 280 * exp(0.0014876*(myDate.year -1840));//exponential change in CO2 concentration (ppm) } else { - atmCO2= 350 * exp(0.00630*(myDate.year - 1990)); + atmCO2= 353 * exp(0.00630*(myDate.year - 1990)); } atmCO2 += 3*cos(2*PI*getDoyFromDate(myDate)/365.0); // to consider the seasonal effects return atmCO2*getPressureFromElevation(myTemperature, myElevation)/1000000 ; // [Pa] in +- ppm/10 @@ -38,3 +54,18 @@ double getPressureFromElevation(double myTemperature, double myElevation) { return SEA_LEVEL_PRESSURE * exp((- GRAVITY * M_AIR * myElevation) / (R_GAS * myTemperature)); } + +double getLAI() +{ + // TODO + return 4; +} +/* +double meanLastMonthTemperature(double previousLastMonthTemp, double simulationStepInSeconds, double myInstantTemp) +{ + double newTemperature; + double monthFraction; + monthFraction = simulationStepInSeconds/(2592000.0); // seconds of 30 days + newTemperature = previousLastMonthTemp * (1 - monthFraction) + myInstantTemp * monthFraction ; + return newTemperature; +}*/ diff --git a/hydrall/hydrall.h b/hydrall/hydrall.h index 2882e175..19726653 100644 --- a/hydrall/hydrall.h +++ b/hydrall/hydrall.h @@ -24,8 +24,10 @@ #define NOT_INITIALIZED_VINE -1 + bool computeHydrall(Crit3DDate myDate, double myTemperature, double myElevation, int secondPerStep); double getCO2(Crit3DDate myDate, double myTemperature, double myElevation); double getPressureFromElevation(double myTemperature, double myElevation); - + double getLAI(); + double meanLastMonthTemperature(double previousLastMonthTemp, double simulationStepInSeconds, double myInstantTemp); #endif // HYDRALL_H From 598646c998275c13ff224c9117d1676a5cb66354 Mon Sep 17 00:00:00 2001 From: Atena Niazi Nasihati Date: Tue, 21 Jan 2025 14:59:35 +0100 Subject: [PATCH 11/12] update statistical summary --- commonDialogs/formSelectionSource.cpp | 8 +++++++- commonDialogs/formSelectionSource.h | 1 + 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/commonDialogs/formSelectionSource.cpp b/commonDialogs/formSelectionSource.cpp index 73130e78..27165bac 100644 --- a/commonDialogs/formSelectionSource.cpp +++ b/commonDialogs/formSelectionSource.cpp @@ -20,6 +20,7 @@ FormSelectionSource::FormSelectionSource() sourceLayout->addWidget(pointButton); sourceLayout->addWidget(interpolationButton); + QGroupBox *sourceGroupBox = new QGroupBox("Source"); sourceGroupBox->setLayout(sourceLayout); @@ -56,6 +57,12 @@ void FormSelectionSource::done(int res) } } +void FormSelectionSource::disableRadioButtons(bool pointDisable, bool gridDisable, bool interpolationDisable) +{ + pointButton->setCheckable(!pointDisable); + gridButton->setCheckable(!gridDisable); + interpolationButton->setCheckable(!interpolationDisable); +} int FormSelectionSource::getSourceSelectionId() { @@ -76,4 +83,3 @@ int FormSelectionSource::getSourceSelectionId() return NODATA; } } - diff --git a/commonDialogs/formSelectionSource.h b/commonDialogs/formSelectionSource.h index efef9e4c..6998f020 100644 --- a/commonDialogs/formSelectionSource.h +++ b/commonDialogs/formSelectionSource.h @@ -11,6 +11,7 @@ FormSelectionSource(); int getSourceSelectionId(); + void disableRadioButtons(bool pointDisable, bool gridDisable, bool interpolationDisable); private: QRadioButton* pointButton; From 32b0f6f76a40749cc94198f67c13f4bbd14cef95 Mon Sep 17 00:00:00 2001 From: ftomei Date: Tue, 21 Jan 2025 20:21:15 +0100 Subject: [PATCH 12/12] update 3D --- gis/color.cpp | 24 +-- graphics/mapGraphicsRasterUtm.cpp | 2 +- mathFunctions/commonConstants.h | 3 +- meteo/meteo.cpp | 2 +- soilFluxes3D/header/parameters.h | 2 +- soilFluxes3D/header/solver.h | 4 +- soilFluxes3D/header/types.h | 2 +- soilFluxes3D/heat.cpp | 2 +- soilFluxes3D/memory.cpp | 42 +++-- soilFluxes3D/soilFluxes3D.cpp | 3 +- soilFluxes3D/solver.cpp | 284 ++++++++++++++++++++++-------- soilFluxes3D/water.cpp | 19 +- 12 files changed, 267 insertions(+), 122 deletions(-) diff --git a/gis/color.cpp b/gis/color.cpp index 2ed7ab2f..c4cd4d25 100644 --- a/gis/color.cpp +++ b/gis/color.cpp @@ -97,6 +97,9 @@ bool Crit3DColorScale::classify() nrIntervals = MAXVALUE(_nrKeyColors -1, 1); nrStep = _nrColors / nrIntervals; + _nrColors = nrStep * nrIntervals; + color.resize(_nrColors); + for (i = 0; i < nrIntervals; i++) { dRed = float(keyColor[i+1].red - keyColor[i].red) / float(nrStep); @@ -120,24 +123,7 @@ bool Crit3DColorScale::classify() Crit3DColor* Crit3DColorScale::getColor(float value) { - unsigned int index = 0; - - if (value <= _minimum) - { - index = 0; - } - else if (value >= _maximum) - { - index = _nrColors-1; - } - else - { - if (_classification == classificationMethod::EqualInterval) - { - index = unsigned(float(_nrColors-1) * ((value - _minimum) / (_maximum - _minimum))); - } - } - + unsigned int index = getColorIndex(value); return &color[index]; } @@ -237,7 +223,7 @@ bool setAnomalyScale(Crit3DColorScale* myScale) bool setPrecipitationScale(Crit3DColorScale* myScale) { - myScale->initialize(6, 252); + myScale->initialize(6, 256); myScale->keyColor[0] = Crit3DColor(255, 255, 255); /*!< white */ myScale->keyColor[1] = Crit3DColor(0, 0, 255); /*!< blue */ diff --git a/graphics/mapGraphicsRasterUtm.cpp b/graphics/mapGraphicsRasterUtm.cpp index 4a38cb3c..99d84ca0 100644 --- a/graphics/mapGraphicsRasterUtm.cpp +++ b/graphics/mapGraphicsRasterUtm.cpp @@ -370,7 +370,7 @@ bool RasterUtmObject::drawRaster(QPainter* painter) // check outliers (transparent) if (_rasterPointer->colorScale->isHideOutliers()) { - if (value <= _rasterPointer->colorScale->minimum() || value > _rasterPointer->colorScale->maximum()) + if (isEqual(value, 0) || value <= _rasterPointer->colorScale->minimum() || value > _rasterPointer->colorScale->maximum()) continue; } diff --git a/mathFunctions/commonConstants.h b/mathFunctions/commonConstants.h index cce50ef9..579ee44c 100644 --- a/mathFunctions/commonConstants.h +++ b/mathFunctions/commonConstants.h @@ -109,7 +109,8 @@ #define BOUNDARY_SOLUTEFLUX 30 #define BOUNDARY_NONE 99 - #define RELAXATION 1 + #define GAUSS_SEIDEL 1 + #define JACOBI 2 // --------------- heat model ----------------- #define SAVE_HEATFLUXES_NONE 0 diff --git a/meteo/meteo.cpp b/meteo/meteo.cpp index 7e3222c4..96cde950 100644 --- a/meteo/meteo.cpp +++ b/meteo/meteo.cpp @@ -792,7 +792,7 @@ bool setColorScale(meteoVariable variable, Crit3DColorScale *colorScale) case snowFall: case snowWaterEquivalent: case snowLiquidWaterContent: case snowMelt: case dailyWaterTableDepth: setPrecipitationScale(colorScale); - if (variable == snowFall || variable == snowWaterEquivalent + if (variable == precipitation || variable == snowFall || variable == snowWaterEquivalent || variable == snowLiquidWaterContent || variable == snowMelt) { colorScale->setHideOutliers(true); diff --git a/soilFluxes3D/header/parameters.h b/soilFluxes3D/header/parameters.h index 17fa0a30..00687832 100644 --- a/soilFluxes3D/header/parameters.h +++ b/soilFluxes3D/header/parameters.h @@ -22,7 +22,7 @@ void initialize() { - numericalSolutionMethod = RELAXATION; + numericalSolutionMethod = GAUSS_SEIDEL; delta_t_min = 1; delta_t_max = 600; current_delta_t = delta_t_max; diff --git a/soilFluxes3D/header/solver.h b/soilFluxes3D/header/solver.h index 9c0129ea..ccf08a66 100644 --- a/soilFluxes3D/header/solver.h +++ b/soilFluxes3D/header/solver.h @@ -11,9 +11,7 @@ double arithmeticMean(double v1, double v2); - bool GaussSeidelRelaxation (int myApproximation, double myResidualTolerance, int myProcess); - - bool JacobiRelaxation (int approximation, double toleranceThreshold, int process); + bool solver(int myApproximation, double myResidualTolerance, int myProcess); #endif // SOLVER_H diff --git a/soilFluxes3D/header/types.h b/soilFluxes3D/header/types.h index 35b61f47..4a689005 100644 --- a/soilFluxes3D/header/types.h +++ b/soilFluxes3D/header/types.h @@ -145,7 +145,7 @@ extern TCrit3Dnode *nodeList; extern TmatrixElement **A; extern Tculvert myCulvert; - extern double *b, *C, *X; + extern double *b, *C, *X, *X1; extern double *invariantFlux; // array accessorio per flussi avvettivi e latenti extern double Courant; diff --git a/soilFluxes3D/heat.cpp b/soilFluxes3D/heat.cpp index bdd77bf0..c3e81da2 100644 --- a/soilFluxes3D/heat.cpp +++ b/soilFluxes3D/heat.cpp @@ -991,7 +991,7 @@ bool HeatComputation(double timeStep, double timeStepWater) return (false); } - GaussSeidelRelaxation(0, myParameters.ResidualTolerance, PROCESS_HEAT); + solver(0, myParameters.ResidualTolerance, PROCESS_HEAT); for (i = 1; i < myStructure.nrNodes; i++) nodeList[i].extra->Heat->T = X[i]; diff --git a/soilFluxes3D/memory.cpp b/soilFluxes3D/memory.cpp index be596131..4f75e9e8 100644 --- a/soilFluxes3D/memory.cpp +++ b/soilFluxes3D/memory.cpp @@ -44,10 +44,11 @@ void cleanArrays() } /*! free arrays */ + if (invariantFlux != nullptr) { free(invariantFlux); invariantFlux = nullptr; } if (b != nullptr) { free(b); b = nullptr; } if (C != nullptr) { free(C); C = nullptr; } - if (invariantFlux != nullptr) { free(invariantFlux); invariantFlux = nullptr; } if (X != nullptr) { free(X); X = nullptr; } + if (X1 != nullptr) { free(X1); X1 = nullptr; } } @@ -72,38 +73,49 @@ void cleanNodes() */ int initializeArrays() { - long i, j, n; - /*! clean previous arrays */ cleanArrays(); - /*! matrix solver: rows */ + /*! matrix A: rows */ A = (TmatrixElement **) calloc(myStructure.nrNodes, sizeof(TmatrixElement *)); - /*! matrix solver: columns */ - for (i = 0; i < myStructure.nrNodes; i++) - A[i] = (TmatrixElement *) calloc(myStructure.maxNrColumns, sizeof(TmatrixElement)); + /*! matrix A: columns */ + for (long i = 0; i < myStructure.nrNodes; i++) + { + A[i] = (TmatrixElement *) calloc(myStructure.maxNrColumns, sizeof(TmatrixElement)); + } - /*! initialize matrix solver */ - for (i = 0; i < myStructure.nrNodes; i++) - for (j = 0; j < (myStructure.nrLateralLinks + 2); j++) + /*! initialize matrix A */ + for (long i = 0; i < myStructure.nrNodes; i++) + { + for (int j = 0; j < myStructure.maxNrColumns; j++) { A[i][j].index = NOLINK; A[i][j].val = 0.; } + } b = (double *) calloc(myStructure.nrNodes, sizeof(double)); - for (n = 0; n < myStructure.nrNodes; n++) b[n] = 0.; + for (long i = 0; i < myStructure.nrNodes; i++) + { + b[i] = 0.; + } X = (double *) calloc(myStructure.nrNodes, sizeof(double)); - /*! mass diagonal matrix */ - C = (double *) calloc(myStructure.nrNodes, sizeof(double)); - for (n = 0; n < myStructure.nrNodes; n++) C[n] = 0.; + if (myParameters.numericalSolutionMethod == JACOBI) + { + X1 = (double *) calloc(myStructure.nrNodes, sizeof(double)); + } /*! mass diagonal matrix */ + C = (double *) calloc(myStructure.nrNodes, sizeof(double)); invariantFlux = (double *) calloc(myStructure.nrNodes, sizeof(double)); - for (n = 0; n < myStructure.nrNodes; n++) invariantFlux[n] = 0.; + for (long n = 0; n < myStructure.nrNodes; n++) + { + C[n] = 0.; + invariantFlux[n] = 0.; + } if (A == nullptr) return MEMORY_ERROR; diff --git a/soilFluxes3D/soilFluxes3D.cpp b/soilFluxes3D/soilFluxes3D.cpp index 3783c290..e1795620 100644 --- a/soilFluxes3D/soilFluxes3D.cpp +++ b/soilFluxes3D/soilFluxes3D.cpp @@ -53,8 +53,9 @@ TmatrixElement **A = nullptr; double *invariantFlux = nullptr; double *C = nullptr; -double *X = nullptr; double *b = nullptr; +double *X = nullptr; +double *X1 = nullptr; std::vector< std::vector> Soil_List; std::vector Surface_List; diff --git a/soilFluxes3D/solver.cpp b/soilFluxes3D/solver.cpp index cc17334f..eaafa9be 100644 --- a/soilFluxes3D/solver.cpp +++ b/soilFluxes3D/solver.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #include "basicMath.h" #include "header/types.h" @@ -108,6 +109,34 @@ int getMaxIterationsNr(int approximationNr) } +double GaussSeidelIterationHeat() +{ + double delta, new_x, norma_inf = 0.; + short j; + + for (long i = 1; i < myStructure.nrNodes; i++) + if (!nodeList[i].isSurface) + { + if (A[i][0].val != 0.) + { + j = 1; + new_x = b[i]; + while ((A[i][j].index != NOLINK) && (j < myStructure.maxNrColumns)) + { + new_x -= A[i][j].val * X[A[i][j].index]; + j++; + } + + delta = fabs(new_x - X[i]); + if (delta > norma_inf) norma_inf = delta; + X[i] = new_x; + } + } + + return norma_inf; +} + + double GaussSeidelIterationWater(short direction) { long firstIndex, lastIndex; @@ -138,20 +167,20 @@ double GaussSeidelIterationWater(short direction) } // surface check (H cannot go below z) - if (nodeList[i].isSurface) + if (nodeList[i].isSurface && newX < nodeList[i].z) { - if (newX < double(nodeList[i].z)) - newX = double(nodeList[i].z); + newX = nodeList[i].z; } - // water potential [m] - double psi = fabs(newX - double(nodeList[i].z)); + currentNorm = fabs(newX - X[i]); - // infinity norm (normalized if psi > 1m) + // water potential [m] + double psi = fabs(newX - nodeList[i].z); if (psi > 1) - currentNorm = (fabs(newX - X[i])) / psi; - else - currentNorm = fabs(newX - X[i]); + { + // normalized if psi > 1m + currentNorm /= psi; + } if (currentNorm > infinityNorm) infinityNorm = currentNorm; @@ -165,35 +194,180 @@ double GaussSeidelIterationWater(short direction) } -double GaussSeidelIterationHeat() +//------------- THREADS ------------- + +void GaussSeidel(int start, int end, double *infinityNorm, short direction) { - double delta, new_x, norma_inf = 0.; - short j; + long firstIndex, lastIndex; - for (long i = 1; i < myStructure.nrNodes; i++) - if (!nodeList[i].isSurface) + if (direction == UP) + { + firstIndex = start; + lastIndex = end + 1; + } + else + { + firstIndex = end; + lastIndex = start - 1; + } + + double currentNorm = 0.; + *infinityNorm = 0.; + + long i = firstIndex; + while (i != lastIndex) + { + double newX = b[i]; + short j = 1; + while ((A[i][j].index != NOLINK) && (j < myStructure.maxNrColumns)) { - if (A[i][0].val != 0.) - { - j = 1; - new_x = b[i]; - while ((A[i][j].index != NOLINK) && (j < myStructure.maxNrColumns)) - { - new_x -= A[i][j].val * X[A[i][j].index]; - j++; - } + newX -= A[i][j].val * X[A[i][j].index]; + j++; + } - delta = fabs(new_x - X[i]); - if (delta > norma_inf) norma_inf = delta; - X[i] = new_x; - } + // surface check (H cannot go below z) + if (nodeList[i].isSurface && newX < nodeList[i].z) + { + newX = nodeList[i].z; } - return norma_inf; - } + currentNorm = fabs(newX - X[i]); + + // water potential [m] + double psi = fabs(newX - nodeList[i].z); + if (psi > 1) + { + // normalized if psi > 1m + currentNorm /= psi; + } + + if (currentNorm > *infinityNorm) + *infinityNorm = currentNorm; + + X[i] = newX; + + (direction == UP)? i++ : i--; + } +} + + +double GaussSeidelThreads(short direction) +{ + int nrThreads = 4; + int step = myStructure.nrNodes / nrThreads; + + std::vector threads(nrThreads); + std::vector norm(nrThreads); + + for (int n = 0; n < nrThreads; n++) + { + int start = n * step; + int end = (n+1) * step - 1; + + if (n == nrThreads-1) + end = myStructure.nrNodes-1; + + threads[n] = std::thread(GaussSeidel, start, end, &(norm[n]), direction); + } + + // wait threads + for (auto& th : threads) { + th.join(); + } + + // compute norm + double infinityNorm = norm[0]; + for (int n = 1; n < nrThreads; n++) + { + if (norm[n] > infinityNorm) + infinityNorm = norm[n]; + } + + return infinityNorm; +} + + +void Jacobi(int start, int end, double *infinityNorm) +{ + *infinityNorm = 0; + double currentNorm = 0; + + for (long i = start; i <= end; i++) + { + X1[i] = b[i]; + short j = 1; + while ((A[i][j].index != NOLINK) && (j < myStructure.maxNrColumns)) + { + X1[i] -= (A[i][j].val * X[A[i][j].index]); + j++; + } + + // surface check (H cannot go below z) + if (nodeList[i].isSurface && X1[i] < nodeList[i].z) + { + X1[i] = nodeList[i].z; + } + + currentNorm = fabs(X1[i] - X[i]); + + // water potential [m] + double psi = fabs(X1[i] - nodeList[i].z); + if (psi > 1) + { + // normalized if psi > 1m + currentNorm /= psi; + } + + if (currentNorm > *infinityNorm) + *infinityNorm = currentNorm; + } +} + + +double JacobiThreads() +{ + int nrThreads = 4; + int step = myStructure.nrNodes / nrThreads; + + std::vector threads(nrThreads); + std::vector norm(nrThreads); + + for (int n = 0; n < nrThreads; n++) + { + int start = n * step; + int end = (n+1) * step - 1; + + if (n == nrThreads-1) + end = myStructure.nrNodes-1; + + threads[n] = std::thread(Jacobi, start, end, &(norm[n])); + } + + // wait threads + for (auto& th : threads) { + th.join(); + } + + // update X + for (long i = 0; i < myStructure.nrNodes; i++) + { + X[i] = X1[i]; + } + + double infinityNorm = norm[0]; + for (int n = 1; n < nrThreads; n++) + { + if (norm[n] > infinityNorm) + infinityNorm = norm[n]; + } + + return infinityNorm; +} + +//------------- SOLVER ------------- -bool GaussSeidelRelaxation (int approximation, double toleranceThreshold, int process) +bool solver (int approximation, double toleranceThreshold, int process) { double currentNorm = 1.0; double bestNorm = currentNorm; @@ -208,13 +382,20 @@ bool GaussSeidelRelaxation (int approximation, double toleranceThreshold, int pr else if (process == PROCESS_WATER) { - if (iteration%2 == 0) + if (myParameters.numericalSolutionMethod == GAUSS_SEIDEL) { - currentNorm = GaussSeidelIterationWater(DOWN); + if (iteration%2 == 0) + { + currentNorm = GaussSeidelThreads(DOWN); + } + else + { + currentNorm = GaussSeidelThreads(UP); + } } - else + else if (myParameters.numericalSolutionMethod == JACOBI) { - currentNorm = GaussSeidelIterationWater(UP); + currentNorm = JacobiThreads(); } if (currentNorm > (bestNorm * 10.0)) @@ -234,38 +415,3 @@ bool GaussSeidelRelaxation (int approximation, double toleranceThreshold, int pr return true; } - -bool JacobiRelaxation (int approximation, double toleranceThreshold, int process) -{ - double currentNorm = 1.0; - double bestNorm = currentNorm; - - int maxIterationsNr = getMaxIterationsNr(approximation); - int iteration = 0; - - while ((currentNorm > toleranceThreshold) && (iteration < maxIterationsNr)) - { - if (process == PROCESS_HEAT) - { - currentNorm = GaussSeidelIterationHeat(); - } - if (process == PROCESS_WATER) - { - // TODO Jacobi method - - if (currentNorm < bestNorm) - { - bestNorm = currentNorm; - } - else if (currentNorm > (bestNorm * 10.0)) - { - // non-convergent system - return false; - } - } - - iteration++; - } - - return true; -} diff --git a/soilFluxes3D/water.cpp b/soilFluxes3D/water.cpp index 709c590b..594b0b99 100644 --- a/soilFluxes3D/water.cpp +++ b/soilFluxes3D/water.cpp @@ -75,7 +75,6 @@ double getWaterExchange(long i, TlinkedNode *link, double deltaT) double runoff(long i, long j, TlinkedNode *link, double deltaT, unsigned long approximationNr) { double Hi, Hj; - double const EPSILON_mm = 0.0001; if (approximationNr == 0) { @@ -86,23 +85,25 @@ double runoff(long i, long j, TlinkedNode *link, double deltaT, unsigned long ap } else { - - Hi = nodeList[i].H; - Hj = nodeList[j].H; - /* + Hi = nodeList[i].H; + Hj = nodeList[j].H; + + /* Hi = (nodeList[i].H + nodeList[i].oldH) / 2.0; Hj = (nodeList[j].H + nodeList[j].oldH) / 2.0; - */ + */ } double H = MAXVALUE(Hi, Hj); double z = MAXVALUE(nodeList[i].z + nodeList[i].pond, nodeList[j].z + nodeList[j].pond); double Hs = H - z; - if (Hs <= 0.) return(0.); + if (Hs <= 0.) + return 0.; double dH = fabs(Hi - Hj); double cellDistance = distance2D(i,j); - if ((dH/cellDistance) < EPSILON_mm) return(0.); + if ((dH/cellDistance) < EPSILON) + return 0.; double roughness = (nodeList[i].Soil->Roughness + nodeList[j].Soil->Roughness) / 2.; @@ -334,7 +335,7 @@ bool waterFlowComputation(double deltaT) return false; } - if (! GaussSeidelRelaxation(approximationNr, myParameters.ResidualTolerance, PROCESS_WATER)) + if (! solver(approximationNr, myParameters.ResidualTolerance, PROCESS_WATER)) { if (deltaT > myParameters.delta_t_min) {