diff --git a/src/Model/FluidMechanics.cpp b/src/Model/FluidMechanics.cpp index 0400a8e6..2adc2022 100644 --- a/src/Model/FluidMechanics.cpp +++ b/src/Model/FluidMechanics.cpp @@ -117,7 +117,6 @@ namespace GlobalFlow { //Transmissivity = K * thickness of prism quantity transmissivity_self = deltaV_self * k_self; quantity 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) diff --git a/src/Model/Node.hpp b/src/Model/Node.hpp index d88be8ce..30794236 100644 --- a/src/Model/Node.hpp +++ b/src/Model/Node.hpp @@ -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; @@ -1343,6 +1343,7 @@ Calculate * @return */ t_vol_t getGNCFromNodes(){ + //LOG(debug) << "getGNCFromNodes"; return -calculateGhostNodeCorrection(getPossibleRefinedNeighbours()); } @@ -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 neighbours_LRFB = getPossibleNeighbours_LRFB(); + std::forward_list 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(), unrefNeigPos, + NeighbourPosition refNeigPos = getRefNeigPosToUnrefNeig(get(), neigPos, getRefinedInto()); + //LOG(debug) << "getGNCToRefinedNode, nodeID: " << getID() << ", neig nodeID: " << at(neig)->getID(); - out += at(unrefNeig)->calculateGhostNodeCorrection({refNeigPos}); + out += at(neig)->calculateGhostNodeCorrection({refNeigPos}); } } } @@ -1437,31 +1439,32 @@ Calculate } t_vol_t calculateGhostNodeCorrection(const std::forward_list& 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(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(refinedNeig)); if (neigRefInto == 4) { multiplierNodeInner = (1.0 / 4.0) * si::si_dimensionless; @@ -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() << - " to refined nodeID " << getAt(refinedNeig) << - " with contributor " << getAt(contributor) << - " = " << out.value();*/ + gnc = conductance * (alpha * (getHead() - at(contributor)->getHead())); + out += gnc; + if (getID() == 46UL){ + LOG(debug) << "GNC for nodeID " << get() << + " to refined nodeID " << getAt(refinedNeig) << + " with contributor " << getAt(contributor) << + " = " << gnc.value(); + } } } return out; @@ -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); } @@ -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 - getPossibleNeighbours_LRFB(){ + getNeigPos_LRFB(){ return {NeighbourPosition::BACK, NeighbourPosition::FRONT, NeighbourPosition::LEFT, NeighbourPosition::RIGHT}; } std::forward_list getPossibleNeighbours_horizontal(){ - std::forward_list possiblePositions = getPossibleNeighbours_LRFB(); + std::forward_list possiblePositions = getNeigPos_LRFB(); possiblePositions.merge(getPossibleRefinedNeighbours()); return possiblePositions; } @@ -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 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 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); } /** @@ -2966,7 +2995,8 @@ Calculate * The left hand side of the equation */ std::unordered_map 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 out; out.reserve(numC); t_s_meter_t conduct; @@ -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; } @@ -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; diff --git a/src/Solver/Equation.cpp b/src/Solver/Equation.cpp index b331182c..134419f2 100644 --- a/src/Solver/Equation.cpp +++ b/src/Solver/Equation.cpp @@ -138,6 +138,7 @@ Equation::updateMatrix() { addToA(nodes->at(id)); //---------------------Right b(id) = nodes->at(id)->getRHS().value(); + NANChecker(b[id], "Right hand side"); } diff --git a/tests/LargeModel/config_gnc.json b/tests/LargeModel/config_gnc.json index 85a95908..60e8fa76 100644 --- a/tests/LargeModel/config_gnc.json +++ b/tests/LargeModel/config_gnc.json @@ -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, @@ -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, @@ -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", diff --git a/tests/LargeModel/runner.cpp b/tests/LargeModel/runner.cpp index 1fa26b80..76f4b1c2 100644 --- a/tests/LargeModel/runner.cpp +++ b/tests/LargeModel/runner.cpp @@ -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() { diff --git a/tests/SimpleRefinedModel/SimpleRefinedDataReader.hpp b/tests/SimpleRefinedModel/SimpleRefinedDataReader.hpp index 18c6b377..e27f9b19 100644 --- a/tests/SimpleRefinedModel/SimpleRefinedDataReader.hpp +++ b/tests/SimpleRefinedModel/SimpleRefinedDataReader.hpp @@ -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"; diff --git a/tests/SimpleRefinedModel/config.json b/tests/SimpleRefinedModel/config.json index eef85d38..e6a71ffd 100644 --- a/tests/SimpleRefinedModel/config.json +++ b/tests/SimpleRefinedModel/config.json @@ -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", diff --git a/tests/SimpleRefinedModel/config_coast.json b/tests/SimpleRefinedModel/config_coast.json index 2a79a447..bf7d4487 100644 --- a/tests/SimpleRefinedModel/config_coast.json +++ b/tests/SimpleRefinedModel/config_coast.json @@ -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", @@ -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, @@ -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", @@ -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, diff --git a/tests/SimpleRefinedModel/data/elevation.csv b/tests/SimpleRefinedModel/data/elevation.csv index e4e2b921..da79daea 100644 --- a/tests/SimpleRefinedModel/data/elevation.csv +++ b/tests/SimpleRefinedModel/data/elevation.csv @@ -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 diff --git a/tests/SimpleRefinedModel/data/elevation_coast.csv b/tests/SimpleRefinedModel/data/elevation_coast.csv index 5e991d60..034df312 100644 --- a/tests/SimpleRefinedModel/data/elevation_coast.csv +++ b/tests/SimpleRefinedModel/data/elevation_coast.csv @@ -1,17 +1,17 @@ spatID,data,refID -0,0,1 +0,0.1,1 0,0.3,2 0,0.6,3 0,1,4 -0,0,5 +0,0.1,5 0,0.3,6 0,0.6,7 0,1,8 -0,0,9 +0,0.1,9 0,0.3,10 0,0.6,11 0,1,12 -0,0,13 +0,0.1,13 0,0.3,14 0,0.6,15 0,1,16 @@ -24,19 +24,19 @@ spatID,data,refID 7,5,0 8,5,0 9,5,0 -10,0,1 +10,0.1,1 10,0.3,2 10,0.6,3 10,1,4 -10,0,5 +10,0.1,5 10,0.3,6 10,0.6,7 10,1,8 -10,0,9 +10,0.1,9 10,0.3,10 10,0.6,11 10,1,12 -10,0,13 +10,0.1,13 10,0.3,14 10,0.6,15 10,1,16 @@ -49,19 +49,19 @@ spatID,data,refID 17,5,0 18,5,0 19,5,0 -20,0,1 +20,0.1,1 20,0.3,2 20,0.6,3 20,1,4 -20,0,5 +20,0.1,5 20,0.3,6 20,0.6,7 20,1,8 -20,0,9 +20,0.1,9 20,0.3,10 20,0.6,11 20,1,12 -20,0,13 +20,0.1,13 20,0.3,14 20,0.6,15 20,1,16 @@ -74,19 +74,19 @@ spatID,data,refID 27,5,0 28,5,0 29,5,0 -30,0,1 +30,0.1,1 30,0.3,2 30,0.6,3 30,1,4 -30,0,5 +30,0.1,5 30,0.3,6 30,0.6,7 30,1,8 -30,0,9 +30,0.1,9 30,0.3,10 30,0.6,11 30,1,12 -30,0,13 +30,0.1,13 30,0.3,14 30,0.6,15 30,1,16 @@ -99,19 +99,19 @@ spatID,data,refID 37,5,0 38,5,0 39,5,0 -40,0,1 +40,0.1,1 40,0.3,2 40,0.6,3 40,1,4 -40,0,5 +40,0.1,5 40,0.3,6 40,0.6,7 40,1,8 -40,0,9 +40,0.1,9 40,0.3,10 40,0.6,11 40,1,12 -40,0,13 +40,0.1,13 40,0.3,14 40,0.6,15 40,1,16 diff --git a/tests/SimpleRefinedModel/data/grid.csv b/tests/SimpleRefinedModel/data/grid.csv index 262c6101..908700de 100644 --- a/tests/SimpleRefinedModel/data/grid.csv +++ b/tests/SimpleRefinedModel/data/grid.csv @@ -45,10 +45,7 @@ spatID,lon,lat,area,refID 33,3.25,2.75,0.25,2 33,2.75,3.25,0.25,3 33,3.25,3.25,0.25,4 -34,3.75,2.75,0.25,1 -34,4.25,2.75,0.25,2 -34,3.75,3.25,0.25,3 -34,4.25,3.25,0.25,4 +34,4,3,1,0 35,4.75,2.75,0.25,1 35,5.25,2.75,0.25,2 35,4.75,3.25,0.25,3 diff --git a/tests/SimpleRefinedModel/data/heads.csv b/tests/SimpleRefinedModel/data/heads.csv index b2fcbce2..e698c0aa 100644 --- a/tests/SimpleRefinedModel/data/heads.csv +++ b/tests/SimpleRefinedModel/data/heads.csv @@ -45,10 +45,7 @@ spatID,data,refID 33,5,2 33,5,3 33,5,4 -34,5,1 -34,5,2 -34,5,3 -34,5,4 +34,5,0 35,5,1 35,5,2 35,5,3 diff --git a/tests/SimpleRefinedModel/data/lithology.csv b/tests/SimpleRefinedModel/data/lithology.csv index 9c4bd816..db97bfc1 100644 --- a/tests/SimpleRefinedModel/data/lithology.csv +++ b/tests/SimpleRefinedModel/data/lithology.csv @@ -45,10 +45,7 @@ spatID,data,refID 33,0.8,2 33,0.8,3 33,0.8,4 -34,0.8,1 -34,0.8,2 -34,0.8,3 -34,0.8,4 +34,0.8,0 35,0.8,1 35,0.8,2 35,0.8,3 diff --git a/tests/SimpleRefinedModel/data/rivers_coast.csv b/tests/SimpleRefinedModel/data/rivers_coast.csv index d7278f93..b7453c9e 100644 --- a/tests/SimpleRefinedModel/data/rivers_coast.csv +++ b/tests/SimpleRefinedModel/data/rivers_coast.csv @@ -1,17 +1,13 @@ spatID,Head,Bottom,Conduct,refID -30,0,0,480,1 30,0.2,0,480,2 30,0.3,0,480,3 30,0.4,0,480,4 -30,0,0,480,5 30,0.2,0,480,6 30,0.3,0,480,7 30,0.4,0,480,8 -30,0,0,480,9 30,0.2,0,480,10 30,0.3,0,480,11 30,0.4,0,480,12 -30,0,0,480,13 30,0.2,0,480,14 30,0.3,0,480,15 30,0.4,0,480,16 diff --git a/tests/SimpleRefinedModel/simpleRefined.cpp b/tests/SimpleRefinedModel/simpleRefined.cpp index 6b99933c..0044579b 100644 --- a/tests/SimpleRefinedModel/simpleRefined.cpp +++ b/tests/SimpleRefinedModel/simpleRefined.cpp @@ -5,7 +5,7 @@ namespace GlobalFlow { void StandaloneRunner::loadSettings() { op = Simulation::Options(); - op.load("data/config_coast_simpleRefined.json"); + op.load("data/config_simpleRefined.json"); } void StandaloneRunner::setupSimulation() {