Skip to content

Commit

Permalink
implementing refinement into 16 (3/3)
Browse files Browse the repository at this point in the history
  • Loading branch information
@daniel-v-k
daniel-v-k committed Oct 12, 2023
1 parent 17d4b59 commit 8bca820
Show file tree
Hide file tree
Showing 15 changed files with 114 additions and 102 deletions.
1 change: 0 additions & 1 deletion src/Model/FluidMechanics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,6 @@ namespace GlobalFlow {
//Transmissivity = K * thickness of prism
quantity<MeterSquaredPerTime> transmissivity_self = deltaV_self * k_self;
quantity<MeterSquaredPerTime> transmissivity_neig = deltaV_neig * k_neig;

if (transmissivity_neig != 0 * si::square_meter / day and
transmissivity_self != 0 * si::square_meter / day) {
out = nodeWidth * ((transmissivity_self * transmissivity_neig)
Expand Down
96 changes: 63 additions & 33 deletions src/Model/Node.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1120,7 +1120,7 @@ Calculate
mapToNeig = getMapNeighboursRefinedToFour(neighbourPosition);
} else if (nodeIDs.size() == 9){
mapToNeig = getMapNeighboursRefinedToNine(neighbourPosition);
} if (nodeIDs.size() == 16){
} else if (nodeIDs.size() == 16){
mapToNeig = getMapNeighboursRefinedToSixteen(neighbourPosition);
}else {
return;
Expand Down Expand Up @@ -1343,6 +1343,7 @@ Calculate
* @return
*/
t_vol_t getGNCFromNodes(){
//LOG(debug) << "getGNCFromNodes";
return -calculateGhostNodeCorrection(getPossibleRefinedNeighbours());
}

Expand All @@ -1356,18 +1357,19 @@ Calculate
if (getRefinedInto() == 1) { return out;} // if this node is not refined -> return

// get neighbours at left, right, front or back
std::forward_list<NeighbourPosition> neighbours_LRFB = getPossibleNeighbours_LRFB();
std::forward_list<NeighbourPosition> neigPos_LRFB = getNeigPos_LRFB();

for (const auto &unrefNeigPos: neighbours_LRFB) {
auto unrefNeig = neighbours.find(unrefNeigPos);
if (unrefNeig != neighbours.end()) {
for (const auto &neigPos: neigPos_LRFB) {
auto neig = neighbours.find(neigPos);
if (neig != neighbours.end()) {
// if neighbour is unrefined
if (at(unrefNeig)->getRefinedInto() == 1) { // if neig is not refined
if (at(neig)->getRefinedInto() == 1) { // if neig is not refined
// get the refined neighbour position this node has relative to that unrefined neighbour
NeighbourPosition refNeigPos = getRefNeigPosToUnrefNeig(get<large_num, RefID>(), unrefNeigPos,
NeighbourPosition refNeigPos = getRefNeigPosToUnrefNeig(get<large_num, RefID>(), neigPos,
getRefinedInto());
//LOG(debug) << "getGNCToRefinedNode, nodeID: " << getID() << ", neig nodeID: " << at(neig)->getID();

out += at(unrefNeig)->calculateGhostNodeCorrection({refNeigPos});
out += at(neig)->calculateGhostNodeCorrection({refNeigPos});
}
}
}
Expand Down Expand Up @@ -1437,31 +1439,32 @@ Calculate
}

t_vol_t calculateGhostNodeCorrection(const std::forward_list<NeighbourPosition>& possibleRefNeigPos){
t_vol_t gnc;
t_vol_t out = 0.0 * (si::cubic_meter / day);
t_meter head = getHead();
t_dim multiplierContributor{};
t_dim multiplierNodeInner{};
t_dim multiplierNodeOuter{};
int neigRefInto{};

for (const auto &position: possibleRefNeigPos) {
auto refinedNeig = neighbours.find(position);
if (refinedNeig != neighbours.end()) {
// todo: if contributor is refined: loop starts here
auto contributor = neighbours.find(getPotentialContributor(refinedNeig)); // contributor is named "j" in USG documentation
if (contributor == neighbours.end()) { // at model boundary: no contributor
continue; // todo implement impact of GHB (not required if all nodes at coast/GHB are refined)
}

neigRefInto = (int) at(refinedNeig)->getRefinedInto();
// conductance to the refined neighbour node
t_s_meter_t conductance =
mechanics.calculateHarmonicMeanConductance(createDataTuple<Head>(refinedNeig));

// todo: if contributor is refined: loop starts here
multiplierContributor = 0.5 * si::si_dimensionless; // if contributor refined: 1 / refinedInto
multiplierContributor = 0.5 * si::si_dimensionless; // if contributor refined: 1 / (2*sqrt(refinedInto))

t_s_meter_t contributorConductance = 0.0 * (si::square_meter / day);
t_s_meter_t nodeConductance = 0.0 * (si::square_meter / day);

auto contributor = neighbours.find(getPotentialContributor(refinedNeig)); // contributor is named "j" in USG documentation
if (contributor == neighbours.end()) { // at model boundary: no contributor
continue; // todo implement impact of GHB here
}
// conductance to the refined neighbour node
t_s_meter_t conductance =
mechanics.calculateHarmonicMeanConductance(createDataTuple<Head>(refinedNeig));

if (neigRefInto == 4) {
multiplierNodeInner = (1.0 / 4.0) * si::si_dimensionless;
Expand Down Expand Up @@ -1505,11 +1508,14 @@ Calculate
//LOG(debug) << "alpha = " << alpha << ", contributorConductance = " << contributorConductance.value();

// calculate ghost node correction
out += conductance * (alpha * (getHead() - at(contributor)->getHead()));
/*LOG(debug) << "GNC for nodeID " << get<large_num, ID>() <<
" to refined nodeID " << getAt<large_num, ID>(refinedNeig) <<
" with contributor " << getAt<large_num, ID>(contributor) <<
" = " << out.value();*/
gnc = conductance * (alpha * (getHead() - at(contributor)->getHead()));
out += gnc;
if (getID() == 46UL){
LOG(debug) << "GNC for nodeID " << get<large_num, ID>() <<
" to refined nodeID " << getAt<large_num, ID>(refinedNeig) <<
" with contributor " << getAt<large_num, ID>(contributor) <<
" = " << gnc.value();
}
}
}
return out;
Expand All @@ -1523,18 +1529,22 @@ Calculate
{NeighbourPosition::BACKLEFT, NeighbourPosition::LEFT},
{NeighbourPosition::FRONTFRONTLEFT, NeighbourPosition::LEFT},
{NeighbourPosition::BACKBACKLEFT, NeighbourPosition::LEFT},
{NeighbourPosition::LEFTLEFT, NeighbourPosition::LEFT}, // -> no contributor
{NeighbourPosition::FRONTRIGHT, NeighbourPosition::RIGHT},
{NeighbourPosition::BACKRIGHT, NeighbourPosition::RIGHT},
{NeighbourPosition::FRONTFRONTRIGHT, NeighbourPosition::RIGHT},
{NeighbourPosition::BACKBACKRIGHT, NeighbourPosition::RIGHT},
{NeighbourPosition::RIGHTRIGHT, NeighbourPosition::RIGHT}, // -> no contributor
{NeighbourPosition::LEFTFRONT, NeighbourPosition::FRONT},
{NeighbourPosition::RIGHTFRONT, NeighbourPosition::FRONT},
{NeighbourPosition::LEFTLEFTFRONT, NeighbourPosition::FRONT},
{NeighbourPosition::RIGHTRIGHTFRONT, NeighbourPosition::FRONT},
{NeighbourPosition::FRONTFRONT, NeighbourPosition::FRONT}, // -> no contributor
{NeighbourPosition::LEFTBACK, NeighbourPosition::BACK},
{NeighbourPosition::RIGHTBACK, NeighbourPosition::BACK},
{NeighbourPosition::LEFTLEFTBACK, NeighbourPosition::BACK},
{NeighbourPosition::RIGHTRIGHTBACK, NeighbourPosition::BACK},
{NeighbourPosition::BACKBACK, NeighbourPosition::BACK}, // -> no contributor
};
return mapPotentialContributor.at(refinedNeig->first);
}
Expand All @@ -1551,20 +1561,20 @@ Calculate
NeighbourPosition::BACKLEFT, NeighbourPosition::BACKBACK, NeighbourPosition::BACKRIGHT,
NeighbourPosition::BACKBACKLEFT, NeighbourPosition::BACKBACKRIGHT,
NeighbourPosition::LEFTFRONT, NeighbourPosition::LEFTLEFT, NeighbourPosition::LEFTBACK,
NeighbourPosition::LEFTLEFTFRONT, NeighbourPosition::LEFTLEFTBACK, NeighbourPosition::LEFTBACK,
NeighbourPosition::LEFTLEFTFRONT, NeighbourPosition::LEFTLEFTBACK,
NeighbourPosition::RIGHTFRONT, NeighbourPosition::RIGHTRIGHT, NeighbourPosition::RIGHTBACK,
NeighbourPosition::RIGHTRIGHTFRONT, NeighbourPosition::RIGHTRIGHTBACK};
}

static std::forward_list<NeighbourPosition>
getPossibleNeighbours_LRFB(){
getNeigPos_LRFB(){
return {NeighbourPosition::BACK, NeighbourPosition::FRONT,
NeighbourPosition::LEFT, NeighbourPosition::RIGHT};
}

std::forward_list<NeighbourPosition>
getPossibleNeighbours_horizontal(){
std::forward_list<NeighbourPosition> possiblePositions = getPossibleNeighbours_LRFB();
std::forward_list<NeighbourPosition> possiblePositions = getNeigPos_LRFB();
possiblePositions.merge(getPossibleRefinedNeighbours());
return possiblePositions;
}
Expand Down Expand Up @@ -2897,13 +2907,32 @@ Calculate
}

bool isLeftOrRight(NeighbourPosition neig) {
return (neig == LEFT or neig == RIGHT or
neig == LEFTFRONT or neig == LEFTBACK or neig == RIGHTFRONT or neig == RIGHTBACK);
std::unordered_map<NeighbourPosition, bool> leftRightMap = {
{LEFT, true}, {RIGHT, true},
{LEFTFRONT, true}, {LEFTBACK, true}, {RIGHTFRONT, true}, {RIGHTBACK, true},
{LEFTLEFTFRONT, true}, {LEFTLEFTBACK, true}, {RIGHTRIGHTFRONT, true}, {RIGHTRIGHTBACK, true},
{LEFTLEFT, true}, {RIGHTRIGHT, true}
};
try{
return leftRightMap.at(neig);
} catch (const std::out_of_range &e) { return false;}
//return (neig == LEFT or neig == RIGHT or
// neig == LEFTFRONT or neig == LEFTBACK or neig == RIGHTFRONT or neig == RIGHTBACK or);
}

bool isFrontOrBack(NeighbourPosition neig) {
return (neig == FRONT or neig == BACK or
neig == FRONTLEFT or neig == FRONTRIGHT or neig == BACKLEFT or neig == BACKRIGHT);
std::unordered_map<NeighbourPosition, bool> frontBackMap = {
{FRONT, true}, {BACK, true},
{FRONTLEFT, true}, {FRONTRIGHT, true}, {BACKLEFT, true}, {BACKRIGHT, true},
{FRONTFRONTLEFT, true}, {FRONTFRONTRIGHT, true}, {BACKBACKLEFT, true}, {BACKBACKRIGHT, true},
{FRONTFRONT, true}, {BACKBACK, true}
};

try{
return frontBackMap.at(neig);
} catch (const std::out_of_range &e) { return false;}
//return (neig == FRONT or neig == BACK or
// neig == FRONTLEFT or neig == FRONTRIGHT or neig == BACKLEFT or neig == BACKRIGHT);
}

/**
Expand Down Expand Up @@ -2966,7 +2995,8 @@ Calculate
* The left hand side of the equation
*/
std::unordered_map<large_num, t_s_meter_t> getMatrixEntries() {
size_t numC = 11;
size_t numC = getNumofNeighbours() + 1; // matrix needs 1 entry per neighbour + 1 entry for this node
//LOG(debug) << "numC: " << numC;
std::unordered_map<large_num, t_s_meter_t> out;
out.reserve(numC);
t_s_meter_t conduct;
Expand Down Expand Up @@ -3089,8 +3119,8 @@ Calculate

if(useGhostNodeCorrection) {
gncFromNodes = getGNCFromNodes();
gncToRefined = getGNCToRefinedNode();
//LOG(debug) << "gncFromNodes: " << gncFromNodes.value() << std::endl;
gncToRefined = getGNCToRefinedNode();
//LOG(debug) << "gncToRefined: " << gncToRefined.value() << std::endl;
}

Expand All @@ -3111,7 +3141,7 @@ Calculate

out += pseudoSourceNode + verticalFluxCorrections;
}
//LOG(debug) << "getRHS: " << out.value() << std::endl;
//LOG(debug) << "getRHS (nodeID: " << getID() << "): " << out.value() << std::endl;

NANChecker(out.value(), "RHS");
return out;
Expand Down
1 change: 1 addition & 0 deletions src/Solver/Equation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,7 @@ Equation::updateMatrix() {
addToA(nodes->at(id));
//---------------------Right
b(id) = nodes->at(id)->getRHS().value();

NANChecker(b[id], "Right hand side");
}

Expand Down
22 changes: 11 additions & 11 deletions tests/LargeModel/config_gnc.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,19 +3,22 @@
"model_config": {
"nodes": "grid_nz_refined.csv",
"number_of_nodes_per_layer": 6322,
"lon_range": 360,
"lat_range": 180,
"x_range": 360,
"y_range": 180,
"is_global": "true",
"resolution": 0.08333,
"resolution_in_degree": 0.08333,
"edge_length_left_right": -0,
"edge_length_front_back": -0,
"threads": 50,
"layers": 1,
"use_efolding": "false",
"confinement": [
"true"
],
"grid_refined": "true",
"boundary_condition": "GeneralHeadBoundary",
"sensitivity": "false"
},
"vdf_config": {
"density_variable": "true",
"density_zones": [
1000,
Expand All @@ -28,17 +31,14 @@
"max_toe_slope": 0.2,
"min_depth_factor": 0.1,
"slope_adj_factor": 0.1,
"vdf_lock": 0.001,
"cache": "false",
"adaptive_step_size": "false",
"boundary_condition": "GeneralHeadBoundary",
"sensitivity": "false"
"vdf_lock": 0.001
},
"numerics": {
"threads": 50,
"solver": "PCG",
"iterations": 10,
"inner_itter": 10,
"closing_crit_head": 1e-50,
"closing_crit_head": 1e-150,
"closing_crit_zeta": 1e-8,
"head_change": 0.1,
"zeta_change": 0.1,
Expand All @@ -53,7 +53,7 @@
"efold_as_array": "true",
"eq_wtd_from_file": "true",
"initial_head_from_file": "false",
"initial_zetas_as_array": "true",
"initial_zetas_as_array": "false",
"k_from_file": "true",
"k_ghb_from_file": "true",
"k_river_from_file": "false",
Expand Down
2 changes: 1 addition & 1 deletion tests/LargeModel/runner.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ namespace GlobalFlow {

void Runner::loadSettings() {
op = Simulation::Options();
op.load("data/config.json");
op.load("data/config_gnc.json");
}

void Runner::setupSimulation() {
Expand Down
6 changes: 4 additions & 2 deletions tests/SimpleRefinedModel/SimpleRefinedDataReader.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,8 +75,10 @@ class SimpleRefinedDataReader : public DataReader {
LOG(userinfo) << "Initializing head";
readInitialHeads((buildDir(op.getInitialHeadsDir())));

LOG(userinfo) << "Defining rivers";
readRiverConductance(buildDir(op.getKRiver()));
if(op.isKRiverFromFile()) {
LOG(userinfo) << "Defining rivers";
readRiverConductance(buildDir(op.getKRiver()));
}

if(op.isKGHBFromFile()) {
LOG(userinfo) << "Reading the boundary condition";
Expand Down
2 changes: 1 addition & 1 deletion tests/SimpleRefinedModel/config.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
"config": {
"model_config": {
"nodes": "grid_simpleRefined.csv",
"number_of_nodes_per_layer": 127,
"number_of_nodes_per_layer": 124,
"x_range": 10,
"y_range": 10,
"is_global": "false",
Expand Down
20 changes: 8 additions & 12 deletions tests/SimpleRefinedModel/config_coast.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,15 +4,14 @@
"nodes": "grid_coast_simpleRefined.csv",
"number_of_nodes_per_layer": 125,
"x_range": 10,
"y_range": 4,
"y_range": 5,
"is_global": "false",
"resolution_in_degree": 1,
"edge_length_left_right": 0,
"edge_length_front_back": 0,
"layers": 2,
"layers": 1,
"use_efolding": "false",
"confinement": [
"true",
"true"
],
"grid_refined": "true",
Expand All @@ -33,9 +32,9 @@
"numerics": {
"threads": 1,
"solver": "PCG",
"iterations": 100,
"inner_itter": 30,
"closing_crit_head": 1e-90,
"iterations": 1000,
"inner_itter": 100,
"closing_crit_head": 1e-20,
"closing_crit_zeta": 0,
"head_change": 0.01,
"zeta_change": 0,
Expand All @@ -49,7 +48,7 @@
"k_ghb_from_file": "true",
"specific_storage_from_file": "false",
"specific_yield_from_file": "false",
"k_river_from_file": "false",
"k_river_from_file": "true",
"aquifer_depth_from_file": "false",
"initial_head_from_file": "false",
"initial_zetas_as_array": "false",
Expand All @@ -62,17 +61,14 @@
"effective_porosity": 0.2,
"initial_head": 100,
"K": [
0.8,
0.8
],
"ghb_K": 800,
"aquifer_thickness": [
10,
10
100
],
"anisotropy": [
10,
10
1
],
"specific_yield": 0.15,
"specific_storage": 0.000015,
Expand Down
5 changes: 1 addition & 4 deletions tests/SimpleRefinedModel/data/elevation.csv
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,10 +45,7 @@ spatID,data,refID
33,3,2
33,3,3
33,3,4
34,3,1
34,3,2
34,3,3
34,3,4
34,3,0
35,3,1
35,3,2
35,3,3
Expand Down
Loading

0 comments on commit 8bca820

Please sign in to comment.