diff --git a/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp b/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp index 9323d44..87ef100 100644 --- a/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp +++ b/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp @@ -1,168 +1,183 @@ #include "catch.hpp" #include "reactmicp/systems/unsaturated/transport_constraints.hpp" #include "reactmicp/systems/unsaturated/saturation_equation.hpp" #include "reactmicp/systems/unsaturated/aqueous_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; static scalar_t one(index_t node, scalar_t saturation) {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("Aqueous equation") { +TEST_CASE("Aqueous equation", "[transport],[aqueous]") { // begin initialisation 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}); + mesh::Uniform1DMeshGeometry geom; + geom.dx = dx; + geom.nb_nodes = nb_nodes; + geom.section = cross_section; + + + const scalar_t D = 1e-3; + const scalar_t v = 1e-2; + + auto the_mesh = uniform_mesh1d(geom); auto raw_data = get_database(); VariablesInterface vars_interface(the_mesh, raw_data, {0}); Vector aq_conc(nb_nodes); aq_conc.setConstant(1.0); aq_conc(2) = 2.0; vars_interface.set_aqueous_concentration(2, aq_conc); + vars_interface.set_liquid_saturation(1.0); vars_interface.set_porosity(0.5); - vars_interface.set_liquid_diffusivity(1.0); + vars_interface.set_liquid_diffusivity(D); vars_interface.set_relative_liquid_diffusivity_model(one); - vars_interface.set_advection_flux(0.3); + vars_interface.set_advection_flux(v); TransportConstraints constraints; constraints.add_fixed_node(0); // end initialisation SECTION("Simple Residuals") { UnsaturatedVariablesPtr vars = vars_interface.get_variables(); LiquidAqueousComponentVariableBox aq_vars = vars->get_liquid_aqueous_component_variables(2); MainVariable& aq_concentration = aq_vars.aqueous_concentration; AqueousTransportEquation equation(the_mesh, aq_vars, constraints); REQUIRE(equation.get_neq() == 9); Vector residuals; equation.compute_residuals(aq_concentration.variable, aq_concentration.velocity, residuals); REQUIRE(residuals.rows() == equation.get_neq()); CHECK(residuals(4) == Approx(0.0).epsilon(1e-10)); - CHECK(residuals(0) == Approx(-(2.0-1.0)/0.1).epsilon(1e-10)); - CHECK(residuals(1) == Approx(+2*(2.0-1.0)/0.1-0.3*(1.0-2.0)).epsilon(1e-10)); - CHECK(residuals(2) == Approx(-(2.0-1.0)/0.1-0.3*(2.0-1.0)).epsilon(1e-10)); + CHECK(residuals(0) == Approx(-(2.0-1.0)/0.1*D).epsilon(1e-10)); + CHECK(residuals(1) == Approx(+2*(2.0-1.0)/0.1*D-v*(1.0-2.0)).epsilon(1e-10)); + CHECK(residuals(2) == Approx(-(2.0-1.0)/0.1*D-v*(2.0-1.0)).epsilon(1e-10)); Eigen::SparseMatrix jacobian; equation.compute_jacobian(aq_concentration.variable, aq_concentration.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.1-0.3)); - CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1)); + CHECK(jacobian.coeff(1, 0) == Approx(-1/0.1*D-v)); + CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1*D)); - vars_interface.set_advection_flux(-0.3); + vars_interface.set_advection_flux(-v); aq_concentration.transport_fluxes.setZero(); equation.compute_residuals(aq_concentration.variable, aq_concentration.velocity, residuals); REQUIRE(residuals.rows() == equation.get_neq()); CHECK(residuals(4) == Approx(0.0).epsilon(1e-10)); - CHECK(residuals(0) == Approx(-(2.0-1.0)/0.1-0.3*(2.0-1.0)).epsilon(1e-10)); - CHECK(residuals(1) == Approx(+2*(2.0-1.0)/0.1+0.3*(2.0-1.0)).epsilon(1e-10)); - CHECK(residuals(2) == Approx(-(2.0-1.0)/0.1).epsilon(1e-10)); + CHECK(residuals(0) == Approx(-(2.0-1.0)/0.1*D-v*(2.0-1.0)).epsilon(1e-10)); + CHECK(residuals(1) == Approx(+2*(2.0-1.0)/0.1*D+v*(2.0-1.0)).epsilon(1e-10)); + CHECK(residuals(2) == Approx(-(2.0-1.0)/0.1*D).epsilon(1e-10)); equation.compute_jacobian(aq_concentration.variable, aq_concentration.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.1)); - CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1-0.3)); + CHECK(jacobian.coeff(1, 0) == Approx(-1/0.1*D)); + CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1*D-v)); + std::cout << "jacobian : \n " << jacobian.toDense() << std::endl; } SECTION("Solving") { std::cout << "---------------- \n Aqueous component transport \n --------------\n"; UnsaturatedVariablesPtr vars = vars_interface.get_variables(); LiquidAqueousComponentVariableBox aq_vars = vars->get_liquid_aqueous_component_variables(2); MainVariable& aq_concentration = aq_vars.aqueous_concentration; AqueousTransportEquation equation(the_mesh, aq_vars, constraints); dfpmsolver::ParabolicDriver solver(equation); - solver.get_options().sparse_solver = sparse_solvers::SparseSolver::GMRES; - - dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(0.01, aq_concentration.variable); + solver.get_options().sparse_solver = sparse_solvers::SparseSolver::SparseLU; + std::cout << aq_concentration.variable << "\n ----- " << std::endl; + dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(1.0, aq_concentration.variable); REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); - std::cout << aq_concentration.variable << std::endl; + std::cout << aq_concentration.variable << "\n ----- " << std::endl; Vector save = aq_concentration.variable; aq_concentration.set_constant(1.0); - aq_concentration(6) = 2.0; + //aq_concentration(6) = 2.0; aq_concentration(7) = 2.0; aq_concentration.velocity.setZero(); aq_concentration.transport_fluxes.setZero(); - vars_interface.set_advection_flux(-3.0); + std::cout << aq_concentration.variable << "\n ----- " << std::endl; + + vars_interface.set_advection_flux(-v); constraints = TransportConstraints(); constraints.add_fixed_node(9); AqueousTransportEquation equation2(the_mesh, aq_vars, constraints); dfpmsolver::ParabolicDriver solver2(equation2); + solver2.get_options().sparse_solver = sparse_solvers::SparseSolver::SparseLU; - retcode = solver2.solve_timestep(0.01, aq_concentration.variable); + retcode = solver2.solve_timestep(1.0, aq_concentration.variable); REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); + std::cout << solver2.velocity() << "\n ----- " << std::endl; std::cout << aq_concentration.variable << std::endl; for (int i=0; i