Page MenuHomec4science
Files F80128314

momas_benchmark.cpp
No OneTemporary
Actions

File Metadata

Created
Fri, Aug 30, 01:14

momas_benchmark.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. *
============================================================================= */
#include "reactmicp/reactmicp.hpp"
#include "specmicp_database/aqueous_selector.hpp"
#include <iostream>
#include <fstream>
#include <chrono>
#include <ctime>
// parameters
#define DX 0.05
#define IS_ADVECTIVE true
#define MIN_DT 0.001
#define MAX_DT 100.0
#define RESTART_DT MIN_DT
using namespace specmicp;
RawDatabasePtr get_momas_db()
{
database::Database thedatabase("../data/momas_benchmark.yaml");
return thedatabase.get_database();
}
void prepare_db_for_easy_case(RawDatabasePtr raw_data)
{
database::Database db_handler(raw_data);
db_handler.remove_components({"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.initialize_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) {
return initialize(var.get());
}
virtual void initialize(VariablesBase * const var) override
{
SaturatedVariables * const true_var = static_cast<SaturatedVariables * const>(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
virtual void initialize_timestep(scalar_t _, VariablesBase * const __) override
{
}
//! \brief Solve the equation for the timestep
virtual StaggerReturnCode restart_timestep(VariablesBase * const _) override
{
return StaggerReturnCode::ResidualMinimized;
}
private:
RawDatabasePtr m_data;
units::UnitsSet m_units_set;
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)
{
mesh::Uniform1DMeshGeometry geom;
scalar_t tot_length = 2.1+dx;
geom.nb_nodes = tot_length/dx + 2;
geom.dx = dx;
geom.section = 1.0;
return mesh::uniform_mesh1d(geom);
}
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();
io::OutFile output_gen("out_momas_out.txt");
io::print_reactmicp_header(output_gen);
specmicp::logger::ErrFile::stream() = &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
io::print_mesh("out_momas_mesh.dat", the_mesh, the_units.length);
//std::cout << variables->displacement() << std::endl;
io::OutFile 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;
}

Event Timeline

c4science · Help