Homec4science
Diffusion SpecMiCP / ReactMiCP (single) Browse

specmicp/tests/reactmicp/systems/chloride/transport_program.cpp
c90630f67f84single

Last Change
Raw File

transport_program.cpp
View Options

#include "catch.hpp"
#include "reactmicp/systems/chloride/variables.hpp"
#include "reactmicp/systems/chloride/variables_interface.hpp"
#include "reactmicp/systems/chloride/transport_program.hpp"
#include "specmicp_common/log.hpp"
#include "specmicp_database/io/configuration.hpp"
#include "specmicp/problem_solver/smart_solver.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp/adimensional/adimensional_system_solver_structs.hpp"
#include "reactmicp/systems/chloride/transport_program_struct.hpp"
#include "reactmicp/systems/chloride/equilibrium_stagger.hpp"
#include "reactmicp/systems/chloride/boundary_conditions.hpp"
#include "reactmicp/systems/common/equilibrium_options_vector.hpp"
#include "specmicp/adimensional/adimensional_system_structs.hpp"
#include <vector>
#include <iostream>
using namespace specmicp;
using namespace specmicp::reactmicp::systems::chloride;
static mesh::Mesh1DPtr get_mesh(index_t nb_nodes)
{
mesh::Uniform1DMeshGeometry amesh;
amesh.dx = 1.0;
amesh.section = 1.0;
amesh.nb_nodes = nb_nodes;
return mesh::uniform_mesh1d(amesh);
}
static specmicp::database::RawDatabasePtr get_database()
{
static database::RawDatabasePtr raw_data {nullptr};
if (raw_data == nullptr)
{
specmicp::database::Database thedatabase(TEST_CEMDATA_PATH);
thedatabase.keep_only_components({"H2O","Na[+]", "Cl[-]"});
thedatabase.remove_half_cell_reactions();
raw_data = thedatabase.get_database();
raw_data->freeze_db();
}
return raw_data;
}
static specmicp::AdimensionalSystemSolution low_conc(database::RawDatabasePtr data, units::UnitsSet tunits)
{
ReactantBox box(data, tunits);
box.set_solution(0.3, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.7);
box.add_aqueous_species("NaCl", 0.001, "mol/kg");
box.set_charge_keeper("Cl[-]");
SmartAdimSolver solver(data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve init cond");
}
return solver.get_solution();
}
static specmicp::AdimensionalSystemSolution high_conc(database::RawDatabasePtr data, units::UnitsSet tunits)
{
ReactantBox box(data, tunits);
box.set_solution(0.3, "L/L");
box.set_saturated_system();
box.set_inert_volume_fraction(0.7);
box.add_aqueous_species("NaCl", 0.5, "mol/kg");
SmartAdimSolver solver(data, box, false);
if (! solver.solve()) {
throw std::runtime_error("Cannot solve bound cond");
}
return solver.get_solution();
}
TEST_CASE("Transport program", "[Transport,dfpm,program]") {
SECTION ("Init") {
const index_t nb_nodes {50};
mesh::Mesh1DPtr amesh = get_mesh(nb_nodes);
auto tunits = units::SI_units;
auto data = get_database();
ChlorideVariablesInterface interface(data, amesh, tunits);
std::vector<AdimensionalSystemSolution> sols = {low_conc(data, tunits), high_conc(data, tunits)};
std::vector<index_t> indices(nb_nodes, 0);
indices[0] = 1;
interface.set_adim_solutions(indices, sols);
index_t i_na = data->get_id_component("Na[+]");
index_t i_cl = data->get_id_component("Cl[-]");
index_t i_ho = data->get_id_component("HO[-]");
interface.set_diffusion_coefficient_component(i_na, 0.3);
interface.set_default_diffusion_coefficient_component(0.5);
interface.set_diffusion_coefficient_component("Cl[-]", 0.2);
interface.set_default_diffusion_resistance_factor(1.0);
auto variables = interface.get_variables();
std::shared_ptr<BoundaryConditions> bcs = BoundaryConditions::make(amesh, data);
bcs->fix_node(0);
AdimensionalSystemConstraints& constraints = bcs->get_constraint("default");
constraints.set_saturated_system();
std::shared_ptr<reactmicp::systems::EquilibriumOptionsVector> opts =
reactmicp::systems::EquilibriumOptionsVector::make(nb_nodes);
auto& opt = opts->get("default");
opt.units_set = tunits;
opt.system_options.cutoff_total_concentration = 1e-14;
opt.solver_options.set_tolerance(1e-12,1e-14);
opt.solver_options.set_maximum_step_length(1.0, 200);
std::shared_ptr<ChlorideChemistryStaggerBase> chem_stagger =
std::make_shared<ChlorideEquilibriumStagger>(variables, bcs, opts);
ChlorideProgram program(variables, bcs->get_fixed_dofs());
program.register_chem_stagger(chem_stagger);
REQUIRE(program.get_ndf() == 2+1);
REQUIRE(program.get_tot_ndf() == (2+1)*nb_nodes);
REQUIRE(program.get_neq() == (2+1)*(nb_nodes-1));
Vector displacement = variables->main_variables();
Vector velocity = variables->velocities();
Vector residuals = Vector::Zero(program.get_neq());
program.compute_residuals(displacement, velocity, residuals);
//std::cout << residuals << std::endl;
Vector eresiduals = Vector::Zero(2*3);
program.compute_element_residuals(0, displacement, velocity, variables->adim_solution(0), variables->adim_solution(1), eresiduals, false, true);
}
SECTION ("Jacobian") {
const index_t nb_nodes {10};
mesh::Mesh1DPtr amesh = get_mesh(nb_nodes);
auto tunits = units::SI_units;
auto data = get_database();
ChlorideVariablesInterface interface(data, amesh, tunits);
std::vector<AdimensionalSystemSolution> sols = {low_conc(data, tunits), high_conc(data, tunits)};
std::vector<index_t> indices(nb_nodes, 0);
indices[0] = 1;
interface.set_adim_solutions(indices, sols);
index_t i_na = data->get_id_component("Na[+]");
index_t i_cl = data->get_id_component("Cl[-]");
index_t i_ho = data->get_id_component("HO[-]");
interface.set_diffusion_coefficient_component(i_na, 0.3);
interface.set_default_diffusion_coefficient_component(0.5);
interface.set_diffusion_coefficient_component("Cl[-]", 0.2);
interface.set_default_diffusion_resistance_factor(1.0);
auto variables = interface.get_variables();
std::shared_ptr<BoundaryConditions> bcs = BoundaryConditions::make(amesh, data);
bcs->fix_node(0);
AdimensionalSystemConstraints& constraints = bcs->get_constraint("default");
constraints.set_saturated_system();
std::shared_ptr<reactmicp::systems::EquilibriumOptionsVector> opts =
reactmicp::systems::EquilibriumOptionsVector::make(nb_nodes);
auto& opt = opts->get("default");
opt.units_set = tunits;
opt.system_options.cutoff_total_concentration = 1e-14;
opt.solver_options.set_tolerance(1e-12,1e-14);
opt.solver_options.set_maximum_step_length(1.0, 200);
std::shared_ptr<ChlorideChemistryStaggerBase> chem_stagger =
std::make_shared<ChlorideEquilibriumStagger>(variables, bcs, opts);
ChlorideProgram program(variables, bcs->get_fixed_dofs());
program.register_chem_stagger(chem_stagger);
REQUIRE(program.get_ndf() == 2+1);
REQUIRE(program.get_tot_ndf() == (2+1)*nb_nodes);
REQUIRE(program.get_neq() == (2+1)*(nb_nodes-1));
Vector displacement = variables->main_variables();
Vector velocity = variables->velocities();
Vector residuals = Vector::Zero(program.get_neq());
program.compute_residuals(displacement, velocity, residuals);
Vector eresiduals = Vector::Zero(2*3);
program.compute_element_residuals(0, displacement, velocity, variables->adim_solution(0), variables->adim_solution(1), eresiduals, false, true);
program.set_approx_method(AdimSolutionPerturbationMethod::component);
Eigen::SparseMatrix<scalar_t> jacob(program.get_neq(), program.get_neq());
program.compute_jacobian(displacement, velocity, jacob, 1.0);
std::cout << jacob.toDense() << std::endl;
}
}
c4science · Help