Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/tests/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp
febe8f3e40e1master

Last Change
Raw File

saturation_pressure_equation.cpp
View Options

#include "catch.hpp"
#include "reactmicp/systems/unsaturated/transport_constraints.hpp"
#include "reactmicp/systems/unsaturated/saturation_pressure_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;
extern template class dfpmsolver::ParabolicDriver<SaturationPressureEquation>;
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 vapor_pressure(index_t _, scalar_t saturation) {return 300*(saturation);}
static scalar_t one(index_t _, scalar_t __) {return 1.0;}
static scalar_t zero(index_t _, scalar_t __) {return 0.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 pressure 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);
vars_interface.set_relative_gas_diffusivity_model(one);
vars_interface.set_resistance_gas_diffusivity(1.0);
vars_interface.set_binary_gas_diffusivity(0, 1.0);
SECTION("Simple Residuals") {
vars_interface.set_vapor_pressure_model(zero);
// 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();
SaturationPressureVariableBox sat_vars = vars->get_saturation_pressure_variables();
MainVariable& saturation = sat_vars.liquid_saturation;
TransportConstraints constraints;
constraints.add_fixed_node(0);
SaturationPressureEquation equation(the_mesh, sat_vars, constraints);
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") {
std::cout << " ---------- \n "
<< "Pressure saturation equation "
<< "\n ----------- "
<< std::endl;
vars_interface.set_vapor_pressure_model(vapor_pressure);
// 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();
SaturationPressureVariableBox sat_vars = vars->get_saturation_pressure_variables();
MainVariable& saturation = sat_vars.liquid_saturation;
MainVariable& vap_pressure = sat_vars.partial_pressure;
TransportConstraints constraints;
constraints.add_fixed_node(0);
SaturationPressureEquation equation(the_mesh, sat_vars, constraints);
dfpmsolver::ParabolicDriver<SaturationPressureEquation> 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())
{
vars->set_vapor_pressure(node);
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);
}
CHECK(vap_pressure.variable(node) >= 150.0);
CHECK(vap_pressure.variable(node) < 240.0);
}
std::cout << saturation.variable << "\n ---------" << std::endl;
std::cout << sat_vars.advection_flux.variable << "\n ---------" << std::endl;
std::cout << sat_vars.partial_pressure.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())
{
vars->set_vapor_pressure(node);
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)
{
for (index_t node: the_mesh->range_nodes())
{
vars->set_vapor_pressure(node);
}
retcode = solver.solve_timestep(0.05, saturation.variable);
REQUIRE((int) retcode > (int) dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet);
}
for (index_t node: the_mesh->range_nodes())
{
vars->set_vapor_pressure(node);
CHECK(saturation.variable(node) == Approx(0.5).epsilon(1e-4));
CHECK(sat_vars.partial_pressure(node) == Approx(150).epsilon(1e-4));
}
std::cout << saturation.variable << std::endl;
std::cout << sat_vars.partial_pressure.variable << "\n ------- " << std::endl;
}
SECTION("Solving") {
std::cout << " ---------- \n "
<< "Pressure saturation equation - pressure driven"
<< "\n ----------- "
<< std::endl;
vars_interface.set_vapor_pressure_model(vapor_pressure);
vars_interface.set_liquid_permeability(1.0);
// Initialisation
// --------------
Vector saturation_init(nb_nodes);
saturation_init.setConstant(0.8);
saturation_init(0) = 0.0;
vars_interface.set_liquid_saturation(saturation_init);
Vector vappressure_init(nb_nodes);
vappressure_init.setConstant(240);
vappressure_init(0) = 10;
vars_interface.set_partial_pressure(0, vappressure_init);
vars_interface.set_capillary_pressure_model(decreasing_linear);
UnsaturatedVariablesPtr vars = vars_interface.get_variables();
SaturationPressureVariableBox sat_vars = vars->get_saturation_pressure_variables();
MainVariable& saturation = sat_vars.liquid_saturation;
MainVariable& vap_pressure = sat_vars.partial_pressure;
TransportConstraints constraints;
constraints.add_fixed_node(0);
constraints.add_gas_node(0);
SaturationPressureEquation equation(the_mesh, sat_vars, constraints);
dfpmsolver::ParabolicDriver<SaturationPressureEquation> solver(equation);
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(0.1, saturation.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
for (index_t node: the_mesh->range_nodes())
{
if (node != 0) {
vars->set_vapor_pressure(node);
CHECK(saturation.variable(node) < 0.8);
CHECK(vap_pressure.variable(node) < 240.0);
}
}
std::cout << saturation.variable << "\n ---------" << std::endl;
std::cout << sat_vars.advection_flux.variable << "\n ---------" << std::endl;
std::cout << sat_vars.partial_pressure.variable << "\n ------- " << std::endl;
retcode = solver.solve_timestep(0.1, saturation.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
for (index_t node: the_mesh->range_nodes())
{
if (node != 0) {vars->set_vapor_pressure(node);}
CHECK(saturation.variable(node) < 0.8);
}
std::cout << saturation.variable << "\n ---------" << std::endl;
//for (int i=0; i<500; ++i)
//{
// retcode = solver.solve_timestep(1.0, 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));
// CHECK(sat_vars.partial_pressure(node) == Approx(150).epsilon(1e-4));
//}
std::cout << sat_vars.partial_pressure.variable << "\n ------- " << std::endl;
}
}
c4science · Help