Skip to content

Commit

Permalink
fixing eqhead issue
Browse files Browse the repository at this point in the history
  • Loading branch information
@daniel-v-k
daniel-v-k committed Sep 2, 2024
1 parent 4e4e9c4 commit 11e1804
Show file tree
Hide file tree
Showing 4 changed files with 111 additions and 33 deletions.
70 changes: 57 additions & 13 deletions src/DataProcessing/DataReader.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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";
};

/**
Expand Down Expand Up @@ -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";
};

/**
Expand Down
61 changes: 46 additions & 15 deletions src/Model/Node.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -341,9 +341,22 @@ Set Properties
});
};

void setHead(const t_meter& head) noexcept {
NANChecker(head.value(), "Set Head");
set<t_meter, Head>(head);
}

void setHead_allLayers(t_meter head) noexcept {
setHead(head);
set<t_meter, EQHead>(head);
applyToAllLayers([&head](NodeInterface *nodeInterface) {
nodeInterface->setHead(head);
nodeInterface->set<t_meter, EQHead>(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) {
Expand All @@ -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<t_meter, EQHead>();
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<t_meter, EQHead>(head);
applyToAllLayers([head](NodeInterface *nodeInterface) {
try {
nodeInterface->set<t_meter, EQHead>(head);
}
catch (...) {}
});
}

/**
* @brief Update the current head change (in comparison to last time step)
* @note Should only be called at end of time step
Expand All @@ -378,20 +423,6 @@ Set Properties

void setSourceZoneGHB(int sourceZoneGHB) { set<int, SourceZoneGHB>(sourceZoneGHB); }

void setHead(const t_meter& head) noexcept {
NANChecker(head.value(), "Set Head");
set<t_meter, Head>(head);
}

void setHead_allLayers(t_meter head) noexcept {
setHead(head);
set<t_meter, EQHead>(head);
applyToAllLayers([&head](NodeInterface *nodeInterface) {
nodeInterface->setHead(head);
nodeInterface->set<t_meter, EQHead>(head);
});
}

/*****************************************************************
Helpers
******************************************************************/
Expand Down
11 changes: 7 additions & 4 deletions tests/LargeModel/GlobalDataReader.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
}

/*
Expand Down
2 changes: 1 addition & 1 deletion tests/LargeTransientModel/transient_runner.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ namespace GlobalFlow {
std::string pathToRecharge;
std::string pathToGHB;

int year = 1891;
int year = 2048; // 1891;

std::vector<bool> isSteadyState = op.getStressPeriodSteadyState();
std::vector<int> numberOfSteps = op.getStressPeriodSteps();
Expand Down

0 comments on commit 11e1804

Please sign in to comment.