Skip to content

Commit

Permalink
debugging large model
Browse files Browse the repository at this point in the history
  • Loading branch information
@daniel-v-k
daniel-v-k committed Dec 11, 2023
1 parent d03821b commit a07ccfb
Show file tree
Hide file tree
Showing 15 changed files with 397,046 additions and 452,887 deletions.
150 changes: 108 additions & 42 deletions src/DataProcessing/DataReader.hpp
View file

Large diffs are not rendered by default.

39 changes: 25 additions & 14 deletions src/DataProcessing/Neighbouring.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,16 +5,17 @@ namespace GlobalFlow {

/**
* @brief function to add boundaries to nodes with no neighbours
* @return total sum of boundaries added
*
*/
void addBoundary(NodeVector const& nodes,
int addBoundary(NodeVector const& nodes,
double boundaryConduct,
Simulation::Options::BoundaryCondition boundaryCondition,
large_num nodeID,
int layer, bool isGlobal) {
int layer, bool isGlobal, int sumBoundaries) {

if (layer > 0) {
return;
if (layer > 0 or nodes->at(nodeID)->hasGHB()) {
return sumBoundaries;
}

if (isGlobal){
Expand All @@ -32,6 +33,7 @@ void addBoundary(NodeVector const& nodes,
break;
}
}
return sumBoundaries + 1;
}

/**
Expand Down Expand Up @@ -61,6 +63,7 @@ void buildBySpatID(NodeVector nodes,
large_num refID;
large_num nodeID;
Model::NeighbourPosition neigPos;
int sumBoundaries{0};

for (large_num nodeIDTopLayer = 0; nodeIDTopLayer < numberOfNodesPerLayer; ++nodeIDTopLayer) {
spatID = nodes->at(nodeIDTopLayer)->getSpatID();
Expand All @@ -78,22 +81,26 @@ void buildBySpatID(NodeVector nodes,
nodeIDs_neig = spatIDtoNodeIDs.at(spatID_neig).at(layer);
nodes->at(nodeID)->setNeighbours(nodeIDs_neig, neigPositions[neigPosID]);
} else {
addBoundary(nodes, boundaryConduct, boundaryCondition, nodeID, layer, isGlobal);
sumBoundaries = addBoundary(nodes, boundaryConduct, boundaryCondition, nodeID, layer, isGlobal,
sumBoundaries);
}
} else { // #### set neighbour of refined node ####
setNeigOfRefinedNode(nodes, spatID, neigPos, resolution, xRange, yRange, isGlobal, refID, nodeID,
layer,spatIDtoNodeIDs, boundaryConduct, boundaryCondition);
sumBoundaries = setNeigOfRefinedNode(nodes, spatID, neigPos, resolution, xRange, yRange, isGlobal,
refID, nodeID, layer,spatIDtoNodeIDs, boundaryConduct,
boundaryCondition, sumBoundaries);
}
}
}
}
LOG(debug) << " Added " << sumBoundaries << " boundaries";
}


void setNeigOfRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPosition neigPos, double resolution,
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) {
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 All @@ -103,7 +110,8 @@ void setNeigOfRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPo
try {
refID_neig = mapOutside.at(neigPos).at(refID);
setNeighbourOutsideRefinedNode(nodes, spatID, neigPos, resolution, xRange, yRange, isGlobal, nodeID, layer,
spatIDtoNodeIDs, boundaryConduct, boundaryCondition, refID_neig, neigPos);
spatIDtoNodeIDs, boundaryConduct, boundaryCondition, refID_neig, neigPos,
sumBoundaries);
} catch (const std::out_of_range &ex) {}

// define map to set neighbour WITHIN refined node: neigPos >maps to> refID >maps to> refID_neig
Expand All @@ -115,6 +123,7 @@ void setNeigOfRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPo
large_num nodeID_ref_neig = spatIDtoNodeIDs.at(spatID).at(layer).at(refID_neig); // layer = 0
nodes->at(nodeID)->setNeighbour(nodeID_ref_neig, neigPos);
} catch (const std::out_of_range &ex) {}
return sumBoundaries;
}

/**
Expand Down Expand Up @@ -238,11 +247,12 @@ defineMapInside(large_num refinedInto) {
return mapInside;
}

void setNeighbourOutsideRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPosition neigPos, double resolution,
int setNeighbourOutsideRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPosition neigPos, double resolution,
large_num xRange, large_num yRange, bool isGlobal, 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,
large_num ref_id_neig, Model::NeighbourPosition neighbourPosition) {
large_num ref_id_neig, Model::NeighbourPosition neighbourPosition,
int sumBoundaries) {
int spatID_neig = getNeighbourSpatID((int) spatID, neigPos, resolution, xRange, yRange, isGlobal);
if (spatIDtoNodeIDs.contains(spatID_neig)) {
std::unordered_map<large_num, large_num> nodeIDs_neig = spatIDtoNodeIDs.at(spatID_neig).at(layer); // layer = 0
Expand All @@ -252,8 +262,9 @@ void setNeighbourOutsideRefinedNode(NodeVector nodes, large_num spatID, Model::N
nodes->at(nodeID)->setNeighbour(nodeIDs_neig.at(ref_id_neig), neighbourPosition);
}
} else {
addBoundary(nodes, boundaryConduct, boundaryCondition, nodeID, layer, isGlobal); // layer = 0
sumBoundaries = addBoundary(nodes, boundaryConduct, boundaryCondition, nodeID, layer, isGlobal, sumBoundaries); // layer = 0
}
return sumBoundaries;
};


Expand Down
14 changes: 8 additions & 6 deletions src/DataProcessing/Neighbouring.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,25 +61,27 @@ void copyNeighboursToBottomLayers(NodeVector nodes, int layers);
int getNeighbourSpatID(int spatID, Model::NeighbourPosition neigPos, double res, large_num xRange, large_num yRange,
bool isGlobal);

void addBoundary(NodeVector const& nodes, double boundaryConduct, Simulation::Options::BoundaryCondition boundaryCondition,
large_num nodeID, int layer, bool isGlobal);
int addBoundary(NodeVector const& nodes, double boundaryConduct, Simulation::Options::BoundaryCondition boundaryCondition,
large_num nodeID, int layer, bool isGlobal, int sumBoundaries);

void setNeigOfRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPosition neigPos, double resolution,
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);
double boundaryConduct, Simulation::Options::BoundaryCondition boundaryCondition,
int sumBoundaries);

std::unordered_map<Model::NeighbourPosition, std::unordered_map<large_num, large_num>>
defineMapOutside(large_num refinedInto);

std::unordered_map<Model::NeighbourPosition, std::unordered_map<large_num, large_num>>
defineMapInside(large_num refinedInto);

void setNeighbourOutsideRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPosition neigPos, double resolution,
int setNeighbourOutsideRefinedNode(NodeVector nodes, large_num spatID, Model::NeighbourPosition neigPos, double resolution,
large_num xRange, large_num yRange, bool isGlobal, 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,
large_num ref_id_neig, Model::NeighbourPosition neighbourPosition);
large_num ref_id_neig, Model::NeighbourPosition neighbourPosition,
int sumBoundaries);

/**
* Builds a map of neighbouring nodes based on x and y coordinates
Expand Down
48 changes: 24 additions & 24 deletions src/Model/ExternalFlows.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
namespace GlobalFlow {
namespace Model {

t_s_meter_t ExternalFlow::getP(t_meter eq_head, t_meter head,
t_s_meter_t ExternalFlow::getP(t_meter eq_gw_head, t_meter gw_head,
t_vol_t recharge,
t_vol_t eqFlow) const noexcept {

Expand All @@ -17,7 +17,7 @@ t_s_meter_t ExternalFlow::getP(t_meter eq_head, t_meter head,
return out;
case EVAPOTRANSPIRATION:
//flowHead = surface, bottom = extinction depth
if ((head < flowHead - bottom) xor (head > flowHead)) {
if ((gw_head < flowHead - bottom) xor (gw_head > flowHead)) {
return out;
} else {
return -special_flow / bottom;
Expand All @@ -30,19 +30,19 @@ t_s_meter_t ExternalFlow::getP(t_meter eq_head, t_meter head,
// River head <= river bottom (may happen in transient coupling)
if (flowHead <= bottom){
// Groundwater head >= river bottom --> allow gaining conditions of river!
if(head >= bottom) {
return -calcERC(recharge, eq_head, head, eqFlow);
if(gw_head >= bottom) {
return -calcERC(recharge, eq_gw_head, gw_head, eqFlow);
} else {
return out;
}
} else {
// River head > river bottom
return -calcERC(recharge, eq_head, head, eqFlow);
return -calcERC(recharge, eq_gw_head, gw_head, eqFlow);
}
case WETLAND:
//Can happen in transient coupling
if (flowHead <= bottom) {
if (head >= bottom) {
if (gw_head >= bottom) {
return -conductance;
} else {
return out;
Expand All @@ -53,7 +53,7 @@ t_s_meter_t ExternalFlow::getP(t_meter eq_head, t_meter head,
case GLOBAL_WETLAND:
//Can happen in transient coupling
if (flowHead <= bottom) {
if (head >= bottom) {
if (gw_head >= bottom) {
return -conductance;
} else {
return out;
Expand All @@ -64,7 +64,7 @@ t_s_meter_t ExternalFlow::getP(t_meter eq_head, t_meter head,
case LAKE:
//Can happen in transient coupling
if (flowHead <= bottom) {
if (head >= bottom) {
if (gw_head >= bottom) {
return -conductance;
} else {
return out;
Expand All @@ -75,7 +75,7 @@ t_s_meter_t ExternalFlow::getP(t_meter eq_head, t_meter head,
case GLOBAL_LAKE:
//Can happen in transient coupling
if (flowHead <= bottom) {
if (head >= bottom) {
if (gw_head >= bottom) {
return -conductance;
} else {
return out;
Expand All @@ -84,8 +84,8 @@ t_s_meter_t ExternalFlow::getP(t_meter eq_head, t_meter head,
return -conductance;
}
case DRAIN:
if (head > flowHead) {
return -calcERC(recharge, eq_head, head, eqFlow);
if (gw_head > flowHead) {
return -calcERC(recharge, eq_gw_head, gw_head, eqFlow);
} else {
return out;
}
Expand All @@ -95,7 +95,7 @@ t_s_meter_t ExternalFlow::getP(t_meter eq_head, t_meter head,
return out;
}

t_vol_t ExternalFlow::getQ(t_meter eq_head, t_meter head,
t_vol_t ExternalFlow::getQ(t_meter eq_gw_head, t_meter gw_head,
t_vol_t recharge,
t_vol_t eqFlow) const noexcept {
quantity<VolumePerTime, double> out = 0.0 * (si::cubic_meter / day);
Expand All @@ -107,32 +107,32 @@ t_vol_t ExternalFlow::getQ(t_meter eq_head, t_meter head,
case NET_ABSTRACTION:
return this->special_flow;
case EVAPOTRANSPIRATION:
if (head < flowHead - bottom) {
if (gw_head < flowHead - bottom) {
return out;
} else if (flowHead - bottom <= head and head <= flowHead) {
} else if (flowHead - bottom <= gw_head and gw_head <= flowHead) {
return -this->special_flow + (this->special_flow * flowHead) / bottom;
} else {
return -this->special_flow;
}
case FLOODPLAIN_DRAIN:
return -calculateFloodplainDrainage(head);
return -calculateFloodplainDrainage(gw_head);
case RIVER:
return conductance * flowHead;
case RIVER_MM:
//Can happen in transient coupling
if (flowHead <= bottom) {
if (head >= bottom) {
return calcERC(recharge, eq_head, head, eqFlow) * flowHead;
if (gw_head >= bottom) {
return calcERC(recharge, eq_gw_head, gw_head, eqFlow) * flowHead;
} else {
return out;
}
} else {
return calcERC(recharge, eq_head, head, eqFlow) * flowHead;
return calcERC(recharge, eq_gw_head, gw_head, eqFlow) * flowHead;
}
case WETLAND:
//Can happen in transient coupling
if (flowHead <= bottom) {
if (head >= bottom) {
if (gw_head >= bottom) {
return conductance * flowHead;
} else {
return out;
Expand All @@ -142,7 +142,7 @@ t_vol_t ExternalFlow::getQ(t_meter eq_head, t_meter head,
case GLOBAL_WETLAND:
//Can happen in transient coupling
if (flowHead <= bottom) {
if (head >= bottom) {
if (gw_head >= bottom) {
return conductance * flowHead;
} else {
return out;
Expand All @@ -152,7 +152,7 @@ t_vol_t ExternalFlow::getQ(t_meter eq_head, t_meter head,
case LAKE:
//Can happen in transient coupling
if (flowHead <= bottom) {
if (head >= bottom) {
if (gw_head >= bottom) {
return conductance * flowHead;
} else {
return out;
Expand All @@ -162,16 +162,16 @@ t_vol_t ExternalFlow::getQ(t_meter eq_head, t_meter head,
case GLOBAL_LAKE:
//Can happen in transient coupling
if (flowHead <= bottom) {
if (head >= bottom) {
if (gw_head >= bottom) {
return conductance * flowHead;
} else {
return out;
}
}
return conductance * flowHead;
case DRAIN:
if (head > flowHead) {
return calcERC(recharge, eq_head, head, eqFlow) * flowHead;
if (gw_head > flowHead) {
return calcERC(recharge, eq_gw_head, gw_head, eqFlow) * flowHead;
} else {
return out;
}
Expand Down
28 changes: 19 additions & 9 deletions src/Model/Node.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,15 +155,15 @@ namespace GlobalFlow {
* @return A dimensionless factor that can be used to modify
* hydraulic conductance depending on depth.
*
* E-Folding function as defined by Ying Fan et al. e^(-(Depth - Factor)).
* E-Folding function as defined by Ying Fan et al. e^(-Depth / E-Folding Depth).
*/
t_dim efoldingFromData(t_meter z) {
t_dim efoldingFromData(t_meter depth) {
t_meter folding = get<t_meter, EFolding>();
if (folding == 0.0 * si::meter)
return 1 * si::si_dimensionless;
//Alter if a different size should be used and not full vertical size
//z = (z / (2 * si::si_dimensionless));
t_dim out = exp(-z / folding);
//depth = (depth / (2 * si::si_dimensionless));
t_dim out = exp(-depth / folding);
if (out == 0 * si::si_dimensionless)
return 1e-7 * si::si_dimensionless;
return out;
Expand Down Expand Up @@ -876,8 +876,11 @@ Calculate
flow.getQ(eq_head, head, recharge, eqFlow)) * get<t_dim, StepModifier>();
}
} else { // GENERAL_HEAD_BOUNDARY (Question: what about FLOODPLAIN_DRAIN, EVAPOTRANSPIRATION, FAST_SURFACE_RUNOFF)
ex = (flow.getP(eq_head, head, recharge, eqFlow) * head + // = conductance * gw_head
flow.getQ(eq_head, head, recharge, eqFlow)) * get<t_dim, StepModifier>(); // = conductance * ghb_elevation (in examples = 0) * timestep
//LOG(debug) << "getting GHB flow";
ex = (flow.getP(eq_head, head, recharge, eqFlow) * head + // = ghb_conductance * gw_head
flow.getQ(eq_head, head, recharge, eqFlow)) // = ghb_conductance * ghb_elevation (often = 0)
* get<t_dim, StepModifier>();
//LOG(debug) << "GHB flow:" << ex.value();
}
return ex;
}
Expand Down Expand Up @@ -1662,16 +1665,20 @@ Calculate
zeta = getBottom();
}

// add initial values to vectors "zetas" and "zetasChange"
// add initial values to vectors "zetas", "zetasTZero" and "zetasChange"
std::vector<t_meter> zetas;
std::vector<t_meter> zetasTZero;
std::vector<t_meter> zetasChange;
if (localZetaID == 0) {
zetas.push_back(zeta);
zetasTZero.push_back(zeta);
zetasChange.push_back(0 * si::meter);
} else {
zetas = getZetas();
zetasTZero = getZetasTZero();
zetasChange = getZetasChange();
zetas.insert(zetas.begin() + localZetaID, zeta);
zetasTZero.insert(zetasTZero.begin() + localZetaID, zeta);
zetasChange.insert(zetasChange.begin() + localZetaID, 0 * si::meter);
}

Expand All @@ -1685,7 +1692,8 @@ Calculate
}

setZetas(zetas);
set<std::vector<t_meter>,ZetasChange>(zetas);
set<std::vector<t_meter>,Zetas_TZero>(zetasTZero);
set<std::vector<t_meter>,ZetasChange>(zetasChange);
//LOG(debug) << "nodeID: " << getID() << ", localZetaID: " << localZetaID << ", zeta: " << zeta.value() << ", getZeta: " << getZeta(localZetaID).value();
}

Expand Down Expand Up @@ -3316,7 +3324,9 @@ Calculate
* @param delta
*/
virtual void __setHeadChange(t_meter change) {
NANChecker(change.value(), "Set Head Change");
// todo pass a string including the nodeID
//std::string message = "Set Head Change at nodeID = " + std::to_string(getID());
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);
Expand Down
Loading

0 comments on commit a07ccfb

Please sign in to comment.