Page MenuHomec4science
Files F91980646

solver.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Nov 16, 08:38

solver.cpp
View Options

#include "catch.hpp"
#include "utils.hpp"
#include "dfpmsolver/parabolic_driver.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_program.hpp"
#include "reactmicp/systems/saturated_diffusion/variables.hpp"
#include "specmicp/reduced_system_solver.hpp"
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems;
using namespace specmicp::reactmicp::systems::siasaturated;
TEST_CASE("dfpm solver", "[dfpm, solver]") {
SECTION("transport program") {
std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
int nb_nodes = 5;
EquilibriumState initial_state = sample_carbo_composition(database);
std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(nb_nodes-1, 1e-3, 0.001);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(1e-8, 0.2);
SIASaturatedVariables variables(database, themesh, parameters);
SaturatedDiffusionProgram the_program(themesh, database, parameters);
Eigen::VectorXd& concentrations = variables.total_mobile_concentrations();
concentrations.setOnes();
concentrations.block(0, 0, the_program.get_ndf(), 1).setZero();
the_program.number_equations({0,});
dfpmsolver::ParabolicDriver<SaturatedDiffusionProgram> solver(the_program);
solver.get_velocity() = Eigen::VectorXd::Zero(concentrations.rows());
double dt = 100;
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(dt, concentrations);
REQUIRE(retcode == dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
for (int node=1; node<themesh->nnodes(); ++node)
{
for (int component: database->range_aqueous_component())
{
the_program.external_flux(node, component) =
parameters->porosity(0)*(1.0 - concentrations(the_program.get_dof_component(node, component)))
/ (parameters->density_water()*dt*themesh->cell_volume(node));
}
}
retcode = solver.restart_timestep(concentrations);
REQUIRE(retcode == dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
}
SECTION("micpsolver") {
std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
database::Database database_handler(database);
int nb_nodes = 5;
std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(nb_nodes-1, 1e-3, 0.001);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(1e-8, 0.2);
SIASaturatedVariables variables(database, themesh, parameters);
EquilibriumState initial_state = sample_carbo_composition(database);
variables.update_composition(1, initial_state);
Eigen::VectorXd totconc;
variables.nodal_update_total_amount(1, totconc);
totconc(database_handler.component_label_to_id("HCO3[-]")) += 0.05;
totconc(database_handler.component_label_to_id("HO[-]")) -= 0.05;
ReducedSystemSolverOptions options;
options.solver_options.max_factorization_step = 1;
ReducedSystemSolver solver(database, totconc, options);
Eigen::VectorXd spec_var(variables.speciation_variables().col(1));
specmicp::micpsolver::MiCPPerformance perf = solver.solve(spec_var);
REQUIRE((int) perf.return_code == (int) specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized);
EquilibriumState final_compo = solver.get_solution(spec_var);
variables.update_composition(1, final_compo);
}
}

Event Timeline

c4science · Help