Page MenuHomec4science
Files F121736719

saturation_equation.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Jul 13, 11:50

saturation_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/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 identity(index_t _, scalar_t saturation) {return saturation;}
static scalar_t decreasing_linear(index_t _, scalar_t saturation) {return 10*(1-saturation);}
static scalar_t decreasing_linear_2(index_t node, scalar_t saturation) {
if (node < 5)
return 10*(1-saturation);
else
return 30*(1-saturation);
}
static scalar_t one(index_t _, scalar_t __) {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("Saturation equation") {
index_t nb_nodes = 10;
scalar_t dx = 0.1;
scalar_t cross_section = 1.0;
auto the_mesh = uniform_mesh1d({dx, nb_nodes, cross_section});
auto raw_data = get_database();
VariablesInterface vars_interface(the_mesh, raw_data, {0, });
vars_interface.set_porosity(0.5);
vars_interface.set_water_aqueous_concentration(1.0);
vars_interface.set_liquid_permeability(1.0e-4);
vars_interface.set_liquid_diffusivity(1.0);
vars_interface.set_relative_liquid_diffusivity_model(one);
vars_interface.set_relative_liquid_permeability_model(one);
SECTION("Simple Residuals") {
// Initialisation
// --------------
vars_interface.set_liquid_saturation(1.0);
vars_interface.set_liquid_saturation(2.0, 0.5);
vars_interface.set_capillary_pressure_model(identity);
// Initialisation of equation
// --------------------------
UnsaturatedVariablesPtr vars = vars_interface.get_variables();
SaturationVariableBox sat_vars = vars->get_saturation_variables();
MainVariable& saturation = sat_vars.liquid_saturation;
auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
bcs->add_fixed_node(0);
SaturationEquation equation(0, the_mesh, sat_vars, bcs);
REQUIRE(equation.get_neq() == 9);
Vector residuals;
CHECK(sat_vars.capillary_pressure(1) == 1.0);
CHECK(sat_vars.capillary_pressure(2) == 0.5);
CHECK(sat_vars.capillary_pressure(3) == 1.0);
CHECK(sat_vars.liquid_permeability(1) == 1e-4);
CHECK(sat_vars.liquid_permeability(2) == 1e-4);
CHECK(sat_vars.liquid_permeability(3) == 1e-4);
CHECK(sat_vars.relative_liquid_diffusivity(1) == 1.0);
CHECK(sat_vars.relative_liquid_diffusivity(2) == 1.0);
CHECK(sat_vars.relative_liquid_permeability(1) == 1.0);
CHECK(sat_vars.relative_liquid_permeability(2) == 1.0);
equation.compute_residuals(saturation.variable, saturation.velocity, residuals);
REQUIRE(residuals.rows() == equation.get_neq());
CHECK(residuals(4) == Approx(0.0).epsilon(1e-10));
CHECK(residuals(0) == Approx(1.0/8.937*(0.5-1.0)/0.1).epsilon(1e-10));
CHECK(residuals(1) == Approx(-2.0*1.0/8.937*(0.5-1.0)/0.1).epsilon(1e-10));
CHECK(residuals(2) == Approx(1.0/8.937*(0.5-1.0)/0.1).epsilon(1e-10));
CHECK(residuals(0) == - saturation.transport_fluxes(1));
CHECK(saturation.transport_fluxes(0) == 0.0);
Eigen::SparseMatrix<scalar_t> jacobian;
equation.compute_jacobian(saturation.variable, saturation.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/8.937/0.1));
std::cout << jacobian.toDense() << std::endl;
}
SECTION("Solving") {
// Initialisation
// --------------
Vector saturation_init(nb_nodes);
saturation_init.setConstant(0.5);
saturation_init(2) = 0.8;
saturation_init(3) = 0.8;
saturation_init(4) = 0.8;
vars_interface.set_liquid_saturation(saturation_init);
vars_interface.set_capillary_pressure_model(decreasing_linear);
UnsaturatedVariablesPtr vars = vars_interface.get_variables();
SaturationVariableBox sat_vars = vars->get_saturation_variables();
MainVariable& saturation = sat_vars.liquid_saturation;
auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
bcs->add_fixed_node(0);
SaturationEquation equation(0, the_mesh, sat_vars, bcs);
dfpmsolver::ParabolicDriver<SaturationEquation> solver(equation);
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(0.01, saturation.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
for (index_t node: the_mesh->range_nodes())
{
CHECK(saturation.variable(node) >= 0.5);
CHECK(saturation.variable(node) < 0.8);
CHECK(saturation.velocity(node) == 0.0);
if (node >0)
{
CHECK(solver.velocity()(node) != 0.0);
}
}
std::cout << saturation.variable << "\n ---------" << std::endl;
std::cout << sat_vars.advection_flux.variable << "\n ---------" << std::endl;
retcode = solver.solve_timestep(0.01, saturation.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
for (index_t node: the_mesh->range_nodes())
{
CHECK(saturation.variable(node) >= 0.5);
CHECK(saturation.variable(node) < 0.8);
}
std::cout << saturation.variable << "\n ---------" << std::endl;
for (int i=0; i<500; ++i)
{
retcode = solver.solve_timestep(0.01, saturation.variable);
REQUIRE((int) retcode > (int) dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet);
}
for (index_t node: the_mesh->range_nodes())
{
CHECK(saturation.variable(node) == Approx(0.5).epsilon(1e-4));
}
std::cout << saturation.variable << std::endl;
}
SECTION("Solving - layered") {
std::cout << "-------- \n Layered system \n ---------" << std::endl;
Vector saturation_init(nb_nodes);
saturation_init.setConstant(0.5);
saturation_init(2) = 0.8;
saturation_init(3) = 0.8;
saturation_init(4) = 0.8;
vars_interface.set_liquid_saturation(saturation_init);
vars_interface.set_capillary_pressure_model(decreasing_linear_2);
UnsaturatedVariablesPtr vars = vars_interface.get_variables();
SaturationVariableBox sat_vars = vars->get_saturation_variables();
MainVariable& saturation = sat_vars.liquid_saturation;
auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
bcs->add_fixed_node(0);
SaturationEquation equation(0, the_mesh, sat_vars, bcs);
dfpmsolver::ParabolicDriver<SaturationEquation> solver(equation);
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(0.01, saturation.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
std::cout << saturation.variable << "\n ---------" << std::endl;
std::cout << sat_vars.advection_flux.variable << "\n ---------" << std::endl;
retcode = solver.solve_timestep(0.01, saturation.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
std::cout << saturation.variable << "\n ---------" << std::endl;
for (int i=0; i<500; ++i)
{
retcode = solver.solve_timestep(0.01, saturation.variable);
REQUIRE((int) retcode > (int) dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet);
}
REQUIRE(saturation.variable(3) == Approx(saturation.variable(4)).epsilon(1e-7));
REQUIRE(saturation.variable(4) != saturation.variable(5));
REQUIRE(saturation.variable(5) == Approx(saturation.variable(6)).epsilon(1e-7));
std::cout << saturation.variable << std::endl;
}
}

Event Timeline

c4science · Help