diff --git a/tests/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp b/tests/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp index fddaee3..05972a1 100644 --- a/tests/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp +++ b/tests/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp @@ -1,304 +1,317 @@ #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 "database/database.hpp" #include "dfpm/meshes/generic_mesh1d.hpp" #include "dfpmsolver/parabolic_driver.hpp" #include using namespace specmicp; using namespace specmicp::mesh; using namespace specmicp::reactmicp::systems::unsaturated; extern template class dfpmsolver::ParabolicDriver; 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 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 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) { - retcode = solver.solve_timestep(0.01, saturation.variable); + + 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 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()) { - CHECK(saturation.variable(node) < 0.8); - - CHECK(vap_pressure.variable(node) < 240.0); + 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 << saturation.variable << std::endl; std::cout << sat_vars.partial_pressure.variable << "\n ------- " << std::endl; } }