Page MenuHomec4science
Files F86222427

unsaturated.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Oct 5, 01:18

unsaturated.cpp
View Options

/* =============================================================================
Copyright (c) 2014 - 2016
F. Georget <fabieng@princeton.edu> Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
3. Neither the name of the copyright holder nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
============================================================================= */
//! \file bin/reactmicp/unsaturated.cpp
//! \brief Unsaturated
#include "specmicp_common/types.hpp"
#include "specmicp_common/cli/parser.hpp"
#include "specmicp_common/plugins/plugin_manager.hpp"
#include "specmicp_common/io/yaml.hpp"
#include "specmicp_common/io/safe_config.hpp"
#include "specmicp_common/compat.hpp"
#include "specmicp_common/log.hpp"
#include "specmicp_common/filesystem.hpp"
#include "specmicp_common/string_algorithms.hpp"
#include "specmicp_common/physics/io/configuration.hpp"
#include "specmicp_database/io/configuration.hpp"
#include "dfpm/io/configuration.hpp"
#include "specmicp_common/io/configuration.hpp"
#include "reactmicp/reactmicp_unsaturated.hpp"
#include "reactmicp/systems/unsaturated/variables.hpp"
#include "dfpm/solver/parabolic_structs.hpp"
#include "reactmicp/io/hdf5_unsaturated.hpp"
#include <string>
#include <memory>
#include <stdexcept>
#include <iostream>
#include "specmicp_common/io/config_yaml_sections.h"
#define PLUGIN_MODULE_UPSCALING "upscaling_factory"
#define PLUGIN_MODULE_VARIABLES "variables_factory"
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems::unsaturated;
const char* welcome =
R"(===============================
= Reactmicp - unsaturated =
===============================
(C) copyright 2014-2016 F.Georget, Princeton University
version: not good enough yet
)";
int main (int argc, char* argv[]);
plugins::PluginManager& initialize_plugin_manager(
io::YAMLConfigHandle& config
);
SimulationData initialize_simul_data(
io::YAMLConfigHandle& safe_config,
const std::vector<std::string>& db_dirs
);
std::shared_ptr<UnsaturatedVariables> initialize_variables(
std::shared_ptr<database::DataContainer> raw_data,
std::shared_ptr<mesh::Mesh1D> the_mesh,
const units::UnitsSet& the_units,
io::YAMLConfigHandle&& configuration
);
struct Staggers
{
std::shared_ptr<solver::UpscalingStaggerBase> upscaling;
std::shared_ptr<solver::ChemistryStaggerBase> chemistry;
std::shared_ptr<solver::TransportStaggerBase> transport;
};
Staggers initialize_staggers(
SimulationData& simul_data,
io::YAMLConfigHandle&& configuration
);
std::shared_ptr<solver::UpscalingStaggerBase>
initialize_upscaling_stagger(
SimulationData& simul_data,
io::YAMLConfigHandle& conf_staggers
);
std::shared_ptr<solver::ChemistryStaggerBase>
initialize_chemistry_stagger(
SimulationData& simul_data,
io::YAMLConfigHandle& conf_staggers
);
std::shared_ptr<solver::TransportStaggerBase>
initialize_transport_stagger(
SimulationData& simul_data,
io::YAMLConfigHandle& conf_staggers
);
std::string initialize_output(
io::YAMLConfigHandle&& conf
);
int main(int argc, char *argv[])
{
std::cout << welcome << std::endl;
init_logger(&std::cout, logger::Warning);
// CLI options
cli::CommandLineParser parser;
parser.add_option('i', "input",
cli::ValueType::string,
"input_file (yaml format)");
parser.add_option('d', "database_dir",
std::string(""),
"directory where database are stored");
parser.set_help_message("ReactMiCP Unsaturated System");
const auto ret_parse = parser.parse(argc, argv);
if (ret_parse > 0) {
return EXIT_SUCCESS;
}
std::string input_filepath =
parser.get_option<cli::ValueType::string>("input");
std::string database_dir =
parser.get_option<cli::ValueType::string>("database_dir");
std::vector<std::string> db_dirs;
const auto env = utils::get_env("SPECMICP_DATABASE_PATH");
if (not env.empty()) {
db_dirs = utils::split(env, ';');
}
if (not database_dir.empty()) {db_dirs.push_back(database_dir);}
// Configuration File
// ------------------
auto config = io::YAMLConfigFile::make(input_filepath);
initialize_plugin_manager(*(config.get()));
// Eigen initialization
// --------------------
#if SPECMICP_USE_OPENMP
Eigen::initParallel();
#endif
// Simulation Data
// ---------------
SimulationData simul_data = initialize_simul_data(*config.get(), db_dirs);
// Staggers
// -------
Staggers staggers = initialize_staggers(
simul_data,
config->get_section(SPC_CF_S_STAGGERS)
);
// Check the variables
// -------------------
auto res = VariablesInterface(simul_data.variables).check_variables();
if (res != VariablesValidity::good) {
const auto msg = "Variables were not initialized correctly."
"Please see log for errors.";
ERROR << msg;
throw std::invalid_argument(msg);
}
// Solver
// ------
solver::ReactiveTransportSolver react_solver(
staggers.transport, staggers.chemistry, staggers.upscaling);
io::configure_reactmicp_options(
react_solver.get_options(),
config->get_section(SPC_CF_S_REACTMICP)
);
// Initialize output
std::unique_ptr<io::UnsaturatedHDF5Saver> saver;
solver::SimulationInformation simul_info =
io::configure_simulation_information(
config->get_section(SPC_CF_S_SIMULINFO));
solver::ReactiveTransportRunner runner(react_solver,
1.0, 10.0, // dummy value, set later
simul_info);
io::configure_reactmicp_timestepper(
runner.get_timestepper_options(),
config->get_section(SPC_CF_S_TIMESTEPPER)
);
//
if (config->has_section(SPC_CF_S_REACTOUTPUT)) {
auto filepath = initialize_output(
config->get_section(SPC_CF_S_REACTOUTPUT));
saver = make_unique<io::UnsaturatedHDF5Saver>(
filepath, simul_data.variables, simul_data.units);
auto* saver_ptr = saver.get();
auto out_pol = [saver_ptr](
scalar_t timestep,
solver::VariablesBasePtr _) {
saver_ptr->save_timestep(timestep);
};
saver->save_timestep(0.0);
runner.set_output_policy(out_pol);
}
// run info
auto run_section = config->get_section(SPC_CF_S_RUN);
auto run_until = run_section.get_required_attribute<scalar_t>(
SPC_CF_S_RUN_A_RUNUTNIL);
// we destroy config to get information
config.reset(nullptr);
runner.run_until(run_until, simul_data.variables);
return EXIT_SUCCESS;
}
// ==================== //
// Simulation data //
// ==================== //
SimulationData initialize_simul_data(
io::YAMLConfigHandle& safe_config,
const std::vector<std::string>& db_dirs
)
{
// Package the information
SimulationData simul_data;
// The units
// ---------
if (safe_config.has_section(SPC_CF_S_UNITS))
{
simul_data.units = io::configure_units(
safe_config.get_section(SPC_CF_S_UNITS));
}
// The database
// ------------
simul_data.raw_data = io::configure_database(
safe_config.get_section(SPC_CF_S_DATABASE),
db_dirs
);
if (not simul_data.raw_data->is_valid()) {
CRITICAL << "Invalid database from the configuration.";
throw std::runtime_error("Database parsed from the configuration has"
" been detected to be invalid. Abort.");
}
// The mesh
// --------
simul_data.mesh1d =
io::configure_mesh(safe_config.get_section(SPC_CF_S_MESH));
// The variables
// -------------
simul_data.variables = initialize_variables(
simul_data.raw_data,
simul_data.mesh1d,
simul_data.units,
safe_config.get_section(SPC_CF_S_REACTMICPVARIABLES)
);
if (not simul_data.raw_data->is_valid()) {
CRITICAL << "Invalid database after variable initialization.";
throw std::runtime_error("Database has been detected to be invalid"
" after the variables initialization. Abort.");
}
// The database is supposed set after this point !
simul_data.raw_data->freeze_db();
// The boundary conditions
// -----------------------
simul_data.bcs =
io::configure_unsaturated_boundary_conditions(
simul_data.mesh1d->nb_nodes(),
simul_data.raw_data.get(),
safe_config.get_section(SPC_CF_S_BOUNDARYCONDITIONS)
);
// Simulation data ready
return simul_data;
}
std::shared_ptr<UnsaturatedVariables> initialize_variables(
std::shared_ptr<database::DataContainer> raw_data,
std::shared_ptr<mesh::Mesh1D> the_mesh,
const units::UnitsSet& the_units,
io::YAMLConfigHandle&& configuration
)
{
auto type = configuration.get_optional_attribute<std::string>(
SPC_CF_S_REACTMICPVARIABLES_A_TYPE, SPC_CF_V_PLUGIN);
if (type == SPC_CF_V_DEFAULT) type = SPC_CF_V_PLUGIN;
if (type != SPC_CF_V_PLUGIN) {
configuration.report_error(
io::YAMLConfigError::InvalidArgument,
"Only accepted value for attribute '"
SPC_CF_S_REACTMICPVARIABLES_A_TYPE
"' is '" SPC_CF_V_PLUGIN "'."
);
}
auto plugin_name = configuration.get_required_attribute<std::string>(
SPC_CF_A_PLUGIN_NAME
);
plugins::PluginManager& plugin_manager = plugins::get_plugin_manager();
auto extra_plugin_file = configuration.get_optional_attribute<std::string>(
SPC_CF_A_PLUGIN_FILE, "");
if (extra_plugin_file != "") {
auto retcode = plugin_manager.load_plugin(extra_plugin_file);
if (not retcode) {
CRITICAL << "Failed to load plugin file : "
<< extra_plugin_file << std::endl;
std::runtime_error("Fail to load plugin file : "
+ extra_plugin_file + ".");
}
}
auto initializer = plugin_manager.get_object<InitializeVariables>(
PLUGIN_MODULE_VARIABLES, plugin_name);
if (initializer == nullptr) {
CRITICAL << "No variable initializer found";
throw std::runtime_error("No variable initializer found");
}
return initializer->initialize_variables(
raw_data, the_mesh, the_units, configuration);
}
// ================ //
// Plugin Manager //
// ================ //
plugins::PluginManager& initialize_plugin_manager(
io::YAMLConfigHandle& config
)
{
plugins::PluginManager& plugin_manager = plugins::get_plugin_manager();
plugin_manager.register_module(PLUGIN_MODULE_VARIABLES,
make_unique<plugins::ModuleBase>());
plugin_manager.register_module(PLUGIN_MODULE_UPSCALING,
make_unique<plugins::ModuleBase>());
if (config.has_section(SPC_CF_S_PLUGINS)) {
io::configure_plugin_manager(config.get_section(SPC_CF_S_PLUGINS));
}
return plugin_manager;
}
// ================ //
// Staggers //
// ================ //
Staggers initialize_staggers(
SimulationData& simul_data,
io::YAMLConfigHandle&& conf_staggers
)
{
Staggers staggers;
auto* var_ptr = simul_data.variables.get();
staggers.upscaling =
initialize_upscaling_stagger(simul_data, conf_staggers);
staggers.upscaling->initialize(var_ptr);
staggers.chemistry =
initialize_chemistry_stagger(simul_data, conf_staggers);
staggers.chemistry->initialize(var_ptr);
staggers.transport =
initialize_transport_stagger(simul_data, conf_staggers);
staggers.transport->initialize(var_ptr);
// check the database
if (not simul_data.raw_data->is_valid()) {
CRITICAL << "Change in database during stagger initialization !";
throw std::logic_error("Change in the database has been detected"
" during the staggers initialization. Abort.");
}
return staggers;
}
std::shared_ptr<solver::UpscalingStaggerBase>
initialize_upscaling_stagger(
SimulationData& simul_data,
io::YAMLConfigHandle& config_staggers
)
{
auto conf = config_staggers.get_section(SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER);
auto type = conf.get_optional_attribute<std::string>(
SPC_CF_S_STAGGERS_SS_UPSCALINGSTAGGER_A_TYPE,
SPC_CF_V_PLUGIN);
// only plugin accepted for now
if (type == SPC_CF_V_DEFAULT) type = SPC_CF_V_PLUGIN;
if (type != SPC_CF_V_PLUGIN) {
conf.report_error(
io::YAMLConfigError::InvalidArgument,
"Invalid argument for attribute '"
SPC_CF_A_TYPE
"' only value accepted is '" SPC_CF_V_PLUGIN "'.");
}
auto plugin_name = conf.get_required_attribute<std::string>(
SPC_CF_A_PLUGIN_NAME);
plugins::PluginManager& plugin_manager = plugins::get_plugin_manager();
auto extra_plugin_file = conf.get_optional_attribute<std::string>(
SPC_CF_A_PLUGIN_FILE, "");
if (extra_plugin_file != "") {
auto retcode = plugin_manager.load_plugin(extra_plugin_file);
if (not retcode) {
CRITICAL << "Failed to load plugin file : " << extra_plugin_file << std::endl;
std::runtime_error("Fail to load plugin file : "
+ extra_plugin_file + ".");
}
}
auto upscaling_factory = plugin_manager.get_object<UpscalingStaggerFactory>(
PLUGIN_MODULE_UPSCALING, plugin_name);
if (upscaling_factory == nullptr) {
CRITICAL << "No upscaling stagger factory found";
throw std::runtime_error("No upscaling stagger factory found");
}
return upscaling_factory->get_upscaling_stagger(
simul_data, std::move(conf));
}
std::shared_ptr<solver::ChemistryStaggerBase>
initialize_chemistry_stagger(
SimulationData& simul_data,
io::YAMLConfigHandle& conf_staggers
)
{
// no conf
// -------
if (not conf_staggers.has_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER)) {
auto opts = EquilibriumOptionsVector::make(
simul_data.mesh1d->nb_nodes());
return EquilibriumStagger::make(
simul_data.variables,
simul_data.bcs,
opts
);
}
// conf is given
// --------------
auto conf_chem = conf_staggers.get_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER);
auto type = conf_chem.get_optional_attribute<std::string>(
SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE,
SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM);
if (type == SPC_CF_V_DEFAULT)
type = SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM;
if (type !=SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_A_TYPE_V_EQUILIBRIUM) {
CRITICAL << "Unknown chemistry stagger";
throw std::invalid_argument("Only equilibrium stagger supported.");
}
std::shared_ptr<EquilibriumOptionsVector> opts {nullptr};
if (conf_chem.has_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS)) {
opts = io::configure_unsaturated_equilibrium_options(
simul_data.mesh1d->nb_nodes(),
simul_data.units,
conf_chem.get_section(SPC_CF_S_STAGGERS_SS_CHEMISTRYSTAGGER_SS_EQUILIBRIUM_OPTS)
);
} else {
opts = EquilibriumOptionsVector::make(simul_data.mesh1d->nb_nodes());
}
return EquilibriumStagger::make(
simul_data.variables,
simul_data.bcs,
opts
);
}
std::shared_ptr<solver::TransportStaggerBase> initialize_transport_stagger(
SimulationData& simul_data,
io::YAMLConfigHandle& conf_staggers
)
{
const UnsaturatedVariables* const variables = simul_data.variables.get();
const database::DataContainer* const raw_data = simul_data.raw_data.get();
// No conf
// -------
if (not conf_staggers.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER)) {
return UnsaturatedTransportStagger::make(
simul_data.variables,
simul_data.bcs,
true, true);
}
// Conf is given
// -------------
auto transport_conf = conf_staggers.get_section(
SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER);
const bool merge_sat = transport_conf.get_optional_attribute<bool>(
SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_SATURATION,
true);
const bool merge_aq = transport_conf.get_optional_attribute<bool>(
SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_A_MERGE_AQUEOUS,
true);
auto stagger = UnsaturatedTransportStagger::make(
simul_data.variables,
simul_data.bcs,
merge_sat, merge_aq);
// ! scaling => done in variables
// default options
// ---------------
dfpmsolver::ParabolicDriverOptions default_opts;
bool rewrite_default = false;
if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS))
{
io::configure_transport_options(
default_opts,
transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_DEFAULT_OPTS)
);
rewrite_default = true;
}
// Saturation option
// -----------------
auto& sat_opts = *(stagger->get_saturation_options());
sat_opts = default_opts;
if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS))
{
io::configure_transport_options(
sat_opts,
transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_SATURATION_OPTS)
);
}
// Aqueous options
// ---------------
if (transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS))
{
auto conf = transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_AQUEOUS_OPTS);
if (not conf.is_sequence()) {
// config for every aqueous equations
dfpmsolver::ParabolicDriverOptions opts;
if (rewrite_default) opts = default_opts;
io::configure_transport_options(
opts, std::move(conf));
for (index_t aq: raw_data->range_aqueous_component()) {
*(stagger->get_aqueous_options(aq)) = opts;
}
}
else {
// First set default is needed
if (rewrite_default) {
for (index_t aq: simul_data.raw_data->range_aqueous_component()) {
*(stagger->get_aqueous_options(aq)) = default_opts;
}
}
uindex_t size = conf.size();
for (uindex_t ind=0; ind<size; ++ind)
{
auto subconf = conf.get_section(ind);
// personalized options per component
auto comp_label = subconf.get_required_attribute<std::string>(
SPC_CF_A_COMPONENT);
auto comp_id = raw_data->get_id_component(comp_label);
io::configure_transport_options(
*(stagger->get_aqueous_options(comp_id)),
std::move(subconf)
);
}
}
}
// Gas options
// -----------
if ( transport_conf.has_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS)
and not (merge_aq and merge_sat)) // skip conf if not needed
{
auto conf = transport_conf.get_section(SPC_CF_S_STAGGERS_SS_TRANSPORTSTAGGER_SS_GAS_OPTS);
if (not conf.is_sequence()) {
// config for every gas equations
dfpmsolver::ParabolicDriverOptions opts;
if (rewrite_default) opts = default_opts;
io::configure_transport_options(
opts, std::move(conf));
if ( not merge_sat and variables->component_has_gas(0)) {
*(stagger->get_gas_options(0)) = opts;
}
if ( not merge_aq) {
for (index_t aq: raw_data->range_aqueous_component()) {
if (simul_data.variables->component_has_gas(aq) )
*(stagger->get_aqueous_options(aq)) = opts;
}
}
}
else
{
// First rewrite default
if (rewrite_default) {
// H2O
if (not merge_sat and variables->component_has_gas(0)) {
*(stagger->get_gas_options(0)) = default_opts;
}
// aqueous components
if (not merge_aq)
{
for (index_t aq: raw_data->range_aqueous_component()) {
if (simul_data.variables->component_has_gas(aq))
*(stagger->get_aqueous_options(aq)) = default_opts;
}
}
}
// Then go over all conf
uindex_t size = conf.size();
for (uindex_t ind=0; ind<size; ++ind)
{
auto subconf = conf.get_section(ind);
auto comp_label = subconf.get_required_attribute<std::string>(
SPC_CF_A_COMPONENT);
auto comp_id = raw_data->get_id_component(comp_label);
// H2O
if (comp_id == 0) {
if (merge_sat) continue;
io::configure_transport_options(
*(stagger->get_gas_options(0)),
std::move(subconf)
);
}
// aqueous component
else {
if (merge_aq) continue;
io::configure_transport_options(
*(stagger->get_gas_options(comp_id)),
std::move(subconf)
);
}
}
}
}
return stagger;
}
std::string initialize_output(
io::YAMLConfigHandle&& conf
)
{
auto type = conf.get_optional_attribute<std::string>(
SPC_CF_S_REACTOUTPUT_A_TYPE, SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5);
if (type == SPC_CF_V_DEFAULT) type = SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5;
if (type != SPC_CF_S_REACTOUTPUT_A_TYPE_V_HDF5) {
conf.report_error(
io::YAMLConfigError::InvalidArgument,
"Only hdf5 is accepted for attribute "
SPC_CF_S_REACTOUTPUT_A_TYPE " for now.");
}
auto filepath = conf.get_required_attribute<std::string>(
SPC_CF_S_REACTOUTPUT_A_FILEPATH);
return filepath;
}

Event Timeline

c4science · Help