Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/tests/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
78be35d16fbereactivetransport

Last Change
Raw File

equilibrium_stagger.cpp
View Options

#include "catch.hpp"
#include "speciation_system.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "specmicp/adimensional/adimensional_system_solver_structs.hpp"
#include "dfpm/mesh.hpp"
#include "reactmicp/systems/saturated_react/variables.hpp"
#include "reactmicp/systems/saturated_react/init_variables.hpp"
#include "reactmicp/systems/saturated_react/equilibrium_stagger.hpp"
#include "dfpmsolver/parabolic_driver.hpp"
#include "reactmicp/systems/saturated_react/transport_stagger.hpp"
#include <iostream>
#include "utils/log.hpp"
TEST_CASE("Equilibrium stagger", "[SaturatedReact, chemistry, equilibrium, stagger]") {
specmicp::logger::ErrFile::stream() = &std::cerr;
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::units::UnitsSet the_units;
the_units.mass = specmicp::units::MassUnit::kilogram;
the_units.length = specmicp::units::LengthUnit::centimeter;
specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(5, 0.1, 5.0);
specmicp::RawDatabasePtr the_database = leaching_database();
std::vector<specmicp::AdimensionalSystemSolution> list_initial_composition;
list_initial_composition.push_back(initial_leaching_pb(the_database, 0.005, the_units));
list_initial_composition.push_back(initial_blank_leaching_pb(the_database, 0.005, the_units));
std::vector<specmicp::index_t> list_fixed_nodes = {0, };
std::vector<int> index_initial_state = {1, 0, 0, 0, 0};
specmicp::reactmicp::systems::satdiff::SaturatedVariablesPtr variables =
specmicp::reactmicp::systems::satdiff::init_variables(the_mesh, the_database, the_units,
list_fixed_nodes,
list_initial_composition,
index_initial_state);
for (specmicp::index_t node: the_mesh->range_nodes())
{
variables->porosity(node) = 0.2;
variables->diffusion_coefficient(node) = 1.0e-8;
}
SECTION("Initialisation") {
specmicp::AdimensionalSystemConstraints constraints;
constraints.set_saturated_system();
constraints.charge_keeper = 1;
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 10.0;
options.solver_options.max_iter = 100;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = 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-6;
options.solver_options.steptol = 1e-14;
options.system_options.non_ideality_tolerance = 1e-10;
options.units_set = the_units;
variables->predictor() = variables->displacement();
specmicp::reactmicp::systems::satdiff::EquilibriumStagger stagger(constraints, options);
stagger.initialize_timestep(10, variables);
stagger.restart_timestep(variables);
REQUIRE(variables->displacement() == variables->predictor());
}
SECTION("Transport + equilibrium") {
specmicp::AdimensionalSystemConstraints constraints;
constraints.set_saturated_system();
constraints.charge_keeper = 1;
specmicp::AdimensionalSystemSolverOptions options;
options.solver_options.maxstep = 10.0;
options.solver_options.max_iter = 100;
options.solver_options.maxiter_maxstep = 100;
options.solver_options.use_crashing = false;
options.solver_options.use_scaling = 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-6;
options.solver_options.steptol = 1e-14;
options.system_options.non_ideality_tolerance = 1e-10;
options.units_set = the_units;
variables->predictor() = variables->displacement();
specmicp::reactmicp::systems::satdiff::EquilibriumStagger equilibrium_stagger(constraints, options);
specmicp::reactmicp::systems::satdiff::SaturatedTransportStagger transport_stagger(variables, list_fixed_nodes);
transport_stagger.initialize_timestep(10, variables);
equilibrium_stagger.initialize_timestep(10, variables);
int retcode;
retcode = (int) transport_stagger.restart_timestep(variables);
REQUIRE(retcode > 0);
retcode = (int) equilibrium_stagger.restart_timestep(variables);
REQUIRE(retcode > 0);
REQUIRE(variables->displacement() != variables->predictor());
}
}
c4science · Help