Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88892991
saturation_equation.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Oct 21, 05:24
Size
8 KB
Mime Type
text/x-c
Expires
Wed, Oct 23, 05:24 (2 d)
Engine
blob
Format
Raw Data
Handle
21832676
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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
Log In to Comment