Page MenuHomec4science
Files F79798852

pressure_equation.cpp
No OneTemporary
Actions

File Metadata

Created
Wed, Aug 28, 07:25

pressure_equation.cpp
View Options

#include "catch.hpp"
#include "reactmicp/systems/unsaturated/pressure_equation.hpp"
#include "reactmicp/systems/unsaturated/variables.hpp"
#include "reactmicp/systems/unsaturated/variables_box.hpp"
#include "reactmicp/systems/unsaturated/variables_interface.hpp"
#include "reactmicp/systems/unsaturated/boundary_conditions.hpp"
#include "specmicp_database/database.hpp"
#include "dfpm/meshes/generic_mesh1d.hpp"
#include "dfpm/solver/parabolic_driver.hpp"
#include "specmicp_common/physics/constants.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 constexpr scalar_t rt = constants::gas_constant*(273.15+25);
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("Pressure equation", "[pressure]") {
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_liquid_saturation(0.0);
vars_interface.set_partial_pressure(0, 0.5);
vars_interface.set_partial_pressure(0, 2, 0.8);
vars_interface.set_porosity(0.5);
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);
auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component());
bcs->add_fixed_node(0);
SECTION("Simple Residuals") {
UnsaturatedVariablesPtr vars = vars_interface.get_variables();
PressureVariableBox pressure_vars = vars->get_vapor_pressure_variables();
MainVariable& pressure = pressure_vars.partial_pressure;
PressureEquation equation(0, the_mesh, pressure_vars, bcs);
REQUIRE(equation.get_neq() == 9);
Vector residuals;
equation.compute_residuals(pressure.variable, pressure.velocity, residuals);
std::cout << residuals << std::endl;
REQUIRE(residuals(0) == Approx(-(0.8-0.5)/0.1/rt).epsilon(1e-6));
REQUIRE(residuals(1) == Approx(2*(0.8-0.5)/0.1/rt).epsilon(1e-6));
Eigen::SparseMatrix<scalar_t> jacobian;
equation.compute_jacobian(pressure.variable, pressure.velocity, jacobian, 1.0);
auto mat = jacobian.toDense();
REQUIRE(mat(1, 0) == Approx(- 1/0.1/rt).epsilon(1e-6));
REQUIRE(mat(1, 1) == Approx(0.1/2.0*0.5*0.5/rt + 2*1/0.1/rt).epsilon(1e-4));
}
SECTION("Solving the system") {
std::cout << "========= \n Pressure Equation \n ======= \n";
UnsaturatedVariablesPtr vars = vars_interface.get_variables();
PressureVariableBox pressure_vars = vars->get_vapor_pressure_variables();
MainVariable& pressure = pressure_vars.partial_pressure;
PressureEquation equation(0, the_mesh, pressure_vars, bcs);
dfpmsolver::ParabolicDriver<PressureEquation> solver(equation);
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(0.01, pressure.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
CHECK(pressure.variable(0) == 0.5);
for (index_t node=1; node<10; ++node)
{
CHECK(pressure(node) >= 0.5);
CHECK(pressure(node) < 0.8);
}
std::cout << pressure.variable << std::endl;
for (int i=0; i<500; ++i)
{
dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(0.2, pressure.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
}
for (index_t node=0; node<10; ++node)
{
CHECK(pressure(node) == Approx(0.5).epsilon(1e-6));
}
std::cout << pressure.variable << std::endl;
}
SECTION("Solving with chemistry rate") {
std::cout << "========= \n Pressure Equation - chem rate \n ======= \n";
UnsaturatedVariablesPtr vars = vars_interface.get_variables();
PressureVariableBox pressure_vars = vars->get_vapor_pressure_variables();
MainVariable& pressure = pressure_vars.partial_pressure;
PressureEquation equation(0, the_mesh, pressure_vars, bcs);
dfpmsolver::ParabolicDriver<PressureEquation> solver(equation);
dfpmsolver::ParabolicDriverReturnCode retcode;
solver.initialize_timestep(0.01, pressure.variable);
scalar_t res = 10;
int cnt = 0;
while(res > 1e-4 and cnt < 100)
{
retcode = solver.restart_timestep(pressure.variable);
REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
pressure.chemistry_rate = pressure.transport_fluxes;
solver.velocity().setZero();
pressure.set_constant(0.5);
pressure.variable(2) = 0.8;
Vector residual;
solver.compute_residuals(pressure.variable, residual);
res = residual.norm();
++cnt;
}
CHECK(cnt < 100);
CHECK(pressure.variable(0) == 0.5);
for (index_t node=1; node<10; ++node)
{
CHECK(pressure(node) >= 0.5);
CHECK(pressure(node) <= 0.8);
CHECK(pressure.velocity(node) == Approx(0.0));
}
std::cout << pressure.variable << std::endl;
}
}

Event Timeline

c4science · Help