Skip to content

Commit

Permalink
general debug early 2024
Browse files Browse the repository at this point in the history
  • Loading branch information
@daniel-v-k
daniel-v-k committed Feb 9, 2024
1 parent 108c9ef commit 62ee6e5
Show file tree
Hide file tree
Showing 25 changed files with 312 additions and 317 deletions.
67 changes: 37 additions & 30 deletions src/DataProcessing/DataReader.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -616,7 +616,7 @@ namespace GlobalFlow {
//if Node does not exist ignore entry
continue;
}
nodes->at(nodeID)->setElevation(elevation * Model::si::meter);
nodes->at(nodeID)->setElevation_allLayers(elevation * Model::si::meter);
i++;
}
LOG(debug) << " ... for " << i << " nodes";
Expand All @@ -636,7 +636,7 @@ namespace GlobalFlow {
int layer{0};
int refID{0};
double elevation{0};
large_num nodeID;
large_num nodeID{0};

int i{0};
while (in.read_row(spatID, refID, elevation)) {
Expand All @@ -647,7 +647,7 @@ namespace GlobalFlow {
//if Node does not exist ignore entry
continue;
}
nodes->at(nodeID)->setElevation(elevation * Model::si::meter);
nodes->at(nodeID)->setElevation_allLayers(elevation * Model::si::meter);
i++;
}
LOG(debug) << " ... for " << i << " nodes";
Expand Down Expand Up @@ -694,7 +694,7 @@ namespace GlobalFlow {
int arcID = -1;
std::vector<large_num> spatIDs;
std::unordered_map<large_num, large_num> nodeIDs;
double recharge = 0;
double recharge{0};
int layer{0};
large_num refID{0};

Expand Down Expand Up @@ -758,8 +758,8 @@ namespace GlobalFlow {
int layer{0};
int refID{0};
double conductivity{0};
large_num nodeID;
double threshold{1e-11}; // todo move to config
large_num nodeID{0};
//double threshold{1e-11}; // todo debug Equation: make preconditioner run with threshold (or delete such nodes)
int i{0};
int j{0};
while (in.read_row(spatID, layer, refID, conductivity)) {
Expand All @@ -770,16 +770,19 @@ namespace GlobalFlow {
//if Node does not exist ignore entry
continue;
}
if (conductivity < 1e-11){
conductivity = 1e-11;
}
nodes->at(nodeID)->setK_allLayers(conductivity * (Model::si::meter / Model::day));

i++;
if (conductivity < threshold){
/*if (conductivity < threshold){
nodes->at(nodeID)->setHeadActive_allLayers(false);
j++;
}
}*/
}
LOG(debug) << " ... for " << i << " nodes";
LOG(debug) << " ... < " << threshold << " at " << j << " nodes";
//LOG(debug) << " ... < " << threshold << " at " << j << " nodes";
};

/**
Expand All @@ -792,27 +795,24 @@ namespace GlobalFlow {
in.read_header(io::ignore_no_column, "spatID", "length", "width", "bankfull");

large_num spatID{0};
double length{0};
double width{0};
double riverLength{0};
double riverWidth{0};
double Q_bankfull{0};

std::unordered_map<large_num, large_num> nodeIDs;
std::unordered_map<large_num, large_num> refID_to_nodeID;
int layer{0};
large_num refID{0};

std::unordered_map<large_num, std::array<double, 3>> out;

while (in.read_row(spatID, length, width, Q_bankfull)) {
while (in.read_row(spatID, riverLength, riverWidth, Q_bankfull)) {
try {
nodeIDs = lookupSpatIDtoNodeIDs.at(spatID).at(layer); // layer is 0 since rivers are at surface
refID_to_nodeID = lookupSpatIDtoNodeIDs.at(spatID).at(layer); // layer is 0 since rivers are at surface
}
catch (const std::out_of_range &ex) { //if Node does not exist ignore entry
continue;
}
double bankfull = 0.349 * std::pow(Q_bankfull, 0.341);
double bankfull_depth = 0.349 * std::pow(Q_bankfull, 0.341);

for (auto nodeID : nodeIDs) { // in case the grid is refined: loop over all nodes at refIDs
out[nodeID.second] = {{bankfull, width, length * 1000}};
for (const auto &[refID, nodeID] : refID_to_nodeID) { // in case the grid is refined: loop over all nodes at refIDs
out[nodeID] = {{bankfull_depth, riverWidth, riverLength * 1000}};
}
}
return out;
Expand All @@ -830,37 +830,44 @@ namespace GlobalFlow {
double riverElevation{0};
double riverWidth{0};
double riverLength{0};
double bankfull_depth{0};
double riverConductance{0};
double riverDepth{0};
double riverBottom{0};
std::unordered_map<large_num, large_num> nodeIDs;
std::unordered_map<large_num, large_num> refID_to_nodeID;
int layer{0};
large_num refID{0};

int i{0};
while (in.read_row(spatID, riverElevation)) {
try {
nodeIDs = lookupSpatIDtoNodeIDs.at(spatID).at(layer);
refID_to_nodeID = lookupSpatIDtoNodeIDs.at(spatID).at(layer);
}
catch (const std::out_of_range &ex) { //if Node does not exist ignore entry
continue;
}

for (auto nodeID : nodeIDs) { // in case the grid is refined: loop over all nodes at refIDs
double bankfull_depth = riverStage[nodeID.second][0];
riverWidth = riverStage[nodeID.second][1];
riverLength = riverStage[nodeID.second][2];
for (const auto &[refID, nodeID] : refID_to_nodeID) { // in case the grid is refined: loop over all nodes at refIDs
bankfull_depth = riverStage[nodeID][0];
riverWidth = riverStage[nodeID][1];
riverLength = riverStage[nodeID][2];

if (bankfull_depth <= 1) {
bankfull_depth = 1; // river depth is at least 1 meter
}

if (nodes->at(nodeID)->hasGHB()){
double ghbElevation = nodes->at(nodeID)->getExternalFlowElevation(Model::GENERAL_HEAD_BOUNDARY);
if (riverElevation < ghbElevation){
riverElevation = ghbElevation; // river elevation is at least GHB elevation
}
}
riverBottom = riverElevation - bankfull_depth;
double K = nodes->at(nodeID.second)->getK().value();
double K = nodes->at(nodeID)->getK().value();
// Conductance estimation following Harbaugh (2005)
riverConductance = K * riverLength * riverWidth / (riverElevation - riverBottom);
if (riverConductance <= 1) { riverConductance = 1; }
if (riverConductance <= 1) { riverConductance = 1; } // river conductance is at least 1

nodes->at(nodeID.second)->addExternalFlow(Model::RIVER_MM,
nodes->at(nodeID)->addExternalFlow(Model::RIVER_MM,
riverElevation * Model::si::meter,
riverConductance,
riverBottom * Model::si::meter);
Expand Down
11 changes: 6 additions & 5 deletions src/DataProcessing/Neighbouring.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ void buildBySpatID(NodeVector nodes,
nodes->at(nodeID)->setNeighbours(nodeIDs_neig, neigPositions[neigPosID]);
} else {
sumBoundaries = addBoundary(nodes, boundaryConduct, boundaryCondition, nodeID, layer, isGlobal,
sumBoundaries);
sumBoundaries);
}
} else { // #### set neighbour of refined node ####
sumBoundaries = setNeigOfRefinedNode(nodes, spatID, neigPos, resolution, xRange, yRange, isGlobal,
Expand All @@ -99,10 +99,11 @@ void buildBySpatID(NodeVector nodes,


int setNeigOfRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPosition neigPos, double resolution,
large_num xRange, large_num yRange, bool isGlobal, large_num refID, large_num nodeID, int layer,
std::unordered_map<large_num, std::unordered_map<int, std::unordered_map<large_num, large_num>>> spatIDtoNodeIDs,
double boundaryConduct, Simulation::Options::BoundaryCondition boundaryCondition,
int sumBoundaries) {
large_num xRange, large_num yRange, bool isGlobal, large_num refID, large_num nodeID, int layer,
std::unordered_map<large_num, std::unordered_map<int,
std::unordered_map<large_num, large_num>>> spatIDtoNodeIDs,
double boundaryConduct,
Simulation::Options::BoundaryCondition boundaryCondition, int sumBoundaries) {
large_num refID_neig{};
std::unordered_map<large_num, large_num> nodeIDs = spatIDtoNodeIDs.at(spatID).at(layer);
// define map to set neighbour OUTSIDE refined node: neigPos >maps to> refID >maps to> refID_neig
Expand Down
4 changes: 2 additions & 2 deletions src/DataProcessing/Neighbouring.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,8 @@ int addBoundary(NodeVector const& nodes, double boundaryConduct, Simulation::Opt
int setNeigOfRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPosition neigPos, double resolution,
large_num xRange, large_num yRange, bool isGlobal, large_num refID, large_num nodeID, int layer,
std::unordered_map<large_num, std::unordered_map<int, std::unordered_map<large_num, large_num>>> spatIDtoNodeIDs,
double boundaryConduct, Simulation::Options::BoundaryCondition boundaryCondition,
int sumBoundaries);
double boundaryConduct,
Simulation::Options::BoundaryCondition boundaryCondition, int sumBoundaries);

std::unordered_map<Model::NeighbourPosition, std::unordered_map<large_num, large_num>>
defineMapOutside(large_num refinedInto);
Expand Down
65 changes: 27 additions & 38 deletions src/Model/Node.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,9 +140,6 @@ namespace GlobalFlow {
*/
class NodeInterface {
private:
virtual void
__setHeadChange(t_meter change) = 0;

virtual t_meter
__calcInitialHead(t_meter initialParam) = 0;

Expand Down Expand Up @@ -450,7 +447,7 @@ Set Properties
* @brief Set elevation on top layer and propagate to lower layers
* @param elevation The top elevation (e.g. from DEM)
*/
void setElevation(t_meter elevation) {
void setElevation_allLayers(t_meter elevation) {
set < t_meter, Elevation > (elevation);
set < t_meter, TopElevation > (elevation);
applyToAllLayers([this](NodeInterface *node) {
Expand Down Expand Up @@ -695,7 +692,7 @@ Calculate
void setK_allLayers(t_vel conduct) {
setK(conduct);
applyToAllLayers([&conduct](NodeInterface *nodeInterface) {
nodeInterface->getProperties().set<t_vel, K>(conduct);
nodeInterface->setK(conduct);
});
}

Expand Down Expand Up @@ -3047,33 +3044,37 @@ Calculate
t_s_meter_t conduct;

//Get all conductances from neighbouring cells
for (const auto &[neigPos, neigNodeID]: neighbours) {
// if head of neighbour is not active (K very close to 0): skip
if (!at(neigNodeID)->getHeadActive()) { continue; }
for (const auto &neighbour: neighbours) {
auto neigPos = neighbour.first;
auto neigNodeID = neighbour.second;
//LOG(debug) << "neigNodeID: " << neigNodeID;
conduct = 0 * si::square_meter / day;
if (neigPos == TOP or neigPos == DOWN) {
conduct = mechanics.calculateVerticalConductance(createDataTuple(neigNodeID));
//LOG(debug) << "vertical conductance: " << conduct.value();
} else {
if (get<int, Layer>() > 0 and get<bool, UseEfolding>()) {
conduct = mechanics.calculateEFoldingConductance(createDataTuple<Head>(neigPos, neigNodeID),
get<t_meter, EFolding>(),
getAt<t_meter, EFolding>(neigNodeID));
//LOG(debug) << "horizontal conductance (using e-folding): " << conduct.value();
// if head is active at this and neighboring node: change conductance to neig to value other than 0
if (at(neigNodeID)->getHeadActive() and getHeadActive()) {
if (neigPos == TOP or neigPos == DOWN) {
conduct = mechanics.calculateVerticalConductance(createDataTuple(neigNodeID));
//LOG(debug) << "vertical conductance: " << conduct.value();
} else {
conduct = mechanics.calculateHarmonicMeanConductance(createDataTuple<Head>(neigPos, neigNodeID));
//LOG(debug) << "horizontal conductance: " << conduct.value();
if (get<int, Layer>() > 0 and get<bool, UseEfolding>()) {
conduct = mechanics.calculateEFoldingConductance(createDataTuple<Head>(neigPos, neigNodeID),
get<t_meter, EFolding>(),
getAt<t_meter, EFolding>(neigNodeID));
//LOG(debug) << "horizontal conductance (using e-folding): " << conduct.value();
} else {
conduct = mechanics.calculateHarmonicMeanConductance(createDataTuple<Head>(neigPos, neigNodeID));
//LOG(debug) << "horizontal conductance: " << conduct.value();
}
}
NANChecker(conduct.value(), "Conductances");
}
NANChecker(conduct.value(), "Conductances");

// add conductance to out, the key in the unordered map is the ID of the neighbouring node
// (used to solve for the head at the neighbouring node)
out[neigNodeID] = conduct;
out[neigNodeID] = std::move(conduct);
}

// To solve for the head at this node, the conductances to neighbours and HCOF are used
t_s_meter_t conductNode = 0 * si::square_meter / day;

// To solve for the head at this node, the conductances to neighbours and HCOF are used
// subtract the conductances to neighbours (which were calculated above)
for (const auto &[neigNodeID, conductNeig]: out) { conductNode -= conductNeig; }
// add HCOF
Expand All @@ -3083,12 +3084,10 @@ Calculate
getP());
conductNode += hcof;
//LOG(debug) << "conductNode: " << conductNode.value();

// check for nan
NANChecker(conductNode.value(), "conductNode");

// add resulting conductance to solve for the head at this node to out
out[get<large_num, ID>()] = conductNode;
out[getID()] = std::move(conductNode);
return out;
};

Expand Down Expand Up @@ -3210,7 +3209,8 @@ Calculate
}

void setHeadChange(t_meter change) noexcept {
__setHeadChange(change);
NANChecker(change.value(), "Set Head Change at nodeID = " + std::to_string(getID()));
set<t_meter, HeadChange>(change);
}

t_meter calcInitialHead(t_meter initialParam) noexcept { return __calcInitialHead(initialParam); }
Expand Down Expand Up @@ -3337,17 +3337,6 @@ Calculate
//Learning weight
t_dim weight = 0.1;

/**
* @brief Update heads after one or multiple inner iterations
* @param delta
*/
virtual void __setHeadChange(t_meter change) {
NANChecker(change.value(), "Set Head Change at nodeID = " + std::to_string(getID()));
t_meter current_head = get<t_meter, Head>();
set<t_meter, HeadChange>(change);
set<t_meter, Head>(current_head + change);
};

virtual t_meter
__calcInitialHead(t_meter initialParam) {
t_meter elevation = get<t_meter, TopElevation>();
Expand Down
2 changes: 1 addition & 1 deletion src/Simulation/Options.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ namespace GlobalFlow {
exit(3);
}

BOUNDARY_CONDITION = config.get<std::string>("boundary_condition");
DEFAULT_BOUNDARY_CONDITION = config.get<std::string>("default_boundary_condition");
SENSITIVITY = config.get<bool>("sensitivity");

pt::ptree numerics = tree.get_child("numerics");
Expand Down
14 changes: 9 additions & 5 deletions src/Simulation/Options.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ namespace GlobalFlow {
std::vector<double> ANISOTROPY{10};
double SPECIFIC_YIELD{0.15};
double SPECIFIC_STORAGE{0.000015};
std::string BOUNDARY_CONDITION{"GeneralHeadBoundary"};
std::string DEFAULT_BOUNDARY_CONDITION{"GeneralHeadBoundary"};
bool SENSITIVITY{false};
std::vector<bool> CONFINED{};
// refinement information
Expand Down Expand Up @@ -138,7 +138,8 @@ namespace GlobalFlow {
enum BoundaryCondition {
GENERAL_HEAD_BOUNDARY,
GENERAL_HEAD_NEIGHBOUR,
STATIC_HEAD_SEA_LEVEL
STATIC_HEAD_SEA_LEVEL,
NONE
};

std::vector<int> getSteadyStateStressPeriodSteps() { return STEADY_STATE_STRESS_PERIOD_STEPS; }
Expand Down Expand Up @@ -178,13 +179,16 @@ namespace GlobalFlow {
std::vector<bool> getConfinements() { return CONFINED; }

BoundaryCondition getBoundaryCondition() {
if (BOUNDARY_CONDITION == "GeneralHeadBoundary") {
if (DEFAULT_BOUNDARY_CONDITION == "GeneralHeadBoundary") {
return BoundaryCondition::GENERAL_HEAD_BOUNDARY;
}
if (BOUNDARY_CONDITION == "GeneralHeadNeighbour") {
if (DEFAULT_BOUNDARY_CONDITION == "GeneralHeadNeighbour") {
return BoundaryCondition::GENERAL_HEAD_NEIGHBOUR;
}
return BoundaryCondition::STATIC_HEAD_SEA_LEVEL;
if (DEFAULT_BOUNDARY_CONDITION == "StaticSeaLevel"){
return BoundaryCondition::STATIC_HEAD_SEA_LEVEL;
}
return BoundaryCondition::NONE;
}

bool isSensitivity() { return SENSITIVITY; }
Expand Down
Loading

0 comments on commit 62ee6e5

Please sign in to comment.