Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90849280
momas_benchmark.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Nov 5, 08:01
Size
22 KB
Mime Type
text/x-c++
Expires
Thu, Nov 7, 08:01 (2 d)
Engine
blob
Format
Raw Data
Handle
22143269
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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/selector.hpp"
#include "reactmicp/solver/timestepper.hpp"
#include "utils/io/reactive_transport.hpp"
#include "utils/timer.hpp"
// parameters
#define DX 0.005
#define IS_ADVECTIVE true
#define MIN_DT 0.0001
#define MAX_DT 10.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::DatabaseSelector selector(raw_data);
selector.remove_component({selector.component_label_to_id("X5"),});
specmicp_assert(raw_data->nb_component == 5);
selector.remove_all_minerals();
//std::cout << selector.aqueous_label_to_id("C6") << " - " << selector.aqueous_label_to_id("C7") << std::endl;
selector.remove_aqueous({selector.aqueous_label_to_id("C6"), selector.aqueous_label_to_id("C7")});
//specmicp_assert(raw_data->nb_aqueous == 5);
}
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()
{
Vector total_concentrations(5);
total_concentrations << 1, 0.0, -2.0, 0.0, 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()
{
Vector total_concentrations(5);
total_concentrations << 1, 0.0, -2.0, 0.0, 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()
{
Vector total_concentrations_bc(5);
total_concentrations_bc << 1, 0.3, 0.3, 0.3, 0.0;
AdimensionalSystemConstraints constraints_bc(total_concentrations_bc);
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-12;
options.solver_options.steptol = 1e-16;
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-16;
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()
{
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();
AdimensionalSystemConstraints constraints_b = get_constraints_easy_b();
AdimensionalSystemConstraints constraints_bc = get_constraints_easy_bc();
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};
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);
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()(1) = 0;
variables->displacement()(2) = -2.0;
variables->displacement()(3) = 0;
variables->displacement()(4) = 2.0;
variables->displacement()(6) = 0.0;
variables->displacement()(7) = 0.0;
variables->displacement()(8) = 0.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;
}
Event Timeline
Log In to Comment