Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79796580
saturation_pressure_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
Wed, Aug 28, 07:08
Size
11 KB
Mime Type
text/x-c
Expires
Fri, Aug 30, 07:08 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20234412
Attached To
rSPECMICP SpecMiCP / ReactMiCP
saturation_pressure_equation.cpp
View Options
#include "catch.hpp"
#include "reactmicp/systems/unsaturated/boundary_conditions.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;
auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
bcs->add_fixed_node(0);
SaturationPressureEquation 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") {
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;
auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
bcs->add_fixed_node(0);
SaturationPressureEquation equation(0, the_mesh, sat_vars, bcs);
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;
auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
bcs->add_fixed_node(0);
bcs->add_gas_node(0);
SaturationPressureEquation equation(0, the_mesh, sat_vars, bcs);
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;
}
}
Event Timeline
Log In to Comment