Homec4science
Diffusion SpecMiCP / ReactMiCP (bin/unsaturated) Browse

specmicp/examples/reactmicp/saturated_react/momas_benchmark.cpp
2553b8ca7134bin/unsaturated

Last Change
Raw File

momas_benchmark.cpp
View Options

#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/problem_solver/dissolver.hpp"
#include "specmicp/problem_solver/formulation.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp/adimensional/equilibrium_curve.hpp"
#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
#include "reactmicp/solver/reactive_transport_solver.hpp"
#include "reactmicp/systems/saturated_react/variables.hpp"
#include "reactmicp/systems/saturated_react/transport_stagger.hpp"
#include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp"
#include "reactmicp/systems/saturated_react/init_variables.hpp"
#include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp"
#include "dfpm/meshes/axisymmetric_mesh1d.hpp"
#include "dfpm/meshes/uniform_mesh1d.hpp"
#include "utils/log.hpp"
#include <iostream>
#include <fstream>
#include <chrono>
#include <ctime>
#include "database/database.hpp"
#include "database/aqueous_selector.hpp"
#include "reactmicp/solver/timestepper.hpp"
#include "utils/io/reactive_transport.hpp"
#include "utils/timer.hpp"
// parameters
#define DX 0.01
#define IS_ADVECTIVE true
#define MIN_DT 0.01
#define MAX_DT 100.0
#define RESTART_DT MIN_DT
using namespace specmicp;
RawDatabasePtr get_momas_db()
{
database::Database thedatabase("../data/momas_benchmark.js");
return thedatabase.get_database();
}
void prepare_db_for_easy_case(RawDatabasePtr raw_data)
{
database::Database db_handler(raw_data);
db_handler.remove_components({raw_data->get_id_component("X5"),});
if (raw_data->nb_component() != 6) {
throw std::runtime_error("Not enough or Too many component");
}
db_handler.remove_solid_phases();
if (not raw_data->is_valid()) {
throw std::runtime_error("Database is not valid.");
}
//std::cout << selector.aqueous_label_to_id("C6") << " - " << selector.aqueous_label_to_id("C7") << std::endl;
database::AqueousSelector selector(raw_data);
selector.remove_aqueous({selector.aqueous_label_to_id("C6"), selector.aqueous_label_to_id("C7")});
std::cout << raw_data->aqueous.get_nu_matrix() << std::endl;
if (raw_data->nb_aqueous() != 5) {
throw std::runtime_error("Not enough or Too many aqueous.");
}
if (not raw_data->is_valid()) {
throw std::runtime_error("Database is not valid.");
}
}
AdimensionalSystemSolution easy_material_a(
RawDatabasePtr raw_data,
const AdimensionalSystemConstraints& constraints,
// const AdimensionalSystemSolution& initial_solution,
const AdimensionalSystemSolverOptions& options
)
{
specmicp::Vector variables ;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, true);
if (current_perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized)
throw std::runtime_error("Failed to solve initial solution for material A.");
return solver.get_raw_solution(variables);
}
AdimensionalSystemSolution easy_material_b(
RawDatabasePtr raw_data,
const AdimensionalSystemConstraints& constraints,
const specmicp::AdimensionalSystemSolverOptions& options
)
{
specmicp::Vector variables;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, true);
if (current_perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized)
throw std::runtime_error("Failed to solve initial solution for material B.");
return solver.get_raw_solution(variables);
}
AdimensionalSystemSolution easy_bc(
RawDatabasePtr raw_data,
const AdimensionalSystemConstraints& constraints,
const specmicp::AdimensionalSystemSolverOptions& options
)
{
specmicp::Vector variables;
specmicp::AdimensionalSystemSolver solver(raw_data, constraints, options);
solver.initialise_variables(variables, 1.0, -4.0);
solver.run_pcfm(variables);
specmicp::micpsolver::MiCPPerformance current_perf = solver.solve(variables, false);
if (current_perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized)
throw std::runtime_error("Failed to solve initial solution for BC.");
return solver.get_raw_solution(variables);
}
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::solver;
using namespace specmicp::reactmicp::systems::satdiff;
class MoMasUspcaling: public solver::UpscalingStaggerBase
{
public:
MoMasUspcaling(
RawDatabasePtr the_database,
units::UnitsSet the_units,
bool advective
):
m_data(the_database),
m_units_set(the_units),
m_advective(advective)
{
}
//! \brief Initialize the stagger at the beginning of the computation
void initialize(VariablesBasePtr var)
{
SaturatedVariablesPtr true_var = cast_var_from_base(var);
database::Database db_handler(m_data);
true_var->upscaling_variables().setZero();
for (index_t node=0; node<true_var->nb_nodes(); ++node)
{
// if material b
const scalar_t coord = true_var->get_mesh()->get_position(node);
if (coord < 1.1 and coord >= 1.0)
{
true_var->porosity(node) = 0.5;
true_var->permeability(node) = 1e-5;
if (m_advective)
true_var->diffusion_coefficient(node) = 3.3e-4;
else
true_var->diffusion_coefficient(node) = 330e-3;
true_var->fluid_velocity(node) = 5.5e-3;
}
else // material a
{
true_var->porosity(node) = 0.25;
true_var->permeability(node) = 1e-2;
if (m_advective)
true_var->diffusion_coefficient(node) = 5.5e-5;
else
true_var->diffusion_coefficient(node) = 55e-3;
true_var->fluid_velocity(node) = 5.5e-3;
}
std::cout << "node : " << node << " - porosity : " << true_var->porosity(node) << std::endl;
}
}
//! \brief Initialize the stagger at the beginning of an iteration
void initialize_timestep(scalar_t dt, VariablesBasePtr var)
{
}
//! \brief Solve the equation for the timestep
StaggerReturnCode restart_timestep(VariablesBasePtr var)
{
return StaggerReturnCode::ResidualMinimized;
}
private:
RawDatabasePtr m_data;
units::UnitsSet m_units_set;
index_t m_id_cshj;
scalar_t m_initial_sat_cshj;
scalar_t m_dt;
bool m_advective;
};
AdimensionalSystemConstraints get_constraints_easy_a(RawDatabasePtr& raw_data)
{
Vector total_concentrations(Vector::Zero(raw_data->nb_component()));
total_concentrations(raw_data->get_id_component("X2")) = -2.0;
total_concentrations(raw_data->get_id_component("X4")) = 2.0;
AdimensionalSystemConstraints constraints_a(0.25*total_concentrations);
constraints_a.disable_conservation_water();
constraints_a.surface_model.model_type = SurfaceEquationType::Equilibrium;
constraints_a.surface_model.concentration = 0.25*1.0;
constraints_a.inert_volume_fraction = 0.75;
return constraints_a;
}
AdimensionalSystemConstraints get_constraints_easy_b(RawDatabasePtr& raw_data)
{
Vector total_concentrations(Vector::Zero(raw_data->nb_component()));
total_concentrations(raw_data->get_id_component("X2")) = -2.0;
total_concentrations(raw_data->get_id_component("X4")) = 2.0;
AdimensionalSystemConstraints constraints_b(0.5*total_concentrations);
constraints_b.disable_conservation_water();
constraints_b.surface_model.model_type = SurfaceEquationType::Equilibrium;
constraints_b.surface_model.concentration = 0.5*10.0;
constraints_b.inert_volume_fraction = 0.5;
return constraints_b;
}
AdimensionalSystemConstraints get_constraints_easy_bc(RawDatabasePtr& raw_data)
{
Vector total_concentrations(Vector::Zero(raw_data->nb_component()));
total_concentrations(raw_data->get_id_component("X1")) = 0.3;
total_concentrations(raw_data->get_id_component("X2")) = 0.3;
total_concentrations(raw_data->get_id_component("X3")) = 0.3;
AdimensionalSystemConstraints constraints_bc(total_concentrations);
constraints_bc.disable_conservation_water();
constraints_bc.surface_model.model_type = SurfaceEquationType::NoEquation;
constraints_bc.surface_model.concentration = 0.0;
constraints_bc.inert_volume_fraction = 0.0;
return constraints_bc;
}
AdimensionalSystemSolverOptions get_specmicp_options()
{
AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 1.0;
options.solver_options.max_iter = 200;
options.solver_options.maxiter_maxstep = 200;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = false;
options.solver_options.non_monotone_linesearch = false;
options.solver_options.factor_descent_condition = -1;
options.solver_options.factor_gradient_search_direction = 100;
options.solver_options.projection_min_variable = 1e-9;
options.solver_options.fvectol = 1e-10;
options.solver_options.steptol = 1e-14;
options.solver_options.threshold_cycling_linesearch = 1e-8;
options.system_options.non_ideality_tolerance = 1e-10;
options.system_options.non_ideality = false;
options.system_options.restart_concentration = -2;
options.use_pcfm = true;
return options;
}
mesh::Mesh1DPtr get_mesh(scalar_t dx)
{
scalar_t tot_length = 2.1+dx;
index_t nb_nodes = tot_length/dx +2;
return mesh::uniform_mesh1d(nb_nodes, dx, 1.0);
}
void set_transport_options(dfpmsolver::ParabolicDriverOptions& transport_options)
{
transport_options.maximum_iterations = 450;
transport_options.quasi_newton = 3;
transport_options.residuals_tolerance = 1e-8;
transport_options.step_tolerance = 1e-12;
transport_options.absolute_tolerance = 1e-6;
transport_options.alpha = 1.0;
//transport_options.linesearch = dfpmsolver::ParabolicLinesearch::Strang;
transport_options.sparse_solver = SparseSolver::SparseLU;
//transport_options.sparse_solver = SparseSolver::GMRES;
transport_options.threshold_stationary_point = 1e-2;
}
void set_reactive_solver_options(reactmicp::solver::ReactiveTransportOptions& reactive_options)
{
reactive_options.maximum_iterations = 50; //500;
reactive_options.residuals_tolerance = 1e-2;
reactive_options.good_enough_tolerance = 1.0;
reactive_options.absolute_residuals_tolerance = 1e-6;
reactive_options.implicit_upscaling = false;
}
int main()
{
Eigen::initParallel();
Timer global_timer;
global_timer.start();
std::ofstream output_gen("out_momas_out.txt");
io::print_reactmicp_header(&output_gen);
specmicp::logger::ErrFile::stream() = &output_gen; //&std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
RawDatabasePtr database = get_momas_db();
prepare_db_for_easy_case(database);
units::UnitsSet the_units = units::UnitsSet();
the_units.length = units::LengthUnit::decimeter; // so that rho_w ~ 1
scalar_t dx=DX;
mesh::Mesh1DPtr the_mesh = get_mesh(dx);
AdimensionalSystemConstraints constraints_a = get_constraints_easy_a(database);
AdimensionalSystemConstraints constraints_b = get_constraints_easy_b(database);
AdimensionalSystemConstraints constraints_bc = get_constraints_easy_bc(database);
AdimensionalSystemSolverOptions options = get_specmicp_options();
options.units_set = the_units;
AdimensionalSystemSolution mat_bc = easy_bc(database, constraints_bc, options);
AdimensionalSystemSolution mat_a = easy_material_a(database, constraints_a, options);
AdimensionalSystemSolution mat_b = easy_material_b(database, constraints_b, options);
options.solver_options.maxstep = 1.0;
options.solver_options.use_scaling = true;
options.solver_options.non_monotone_linesearch = true;
index_t nb_nodes = the_mesh->nb_nodes();
std::vector<specmicp::AdimensionalSystemConstraints> list_constraints = {constraints_a,
constraints_b,
constraints_bc};
std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition = {mat_a, mat_b, mat_bc};
std::vector<int> index_initial_state(nb_nodes, 0);
std::vector<int> index_constraints(nb_nodes, 0);
index_t start_mat_b = 1.0/dx;
index_t end_mat_b = 1.1/dx;
for (index_t node=start_mat_b; node<end_mat_b; ++node)
{
index_initial_state[node] = 1;
index_constraints[node] = 1;
}
index_initial_state[0] = 2;
index_constraints[0] = 2;
std::vector<specmicp::index_t> list_fixed_nodes = {0, };
std::map<index_t, index_t> list_slaves_nodes = {{nb_nodes-1, nb_nodes-2}};
std::vector<index_t> list_immobile_components = {0, 1};
systems::satdiff::SaturatedVariablesPtr variables =
systems::satdiff::init_variables(the_mesh, database, the_units,
list_fixed_nodes,
list_initial_composition,
index_initial_state);
ChemistryStaggerPtr equilibrium_stagger =
std::make_shared<reactmicp::systems::satdiff::EquilibriumStagger>(list_constraints, index_constraints, options);
std::shared_ptr<reactmicp::systems::satdiff::SaturatedTransportStagger> transport_stagger =
std::make_shared<reactmicp::systems::satdiff::SaturatedTransportStagger>(
variables, list_fixed_nodes, list_slaves_nodes, list_immobile_components);
set_transport_options(transport_stagger->options_solver());
const bool is_advective = IS_ADVECTIVE;
UpscalingStaggerPtr upscaling_stagger =
std::make_shared<MoMasUspcaling>(database, the_units, is_advective);
ReactiveTransportSolver solver(transport_stagger, equilibrium_stagger, upscaling_stagger);
transport_stagger->initialize(variables);
equilibrium_stagger->initialize(variables);
upscaling_stagger->initialize(variables);
set_reactive_solver_options(solver.get_options());
// mesh output
std::ofstream output_mesh;
output_mesh.open("out_momas_mesh.dat");
output_mesh << "Node\tPosition" << std::endl;
for (index_t node: the_mesh->range_nodes())
{
output_mesh << node << "\t" << the_mesh->get_position(node)-dx/2.0 << std::endl;
}
output_mesh.close();
//std::cout << variables->displacement() << std::endl;
std::ofstream output_iter("out_momas_iter.dat");
std::ofstream output_x1, output_x2, output_x3, output_x4;
std::ofstream output_sx1, output_sx2, output_sx3, output_sx4;
std::ofstream output_c1, output_c2, output_s;
output_x1.open("out_momas_x1.dat");
output_x2.open("out_momas_x2.dat");
output_x3.open("out_momas_x3.dat");
output_x4.open("out_momas_x4.dat");
output_sx1.open("out_momas_sx1.dat");
output_sx2.open("out_momas_sx2.dat");
output_sx3.open("out_momas_sx3.dat");
output_sx4.open("out_momas_sx4.dat");
output_c1.open("out_momas_c1.dat");
output_c2.open("out_momas_c2.dat");
output_s.open("out_momas_s.dat");
index_t cnt = 0;
scalar_t dt=RESTART_DT;
reactmicp::solver::Timestepper timestepper(MIN_DT, MAX_DT, 5000, 2.0);
timestepper.get_options().alpha_average = 0.3;
timestepper.get_options().decrease_failure = 0.1;
timestepper.get_options().increase_factor = 1.05;
timestepper.get_options().decrease_factor = 0.50;
timestepper.get_options().iteration_lower_target = 1.001;
io::print_reactmicp_performance_long_header(&output_iter);
while (timestepper.get_total() < timestepper.get_total_target())
{
output_gen << "dt : " << dt << std::endl;
Timer step_timer;
step_timer.start();
reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
step_timer.stop();
// if (retcode > reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
//{
io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs());
//}
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
{
dt = RESTART_DT;
variables->reset_main_variables();
retcode = solver.solve_timestep(dt, variables);
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total(), solver.get_perfs());
}
output_gen << "Time : " << timestepper.get_total() << std::endl;
io::print_reactmicp_performance_short(&output_gen, solver.get_perfs());
scalar_t total = timestepper.get_total();
if (cnt % 10== 0)
{
output_x1 << total;
output_x2 << total;
output_x3 << total;
output_x4 << total;
output_sx1 << total;
output_sx2 << total;
output_sx3 << total;
output_sx4 << total;
output_c1 << total;
output_c2 << total;
output_s << total;
for (index_t node: the_mesh->range_nodes())
{
output_x1 << "\t" << variables->aqueous_concentration(node, 1, variables->displacement());
output_x2 << "\t" << variables->aqueous_concentration(node, 2, variables->displacement());
output_x3 << "\t" << variables->aqueous_concentration(node, 3, variables->displacement());
output_x4 << "\t" << variables->aqueous_concentration(node, 4, variables->displacement());
output_sx1 << "\t" << variables->solid_concentration(node, 1, variables->displacement());
output_sx2 << "\t" << variables->solid_concentration(node, 2, variables->displacement());
output_sx3 << "\t" << variables->solid_concentration(node, 3, variables->displacement());
output_sx4 << "\t" << variables->solid_concentration(node, 4, variables->displacement());
output_c1 << "\t" << variables->equilibrium_solution(node).secondary_molalities(0);
output_c2 << "\t" << variables->equilibrium_solution(node).secondary_molalities(1);
output_s << "\t" << pow10(variables->equilibrium_solution(node).main_variables(5));
}
output_x1 << std::endl;
output_x2 << std::endl;
output_x3 << std::endl;
output_x4 << std::endl;
output_sx1 << std::endl;
output_sx2 << std::endl;
output_sx3 << std::endl;
output_sx4 << std::endl;
output_c1 << std::endl;
output_c2 << std::endl;
output_s << std::endl;
}
++ cnt;
}
variables->displacement()(database->get_id_component("X1")) = 0;
variables->displacement()(database->get_id_component("X2")) = -2.0;
variables->displacement()(database->get_id_component("X3")) = 0;
variables->displacement()(database->get_id_component("X4")) = 2.0;
timestepper.set_total_target(6000);
dt = RESTART_DT;
transport_stagger->options_solver().absolute_tolerance = 1e-10;
solver.get_options().absolute_residuals_tolerance = 1e-6;
while (timestepper.get_total() < timestepper.get_total_target())
{
output_gen << "dt : " << dt << std::endl;
Timer step_timer;
step_timer.start();
reactmicp::solver::ReactiveTransportReturnCode retcode = solver.solve_timestep(dt, variables);
step_timer.stop();
// if (retcode > reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
//{
io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total()+dt, solver.get_perfs());
//}
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
if (retcode <= reactmicp::solver::ReactiveTransportReturnCode::NotConvergedYet)
{
dt = RESTART_DT;
variables->reset_main_variables();
retcode = solver.solve_timestep(dt, variables);
dt = timestepper.next_timestep(dt, retcode, solver.get_perfs().nb_iterations);
io::print_reactmicp_performance_long(&output_iter, cnt, timestepper.get_total(), solver.get_perfs());
}
output_gen << "Time : " << timestepper.get_total() << std::endl;
io::print_reactmicp_performance_short(&output_gen, solver.get_perfs());
scalar_t total = timestepper.get_total();
if (cnt % 10== 0)
{
output_x1 << total;
output_x2 << total;
output_x3 << total;
output_x4 << total;
output_sx1 << total;
output_sx2 << total;
output_sx3 << total;
output_sx4 << total;
output_c1 << total;
output_c2 << total;
output_s << total;
for (index_t node: the_mesh->range_nodes())
{
output_x1 << "\t" << variables->aqueous_concentration(node, 1, variables->displacement());
output_x2 << "\t" << variables->aqueous_concentration(node, 2, variables->displacement());
output_x3 << "\t" << variables->aqueous_concentration(node, 3, variables->displacement());
output_x4 << "\t" << variables->aqueous_concentration(node, 4, variables->displacement());
output_sx1 << "\t" << variables->solid_concentration(node, 1, variables->displacement());
output_sx2 << "\t" << variables->solid_concentration(node, 2, variables->displacement());
output_sx3 << "\t" << variables->solid_concentration(node, 3, variables->displacement());
output_sx4 << "\t" << variables->solid_concentration(node, 4, variables->displacement());
output_c1 << "\t" << variables->equilibrium_solution(node).secondary_molalities(0);
output_c2 << "\t" << variables->equilibrium_solution(node).secondary_molalities(1);
output_s << "\t" << pow10(variables->equilibrium_solution(node).main_variables(5));
}
output_x1 << std::endl;
output_x2 << std::endl;
output_x3 << std::endl;
output_x4 << std::endl;
output_sx1 << std::endl;
output_sx2 << std::endl;
output_sx3 << std::endl;
output_sx4 << std::endl;
output_c1 << std::endl;
output_c2 << std::endl;
output_s << std::endl;
}
++ cnt;
}
reactmicp::solver::ReactiveTransportTimer timer = solver.get_timer();
io::print_reactmicp_timer(&output_gen, timer);
global_timer.stop();
output_gen << "Total computation time : " << global_timer.elapsed_time() << "s." << std::endl;
return EXIT_SUCCESS;
}
c4science · Help