From 11e1804e90ab79810b464b82048cc5bc109e72cb Mon Sep 17 00:00:00 2001 From: daniel-v-k Date: Mon, 2 Sep 2024 10:11:53 +0200 Subject: [PATCH] fixing eqhead issue --- src/DataProcessing/DataReader.hpp | 70 +++++++++++++++---- src/Model/Node.hpp | 61 ++++++++++++---- tests/LargeModel/GlobalDataReader.hpp | 11 +-- .../LargeTransientModel/transient_runner.cpp | 2 +- 4 files changed, 111 insertions(+), 33 deletions(-) diff --git a/src/DataProcessing/DataReader.hpp b/src/DataProcessing/DataReader.hpp index 83a0989..acddc1d 100644 --- a/src/DataProcessing/DataReader.hpp +++ b/src/DataProcessing/DataReader.hpp @@ -421,17 +421,39 @@ namespace GlobalFlow { * @brief Read in a custom definition file for initial heads * @param path Where to read the file from */ - virtual void readInitialHeads(std::string path) { - readTwoColumns(path, [this](double data, int nodeID) { - if (nodes->at(nodeID)->hasGHB()){ + virtual void readInitialHeads(std::string path, bool isEqHeadFromFile) { + io::CSVReader<2, io::trim_chars<' ', '\t'>, io::no_quote_escape<','>> in(path); + in.read_header(io::ignore_no_column, "spatID", "data"); + large_num spatID{0}; + double head{0}; + large_num nodeID; + int layer{0}; + + int i{0}; + while (in.read_row(spatID, head)) { + try { + nodeID = lookupSpatIDtoNodeID.at(spatID).at(layer); + } + catch (const std::out_of_range &ex) { + //if Node does not exist ignore entry + continue; + } + + /*if (nodes->at(nodeID)->hasGHB()){ double ghbElevation = nodes->at(nodeID)->getExternalFlowElevation(Model::GENERAL_HEAD_BOUNDARY); - if (data < ghbElevation) { - data = ghbElevation; + if (head < ghbElevation) { + head = ghbElevation; } + }*/ + + nodes->at(nodeID)->setHead_allLayers(head * Model::si::meter); + nodes->at(nodeID)->setHead_TZero_allLayers(head * Model::si::meter); + if (!isEqHeadFromFile){ + nodes->at(nodeID)->setEqHeadToHead_allLayers(); } - nodes->at(nodeID)->setHead_allLayers(data * Model::si::meter); - nodes->at(nodeID)->setHead_TZero_allLayers(data * Model::si::meter); - }); + i++; + } + LOG(debug) << " ... for " << i << " nodes"; }; /** @@ -546,16 +568,38 @@ namespace GlobalFlow { * @brief Read equilibrium water-table information used for the dynamic river computation * @param path Where to read the file from */ - void readEqWTD(std::string path) { - readTwoColumns(path, [this](double data, int nodeID) { - if (nodes->at(nodeID)->hasGHB()){ + void readEqWTD(std::string path, bool isInitialHeadFromFile) { + io::CSVReader<2, io::trim_chars<' ', '\t'>, io::no_quote_escape<','>> in(path); + in.read_header(io::ignore_no_column, "spatID", "data"); + large_num spatID{0}; + int layer{0}; + double wtd{0}; + large_num nodeID{0}; + + int i{0}; + while (in.read_row(spatID, wtd)) { + try { + nodeID = lookupSpatIDtoNodeID.at(spatID).at(layer); + } + catch (const std::out_of_range &ex) { + //if Node does not exist ignore entry + continue; + } + + /*if (nodes->at(nodeID)->hasGHB()){ double ghbElevation = nodes->at(nodeID)->getExternalFlowElevation(Model::GENERAL_HEAD_BOUNDARY); if (data < ghbElevation) { data = ghbElevation; } + }*/ + + nodes->at(nodeID)->setEqHead_allLayers(wtd * Model::si::meter); + if (!isInitialHeadFromFile){ + nodes->at(nodeID)->setHeadToEqHead_allLayers(); } - nodes->at(nodeID)->setEqHead_allLayers(data * Model::si::meter); - }); + i++; + } + LOG(debug) << " ... for " << i << " nodes"; }; /** diff --git a/src/Model/Node.hpp b/src/Model/Node.hpp index bbd578e..4e76a3b 100644 --- a/src/Model/Node.hpp +++ b/src/Model/Node.hpp @@ -341,9 +341,22 @@ Set Properties }); }; + void setHead(const t_meter& head) noexcept { + NANChecker(head.value(), "Set Head"); + set(head); + } + + void setHead_allLayers(t_meter head) noexcept { + setHead(head); + set(head); + applyToAllLayers([&head](NodeInterface *nodeInterface) { + nodeInterface->setHead(head); + nodeInterface->set(head); + }); + } + /** * @brief Calculated equilibrium groundwater-head from eq_wtd - * Assumes that if initialhead = false that the eq_head is also used as initial head * @param head */ void setEqHead_allLayers(const t_meter& wtd) { @@ -357,6 +370,38 @@ Set Properties }); } + /** + * @brief Calculated equilibrium groundwater-head is used as initial head + * Assumes that if initialhead = false that the eq_head is also used as initial head + * @param head + */ + void setHeadToEqHead_allLayers() { + t_meter eqhead = get(); + setHead(eqhead); + applyToAllLayers([eqhead](NodeInterface *nodeInterface) { + try { + nodeInterface->setHead(eqhead); + } + catch (...) {} + }); + } + + /** + * @brief Initial groundwater-head is used as equilibrium head + * Assumes that if eqhead = false that the initial head is also used as eqhead + * @param head + */ + void setEqHeadToHead_allLayers() { + t_meter head = getHead(); + set(head); + applyToAllLayers([head](NodeInterface *nodeInterface) { + try { + nodeInterface->set(head); + } + catch (...) {} + }); + } + /** * @brief Update the current head change (in comparison to last time step) * @note Should only be called at end of time step @@ -378,20 +423,6 @@ Set Properties void setSourceZoneGHB(int sourceZoneGHB) { set(sourceZoneGHB); } - void setHead(const t_meter& head) noexcept { - NANChecker(head.value(), "Set Head"); - set(head); - } - - void setHead_allLayers(t_meter head) noexcept { - setHead(head); - set(head); - applyToAllLayers([&head](NodeInterface *nodeInterface) { - nodeInterface->setHead(head); - nodeInterface->set(head); - }); - } - /***************************************************************** Helpers ******************************************************************/ diff --git a/tests/LargeModel/GlobalDataReader.hpp b/tests/LargeModel/GlobalDataReader.hpp index dd5a6eb..9a4d342 100644 --- a/tests/LargeModel/GlobalDataReader.hpp +++ b/tests/LargeModel/GlobalDataReader.hpp @@ -117,13 +117,16 @@ class GlobalDataReader : public DataReader { LOG(userinfo) << "Reading elevation"; readElevation(buildDir(op.getElevation())); - // read either initial head (default) or equilibrium water table depth from file, if available + // Read initial head from file (sets EqHead to initial head if EqHead not from file) if (op.isInitialHeadFromFile()) { LOG(userinfo) << "Reading initial head"; - readInitialHeads((buildDir(op.getInitialHeadsDir()))); - } else if (op.isEqWTDFromFile()) { + readInitialHeads(buildDir(op.getInitialHeadsDir()), op.isEqWTDFromFile()); + } + + // Read equilibrium water table depth from file (sets head to EqHead if head not from file) + if (op.isEqWTDFromFile()) { LOG(userinfo) << "Reading equal water table depth"; - readEqWTD(buildDir(op.getEqWTD())); // requires elevation to be set + readEqWTD(buildDir(op.getEqWTD()), op.isInitialHeadFromFile()); // requires elevation to be set } /* diff --git a/tests/LargeTransientModel/transient_runner.cpp b/tests/LargeTransientModel/transient_runner.cpp index 7cadab8..2df3a45 100644 --- a/tests/LargeTransientModel/transient_runner.cpp +++ b/tests/LargeTransientModel/transient_runner.cpp @@ -27,7 +27,7 @@ namespace GlobalFlow { std::string pathToRecharge; std::string pathToGHB; - int year = 1891; + int year = 2048; // 1891; std::vector isSteadyState = op.getStressPeriodSteadyState(); std::vector numberOfSteps = op.getStressPeriodSteps();