Homec4science
Diffusion SpecMiCP / ReactMiCP (dev) Browse

specmicp/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp
f8005ed71a51dev

Last Change
Raw File

aqueous_equation.cpp
View Options

#include "catch.hpp"
#include "reactmicp/systems/unsaturated/boundary_conditions.hpp"
#include "reactmicp/systems/unsaturated/saturation_equation.hpp"
#include "reactmicp/systems/unsaturated/aqueous_equation.hpp"
#include "reactmicp/systems/unsaturated/variables.hpp"
#include "reactmicp/systems/unsaturated/variables_box.hpp"
#include "reactmicp/systems/unsaturated/variables_interface.hpp"
#include "specmicp_database/database.hpp"
#include "dfpm/meshes/generic_mesh1d.hpp"
#include "dfpm/solver/parabolic_driver.hpp"
#include <iostream>
using namespace specmicp;
using namespace specmicp::mesh;
using namespace specmicp::reactmicp::systems::unsaturated;
static scalar_t one(index_t node, scalar_t saturation) {return 1.0;}
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", "H[+]", "Ca[2+]", "Si(OH)4"});
raw_data = thedatabase.get_database();
raw_data->freeze_db();
}
return raw_data;
}
TEST_CASE("Aqueous equation", "[transport],[aqueous]") {
// begin initialisation
index_t nb_nodes = 10;
scalar_t dx = 0.1;
scalar_t cross_section = 1.0;
mesh::Uniform1DMeshGeometry geom;
geom.dx = dx;
geom.nb_nodes = nb_nodes;
geom.section = cross_section;
const scalar_t D = 1e-3;
const scalar_t v = 1e-2;
auto the_mesh = uniform_mesh1d(geom);
auto raw_data = get_database();
VariablesInterface vars_interface(the_mesh, raw_data, {0});
Vector aq_conc(nb_nodes);
aq_conc.setConstant(1.0);
aq_conc(2) = 2.0;
vars_interface.set_aqueous_concentration(2, aq_conc);
vars_interface.set_liquid_saturation(1.0);
vars_interface.set_porosity(0.5);
vars_interface.set_liquid_diffusivity(D);
vars_interface.set_relative_liquid_diffusivity_model(one);
vars_interface.set_advection_flux(v);
auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
bcs->add_fixed_node(0);
// end initialisation
SECTION("Simple Residuals") {
UnsaturatedVariablesPtr vars = vars_interface.get_variables();
LiquidAqueousComponentVariableBox aq_vars = vars->get_liquid_aqueous_component_variables(2);
MainVariable& aq_concentration = aq_vars.aqueous_concentration;
AqueousTransportEquation equation(2, the_mesh, aq_vars, bcs);
REQUIRE(equation.get_neq() == 9);
Vector residuals;
equation.compute_residuals(aq_concentration.variable, aq_concentration.velocity, residuals);
REQUIRE(residuals.rows() == equation.get_neq());
CHECK(residuals(4) == Approx(0.0).epsilon(1e-10));
CHECK(residuals(0) == Approx(-(2.0-1.0)/0.1*D).epsilon(1e-10));
CHECK(residuals(1) == Approx(+2*(2.0-1.0)/0.1*D-v*(1.0-2.0)).epsilon(1e-10));
CHECK(residuals(2) == Approx(-(2.0-1.0)/0.1*D-v*(2.0-1.0)).epsilon(1e-10));
Eigen::SparseMatrix<scalar_t> jacobian;
equation.compute_jacobian(aq_concentration.variable, aq_concentration.velocity, jacobian, 1.0);
CHECK(jacobian.rows() == equation.get_neq());
CHECK(jacobian.cols() == equation.get_neq());
CHECK(jacobian.coeff(1, 8) == Approx(0.0).epsilon(1e-10));
CHECK(jacobian.coeff(1, 0) == Approx(-1/0.1*D-v));
CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1*D));
vars_interface.set_advection_flux(-v);
aq_concentration.transport_fluxes.setZero();
equation.compute_residuals(aq_concentration.variable, aq_concentration.velocity, residuals);
REQUIRE(residuals.rows() == equation.get_neq());
CHECK(residuals(4) == Approx(0.0).epsilon(1e-10));
CHECK(residuals(0) == Approx(-(2.0-1.0)/0.1*D-v*(2.0-1.0)).epsilon(1e-10));
CHECK(residuals(1) == Approx(+2*(2.0-1.0)/0.1*D+v*(2.0-1.0)).epsilon(1e-10));
CHECK(residuals(2) == Approx(-(2.0-1.0)/0.1*D).epsilon(1e-10));
equation.compute_jacobian(aq_concentration.variable, aq_concentration.velocity, jacobian, 1.0);
CHECK(jacobian.rows() == equation.get_neq());
CHECK(jacobian.cols() == equation.get_neq());
CHECK(jacobian.coeff(1, 8) == Approx(0.0).epsilon(1e-10));
CHECK(jacobian.coeff(1, 0) == Approx(-1/0.1*D));
CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1*D-v));
std::cout << "jacobian : \n " << jacobian.toDense() << std::endl;
}
SECTION("Solving") {
std::cout << "---------------- \n Aqueous component transport \n --------------\n";
UnsaturatedVariablesPtr vars = vars_interface.get_variables();
LiquidAqueousComponentVariableBox aq_vars = vars->get_liquid_aqueous_component_variables(2);
MainVariable& aq_concentration = aq_vars.aqueous_concentration;
AqueousTransportEquation equation(2, the_mesh, aq_vars, bcs);
dfpmsolver::ParabolicDriver<AqueousTransportEquation> solver(equation);
solver.get_options().sparse_solver = sparse_solvers::SparseSolver::SparseLU;
std::cout << aq_concentration.variable << "\n ----- " << std::endl;
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(1.0, aq_concentration.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
std::cout << aq_concentration.variable << "\n ----- " << std::endl;
Vector save = aq_concentration.variable;
aq_concentration.set_constant(1.0);
//aq_concentration(6) = 2.0;
aq_concentration(7) = 2.0;
aq_concentration.velocity.setZero();
aq_concentration.transport_fluxes.setZero();
std::cout << aq_concentration.variable << "\n ----- " << std::endl;
vars_interface.set_advection_flux(-v);
auto bcs2 = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
bcs2->add_fixed_node(9);
AqueousTransportEquation equation2(2, the_mesh, aq_vars, bcs2);
dfpmsolver::ParabolicDriver<AqueousTransportEquation> solver2(equation2);
solver2.get_options().sparse_solver = sparse_solvers::SparseSolver::SparseLU;
retcode = solver2.solve_timestep(1.0, aq_concentration.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
std::cout << solver2.velocity() << "\n ----- " << std::endl;
std::cout << aq_concentration.variable << std::endl;
for (int i=0; i<nb_nodes; ++i)
{
CHECK(aq_concentration(i) == Approx(save(nb_nodes-1-i)).epsilon(1e-6));
}
}
}
c4science · Help