Skip to content

Commit

Permalink
reading and running variable vertical size
Browse files Browse the repository at this point in the history
  • Loading branch information
@daniel-v-k
daniel-v-k committed Sep 14, 2024
1 parent ae4e452 commit 1040033
Show file tree
Hide file tree
Showing 9 changed files with 154 additions and 149 deletions.
41 changes: 39 additions & 2 deletions src/DataProcessing/DataReader.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -237,7 +237,7 @@ namespace GlobalFlow {
int numberOfLayers,
double defaultK,
double initialHead,
double aquiferDepth,
double verticalSize,
double anisotropy,
double specificYield,
double specificStorage,
Expand Down Expand Up @@ -282,7 +282,7 @@ namespace GlobalFlow {
nodeID,
defaultK * (Model::si::meter / Model::day),
initialHead * Model::si::meter,
aquiferDepth,
verticalSize,
anisotropy,
specificYield,
specificStorage,
Expand Down Expand Up @@ -369,6 +369,7 @@ namespace GlobalFlow {
}
ghbVerticalSize = nodes->at(nodeID)->getVerticalSize().value() -
(nodes->at(nodeID)->getElevation().value() - elevation);
if (ghbVerticalSize < 0) { ghbVerticalSize = 0; }

if (conductivity < conductivity_of_clay) {
conductivity = conductivity_of_clay;
Expand Down Expand Up @@ -970,6 +971,42 @@ namespace GlobalFlow {
LOG(debug) << " ... for " << i << " nodes";
};

/**
* @brief Read elevation data from a specified path
* Used for multiple files
* @note !Uses setElevation() function. Should only be called after all layers are build
* as it affects layers below
* @param path Where to read the file from
* @param files If different files for different regions are given
*/
void readVerticalSize(std::string path, std::vector<std::string> files) {
loopFiles(path, files, [this](std::string path) {
io::CSVReader<3, io::trim_chars<' ', '\t'>, io::no_quote_escape<','>> in(path);
in.read_header(io::ignore_no_column, "spatID", "layer", "vertical_size");
large_num spatID{0};
int layer{0};
double vertical_size{0};
large_num nodeID;

int i{0};
while (in.read_row(spatID, layer, vertical_size)) {
try {
nodeID = lookupSpatIDtoNodeID.at(spatID).at(layer);
}
catch (const std::out_of_range &ex) {
//if Node does not exist ignore entry
continue;
}
if (vertical_size < 1) {
vertical_size = 1;
}
nodes->at(nodeID)->setVerticalSize(vertical_size * Model::si::meter);
i++;
}
LOG(debug) << " ... for " << i << " nodes";
});
};

void setDefaultZetas(large_num numberOfZones) {
for (const auto &node : *nodes) {
// initialize zeta surface at top and bottom
Expand Down
10 changes: 5 additions & 5 deletions src/DataProcessing/Neighbouring.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,7 +167,7 @@ void copyNeighboursToBottomLayers(NodeVector nodes, int numberOfLayers) {
void buildBottomLayers(NodeVector nodes,
int numberOfLayers,
std::vector<bool> confined,
std::vector<int> aquifer_thickness,
std::vector<double> verticalSizes,
std::vector<double> conductances,
std::vector<double> anisotropies) {
assert(numberOfLayers && "AsModel::signing 0 layers does not make any sense");
Expand All @@ -189,7 +189,7 @@ void buildBottomLayers(NodeVector nodes,
Model::quantity<Model::Meter> edgeLengthFrontBack;
Model::quantity<Model::Velocity> K;
Model::quantity<Model::Meter> head;
double aquiferDepth;
double verticalSize;
double anisotropy;
double specificYield;
double specificStorage;
Expand Down Expand Up @@ -222,7 +222,7 @@ void buildBottomLayers(NodeVector nodes,
i)->getProperties().get<Model::quantity<Model::Meter>, Model::EdgeLengthFrontBack>();
K = conductances[layer + 1] * Model::si::meter / Model::day; // or: nodes->at(i)->getK__pure();
head = nodes->at(i)->getProperties().get<Model::quantity<Model::Meter>, Model::Head>();
aquiferDepth = aquifer_thickness[layer + 1];
verticalSize = verticalSizes[layer + 1];
anisotropy = anisotropies[layer + 1]; // or nodes->at(i)->getProperties().get<Model::quantity<Model::Dimensionless>, Model::Anisotropy>().value();
specificYield = nodes->at(i)->getProperties().get<Model::quantity<Model::Dimensionless>,
Model::SpecificYield>().value();
Expand Down Expand Up @@ -263,7 +263,7 @@ void buildBottomLayers(NodeVector nodes,
id,
K,
head,
aquiferDepth,
verticalSize,
anisotropy,
specificYield,
specificStorage,
Expand All @@ -284,7 +284,7 @@ void buildBottomLayers(NodeVector nodes,
nodes->at(id)->getProperties().set<int, Model::Layer>(layer + 1);
nodes->at(id)->getProperties().set<Model::quantity<Model::Meter>, Model::Elevation>(
nodes->at(id)->getProperties().get<Model::quantity<Model::Meter>, Model::Elevation>()
- (aquiferDepth * Model::si::meter));
- (verticalSize * Model::si::meter));
}
//2) Neighbouring for top and bottom
nodes->at(id)->setNeighbour(i + (layer * nodesPerLayer), Model::TOP);
Expand Down
2 changes: 1 addition & 1 deletion src/DataProcessing/Neighbouring.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ void
buildBottomLayers(NodeVector nodes,
int layers,
std::vector<bool> confined,
std::vector<int> aquifer_depth,
std::vector<double> aquifer_depth,
std::vector<double> conductances,
std::vector<double> anisotropies);
}
Expand Down
4 changes: 2 additions & 2 deletions src/Model/FluidMechanics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,8 +113,8 @@ 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) {
if (transmissivity_neig > 0 * si::square_meter / day and
transmissivity_self > 0 * si::square_meter / day) {
out = nodeWidth * ((transmissivity_self * transmissivity_neig)
/ (transmissivity_self * nodeLength_neig * 0.5 + // half of neighbour node's length
transmissivity_neig * nodeLength_self * 0.5)); // half of this node's length
Expand Down
4 changes: 2 additions & 2 deletions src/Model/Node.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ NodeInterface::NodeInterface(NodeVector nodes,
unsigned long int identifier,
quantity<Velocity> conduct,
quantity<Meter> head,
double aquiferDepth,
double verticalSize,
double anisotropy,
double specificYield,
double specificStorage,
Expand Down Expand Up @@ -77,7 +77,7 @@ NodeInterface::NodeInterface(NodeVector nodes,
fields.set<bool, Confinement>(confined);
fields.emplace<quantity<Dimensionless>, SpecificYield>(specificYield * si::si_dimensionless);
fields.emplace<quantity<perUnit>, SpecificStorage>(specificStorage * perMeter);
fields.emplace<quantity<Meter>, VerticalSize>(aquiferDepth * si::meter);
fields.emplace<quantity<Meter>, VerticalSize>(verticalSize * si::meter);
fields.emplace<quantity<Dimensionless>, Anisotropy>(anisotropy * si::si_dimensionless);
fields.emplace<quantity<Meter>, EdgeLengthLeftRight>(edgeLengthLeftRight);
fields.emplace<quantity<Meter>, EdgeLengthFrontBack>(edgeLengthFrontBack);
Expand Down
25 changes: 16 additions & 9 deletions src/Model/Node.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,13 +109,13 @@ namespace GlobalFlow {
*
* E-Folding function as defined by Ying Fan et al. e^(-Depth / E-Folding Depth).
*/
t_dim efoldingFromData(t_meter depth) {
t_meter folding = get<t_meter, EFolding>();
t_dim efoldingFromData(t_meter verticalSize) {
auto 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
//depth = (depth / (2 * si::si_dimensionless));
t_dim out = exp(-depth / folding);
t_dim out = exp(-verticalSize / folding);
if (out == 0 * si::si_dimensionless)
return 1e-7 * si::si_dimensionless;
return out;
Expand Down Expand Up @@ -679,6 +679,12 @@ Set Properties
});
};

/**
* @brief Set aquifer depth from data
* @param aquifer_depth aquifer depth value in meters
*/
void setVerticalSize(const t_meter& aquifer_depth) { set < t_meter, VerticalSize > (aquifer_depth); };

/**
* @brief Set e-folding factor from data on all layers
* @param efold efolding value
Expand Down Expand Up @@ -2947,6 +2953,9 @@ Calculate
getAt<t_meter, EFolding>(neigNodeID));
} else {
conduct = mechanics.calculateHarmonicMeanConductance(createDataTuple<Head>(neigPos, neigNodeID));
if (conduct.value() <= 0) {
LOG(debug) << "conduct = " << conduct.value();
}
}
}
NANChecker(conduct.value(), "Conductances");
Expand Down Expand Up @@ -2986,14 +2995,14 @@ Calculate
std::vector<t_meter> zoneThicknesses;

// if this zeta interface is not active: return directly
if (isZetaInactive(zetaID)){
/*if (isZetaInactive(zetaID)){
LOG(userinfo) << "Asking for matrix entry at inactive zeta, causing matrix row of zeros, and unsolvable equation";
throw "Asking for matrix entry at inactive zeta, causing matrix row of zeros, and unsolvable equation";
}
}*/

for (auto const &[neigPos, neigNodeID]: horizontal_neighbours) {
zetaMovementConductance = 0 * (si::square_meter / day);
if (at(neigNodeID)->isZetaActive(zetaID)) {
if (at(neigNodeID)->isZetaTZeroActive(zetaID)) {
// on left hand side of equation -> need updated zeta heights
zoneThicknesses = calculateZoneThicknesses(neigPos, neigNodeID, getZetas(),
at(neigNodeID)->getZetas());
Expand All @@ -3002,9 +3011,7 @@ Calculate
zetaMovementConductance += delnus[zetaID] * zoneConductanceCum; // in SWI2: SWISOLCC/R
NANChecker(zetaMovementConductance.value(), "zetaMovementConductance");
// add conductance to out, the key in the unordered map is the ID of the neighbouring node
//if (std::abs(zetaMovementConductance.value()) > 1e-20) {
out[neigNodeID] = zetaMovementConductance;
//}
out[neigNodeID] = zetaMovementConductance;
}
}

Expand Down
51 changes: 27 additions & 24 deletions src/Simulation/Options.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ namespace GlobalFlow {
specificstorage_from_file = data_config.get<bool>("specific_storage_from_file");
specificyield_from_file = data_config.get<bool>("specific_yield_from_file");
k_river_from_file = data_config.get<bool>("k_river_from_file");
aquifer_depth_from_file = data_config.get<bool>("aquifer_depth_from_file");
vertical_size_as_array = data_config.get<bool>("vertical_size_as_array");
eq_wtd_from_file = data_config.get<bool>("eq_wtd_from_file");
initial_head_from_file = data_config.get<bool>("initial_head_from_file");
effective_porosity_from_file = data_config.get<bool>("effective_porosity_from_file");
Expand All @@ -102,7 +102,7 @@ namespace GlobalFlow {
GHB_K = default_data.get<double>("ghb_K");
RIVER_CONDUCTIVITY = default_data.get<double>("river_conductivity");
SWB_ELEVATION_FACTOR = default_data.get<double>("swb_elevation_factor");
AQUIFER_DEPTH = getTypeArray<int>("aquifer_thickness", default_data);
VERTICAL_SIZES = getTypeArray<double>("vertical_size", default_data);

ANISOTROPY = getTypeArray<double>("anisotropy", default_data);
SPECIFIC_YIELD = default_data.get<double>("specific_yield");
Expand All @@ -122,43 +122,46 @@ namespace GlobalFlow {

pt::ptree data = input.get_child("data");

bool efoldAsArray = data_config.get<bool>("efold_as_array");
if (efoldAsArray){
efold_as_array = data_config.get<bool>("efold_as_array");
if (efold_as_array){
EFOLDING_a = getArray("E-Folding", data.get_child("e-folding"));
}
EFOLDING = getOptional("e-folding", data);
EFOLDING_DIR = getOptional("e-folding", data);

INITIAL_ZETAS_AS_ARRAY = data_config.get<bool>("initial_zetas_as_array");
if (INITIAL_ZETAS_AS_ARRAY){
initial_zetas_as_array = data_config.get<bool>("initial_zetas_as_array");
if (initial_zetas_as_array){
INITIAL_ZETAS_a = getArray("Zetas", data.get_child("zetas"));
}
INITIAL_ZETAS = getOptional("zetas", data);
INITIAL_ZETAS_DIR = getOptional("zetas", data);

ELEVATION = getOptional("elevation", data);
EQUAL_WATER_TABLE_DEPTH = getOptional("equal_water_table_depth", data);
RIVER_ELEVATION = getOptional("river_elevation", data);
vertical_size_as_array = data_config.get<bool>("vertical_size_as_array");
if (vertical_size_as_array){
VERTICAL_SIZE_a = getArray("VerticalSize", data.get_child("vertical_size"));
}
VERTICAL_SIZE_DIR = getOptional("vertical_size", data);

ELEVATION_FILE = getOptional("elevation", data);
EQUAL_WATER_TABLE_DEPTH_FILE = getOptional("equal_water_table_depth", data);
RIVER_ELEVATION_FILE = getOptional("river_elevation", data);

LITHOLOGY = getOptional("lithology", data);
RECHARGE = getOptional("recharge", data);
LITHOLOGY_FILE = getOptional("lithology", data);
RECHARGE_FILE = getOptional("recharge", data);
ZONES_SOURCES_FILE = getOptional("zones_sources", data);
RIVER = getOptional("river_extent", data);
GLOBAL_WETLANDS = getOptional("global_wetlands", data);
GLOBAL_LAKES = getOptional("global_lakes", data);
LOCAL_LAKES = getOptional("local_lakes", data);
LOCAL_WETLANDS = getOptional("local_wetlands", data);
RIVER_FILE = getOptional("river_extent", data);
GLOBAL_WETLANDS_FILE = getOptional("global_wetlands", data);
GLOBAL_LAKES_FILE = getOptional("global_lakes", data);
LOCAL_LAKES_FILE = getOptional("local_lakes", data);
LOCAL_WETLANDS_FILE = getOptional("local_wetlands", data);

//Optional
K_DIR = getOptional("conductance", data);
RIVER_K = getOptional("river_conductance", data);
GHB_DIR = getOptional("ghb", data);
K_FILE = getOptional("conductance", data);
RIVER_K_FILE = getOptional("river_conductance", data);
GHB_FILE = getOptional("ghb", data);
SS_FILE = getOptional("specific_storage", data);
SY_FILE = getOptional("specific_yield", data);
AQ_DEPTH = getOptional("aquifer_depth", data);

INITIAL_HEAD_FILE = getOptional("initial_head", data);

INITIAL_ZONES = getOptional("initial_zones", data);

EFFECTIVE_POROSITY_FILE = getOptional("effective_porosity", data);

SPATID_ARCID = getOptional("spatID-arcID", data);
Expand Down
Loading

0 comments on commit 1040033

Please sign in to comment.