diff --git a/examples/reactmicp/unsaturated/acc_carbo.cpp b/examples/reactmicp/unsaturated/acc_carbo.cpp index f17e375..24500a3 100644 --- a/examples/reactmicp/unsaturated/acc_carbo.cpp +++ b/examples/reactmicp/unsaturated/acc_carbo.cpp @@ -1,590 +1,590 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ //! \file acc_carbo.cpp //! \brief Accelerated carbonation example //! //! This is an example to show the features of the unsaturated //! reactive transport solver //! It is a simple study of the accelerated carbonation of cement paste //! (with drying) //! //! #include "reactmicp/reactmicp_unsaturated.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" #include #include "reactmicp/io/hdf5_unsaturated.hpp" #include "specmicp_common/io/specmicp_hdf5.hpp" #include "dfpm/io/hdf5_mesh.hpp" #include "specmicp_database/io/hdf5_database.hpp" #include "dfpm/solver/parabolic_structs.hpp" #include "specmicp_common/physics/unsaturated_laws.hpp" #include "specmicp_common/cli/parser.hpp" #include #include using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::unsaturated; using VariablesBase = solver::VariablesBase; using StaggerReturnCode = solver::StaggerReturnCode; class UpscalingStagger; database::RawDatabasePtr get_database(); AdimensionalSystemSolution get_initial_condition( database::RawDatabasePtr the_database, const units::UnitsSet& units_set ); void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species); // Upscaling stagger -> user models class UpscalingStagger: public solver::UpscalingStaggerBase { public: UpscalingStagger() {} ~UpscalingStagger() {} //! \brief Initialize the stagger at the beginning of the computation //! //! \param var a shared_ptr to the variables virtual void initialize(VariablesBase * const var) override; //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables virtual void initialize_timestep( scalar_t dt, VariablesBase * const var) override; //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables virtual StaggerReturnCode restart_timestep( VariablesBase * const var) override; // User models scalar_t capillary_pressure(index_t node, scalar_t saturation); scalar_t vapor_pressure(index_t node, scalar_t saturation); scalar_t relative_liquid_diffusivity(index_t node, scalar_t saturation); scalar_t relative_liquid_permeability(index_t node, scalar_t saturation); scalar_t relative_gas_diffusivity(index_t node, scalar_t saturation); scalar_t intrinsic_permeability(scalar_t porosity); scalar_t intrinsic_liquid_diffusivity(scalar_t porosity); scalar_t resistance_gas_diffusivity(scalar_t porosity); private: scalar_t m_a {37.5479e6}; scalar_t m_b {2.1684}; scalar_t m_alpha {7.27}; scalar_t m_beta {0.440}; scalar_t m_pv_sat {3170}; scalar_t m_exp_r_l_d {3.00}; scalar_t m_exp_i_g_d {2.74}; scalar_t m_exp_r_g_d {4.20}; scalar_t m_exp_r_l_k {0.4191}; scalar_t m_k_0 {1e-21*1e4}; scalar_t m_d_0 {2.3e-13*1e4}; scalar_t m_phi_0 {0.2}; }; // implementations int main(int argc, char* argv[]) { cli::CommandLineParser cli_parser; cli_parser.add_option('n', "min_dt", 0.01, "minimum timestep (s)"); cli_parser.add_option('m', "max_dt", 10.0, "maximum timestep (s)"); cli_parser.add_option('d', "dx", 0.05, "space step (cm)"); cli_parser.add_option('r', "run_until", 10.0*24.0*3600.0, "run simulation until s"); cli_parser.add_option('p', "p_h20", 1000.0, "vapor pressure at the boundary"); cli_parser.add_option('c', "p_co2", 500.0, "vapor pressure at the boundary"); cli_parser.set_help_message("Atmospheric carbonation"); cli_parser.parse(argc, argv); init_logger(&std::cout, logger::Warning); index_t nb_nodes = 25; scalar_t dx = cli_parser.get_option("dx"); // cm scalar_t cross_section = 5; // cm^2 units::UnitsSet units_set; units_set.length = units::LengthUnit::centimeter; specmicp::RawDatabasePtr the_database = get_database(); AdimensionalSystemSolution init = get_initial_condition( the_database, units_set); AdimensionalSystemSolutionExtractor extr(init, the_database, units_set); std::cout << extr.porosity() << " - " << extr.saturation_water() << std::endl; //return EXIT_SUCCESS; mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d({dx, nb_nodes, cross_section}); index_t id_co2 = the_database->get_id_component("CO2"); VariablesInterface variable_interface(the_mesh, the_database, {0, id_co2}, units_set); variable_interface.set_binary_gas_diffusivity(0, 0.240); variable_interface.set_binary_gas_diffusivity(id_co2, 0.160); std::shared_ptr upscaling_stagger = std::make_shared(); auto* call_ptr = upscaling_stagger.get(); variable_interface.set_vapor_pressure_model( [call_ptr](index_t node, scalar_t sat) -> scalar_t { return call_ptr->vapor_pressure(node, sat); }); for (index_t node=1; nodenb_nodes(); ++node) { variable_interface.initialize_variables(node, extr); } variable_interface.set_porosity(0, 1.0); variable_interface.set_liquid_saturation(0, 0.0); variable_interface.set_partial_pressure(id_co2, 0, cli_parser.get_option("p_co2")); variable_interface.set_partial_pressure(0, 0, cli_parser.get_option("p_h20")); std::shared_ptr vars = variable_interface.get_variables(); vars->set_aqueous_scaling(0, 1.0e3/18.3e-3*1e-6); for (auto aqc: the_database->range_aqueous_component()) { vars->set_aqueous_scaling(aqc, 100e-6); } //vars->set_gaseous_scaling(0, 1.0/vars->get_rt()); //vars->set_gaseous_scaling(id_co2, 1.0/vars->get_rt()); std::cout << "Youpla : " << vars->get_rt() << std::endl; - auto bcs = BoundaryConditions::make(the_mesh->nb_nodes()); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), the_database->nb_component()); bcs->add_gas_node(0); bcs->fork_constraint(0, "gas_node"); auto chem_options = EquilibriumOptionsVector::make(the_mesh->nb_nodes()); auto& def_chem_opts = chem_options->get("default"); def_chem_opts.units_set = units_set; micpsolver::MiCPSolverOptions& copts = def_chem_opts.solver_options; copts.set_maximum_iterations(200); copts.set_maximum_step_length(10, 200); copts.set_tolerance(1e-8, 1e-12); def_chem_opts.system_options.non_ideality_tolerance = 1e-14; std::shared_ptr chemistry_stagger = EquilibriumStagger::make(vars, bcs, chem_options); upscaling_stagger->initialize(vars.get()); vars->set_relative_variables(0); chemistry_stagger->initialize(vars.get()); std::shared_ptr transport_stagger = UnsaturatedTransportStagger::make(vars, bcs, true, true); transport_stagger->initialize(vars.get()); bcs->get_constraint("gas_node").set_water_partial_pressure_model( std::bind(&UpscalingStagger::vapor_pressure, upscaling_stagger.get(), 1, std::placeholders::_1) ); auto check = variable_interface.check_variables(); if (check >= VariablesValidity::error) { std::cout << variable_interface.get_liquid_saturation() << std::endl; throw std::runtime_error("Invalid variables"); } solver::ReactiveTransportSolver react_solver(transport_stagger, chemistry_stagger, upscaling_stagger); solver::ReactiveTransportOptions& opts = react_solver.get_options(); opts.residuals_tolerance = 1e-2; opts.good_enough_tolerance = 0.99; opts.maximum_iterations = 100; transport_stagger->get_saturation_options()->maximum_iterations = 500; transport_stagger->get_saturation_options()->step_tolerance = 1e-10; transport_stagger->get_saturation_options()->residuals_tolerance = 1e-8; for (auto ind: the_database->range_aqueous_component()) { auto* opts = transport_stagger->get_aqueous_options(ind); opts->step_tolerance = 1e-10; opts->residuals_tolerance = 1e-8; } //transport_stagger->get_gas_options(3)->maximum_iterations = 5000; opts.step_tolerance = 1e-10; io::UnsaturatedHDF5Saver saver("/tmp/out_acc_carbo.hdf5", vars, units_set); saver.save_timestep(0.0); solver::SimulationInformation simul_info("acc_carbo", 500); auto out_pol = [&saver]( scalar_t timestep, solver::VariablesBasePtr _) { saver.save_timestep(timestep); }; database::Database db_handler(the_database); db_handler.save("acc_carbo_data.yml"); //solver::ReactiveTransportReturnCode retcode = react_solver.solve_timestep(0.01, vars); //retcode = react_solver.solve_timestep(1.0, vars); //std::cout << (int) retcode << std::endl; // retcode = react_solver.solve_timestep(5.0, vars); // std::cout << (int) retcode << std::endl; // retcode = react_solver.solve_timestep(10.0, vars); // std::cout << (int) retcode << std::endl; // retcode = react_solver.solve_timestep(10.0, vars); // std::cout << (int) retcode << std::endl; solver::ReactiveTransportRunner runner( react_solver, cli_parser.get_option("min_dt"), cli_parser.get_option("max_dt"), simul_info); runner.set_output_policy(out_pol); runner.run_until(cli_parser.get_option("run_until"), vars); std::cout << vars->get_liquid_saturation().variable << std::endl; std::cout << vars->get_liquid_saturation()(1) - vars->get_liquid_saturation()(2) << std::endl; std::cout << vars->get_pressure_main_variables(id_co2).variable << std::endl; } void UpscalingStagger::initialize(VariablesBase * const var) { UnsaturatedVariables * const true_var = static_cast(var); // set relative models true_var->set_capillary_pressure_model( [this](index_t x, scalar_t sat){ return this->capillary_pressure(x, sat); }); true_var->set_relative_liquid_permeability_model( [this](index_t x, scalar_t sat){ return this->relative_liquid_permeability(x, sat); }); true_var->set_relative_liquid_diffusivity_model( [this](index_t x, scalar_t sat){ return this->relative_liquid_diffusivity(x, sat); }); true_var->set_relative_gas_diffusivity_model( [this](index_t x, scalar_t sat){ return this->relative_gas_diffusivity(x, sat); }); true_var->set_vapor_pressure_model( [this](index_t x, scalar_t sat){ return this->vapor_pressure(x, sat); }); // initialization SecondaryTransientVariable& porosity = true_var->get_porosity(); SecondaryVariable& permeability = true_var->get_liquid_permeability(); SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity(); SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity(); true_var->set_relative_variables(); gas_diffusivity(0) = resistance_gas_diffusivity(1.0); true_var->get_relative_gas_diffusivity()(0) = relative_gas_diffusivity(0, 1.0); for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { scalar_t phi = porosity(node); permeability(node) = intrinsic_permeability(phi); liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi); gas_diffusivity(node) = resistance_gas_diffusivity(phi); //true_var->set_relative_variables(node); } } scalar_t UpscalingStagger::capillary_pressure(index_t node, scalar_t saturation) { if (saturation == 0.0) return 0.0; return laws::capillary_pressure_BaroghelBouny(saturation, m_a, m_b); } scalar_t UpscalingStagger::vapor_pressure(index_t node, scalar_t saturation) { scalar_t pv; if (saturation == 0.0) pv = 0.0; else if (saturation >= 1.0) pv = m_pv_sat; else pv = m_pv_sat*laws::invert_kelvin_equation( capillary_pressure(node, saturation)); return pv; } scalar_t UpscalingStagger::relative_liquid_diffusivity( index_t node, scalar_t saturation ) { if (saturation == 0.0) return 0.0; return std::pow(saturation, m_exp_r_l_d); } scalar_t UpscalingStagger::relative_liquid_permeability( index_t node, scalar_t saturation ) { if (saturation == 0.0) return 0.0; return laws::relative_liquid_permeability_Mualem(saturation, 1.0/m_b); } scalar_t UpscalingStagger::relative_gas_diffusivity(index_t node, scalar_t saturation) { if (saturation == 0.0) return 1.0; return std::pow(1.0 - saturation, m_exp_r_g_d); } scalar_t UpscalingStagger::intrinsic_permeability(scalar_t porosity) { scalar_t tmp_1 = porosity / m_phi_0; scalar_t tmp_2 = (1.0 - m_phi_0) / (1.0 - porosity); tmp_1 = std::pow(tmp_1, 3); tmp_2 = std::pow(tmp_2, 2); return m_k_0 * tmp_1 * tmp_2; } scalar_t UpscalingStagger::intrinsic_liquid_diffusivity(scalar_t porosity) { return m_d_0 * std::exp(-9.95*porosity); } scalar_t UpscalingStagger::resistance_gas_diffusivity(scalar_t porosity) { return std::pow(porosity, m_exp_i_g_d); } void UpscalingStagger::initialize_timestep( scalar_t dt, VariablesBase * const var ) { } StaggerReturnCode UpscalingStagger::restart_timestep(VariablesBase * const var) { UnsaturatedVariables* true_var = static_cast(var); SecondaryTransientVariable& porosity = true_var->get_porosity(); SecondaryVariable& permeability = true_var->get_liquid_permeability(); SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity(); SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity(); for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { scalar_t phi = porosity(node); permeability(node) = intrinsic_permeability(phi); liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi); gas_diffusivity(node) = resistance_gas_diffusivity(phi); true_var->set_relative_variables(node); } return StaggerReturnCode::ResidualMinimized; } specmicp::RawDatabasePtr get_database() { specmicp::database::Database thedatabase(CEMDATA_PATH); std::map swap = {{"HCO3[-]", "CO2"},}; thedatabase.swap_components(swap); thedatabase.remove_half_cell_reactions(); database::RawDatabasePtr raw_data = thedatabase.get_database(); return raw_data; } AdimensionalSystemSolution get_initial_condition( database::RawDatabasePtr the_database, const units::UnitsSet& units_set ) { Vector oxyde_compo(4); oxyde_compo << 67.98, 22.66, 5.19, 2.96; Vector species_compo; compo_from_oxyde(oxyde_compo, species_compo); scalar_t mult = 1e-6; //species_compo *= 0.947e-2; // poro 0.47 //species_compo *= 1.08e-2; // poro 0.4 species_compo *= 1.25e-2; // poro 0.3 //species_compo *= 1.1e-2; Formulation formulation; formulation.mass_solution = 5e-4; formulation.amount_minerals = { {"C3S", species_compo(0)}, {"C2S", species_compo(1)}, {"C3A", species_compo(2)}, {"Gypsum", species_compo(3)}, {"Calcite", species_compo(0)*1e-3} }; //formulation.extra_components_to_keep = {"HCO3[-]"}; formulation.concentration_aqueous = { {"NaOH", mult*0.0298}, {"KOH", mult*0.0801}, {"Cl[-]", mult*0.0001} }; Dissolver dissolver(the_database); UpscalingStagger upsstag_init; // just for vapor pressure model AdimensionalSystemConstraints constraints; constraints.set_total_concentrations( dissolver.dissolve(formulation, true) ); //constraints.set_inert_volume_fraction(0.1); AdimensionalSystemSolver the_solver(the_database, constraints); the_solver.get_options().units_set = units_set; the_solver.get_options().system_options.non_ideality_tolerance = 1e-12; micpsolver::MiCPSolverOptions* solver_options = &the_solver.get_options().solver_options; //solver_options.maxstep = 10.0; solver_options->maxiter_maxstep = 200; solver_options->max_iter = 200; solver_options->set_tolerance(1e-9); Vector x; the_solver.initialise_variables(x, 0.8, -4); micpsolver::MiCPPerformance perf = the_solver.solve(x); if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) { throw std::runtime_error("Error : failed to solve initial conditions"); } auto prev_sol = the_solver.get_raw_solution(x); prev_sol.main_variables(0) *= 0.8; constraints.set_water_partial_pressure_model( std::bind(&UpscalingStagger::vapor_pressure, &upsstag_init, 1, std::placeholders::_1) ); the_solver = AdimensionalSystemSolver (the_database, constraints, prev_sol); the_solver.get_options().units_set = units_set; the_solver.get_options().system_options.non_ideality_tolerance = 1e-12; solver_options = &the_solver.get_options().solver_options; //solver_options.maxstep = 10.0; solver_options->maxiter_maxstep = 200; solver_options->set_tolerance(1e-9); perf = the_solver.solve(x); if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) { throw std::runtime_error("Error : failed to solve initial conditions"); } return the_solver.get_raw_solution(x); } void compo_from_oxyde(Vector& compo_oxyde, Vector& compo_species) { constexpr double M_CaO = 56.08; constexpr double M_SiO2 = 60.09; constexpr double M_Al2O3 = 101.96; constexpr double M_SO3 = 80.06; Eigen::MatrixXd compo_in_oxyde(4, 4); Vector molar_mass(4); molar_mass << 1.0/M_CaO, 1.0/M_SiO2, 1.0/M_Al2O3, 1.0/M_SO3; compo_oxyde = compo_oxyde.cwiseProduct(molar_mass); //C3S C2S C3A Gypsum compo_in_oxyde << 3, 2, 3, 1, // CaO 1, 1, 0, 0, // SiO2 0, 0, 1, 0, // Al2O3 0, 0, 0, 1; // SO3 Eigen::ColPivHouseholderQR solver(compo_in_oxyde); compo_species = solver.solve(compo_oxyde); } diff --git a/examples/reactmicp/unsaturated/drying.cpp b/examples/reactmicp/unsaturated/drying.cpp index fb70c23..97a1279 100644 --- a/examples/reactmicp/unsaturated/drying.cpp +++ b/examples/reactmicp/unsaturated/drying.cpp @@ -1,579 +1,579 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ // ====================================================== // // // // Drying example // // ============== // // // // ====================================================== // // // // This is a simple example to show how the unsaturated // system of ReactMiCP is able to solve a complex drying // equation. // // It implements a special chemistry stagger (DryingStagger) // which only computes the vapor pressure. // Macro to define the problem // =========================== #define NB_NODES 50 #define DX 0.05*1e-2 // cm #define CROSS_SECTION 5.0*1e-4 // cm^2 #define GAMMA_WATER_GLASS 0.1 // N/m #define GAMMA_WATER_AIR 71.97e-3 // N/m = kg/s^2 #define R_SMALL_BEAD 0.015*1e-2 // cm #define R_BIG_BEAD 0.04*1e-2 // cm #define M_v 18.015e-3 // kg/mol // Includes // ======== #include "specmicp_database/database.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" #include "reactmicp/systems/unsaturated/variables_sub.hpp" #include "reactmicp/systems/unsaturated/transport_stagger.hpp" #include "reactmicp/systems/unsaturated/transport_constraints.hpp" #include "reactmicp/solver/staggers_base/chemistry_stagger_base.hpp" #include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "reactmicp/solver/runner.hpp" #include "dfpm/meshes/generic_mesh1d.hpp" #include "specmicp_common/physics/laws.hpp" #include "specmicp_common/physics/units.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/io/csv_formatter.hpp" #include #include // Namespace // ========= using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems; using namespace specmicp::reactmicp::systems::unsaturated; // Declarations // ============ // Constants // ========== static const scalar_t rho_l_si = constants::water_density_25; static const scalar_t rho_l = 1e-6*constants::water_density_25; static const scalar_t cw_tilde_si = rho_l_si / M_v; static const scalar_t cw_tilde = 1e-6*cw_tilde_si; int main(); database::RawDatabasePtr get_database(); mesh::Mesh1DPtr get_mesh(); scalar_t leverett_function(scalar_t s_eff); scalar_t big_bead_cap_pressure(scalar_t saturation); scalar_t small_bead_cap_pressure(scalar_t saturation); scalar_t vapor_pressure(scalar_t pc); // Drying Stagger // ============== // // This stagger computes the exchange between water vapor and liquid water class DryingStagger: public solver::ChemistryStaggerBase { public: DryingStagger() {} virtual void initialize(solver::VariablesBase * const var) override {} //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables virtual void initialize_timestep( scalar_t dt, solver::VariablesBase * const var ) override { m_dt = dt; } //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables virtual solver::StaggerReturnCode restart_timestep( solver::VariablesBase * const var) override; void compute_one_node(index_t node, UnsaturatedVariables * const vars); private: scalar_t m_dt {-1.0}; }; // Upscaling Stagger // ================= // // Define how we compute capillary pressure, vapor pressure, ... class UpscalingDryingStagger: public solver::UpscalingStaggerBase { public: UpscalingDryingStagger(std::vector& is_big_bead): m_is_big_bead(is_big_bead) {} //! \brief Initialize the stagger at the beginning of the computation //! //! \param var a shared_ptr to the variables virtual void initialize(solver::VariablesBase * const var) override; //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables virtual void initialize_timestep( scalar_t dt, solver::VariablesBase * const var ) override {} //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables virtual solver::StaggerReturnCode restart_timestep( solver::VariablesBase * const var) override { return solver::StaggerReturnCode::ResidualMinimized; } scalar_t capillary_pressure_model(index_t node, scalar_t saturation); scalar_t vapor_pressure_model(index_t node, scalar_t saturation); scalar_t relative_liquid_permeability_model(index_t node, scalar_t saturation); scalar_t relative_liquid_diffusivity_model(index_t node, scalar_t saturation); scalar_t relative_gas_diffusivity_model(index_t node, scalar_t saturation); private: std::vector m_is_big_bead; }; // Implementation // ============== scalar_t UpscalingDryingStagger::capillary_pressure_model(index_t node, scalar_t saturation) { scalar_t pc = 0.0; if (m_is_big_bead[node]) { pc = big_bead_cap_pressure(saturation); } else pc = small_bead_cap_pressure(saturation); return pc; } scalar_t UpscalingDryingStagger::vapor_pressure_model(index_t node, scalar_t saturation) { return vapor_pressure(capillary_pressure_model(node, saturation)); } scalar_t UpscalingDryingStagger::relative_liquid_permeability_model(index_t node, scalar_t saturation) { return std::pow(saturation, 3); } scalar_t UpscalingDryingStagger::relative_liquid_diffusivity_model(index_t node, scalar_t saturation) { return saturation; } scalar_t UpscalingDryingStagger::relative_gas_diffusivity_model(index_t node, scalar_t saturation) { return (1.0 - saturation); } void UpscalingDryingStagger::initialize(solver::VariablesBase * const var) { UnsaturatedVariables * const true_vars = static_cast(var); true_vars->set_capillary_pressure_model( [this](index_t node, scalar_t sat) -> scalar_t { return this->capillary_pressure_model(node, sat);}); true_vars->set_vapor_pressure_model(std::bind(std::mem_fn( &UpscalingDryingStagger::vapor_pressure_model ), this, std::placeholders::_1, std::placeholders::_2)); true_vars->set_relative_liquid_permeability_model(std::bind(std::mem_fn( &UpscalingDryingStagger::relative_liquid_permeability_model ), this, std::placeholders::_1, std::placeholders::_2)); true_vars->set_relative_liquid_diffusivity_model(std::bind(std::mem_fn( &UpscalingDryingStagger::relative_liquid_diffusivity_model ), this, std::placeholders::_1, std::placeholders::_2)); true_vars->set_relative_gas_diffusivity_model(std::bind(std::mem_fn( &UpscalingDryingStagger::relative_gas_diffusivity_model ), this, std::placeholders::_1, std::placeholders::_2)); MainVariable& vapor_pressure = true_vars->get_pressure_main_variables(0); MainVariable& saturation = true_vars->get_liquid_saturation(); SecondaryTransientVariable& porosity = true_vars->get_porosity(); SecondaryVariable& liquid_diffusivity = true_vars->get_liquid_diffusivity(); SecondaryVariable& liquid_permeability = true_vars->get_liquid_permeability(); SecondaryVariable& resistance_gas_diffusivity = true_vars->get_resistance_gas_diffusivity(); true_vars->get_binary_gas_diffusivity(0) = 0.282*1e-4; porosity(0) = 1; resistance_gas_diffusivity(0) = 1.0; for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { vapor_pressure(node) = vapor_pressure_model(node, saturation(node)); if (m_is_big_bead[node]) { porosity(node) = 0.371; liquid_diffusivity(node) = 0; liquid_permeability(node) = 3.52e-11; resistance_gas_diffusivity(node) = 2*0.371/(3-0.371); } else { porosity(node) = 0.387; liquid_diffusivity(node) = 0; liquid_permeability(node) = 8.41e-12; resistance_gas_diffusivity(node) = 2*0.387/(3-0.387); } } porosity.predictor = porosity.variable; true_vars->set_relative_variables(); } solver::StaggerReturnCode DryingStagger::restart_timestep( solver::VariablesBase * const var) { UnsaturatedVariables * const true_vars = static_cast(var); for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { compute_one_node(node, true_vars); } return solver::StaggerReturnCode::ResidualMinimized; } void DryingStagger::compute_one_node( index_t node, UnsaturatedVariables * const vars) { MainVariable& satvars = vars->get_liquid_saturation(); MainVariable& presvars = vars->get_pressure_main_variables(0); scalar_t rt = vars->get_rt(); scalar_t saturation = satvars(node); scalar_t pressure = presvars(node); scalar_t porosity = vars->get_porosity()(node); auto vapor_pressure_f = vars->get_vapor_pressure_model(); const scalar_t tot_conc = cw_tilde_si*saturation + (1.0-saturation)*pressure/rt; pressure = vapor_pressure_f(node, saturation); scalar_t cw_hat = pressure/rt; saturation = (tot_conc - cw_hat) / (cw_tilde_si-cw_hat); pressure = vapor_pressure_f(node, saturation); cw_hat = pressure/rt; scalar_t res = tot_conc - cw_tilde_si*saturation - (1.0-saturation)*cw_hat; while (std::abs(res)/tot_conc > 1e-12) { saturation = (tot_conc - cw_hat) / (cw_tilde_si-cw_hat); pressure = vapor_pressure_f(node, saturation); cw_hat = pressure/rt; res = tot_conc - cw_tilde_si*saturation - (1.0-saturation)*cw_hat; } satvars(node) = saturation; satvars.velocity(node) = (saturation - satvars.predictor(node))/m_dt; presvars(node) = pressure; presvars.velocity(node) = (pressure - presvars.predictor(node))/m_dt; presvars.chemistry_rate(node) = - porosity/(rt*m_dt) * ( (pressure*(1.0-saturation)) - (presvars.predictor(node)*(1.0-satvars.predictor(node))) ) + presvars.transport_fluxes(node); } database::RawDatabasePtr get_database() { specmicp::database::Database thedatabase(CEMDATA_PATH); std::vector to_keep = {"H2O", "H[+]"}; thedatabase.keep_only_components(to_keep); thedatabase.remove_gas_phases(); database::RawDatabasePtr raw_data = thedatabase.get_database(); raw_data->freeze_db(); return raw_data; } mesh::Mesh1DPtr get_mesh() { mesh::Uniform1DMeshGeometry geom; geom.dx = DX; geom.nb_nodes = NB_NODES; geom.section = CROSS_SECTION; return mesh::uniform_mesh1d(geom); } scalar_t big_bead_cap_pressure(scalar_t saturation) { return GAMMA_WATER_AIR/std::sqrt(3.52e-11/0.371)*leverett_function(saturation); } scalar_t small_bead_cap_pressure(scalar_t saturation) { return GAMMA_WATER_AIR/std::sqrt(8.41e-12/0.387)*leverett_function(saturation); } scalar_t vapor_pressure(scalar_t pc) { return 3e3*std::exp(-M_v * pc / (rho_l_si * (8.314*(25+273.15)))); } scalar_t leverett_function(scalar_t s_eff) { return 0.325*std::pow((1.0/s_eff - 1), 0.217); } class Printer { public: Printer(UnsaturatedVariables* vars, std::size_t reserve_size); void register_timestep(scalar_t time, solver::VariablesBasePtr base_vars); void print(const std::string& name); private: mesh::Mesh1DPtr m_mesh; std::vector> m_pressure; std::vector> m_saturation; std::vector m_timestep; }; Printer::Printer(UnsaturatedVariables *vars, std::size_t reserve_size): m_mesh(vars->get_mesh()), m_pressure(m_mesh->nb_nodes()), m_saturation(m_mesh->nb_nodes()) { m_timestep.reserve(reserve_size); for (auto node: m_mesh->range_nodes()) { m_pressure[node].reserve(reserve_size); m_saturation[node].reserve(reserve_size); } } void Printer::register_timestep(scalar_t time, solver::VariablesBasePtr base_vars) { UnsaturatedVariables* vars = static_cast(base_vars.get()); m_timestep.push_back(time); MainVariable& saturation = vars->get_liquid_saturation(); MainVariable& pressure = vars->get_pressure_main_variables(0); for (auto node: m_mesh->range_nodes()) { m_pressure[node].push_back(pressure(node)); m_saturation[node].push_back(saturation(node)); } } void Printer::print(const std::string& name) { io::CSVFile pressure_file(name+"_pressure.dat"); io::CSVFile saturation_file(name+"_saturation.dat"); pressure_file << "Node"; saturation_file << "Node"; for (auto& time: m_timestep) { pressure_file.separator(); pressure_file << time; saturation_file.separator(); saturation_file << time; } pressure_file.eol(); saturation_file.eol(); for (auto node: m_mesh->range_nodes()) { pressure_file << m_mesh->get_position(node); saturation_file << m_mesh->get_position(node); for (std::size_t ind=0; ind has_gas = {true, false, false}; std::vector is_big_bead(the_mesh->nb_nodes(), true); //for (int node=25; node( the_mesh, raw_data, has_gas, units::LengthUnit::meter); - auto bcs = BoundaryConditions::make(the_mesh->nb_nodes()); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); bcs->add_gas_node(0); bcs->fork_constraint(0, "gas_node"); // variables initialisation MainVariable& saturation = vars->get_liquid_saturation(); MainVariable& vapor_pressure = vars->get_pressure_main_variables(0); SecondaryTransientVariable& water_aqueous_concentration = vars->get_water_aqueous_concentration(); water_aqueous_concentration.predictor = water_aqueous_concentration.variable; SecondaryVariable& rel_gas_d = vars->get_relative_gas_diffusivity(); rel_gas_d(0) = 1.0; vapor_pressure(0) = 100; for (index_t node = 1; node < the_mesh->nb_nodes(); ++node) { saturation(node) = 0.9; } water_aqueous_concentration.set_constant(cw_tilde_si); std::shared_ptr upscaling_stagger = std::make_shared(is_big_bead); upscaling_stagger->initialize(vars.get()); std::shared_ptr drying_stagger = std::make_shared(); // init vapor pressure drying_stagger->initialize_timestep(1.0, vars.get()); drying_stagger->restart_timestep(vars.get()); std::shared_ptr transport_stagger = std::make_shared(vars, bcs); solver::ReactiveTransportSolver react_solver(transport_stagger, drying_stagger, upscaling_stagger); solver::ReactiveTransportOptions& opts = react_solver.get_options(); opts.residuals_tolerance = 1e-2; opts.good_enough_tolerance = 0.99; vapor_pressure.chemistry_rate.setZero(); //for (auto it=0; it<100; ++it) //{ // auto retcode = react_solver.solve_timestep(0.1, vars); // std::cout << (int) retcode << std::endl; //} //std::cout << saturation.variable << std::endl; //std::cout << vapor_pressure.variable << std::endl; //std::cout << vars->get_gas_diffusivity().variable << std::endl; //std::cout << vars->get_relative_gas_diffusivity().variable << std::endl; int nb_hours = 9; Printer printer(vars.get(), 4*nb_hours); solver::SimulationInformation simul_info("drying", 900); solver::ReactiveTransportRunner runner(react_solver, 0.01, 10, simul_info); solver::output_f out_func = std::bind(&Printer::register_timestep, &printer, std::placeholders::_1, std::placeholders::_2); runner.set_output_policy(out_func); runner.run_until(nb_hours*3600, vars); printer.print("out_drying"); std::cout << saturation.variable << std::endl; std::cout << vapor_pressure.variable << std::endl; } diff --git a/src/reactmicp/io/configuration_unsaturated.cpp b/src/reactmicp/io/configuration_unsaturated.cpp index 2a082cd..45c26b9 100644 --- a/src/reactmicp/io/configuration_unsaturated.cpp +++ b/src/reactmicp/io/configuration_unsaturated.cpp @@ -1,218 +1,218 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "specmicp_common/io/safe_config.hpp" #include "configuration_unsaturated.hpp" #include "specmicp_database/data_container.hpp" #include "specmicp/adimensional/adimensional_system_solver_structs.hpp" #include "specmicp/io/configuration.hpp" #include "reactmicp/systems/unsaturated/boundary_conditions.hpp" #include "reactmicp/systems/unsaturated/equilibrium_stagger.hpp" #include "specmicp_common/string_algorithms.hpp" #define S_BCS_A_GASNODE "gas_nodes" #define S_BCS_A_FIXEDNODE "fixed_nodes" #define S_BCS_SS_CONSTRAINTS "constraints" #define S_BCS_SS_CONSTRAINTS_A_DEFAULT "default" #define S_BCS_SS_CONSTRAINTS_A_OTHERS "others" #define S_BCS_SS_CONSTRAINTS_A_NAME "name" #define S_BCS_SS_CONSTRAINTS_A_FORK_FROM "fork_from" #define S_BCS_SS_CONSTRAINTS_A_NODE "nodes" #define S_OPTIONS_SS_DEFAULT "default" #define S_OPTIONS_SS_OTHERS "others" #define S_OPTIONS_SS_OTHERS_A_NAME "name" #define S_OPTIONS_SS_OTHERS_A_FORK_FROM "fork_from" #define S_OPTIONS_SS_OTHERS_A_NODE "nodes" namespace specmicp { using reactmicp::systems::unsaturated::BoundaryConditions; using reactmicp::systems::unsaturated::EquilibriumOptionsVector; namespace io { static void configure_bcs_constraints( BoundaryConditions* bcs, const database::DataContainer * const raw_data, YAMLConfigHandle&& conf); //! \brief Configure the boundary condition std::shared_ptr configure_unsaturated_boundary_conditions( uindex_t nb_nodes, const database::DataContainer * const raw_data, YAMLConfigHandle&& configuration ) { - auto bcs = BoundaryConditions::make(nb_nodes); + auto bcs = BoundaryConditions::make(nb_nodes, raw_data->nb_component()); if (configuration.has_node(S_BCS_A_FIXEDNODE)) { std::vector fixed_nodes = configuration.list_to_vector(S_BCS_A_FIXEDNODE); bcs->add_fixed_nodes(fixed_nodes); } if (configuration.has_node(S_BCS_A_GASNODE)) { const std::vector gas_nodes = configuration.list_to_vector(S_BCS_A_GASNODE); bcs->add_gas_nodes(gas_nodes); } if (configuration.has_section(S_BCS_SS_CONSTRAINTS)) { configure_bcs_constraints(bcs.get(), raw_data, configuration.get_section(S_BCS_SS_CONSTRAINTS) ); } return bcs; } void configure_bcs_constraints( BoundaryConditions* bcs, const database::DataContainer * const raw_data, YAMLConfigHandle&& conf ) { // default constraints if (conf.has_section(S_BCS_SS_CONSTRAINTS_A_DEFAULT)) { configure_specmicp_constraints( bcs->get_constraint("default"), raw_data, conf.get_section(S_BCS_SS_CONSTRAINTS_A_DEFAULT) ); } // other constraints if (conf.has_section(S_BCS_SS_CONSTRAINTS_A_OTHERS)) { auto others = conf.get_section(S_BCS_SS_CONSTRAINTS_A_OTHERS); for (uindex_t ind=0; ind( S_BCS_SS_CONSTRAINTS_A_NAME ); if (subconf.has_attribute(S_BCS_SS_CONSTRAINTS_A_FORK_FROM)) { auto old_name = subconf.get_attribute( S_BCS_SS_CONSTRAINTS_A_FORK_FROM ); if (not bcs->has_constraint(old_name)) { subconf.report_error( YAMLConfigError::InvalidArgument, "'" + old_name + "' is not an " "existing constraint" ); } bcs->fork_constraint(old_name, name); } else { bcs->add_constraint(name); } if (subconf.has_node(S_BCS_SS_CONSTRAINTS_A_NODE)) { auto nodes = subconf.list_to_vector( S_BCS_SS_CONSTRAINTS_A_NODE); for (const auto& node: nodes) { bcs->set_constraint(node, name); } } configure_specmicp_constraints( bcs->get_constraint(name), raw_data, std::move(subconf)); } } } std::shared_ptr configure_unsaturated_equilibrium_options( uindex_t nb_nodes, const units::UnitsSet& the_units, YAMLConfigHandle&& configuration ) { auto opts = EquilibriumOptionsVector::make(nb_nodes); if (configuration.has_section(S_OPTIONS_SS_DEFAULT)) { auto& def_opts = opts->get("default"); configure_specmicp_options( def_opts, the_units, configuration.get_section(S_OPTIONS_SS_DEFAULT)); } if (configuration.has_section(S_OPTIONS_SS_OTHERS)) { auto subconf = configuration.get_section(S_OPTIONS_SS_OTHERS); const auto size = subconf.size(); for (uindex_t ind=0; ind( S_OPTIONS_SS_OTHERS_A_NAME); if (conf.has_attribute(S_OPTIONS_SS_OTHERS_A_FORK_FROM)) { const auto fork_from = conf.get_attribute( S_OPTIONS_SS_OTHERS_A_FORK_FROM); if (not opts->has_value(fork_from)) { conf.report_error( YAMLConfigError::InvalidArgument, "'" + fork_from + "' is not an " "existing set of options. \n" ); } opts->fork(fork_from, name); } else { opts->push_back_cache(name, AdimensionalSystemSolverOptions()); } if (conf.has_node(S_OPTIONS_SS_OTHERS_A_NODE)) { const auto nodes = conf.list_to_vector( S_OPTIONS_SS_OTHERS_A_NODE); for (const auto& node: nodes) { opts->set(node, name); } } auto& new_opts = opts->get(name); configure_specmicp_options( new_opts, the_units, std::move(conf)); } } return opts; } } //end namespace io } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/CMakeLists.txt b/src/reactmicp/systems/unsaturated/CMakeLists.txt index e791858..6684767 100644 --- a/src/reactmicp/systems/unsaturated/CMakeLists.txt +++ b/src/reactmicp/systems/unsaturated/CMakeLists.txt @@ -1,43 +1,42 @@ set( reactmicp_unsaturated_srcs aqueous_equation.cpp aqueous_pressure_equation.cpp boundary_conditions.cpp equilibrium_stagger.cpp pressure_equation.cpp saturation_equation.cpp saturation_pressure_equation.cpp transport_stagger.cpp variables.cpp variables_interface.cpp ) set( reactmicp_unsaturated_headers fv_1dof_equation.hpp fv_1dof_equation.inl init_variables.hpp simulation_data.hpp - transport_constraints.hpp types_fwd.hpp upscaling_stagger_factory.hpp variables_box.hpp variables_sub.hpp aqueous_equation.hpp aqueous_pressure_equation.hpp boundary_conditions.hpp equilibrium_stagger.hpp pressure_equation.hpp saturation_equation.hpp saturation_pressure_equation.hpp transport_stagger.hpp variables.hpp variables_interface.hpp ) add_to_reactmicp_srcs_list(reactmicp_unsaturated_srcs) add_to_reactmicp_headers_list(reactmicp_unsaturated_headers) INSTALL(FILES ${reactmicp_unsaturated_headers} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/systems/unsaturated ) diff --git a/src/reactmicp/systems/unsaturated/aqueous_equation.cpp b/src/reactmicp/systems/unsaturated/aqueous_equation.cpp index 64b6f41..0579746 100644 --- a/src/reactmicp/systems/unsaturated/aqueous_equation.cpp +++ b/src/reactmicp/systems/unsaturated/aqueous_equation.cpp @@ -1,310 +1,351 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "aqueous_equation.hpp" #include "variables_box.hpp" -#include "transport_constraints.hpp" +#include "boundary_conditions.hpp" #include "dfpm/meshes/mesh1d.hpp" #include "specmicp_common/compat.hpp" #include "specmicp_common/physics/constants.hpp" #include "specmicp_common/physics/maths.hpp" #include "dfpm/solver/parabolic_driver.hpp" namespace specmicp { namespace dfpmsolver { // explicit template instanciation template class ParabolicDriver; } //end namespace dfpmsolver namespace reactmicp { namespace systems { namespace unsaturated { static constexpr index_t no_equation {-1}; static constexpr index_t no_eq_no_var {-2}; static constexpr index_t not_initialized {-5}; struct SPECMICP_DLL_LOCAL AqueousTransportEquation::AqueousTransportEquationImpl { + uindex_t m_id_component; mesh::Mesh1DPtr m_mesh; LiquidAqueousComponentVariableBox m_vars; std::vector m_ideq; + std::shared_ptr m_bcs; mesh::Mesh1D* mesh() {return m_mesh.get();} LiquidAqueousComponentVariableBox& vars() {return m_vars;} bool node_has_equation(index_t node) { return m_ideq[node] > no_equation; } bool node_can_flux(index_t node) { return m_ideq[node] > no_eq_no_var; } index_t id_equation(index_t node) { return m_ideq[node]; } void fix_node(index_t node) { m_ideq[node] = no_equation; } void gas_node(index_t node) { m_ideq[node] = no_eq_no_var; } AqueousTransportEquationImpl( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, - LiquidAqueousComponentVariableBox the_vars + LiquidAqueousComponentVariableBox the_vars, + std::shared_ptr bcs ): + m_id_component(id_component), m_mesh(the_mesh), m_vars(the_vars), - m_ideq(the_mesh->nb_nodes(), not_initialized) + m_ideq(the_mesh->nb_nodes(), not_initialized), + m_bcs(bcs) {} - void compute_transport_rate(scalar_t dt, const Vector& displacment); + void compute_transport_rate(scalar_t dt, const Vector& displacement); + + + scalar_t get_bc_liquid_flux(uindex_t node) { + return m_bcs->get_flux_liquid_dof(node, m_id_component); + } + + scalar_t get_diff( + uindex_t node_0, + uindex_t node_1, + const scalar_t& adv + ) const; + + uindex_t number_equations(); }; AqueousTransportEquation::AqueousTransportEquation( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, LiquidAqueousComponentVariableBox& variables, - const TransportConstraints& constraints): + std::shared_ptr bcs): base(the_mesh->nb_nodes()), - m_impl(make_unique(the_mesh, variables)) + m_impl(make_unique( + id_component, the_mesh, variables, bcs)) { - number_equations(constraints); + number_equations(); } AqueousTransportEquation::~AqueousTransportEquation() = default; index_t AqueousTransportEquation::id_equation_impl(index_t id_dof) { const auto id_eq = m_impl->id_equation(id_dof); return (id_eq>no_equation)?id_eq:no_equation; } mesh::Mesh1D* AqueousTransportEquation::get_mesh_impl() { return m_impl->mesh(); } void AqueousTransportEquation::pre_nodal_residual_hook_impl( index_t node, const Vector& displacement) {} void AqueousTransportEquation::pre_residual_hook_impl(const Vector& displacement) {} void AqueousTransportEquation::post_residual_hook_impl(const Vector& displacement) {} +scalar_t AqueousTransportEquation::AqueousTransportEquationImpl::get_diff( + uindex_t node_0, + uindex_t node_1, + const scalar_t& adv + ) const +{ + scalar_t diff; + if (adv > 0 ) { + diff = m_vars.liquid_diffusivity(node_0) + * m_vars.relative_liquid_diffusivity(node_0); + } else if (adv < 0) { + diff = m_vars.liquid_diffusivity(node_1) + * m_vars.relative_liquid_diffusivity(node_1); + } else { + const scalar_t diff_0 = m_vars.liquid_diffusivity(node_0) + * m_vars.relative_liquid_diffusivity(node_0); + const scalar_t diff_1 = m_vars.liquid_diffusivity(node_1) + * m_vars.relative_liquid_diffusivity(node_1); + diff = average(diff_0, diff_1); + } + return diff; +} + //! \brief Compute the residuals inside 'element' void AqueousTransportEquation::residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ) { element_residual.setZero(); mesh::Mesh1D* m_mesh = m_impl->mesh(); LiquidAqueousComponentVariableBox& vars = m_impl->m_vars; const scalar_t mass_coeff_0 = m_mesh->get_volume_cell_element(element, 0); const scalar_t mass_coeff_1 = m_mesh->get_volume_cell_element(element, 1); const index_t node_0 = m_mesh->get_node(element, 0); const index_t node_1 = m_mesh->get_node(element, 1); scalar_t flux_0 = 0.0; scalar_t flux_1 = 0.0; if (m_impl->node_can_flux(node_0) and m_impl->node_can_flux(node_1)) { + const scalar_t adv_vel = vars.advection_flux(element); // Diffusion Cw - const scalar_t coeff_diff_0 = vars.liquid_diffusivity(node_0) - * vars.relative_liquid_diffusivity(node_0); - const scalar_t coeff_diff_1 = vars.liquid_diffusivity(node_1) - * vars.relative_liquid_diffusivity(node_1); - - const scalar_t coeff_diff = average( - coeff_diff_0, coeff_diff_1); + const scalar_t coeff_diff = m_impl->get_diff(node_0, node_1, adv_vel); const scalar_t diff_aq = (displacement(node_1) - displacement(node_0) ); const scalar_t diff_flux = coeff_diff * diff_aq/ m_mesh->get_dx(element); // advection - if (vars.advection_flux(element) < 0) + if (adv_vel < 0) { - flux_0 = -vars.advection_flux(element)*diff_aq; + flux_0 = -adv_vel*diff_aq; } - else if (vars.advection_flux(element) > 0) + else if (adv_vel > 0) { - flux_1 = -vars.advection_flux(element)*diff_aq; + flux_1 = -adv_vel*diff_aq; } - flux_0 += diff_flux; + flux_0 += diff_flux + m_impl->get_bc_liquid_flux(node_0); flux_0 *= m_mesh->get_face_area(element); - flux_1 -= diff_flux; + flux_1 -= diff_flux + m_impl->get_bc_liquid_flux(node_1); flux_1 *= m_mesh->get_face_area(element); } // transient if (m_impl->node_has_equation(node_0)) { const scalar_t porosity_0 = vars.porosity(node_0); const scalar_t aq_tot_conc_0 = displacement(node_0); const scalar_t saturation_0 = vars.saturation(node_0); const scalar_t transient_0 = porosity_0 * aq_tot_conc_0 * vars.saturation.velocity(node_0) + saturation_0 * aq_tot_conc_0 * vars.porosity.velocity(node_0) + porosity_0 * saturation_0 * velocity(node_0); scalar_t res = mass_coeff_0*transient_0 - flux_0; if (use_chemistry_rate) { const scalar_t chemistry_0 = vars.aqueous_concentration.chemistry_rate(node_0) + vars.solid_concentration.chemistry_rate(node_0) + vars.partial_pressure.chemistry_rate(node_0); res -= mass_coeff_0*chemistry_0; } element_residual(0) = res/get_scaling(); - } if (m_impl->node_has_equation(node_1)) { const scalar_t porosity_1 = vars.porosity(node_1); const scalar_t aq_tot_conc_1 = displacement(node_1); const scalar_t saturation_1 = vars.saturation(node_1); const scalar_t transient_1 = porosity_1 * aq_tot_conc_1 * vars.saturation.velocity(node_1) + saturation_1 * aq_tot_conc_1 * vars.porosity.velocity(node_1) + porosity_1 * saturation_1 * velocity(node_1); scalar_t res = mass_coeff_1*transient_1 - flux_1; if (use_chemistry_rate) { const scalar_t chemistry_1 = vars.aqueous_concentration.chemistry_rate(node_1) + vars.solid_concentration.chemistry_rate(node_1) + vars.partial_pressure.chemistry_rate(node_1); res -= mass_coeff_1*chemistry_1; } element_residual(1) = res/get_scaling(); } } -void AqueousTransportEquation::number_equations(const TransportConstraints& constraints) +uindex_t AqueousTransportEquation::AqueousTransportEquationImpl::number_equations() { - for (int fixed_node: constraints.fixed_nodes()) - { - m_impl->fix_node(fixed_node); - } - for (int gas_node: constraints.gas_nodes()) - { - m_impl->gas_node(gas_node); - } - index_t neq = 0; - for (index_t node: m_impl->mesh()->range_nodes()) - { - if (m_impl->m_ideq[node] == not_initialized) - { - m_impl->m_ideq[node] = neq; + uindex_t neq = 0; + for (index_t node: m_mesh->range_nodes()){ + switch (m_bcs->get_bcs_liquid_dof(node, m_id_component)) { + case IdBCs::NoEquation: + gas_node(node); + break; + case IdBCs::FixedDof: + fix_node(node); + break; + default: + m_ideq[node] = neq; ++neq; } } + return neq; +} + +void AqueousTransportEquation::number_equations() +{ + auto neq = m_impl->number_equations(); register_number_equations(neq); } void AqueousTransportEquation::compute_transport_rate( scalar_t dt, const Vector& displacement ) { m_impl->compute_transport_rate(dt, displacement); } void AqueousTransportEquation::AqueousTransportEquationImpl::compute_transport_rate( scalar_t dt, const Vector& displacement) { MainVariable& aqueous_concentration = m_vars.aqueous_concentration; const MainVariable& saturation = m_vars.saturation; const MainVariable& solid_conc = m_vars.solid_concentration; const MainVariable& pressure = m_vars.partial_pressure; const SecondaryTransientVariable& porosity = m_vars.porosity; for (index_t node: m_mesh->range_nodes()) { if (! node_has_equation(node)) continue; const scalar_t transient = ( ( porosity(node) * saturation(node) * displacement(node)) - ( porosity.predictor(node) * aqueous_concentration.predictor(node) * saturation.predictor(node)) ) / dt; const scalar_t chem_rates = ( saturation.chemistry_rate(node) + solid_conc.chemistry_rate(node) + pressure.chemistry_rate(node) ); aqueous_concentration.transport_fluxes(node) = transient - chem_rates; } } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/aqueous_equation.hpp b/src/reactmicp/systems/unsaturated/aqueous_equation.hpp index 953efb9..a17f7dc 100644 --- a/src/reactmicp/systems/unsaturated/aqueous_equation.hpp +++ b/src/reactmicp/systems/unsaturated/aqueous_equation.hpp @@ -1,104 +1,106 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_UNSATURATED_AQUEOUSEQUATION_HPP #define SPECMICP_REACTMICP_UNSATURATED_AQUEOUSEQUATION_HPP //! \file unsaturated/aqueous_equation.hpp //! \brief The equation for the liquid transport of aqueous component #include "specmicp_common/types.hpp" #include "variables.hpp" #include "fv_1dof_equation.hpp" #include "dfpm/meshes/mesh1dfwd.hpp" #include namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { -class TransportConstraints; +class BoundaryConditions; //! \brief Equation for liquid transport of aqueous component class SPECMICP_DLL_LOCAL AqueousTransportEquation: public FV1DOFEquation { using base = FV1DOFEquation; using base::get_scaling; using base::register_number_equations; public: AqueousTransportEquation( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, LiquidAqueousComponentVariableBox& variables, - const TransportConstraints& constraints); + std::shared_ptr constraints + ); ~AqueousTransportEquation(); //! \brief Compute the residuals inside 'element' void residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ); index_t id_equation_impl(index_t id_dof); mesh::Mesh1D* get_mesh_impl(); void pre_nodal_residual_hook_impl(index_t node, const Vector& displacement); void pre_residual_hook_impl(const Vector& displacement); void post_residual_hook_impl(const Vector& displacement); void compute_transport_rate(scalar_t dt, const Vector& displacement); private: - void number_equations(const TransportConstraints& constraints); + void number_equations(); struct AqueousTransportEquationImpl; std::unique_ptr m_impl; }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_UNSATURATED_AQUEOUSEQUATION_HPP diff --git a/src/reactmicp/systems/unsaturated/aqueous_pressure_equation.cpp b/src/reactmicp/systems/unsaturated/aqueous_pressure_equation.cpp index 571df21..9d2a944 100644 --- a/src/reactmicp/systems/unsaturated/aqueous_pressure_equation.cpp +++ b/src/reactmicp/systems/unsaturated/aqueous_pressure_equation.cpp @@ -1,325 +1,414 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "aqueous_pressure_equation.hpp" #include "variables_box.hpp" -#include "transport_constraints.hpp" +#include "boundary_conditions.hpp" #include "dfpm/meshes/mesh1d.hpp" #include "specmicp_common/compat.hpp" #include "specmicp_common/physics/constants.hpp" #include "specmicp_common/physics/maths.hpp" #include "dfpm/solver/parabolic_driver.hpp" namespace specmicp { namespace dfpmsolver { // explicit template instanciation template class ParabolicDriver; } //end namespace dfpmsolver namespace reactmicp { namespace systems { namespace unsaturated { -static constexpr index_t no_equation {-1}; -static constexpr index_t no_eq_no_var {-2}; +static constexpr index_t no_equation {-1}; +static constexpr index_t no_eq_no_var {-2}; static constexpr index_t not_initialized {-5}; struct SPECMICP_DLL_LOCAL AqueousGasTransportEquation::AqueousGasTransportEquationImpl { + uindex_t m_id_component; mesh::Mesh1DPtr m_mesh; LiquidGasAqueousVariableBox m_vars; std::vector m_ideq; + std::shared_ptr m_bcs; mesh::Mesh1D* mesh() {return m_mesh.get();} LiquidAqueousComponentVariableBox& vars() {return m_vars;} bool node_has_equation(index_t node) { return m_ideq[node] > no_equation; } bool node_can_flux(index_t node) { return m_ideq[node] > no_eq_no_var; } index_t id_equation(index_t node) { return m_ideq[node]; } void fix_node(index_t node) { m_ideq[node] = no_equation; } void gas_node(index_t node) { m_ideq[node] = no_eq_no_var; } AqueousGasTransportEquationImpl( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, - LiquidGasAqueousVariableBox the_vars + LiquidGasAqueousVariableBox the_vars, + std::shared_ptr bcs ): + m_id_component(id_component), m_mesh(the_mesh), m_vars(the_vars), - m_ideq(the_mesh->nb_nodes(), not_initialized) + m_ideq(the_mesh->nb_nodes(), not_initialized), + m_bcs(bcs) {} - void compute_transport_rate(scalar_t dt, const Vector& displacment); + void compute_transport_rate(scalar_t dt, const Vector& displacement); + + scalar_t get_bc_liquid_flux(uindex_t node) { + return m_bcs->get_flux_liquid_dof(node, m_id_component); + } + scalar_t get_bc_gas_flux(uindex_t node) { + return m_bcs->get_flux_gas_dof(node, m_id_component); + } + + scalar_t get_diff( + uindex_t node_0, + uindex_t node_1, + const scalar_t& adv + ) const; + + uindex_t number_equations(); + + scalar_t get_diff_gas( + uindex_t node_0, + uindex_t node_1, + const scalar_t& adv + ) const; + }; -AqueousGasTransportEquation::AqueousGasTransportEquation(mesh::Mesh1DPtr the_mesh, +AqueousGasTransportEquation::AqueousGasTransportEquation( + uindex_t id_component, + mesh::Mesh1DPtr the_mesh, LiquidGasAqueousVariableBox& variables, - const TransportConstraints& constraints): + std::shared_ptr bcs + ): base(the_mesh->nb_nodes()), - m_impl(make_unique(the_mesh, variables)) + m_impl(make_unique( + id_component, the_mesh, variables, bcs)) { - number_equations(constraints); + number_equations(); } AqueousGasTransportEquation::~AqueousGasTransportEquation() = default; index_t AqueousGasTransportEquation::id_equation_impl(index_t id_dof) { const auto id_eq = m_impl->id_equation(id_dof); return (id_eq>no_equation)?id_eq:no_equation; } mesh::Mesh1D* AqueousGasTransportEquation::get_mesh_impl() { return m_impl->mesh(); } void AqueousGasTransportEquation::pre_nodal_residual_hook_impl( index_t node, const Vector& displacement) {} void AqueousGasTransportEquation::pre_residual_hook_impl(const Vector& displacement) {} void AqueousGasTransportEquation::post_residual_hook_impl(const Vector& displacement) {} +scalar_t AqueousGasTransportEquation::AqueousGasTransportEquationImpl::get_diff( + uindex_t node_0, + uindex_t node_1, + const scalar_t& adv + ) const +{ + scalar_t diff; + if (adv > 0 ) { + diff = m_vars.liquid_diffusivity(node_0) + * m_vars.relative_liquid_diffusivity(node_0); + } else if (adv < 0) { + diff = m_vars.liquid_diffusivity(node_1) + * m_vars.relative_liquid_diffusivity(node_1); + } else { + const scalar_t diff_0 = m_vars.liquid_diffusivity(node_0) + * m_vars.relative_liquid_diffusivity(node_0); + const scalar_t diff_1 = m_vars.liquid_diffusivity(node_1) + * m_vars.relative_liquid_diffusivity(node_1); + diff = average(diff_0, diff_1); + } + return diff; +} + + +scalar_t AqueousGasTransportEquation::AqueousGasTransportEquationImpl::get_diff_gas( + uindex_t node_0, + uindex_t node_1, + const scalar_t& adv + ) const +{ + scalar_t diff; + if (adv > 0) { + diff = m_vars.resistance_gas_diffusivity(node_0) + * m_vars.relative_gas_diffusivity(node_0); + } else if (adv < 0) { + diff = m_vars.resistance_gas_diffusivity(node_1) + * m_vars.relative_gas_diffusivity(node_1); + } else { +// if (m_bcs->is_gas_node(node_0)) { +// diff = m_vars.resistance_gas_diffusivity(node_0) +// * m_vars.relative_gas_diffusivity(node_0); +// } else if (m_bcs->is_gas_node(node_1)) { +// diff = m_vars.resistance_gas_diffusivity(node_1) +// * m_vars.relative_gas_diffusivity(node_1); +// } else { + const scalar_t coeff_diff_0 = m_vars.resistance_gas_diffusivity(node_0) + * m_vars.relative_gas_diffusivity(node_0); + const scalar_t coeff_diff_1 = m_vars.resistance_gas_diffusivity(node_1) + * m_vars.relative_gas_diffusivity(node_1); + + diff = average( + coeff_diff_0, coeff_diff_1); +// } + } + return diff * m_vars.binary_diffusion_coefficient; + +} + //! \brief Compute the residuals inside 'element' void AqueousGasTransportEquation::residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ) { element_residual.setZero(); mesh::Mesh1D* m_mesh = m_impl->mesh(); LiquidGasAqueousVariableBox& vars = m_impl->m_vars; const scalar_t& rt = vars.constants.rt; const scalar_t mass_coeff_0 = m_mesh->get_volume_cell_element(element, 0); const scalar_t mass_coeff_1 = m_mesh->get_volume_cell_element(element, 1); const index_t node_0 = m_mesh->get_node(element, 0); const index_t node_1 = m_mesh->get_node(element, 1); scalar_t flux_0 = 0.0; scalar_t flux_1 = 0.0; + const scalar_t adv_vel = vars.advection_flux(element); + if (m_impl->node_can_flux(node_0) and m_impl->node_can_flux(node_1)) { // Diffusion Cw - const scalar_t coeff_diff_0 = vars.liquid_diffusivity(node_0) - * vars.relative_liquid_diffusivity(node_0); - const scalar_t coeff_diff_1 = vars.liquid_diffusivity(node_1) - * vars.relative_liquid_diffusivity(node_1); - - const scalar_t coeff_diff = average( - coeff_diff_0, coeff_diff_1); - - const scalar_t diff_aq = (displacement(node_1) + const scalar_t coeff_diff = m_impl->get_diff(node_0, node_1, adv_vel); + const scalar_t diff_aq = ( displacement(node_1) - displacement(node_0) - ); - const scalar_t diff_flux = coeff_diff * diff_aq/ m_mesh->get_dx(element); + ); + const scalar_t diff_flux = coeff_diff * diff_aq / m_mesh->get_dx(element); // advection - if (vars.advection_flux(element) < 0) + if (adv_vel < 0) { - flux_0 = -vars.advection_flux(element)*diff_aq; + flux_0 = -adv_vel*diff_aq; } - else if (vars.advection_flux(element) > 0) + else if (adv_vel > 0) { - flux_1 = -vars.advection_flux(element)*diff_aq; + flux_1 = -adv_vel*diff_aq; } flux_0 += diff_flux; flux_1 -= diff_flux; } - const scalar_t coeff_diff_gas_0 = vars.resistance_gas_diffusivity(node_0) - * vars.relative_gas_diffusivity(node_0); - const scalar_t coeff_diff_gas_1 = vars.resistance_gas_diffusivity(node_1) - * vars.relative_gas_diffusivity(node_1); - - const scalar_t coeff_diff_gas = vars.binary_diffusion_coefficient * - average(coeff_diff_gas_0, coeff_diff_gas_1); + const scalar_t coeff_diff_gas = m_impl->get_diff_gas(node_0, node_1, adv_vel); const scalar_t diff_flux_gas = coeff_diff_gas*( vars.partial_pressure(node_1) - vars.partial_pressure(node_0)) / m_mesh->get_dx(element) / rt; + const scalar_t section = m_mesh->get_face_area(element); + flux_0 += diff_flux_gas; - flux_0 *= m_mesh->get_face_area(element); + flux_0 += m_impl->get_bc_gas_flux(node_0) / rt; + flux_0 += m_impl->get_bc_liquid_flux(node_0); + flux_0 *= section; + flux_1 -= diff_flux_gas; - flux_1 *= m_mesh->get_face_area(element); + flux_1 += m_impl->get_bc_gas_flux(node_1) / rt; + flux_1 += m_impl->get_bc_liquid_flux(node_1); + flux_1 *= section; // transient if (m_impl->node_has_equation(node_0)) { - const scalar_t porosity_0 = vars.porosity(node_0); + const scalar_t porosity_0 = vars.porosity(node_0); const scalar_t aq_tot_conc_0 = displacement(node_0); - const scalar_t saturation_0 = vars.saturation(node_0); + const scalar_t saturation_0 = vars.saturation(node_0); const scalar_t transient_0 = porosity_0 * aq_tot_conc_0 * vars.saturation.velocity(node_0) + saturation_0 * aq_tot_conc_0 * vars.porosity.velocity(node_0) + porosity_0 * saturation_0 * velocity(node_0); scalar_t res = mass_coeff_0*transient_0 - flux_0; if (use_chemistry_rate) { const scalar_t chemistry_0 = vars.aqueous_concentration.chemistry_rate(node_0) + vars.solid_concentration.chemistry_rate(node_0) + vars.partial_pressure.chemistry_rate(node_0); res -= mass_coeff_0*chemistry_0; } element_residual(0) = res/get_scaling(); - } if (m_impl->node_has_equation(node_1)) { - const scalar_t porosity_1 = vars.porosity(node_1); + const scalar_t porosity_1 = vars.porosity(node_1); const scalar_t aq_tot_conc_1 = displacement(node_1); - const scalar_t saturation_1 = vars.saturation(node_1); + const scalar_t saturation_1 = vars.saturation(node_1); const scalar_t transient_1 = porosity_1 * aq_tot_conc_1 * vars.saturation.velocity(node_1) + saturation_1 * aq_tot_conc_1 * vars.porosity.velocity(node_1) + porosity_1 * saturation_1 * velocity(node_1); scalar_t res = mass_coeff_1*transient_1 - flux_1; if (use_chemistry_rate) { const scalar_t chemistry_1 = vars.aqueous_concentration.chemistry_rate(node_1) + vars.solid_concentration.chemistry_rate(node_1) + vars.partial_pressure.chemistry_rate(node_1); res -= mass_coeff_1*chemistry_1; } element_residual(1) = res/get_scaling(); } - } -void AqueousGasTransportEquation::number_equations(const TransportConstraints& constraints) + +uindex_t AqueousGasTransportEquation::AqueousGasTransportEquationImpl::number_equations() { - for (int fixed_node: constraints.fixed_nodes()) - { - m_impl->fix_node(fixed_node); - } - for (int gas_node: constraints.gas_nodes()) - { - m_impl->gas_node(gas_node); - } - index_t neq = 0; - for (index_t node: m_impl->mesh()->range_nodes()) - { - if (m_impl->m_ideq[node] == not_initialized) - { - m_impl->m_ideq[node] = neq; + uindex_t neq = 0; + for (index_t node: m_mesh->range_nodes()) { + IdBCs bc_l = m_bcs->get_bcs_liquid_dof(node, m_id_component); + switch (bc_l) { + case IdBCs::FixedDof: + m_ideq[node] = no_equation; + break; + case IdBCs::NoEquation: + m_ideq[node] = no_eq_no_var; + break; + default: + m_ideq[node] = neq; ++neq; } } + return neq; +} + +void AqueousGasTransportEquation::number_equations() +{ + uindex_t neq = m_impl->number_equations(); register_number_equations(neq); } void AqueousGasTransportEquation::compute_transport_rate( scalar_t dt, const Vector& displacement ) { m_impl->compute_transport_rate(dt, displacement); } void AqueousGasTransportEquation::AqueousGasTransportEquationImpl::compute_transport_rate( scalar_t dt, const Vector& displacement) { MainVariable& aqueous_concentration = m_vars.aqueous_concentration; const MainVariable& saturation = m_vars.saturation; const MainVariable& solid_conc = m_vars.solid_concentration; const MainVariable& pressure = m_vars.partial_pressure; const SecondaryTransientVariable& porosity = m_vars.porosity; for (index_t node: m_mesh->range_nodes()) { if (! node_has_equation(node)) continue; const scalar_t transient = ( ( porosity(node) * saturation(node) * displacement(node)) - ( porosity.predictor(node) * aqueous_concentration.predictor(node) * saturation.predictor(node)) ) / dt; const scalar_t chem_rates = ( saturation.chemistry_rate(node) + solid_conc.chemistry_rate(node) + pressure.chemistry_rate(node) ); aqueous_concentration.transport_fluxes(node) = transient - chem_rates; } } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/aqueous_pressure_equation.hpp b/src/reactmicp/systems/unsaturated/aqueous_pressure_equation.hpp index 1180c33..1d8d9df 100644 --- a/src/reactmicp/systems/unsaturated/aqueous_pressure_equation.hpp +++ b/src/reactmicp/systems/unsaturated/aqueous_pressure_equation.hpp @@ -1,103 +1,105 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_UNSATURATED_AQUEOUSPRESSUREEQUATION_HPP #define SPECMICP_REACTMICP_UNSATURATED_AQUEOUSPRESSUREEQUATION_HPP //! \file unsaturated/aqueous_pressure_equation.hpp //! \brief The equation for the liquid and gas transport of aqueous component #include "specmicp_common/types.hpp" #include "fv_1dof_equation.hpp" #include "dfpm/meshes/mesh1dfwd.hpp" #include "variables.hpp" #include namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { -class TransportConstraints; +class BoundaryConditions; //! \brief Equation for liquid transport of aqueous component class SPECMICP_DLL_LOCAL AqueousGasTransportEquation: public FV1DOFEquation { using base = FV1DOFEquation; using base::get_scaling; using base::register_number_equations; public: AqueousGasTransportEquation( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, LiquidGasAqueousVariableBox& variables, - const TransportConstraints& constraints); + std::shared_ptr bcs + ); ~AqueousGasTransportEquation(); //! \brief Compute the residuals inside 'element' void residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ); index_t id_equation_impl(index_t id_dof); mesh::Mesh1D* get_mesh_impl(); void pre_nodal_residual_hook_impl(index_t node, const Vector& displacement); void pre_residual_hook_impl(const Vector& displacement); void post_residual_hook_impl(const Vector& displacement); void compute_transport_rate(scalar_t dt, const Vector& displacement); private: - void number_equations(const TransportConstraints& constraints); + void number_equations(); struct AqueousGasTransportEquationImpl; std::unique_ptr m_impl; }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_UNSATURATED_AQUEOUSPRESSUREEQUATION_HPP diff --git a/src/reactmicp/systems/unsaturated/boundary_conditions.cpp b/src/reactmicp/systems/unsaturated/boundary_conditions.cpp index 96e1b60..ffff34f 100644 --- a/src/reactmicp/systems/unsaturated/boundary_conditions.cpp +++ b/src/reactmicp/systems/unsaturated/boundary_conditions.cpp @@ -1,174 +1,284 @@ + /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "boundary_conditions.hpp" #include "specmicp/adimensional/adimensional_system_structs.hpp" #include "specmicp_common/cached_vector.hpp" #include "specmicp_common/compat.hpp" +#include "specmicp_common/eigen/incl_eigen_sparse_core.hpp" + namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { -enum class NodeType { - Normal, - FixedNode, - GasNode +enum DofType: index_t { + Liquid = 0 + , Gas = 1 }; struct BoundaryConditions::BoundaryConditionsImpl { - BoundaryConditionsImpl(uindex_t nb_nodes): - m_types(nb_nodes, NodeType::Normal), + BoundaryConditionsImpl(uindex_t nb_nodes, uindex_t nb_components): + m_nb_nodes(nb_nodes), + m_nb_components(nb_components), + m_id_bcs(2*nb_nodes*nb_components), + m_flux_bcs(2*nb_nodes*nb_components), m_constraints(nb_nodes, "default", AdimensionalSystemConstraints()) - {} + { + m_id_bcs.setZero(); + m_flux_bcs.setZero(); + } + + //index_t cast_from_id(IdBCs the_id) { + // return reinterpret_cast(the_id); + //} + //IdBCs cast_to_id(index_t val) { + // return reinterpret_cast(val); + //} + const index_t& cast_from_id(const IdBCs& the_id) { + return reinterpret_cast(the_id); + } + const IdBCs& cast_to_id(const index_t& val) { + return reinterpret_cast(val); + } + + + index_t id_dof(uindex_t node, uindex_t component, DofType type) const { + return 2*m_nb_components*node+2*component+type; + } + IdBCs& id_bcs(uindex_t node, uindex_t component, DofType type) { + return reinterpret_cast(m_id_bcs(id_dof(node, component, type))); + } + IdBCs id_bcs(uindex_t node, uindex_t component, DofType type) const { + return static_cast(m_id_bcs(id_dof(node, component, type))); + } + scalar_t& flux_bcs(uindex_t node, uindex_t component, DofType type) { + return m_flux_bcs(id_dof(node, component, type)); + } + scalar_t flux_bcs(uindex_t node, uindex_t component, DofType type) const { + return m_flux_bcs(id_dof(node, component, type)); + } + + uindex_t m_nb_nodes; + uindex_t m_nb_components; + + Eigen::Matrix m_id_bcs; + Vector m_flux_bcs; + + //Eigen::SparseVector m_id_bcs; + //Eigen::SparseVector m_flux_bcs; + - std::vector m_types; utils::NameCachedVector m_constraints; }; //BoundaryConditions::BoundaryConditions(): // m_impl(nullptr) //{} -BoundaryConditions::BoundaryConditions(uindex_t nb_nodes): - m_impl(utils::make_pimpl(nb_nodes)) +BoundaryConditions::BoundaryConditions(uindex_t nb_nodes, uindex_t nb_components): + m_impl(utils::make_pimpl(nb_nodes, nb_components)) {} BoundaryConditions::~BoundaryConditions() = default; -void BoundaryConditions::add_gas_node(uindex_t node) { - m_impl->m_types[node] = NodeType::GasNode; +// tranport + +IdBCs BoundaryConditions::get_bcs_liquid_dof( + uindex_t node, + uindex_t component + ) { + return m_impl->id_bcs(node, component, DofType::Liquid); +} +IdBCs BoundaryConditions::get_bcs_gas_dof( + uindex_t node, + uindex_t component + ) { + return m_impl->id_bcs(node, component, DofType::Gas); +} + +scalar_t BoundaryConditions::get_flux_liquid_dof( + uindex_t node, + uindex_t component + ) const { + return m_impl->flux_bcs(node, component, DofType::Liquid); +} +scalar_t BoundaryConditions::get_flux_gas_dof( + uindex_t node, + uindex_t component + ) const { + return m_impl->flux_bcs(node, component, DofType::Gas); +} + +void BoundaryConditions::reset_flux_liquid_dof( + uindex_t node, + uindex_t component + ) { + m_impl->id_bcs( node, component, DofType::Liquid) = IdBCs::NormalNode; + m_impl->flux_bcs(node, component, DofType::Liquid) = 0; +} +void BoundaryConditions::reset_flux_gas_dof( + uindex_t node, + uindex_t component + ) { + m_impl->id_bcs( node, component, DofType::Gas) = IdBCs::NormalNode; + m_impl->flux_bcs(node, component, DofType::Gas) = 0; +} + + +void BoundaryConditions::add_gas_node(uindex_t node) +{ + for (auto component: range(m_impl->m_nb_components)) { + m_impl->id_bcs(node, component, DofType::Liquid) = IdBCs::NoEquation; + } } void BoundaryConditions::add_gas_nodes(const std::vector& nodes) { for (auto node: nodes) add_gas_node(node); } +bool BoundaryConditions::is_gas_node(uindex_t node) const +{ + return (m_impl->id_bcs(node, 0, DofType::Liquid) == IdBCs::NoEquation); +} void BoundaryConditions::add_fixed_node(uindex_t node) { - m_impl->m_types[node] = NodeType::FixedNode; + + for (auto component: RangeIterator(m_impl->m_nb_components)) { + m_impl->id_bcs(node, component, DofType::Liquid) = IdBCs::FixedDof; + m_impl->id_bcs(node, component, DofType::Gas ) = IdBCs::FixedDof; + } } void BoundaryConditions::add_fixed_nodes(const std::vector& nodes) { for (auto node: nodes) add_fixed_node(node); } -bool BoundaryConditions::is_gas_node(uindex_t node) const { - return (m_impl->m_types[node] == NodeType::GasNode); +void BoundaryConditions::add_fixed_liquid_dof( + uindex_t node, + uindex_t component) { + m_impl->id_bcs(node, component, DofType::Liquid) = IdBCs::FixedDof; } -bool BoundaryConditions::is_fixed_node(uindex_t node) const { - return (m_impl->m_types[node] == NodeType::FixedNode); +void BoundaryConditions::add_fixed_gas_dof( + uindex_t node, + uindex_t component) { + m_impl->id_bcs(node, component, DofType::Gas) = IdBCs::FixedDof; } -bool BoundaryConditions::is_normal_node(uindex_t node) const { - return (m_impl->m_types[node] == NodeType::Normal); + +void BoundaryConditions::add_fixed_flux_liquid_dof( + uindex_t node, + uindex_t component) { + m_impl->id_bcs(node, component, DofType::Liquid) = IdBCs::FixedFlux; +} +void BoundaryConditions::add_fixed_flux_gas_dof( + uindex_t node, + uindex_t component) { + m_impl->id_bcs(node, component, DofType::Gas) = IdBCs::FixedFlux; +} +void BoundaryConditions::set_fixed_flux_liquid_dof( + uindex_t node, + uindex_t component, + scalar_t value + ) { + m_impl->id_bcs(node, component, DofType::Liquid) = IdBCs::FixedFlux; + m_impl->flux_bcs(node, component, DofType::Liquid) = value; +} +void BoundaryConditions::set_fixed_flux_gas_dof( + uindex_t node, + uindex_t component, + scalar_t value + ) { + m_impl->id_bcs(node, component, DofType::Gas) = IdBCs::FixedFlux; + m_impl->flux_bcs(node, component, DofType::Gas) = value; } +// chemistry + bool BoundaryConditions::has_constraint(const std::string &name) const { return m_impl->m_constraints.has_value(name); } const AdimensionalSystemConstraints& BoundaryConditions::get_constraint(uindex_t node) const { return m_impl->m_constraints[node]; } AdimensionalSystemConstraints& BoundaryConditions::get_constraint(uindex_t node) { return m_impl->m_constraints[node]; } const AdimensionalSystemConstraints& BoundaryConditions::get_constraint(const std::string& name) const { return m_impl->m_constraints.get(name); } AdimensionalSystemConstraints& BoundaryConditions::get_constraint(const std::string& name) { return m_impl->m_constraints.get(name); } AdimensionalSystemConstraints& BoundaryConditions::fork_constraint( uindex_t node, const std::string& name ) { return m_impl->m_constraints.fork(node, name); } AdimensionalSystemConstraints& BoundaryConditions::fork_constraint( const std::string& old_name, const std::string& new_name ) { m_impl->m_constraints.fork(old_name, new_name); return m_impl->m_constraints.get(new_name); } AdimensionalSystemConstraints& BoundaryConditions::add_constraint( const std::string& name ) { m_impl->m_constraints.emplace_back_cache(name); return m_impl->m_constraints.get(name); } void BoundaryConditions::set_constraint(uindex_t node, const std::string& name) { return m_impl->m_constraints.set(node, name); } -std::vector BoundaryConditions::get_fixed_nodes() { - std::vector fixed_nodes; - fixed_nodes.reserve(5); - for (auto it=m_impl->m_types.cbegin(); it!=m_impl->m_types.cend(); ++it) { - if(*it == NodeType::FixedNode) { - fixed_nodes.push_back(it - m_impl->m_types.cbegin()); - } - } - return fixed_nodes; -} - -std::vector BoundaryConditions::get_gas_nodes() { - std::vector gas_nodes; - gas_nodes.reserve(5); - for (auto it=m_impl->m_types.cbegin(); it!=m_impl->m_types.cend(); ++it) { - if(*it == NodeType::GasNode) { - gas_nodes.push_back(it - m_impl->m_types.cbegin()); - } - } - return gas_nodes; -} - } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/boundary_conditions.hpp b/src/reactmicp/systems/unsaturated/boundary_conditions.hpp index 51c11d4..b6cb243 100644 --- a/src/reactmicp/systems/unsaturated/boundary_conditions.hpp +++ b/src/reactmicp/systems/unsaturated/boundary_conditions.hpp @@ -1,147 +1,230 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_UNSATURATED_BOUNDARYCONDITIONS_HPP #define SPECMICP_REACTMICP_UNSATURATED_BOUNDARYCONDITIONS_HPP //! \file unsaturated/boundary_conditions.hpp //! \brief Boundary conditions for the unsaturated systems #include "specmicp_common/types.hpp" #include "specmicp_common/pimpl_ptr.hpp" #include #include namespace specmicp { struct AdimensionalSystemConstraints; namespace reactmicp { namespace systems { namespace unsaturated { +enum class IdBCs: index_t { + NoEquation = -1 + , NormalNode = 0 + , FixedDof = 1 + , FixedFlux = 2 +}; + + //! \brief Boundary conditions for the unsaturated system //! //! Contains both the transport and chemistry boundary conditions class SPECMICP_DLL_PUBLIC BoundaryConditions { public: //! \brief Return a pointer to a new set of boundary conditions - static std::shared_ptr make(uindex_t nb_nodes) { - return std::make_shared(nb_nodes); + static std::shared_ptr make( + uindex_t nb_nodes, + uindex_t nb_components + ) { + return std::make_shared( + nb_nodes, nb_components); } - //! \brief Build a BC container for the given number of node + //! \brief Build a BC container for the given number of nodes and components //! //! \sa make - BoundaryConditions(uindex_t nb_nodes); + BoundaryConditions(uindex_t nb_nodes, uindex_t nb_components); ~BoundaryConditions(); - //! \brief Set the given node to be a gas node (S_l=0) - void add_gas_node(uindex_t node); - //! \brief Set the given nodes to be a gas node (S_l=0) + //! \brief Return the equation type for a liquid dof + //! + //! \sa get_type_gas_dof + IdBCs get_bcs_liquid_dof(uindex_t node, uindex_t component); + //! \brief Return the equation type for a gas dof + //! + //! \sa get_type_liquid_dof + IdBCs get_bcs_gas_dof( uindex_t node, uindex_t component); + + //! \brief Return the fixed flux value for a liquid dof + //! + //! \return a value or 0 if the value is not set or it's not a fixed flux dof + //! + //! \sa get_flux_gas_dof + scalar_t get_flux_liquid_dof(uindex_t node, uindex_t component) const; + //! \brief Return the fixed flux value for a gas dof + //! + //! \return a value or 0 if the value is not set or it's not a fixed flux dof + //! + //! \sa get_flux_liquid_dof + scalar_t get_flux_gas_dof( uindex_t node, uindex_t component) const; + + //! \brief Remove boundary conditions for a liquid dof + //! + //! \sa reset_flux_gas_dof + void reset_flux_liquid_dof(uindex_t node, uindex_t component); + //! \brief Remove boundary conditions for a liquid dof + //! + //! \sa reset_flux_liquid_dof + void reset_flux_gas_dof( uindex_t node, uindex_t component); + + //! \brief Add a gas node + //! + //! Porosity is 1 and liquid saturation is 0 at this node. + //! This is just the boundary conditions for the liquid transport + //! equations, the user still need to set up the correct chemistry + //! constraints. + //! + //! + //! \sa add_gas_nodes + void add_gas_node( uindex_t node); + //! \brief Add gas nodes + //! + //! \sa add_gas_node void add_gas_nodes(const std::vector& nodes); - //! \brief Set the given node to be of fixed composition - void add_fixed_node(uindex_t node); - //! \brief Set the given node to be of fixed composition + //! \brief Return true if the given node is a gas node + bool is_gas_node( uindex_t node) const; + + //! \brief Add a fixed node + //! + //! All dof of tgis node will be fixed + void add_fixed_node( uindex_t node); + //! \brief Add fixed nodes + //! + //! All dof of these nodes will be fixed void add_fixed_nodes(const std::vector& nodes); - //! \brief Return true if the node is a gas node - bool is_gas_node(uindex_t node) const; - //! \brief Return true if the node has a fixed composition - bool is_fixed_node(uindex_t node) const; - //! \brief Return true if the node is neither a gas nor a fixed node - bool is_normal_node(uindex_t node) const; + //! \brief Set a liquid dof to be an essential BCs + //! + //! \sa add_fixed_gas_dof + void add_fixed_liquid_dof(uindex_t node, uindex_t component); + //! \brief Set a gas dof to be an essential BCs + //! + //! \sa add_fixed_liquid_dof + void add_fixed_gas_dof( uindex_t node, uindex_t component); + //! \brief Add a fixed flux conditions on a liquid dof + //! + //! \sa add_flux_gas_dof set_fixed_flux_liquid_dof + void add_fixed_flux_liquid_dof(uindex_t node, uindex_t component); + //! \brief Add a fixed flux conditions on a liquid dof + //! + //! \sa add_flux_gas_dof set_fixed_flux_liquid_dof + void add_fixed_flux_gas_dof( uindex_t node, uindex_t component); + //! \brief Add a fixed flux conditions on a liquid dof + //! + //! \sa add_flux_gas_dof set_fixed_flux_liquid_dof + void set_fixed_flux_liquid_dof( + uindex_t node, + uindex_t component, + scalar_t value + ); + //! \brief Add a fixed flux conditions on a liquid dof + //! + //! \sa add_flux_gas_dof set_fixed_flux_liquid_dof + void set_fixed_flux_gas_dof( + uindex_t node, + uindex_t component, + scalar_t value + ); + + //! \brief Return true if the constraints exist bool has_constraint(const std::string& name) const; //! \brief Return the chemical constraint of the given node const AdimensionalSystemConstraints& get_constraint(uindex_t node) const; //! \brief Return the chemical constraint of the given node //! //! \warning Changing the constraint changes it for all node sharing the //! same constraint. To modify only this node, use `fork_constraints' first. //! AdimensionalSystemConstraints& get_constraint(uindex_t node); //! \brief Return the constraint identified by 'name' const AdimensionalSystemConstraints& get_constraint( const std::string& name) const; //! \brief Return the constraint identified by 'name' AdimensionalSystemConstraints& get_constraint(const std::string& name); //! \brief Set a new constraints for the given node AdimensionalSystemConstraints& fork_constraint( uindex_t node, const std::string& name ); //! \brief Fork the given constraint AdimensionalSystemConstraints& fork_constraint( const std::string& old_name, const std::string& new_name ); //! \brief Add a constraint AdimensionalSystemConstraints& add_constraint( const std::string& name ); //! \brief Set the constraint of the given node void set_constraint(uindex_t node, const std::string& name); - //! \brief Return the vector of gas nodes - std::vector get_gas_nodes(); - //! \brief Return the vector of fixed nodes - std::vector get_fixed_nodes(); - private: struct SPECMICP_DLL_LOCAL BoundaryConditionsImpl; utils::pimpl_ptr m_impl; BoundaryConditions(const BoundaryConditions& other) = delete; BoundaryConditions& operator=(const BoundaryConditions& other) = delete; }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_UNSATURATED_BOUNDARYCONDITIONS_HPP diff --git a/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp b/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp index 8381137..e252e5c 100644 --- a/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp +++ b/src/reactmicp/systems/unsaturated/equilibrium_stagger.cpp @@ -1,530 +1,530 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "equilibrium_stagger.hpp" #include "variables.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "dfpm/meshes/mesh1d.hpp" #include "specmicp_database/data_container.hpp" #include "specmicp_common/compat.hpp" #include "specmicp/adimensional/adimensional_system_solution.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/config.h" namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { // ============================================ // // Declaration of implementation details // // ============================================ EquilibriumOptionsVector::EquilibriumOptionsVector(): utils::NameCachedVector() {} EquilibriumOptionsVector::EquilibriumOptionsVector( size_type size, std::string name, AdimensionalSystemSolverOptions& value ): utils::NameCachedVector( size, name, std::forward(value)) {} EquilibriumOptionsVector::EquilibriumOptionsVector( size_type size, std::string name ): utils::NameCachedVector(size, name) { static_assert(std::is_default_constructible::value, "AdimensionalSystemOptions must be default constructible."); } struct EquilibriumStagger::EquilibriumStaggerImpl { scalar_t m_dt {-1}; mesh::Mesh1DPtr m_mesh; database::RawDatabasePtr m_database; std::shared_ptr m_bcs; std::shared_ptr m_options; EquilibriumStaggerImpl( std::shared_ptr variables, std::shared_ptr boundary_conditions, std::shared_ptr options ): m_mesh(variables->get_mesh()), m_database(variables->get_database()), m_bcs(boundary_conditions), m_options(options) {} scalar_t dt() {return m_dt;} int compute_one_node(index_t node, UnsaturatedVariables * const vars); void compute_total_concentrations( index_t node, Vector& total_concentrations, UnsaturatedVariables * const vars ); void analyse_solution( index_t node, AdimensionalSystemSolution& solution, UnsaturatedVariables * const vars ); void initialize_timestep_one_node(index_t node, UnsaturatedVariables* vars); void initialize(UnsaturatedVariables* vars); void initialize_timestep(scalar_t dt, UnsaturatedVariables* vars); StaggerReturnCode restart_timestep(UnsaturatedVariables* vars); units::UnitsSet& get_units() { return m_options->get("default").units_set; } }; using TrueConstPtr = UnsaturatedVariables * const; inline TrueConstPtr cast_to_var(solver::VariablesBase * const var) { return static_cast(var); } // ============================================ // // Implementation // // ============================================ EquilibriumStagger::EquilibriumStagger( std::shared_ptr variables, std::shared_ptr boundary_conditions, std::shared_ptr options ): m_impl(utils::make_pimpl( variables, boundary_conditions, options)) { } EquilibriumStagger::~EquilibriumStagger() = default; //! \brief Initialize the stagger at the beginning of the computation //! //! \param var a shared_ptr to the variables void EquilibriumStagger::initialize(VariablesBase * const var) { TrueConstPtr true_vars = cast_to_var(var); m_impl->initialize(true_vars); } //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables void EquilibriumStagger::initialize_timestep( scalar_t dt, VariablesBase * const var ) { TrueConstPtr true_vars = cast_to_var(var); m_impl->initialize_timestep(dt, true_vars); } //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables EquilibriumStagger::StaggerReturnCode EquilibriumStagger::restart_timestep(VariablesBase * const var) { TrueConstPtr true_vars = cast_to_var(var); return m_impl->restart_timestep(true_vars); } int EquilibriumStagger::EquilibriumStaggerImpl::compute_one_node( index_t node, UnsaturatedVariables * const vars ) { AdimensionalSystemConstraints constraints = m_bcs->get_constraint(node); compute_total_concentrations(node, constraints.total_concentrations, vars); // set partial pressure model { user_model_saturation_f callable = vars->get_vapor_pressure_model(); water_partial_pressure_f pv_model = [callable, node](scalar_t sat){ return callable(node, sat); }; constraints.set_water_partial_pressure_model(pv_model); } // We use the current value of the saturation to avoid converging // to the previous solution // This is necessary because the update to the saturation is very small const AdimensionalSystemSolution& previous_solution = vars->get_adim_solution(node); Vector x; AdimensionalSystemSolver the_solver; if (likely(previous_solution.is_valid)) // should always be true { const auto porosity = vars->get_porosity()(node); const auto saturation = vars->get_liquid_saturation()(node); the_solver = AdimensionalSystemSolver( vars->get_database(), constraints, previous_solution, m_options->operator [](node) ); x = previous_solution.main_variables; // we use the new value of the saturation to avoid a change // to small to be detected by the speciaion solver x(0) = porosity*saturation; // x(0) = volume fraction water } else { const auto porosity = vars->get_porosity()(node); const auto saturation = vars->get_liquid_saturation()(node); the_solver = AdimensionalSystemSolver( vars->get_database(), constraints, m_options->operator [](node) ); the_solver.initialise_variables(x, porosity*saturation, -3); } // micpsolver::MiCPPerformance perf = the_solver.solve(x); if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) { ERROR << "Failed to solve equilibrium at node " << node << ". \n" << "Return code : " << (int) perf.return_code; return -1; } AdimensionalSystemSolution solution = the_solver.get_raw_solution(x); analyse_solution(node, solution, vars); vars->set_adim_solution(node, solution); return 0; } void EquilibriumStagger::EquilibriumStaggerImpl::compute_total_concentrations( index_t node, Vector& total_concentrations, UnsaturatedVariables * const vars ) { total_concentrations.resize(m_database->nb_component()); const scalar_t saturation = vars->get_liquid_saturation()(node); const scalar_t porosity = vars->get_porosity()(node); const scalar_t rt = vars->get_rt(); // water { const scalar_t ctilde_w = vars->get_water_aqueous_concentration()(node); const scalar_t cbar_w = vars->get_solid_concentration(0)(node); scalar_t c_w = saturation*porosity*ctilde_w + cbar_w; if (vars->component_has_gas(0)) { const scalar_t pv_w = vars->get_pressure_main_variables(0)(node); c_w += (1.0-saturation)*porosity*pv_w/rt; } total_concentrations(0) = c_w; } total_concentrations(1) = 0.0; // aqueous components for (index_t component: m_database->range_aqueous_component()) { const scalar_t ctilde_i = vars->get_aqueous_concentration(component)(node); const scalar_t cbar_i = vars->get_solid_concentration(component)(node); scalar_t c_i = saturation*porosity*ctilde_i + cbar_i; if (vars->component_has_gas(component)) { const scalar_t pv_i = vars->get_pressure_main_variables(component)(node); c_i += (1.0-saturation)*porosity*pv_i/rt; } total_concentrations(component) = c_i; } } void EquilibriumStagger::EquilibriumStaggerImpl::analyse_solution( index_t node, AdimensionalSystemSolution& solution, UnsaturatedVariables * const vars ) { AdimensionalSystemSolutionExtractor extr(solution, m_database, get_units()); const scalar_t saturation = extr.saturation_water(); const scalar_t sat_g = 1.0 - saturation; const scalar_t rho_l = extr.density_water(); const scalar_t porosity = extr.porosity(); const scalar_t rt = vars->get_rt(); // porosity SecondaryTransientVariable& porosity_vars = vars->get_porosity(); const scalar_t porosity_0 = porosity_vars.predictor(node); const scalar_t porosity_vel = (porosity - porosity_0)/m_dt; porosity_vars.variable(node) = porosity; porosity_vars.velocity(node) = porosity_vel; // water // saturation MainVariable& satvars = vars->get_liquid_saturation(); const scalar_t saturation_0 = satvars.predictor(node); const scalar_t sat_vel = (saturation - saturation_0)/m_dt; satvars.variable(node) = saturation; satvars.velocity(node) = sat_vel; { // total aqueous concentration SecondaryTransientVariable& cwtilde_vars = vars->get_water_aqueous_concentration(); const scalar_t cwtilde = rho_l*extr.total_aqueous_concentration(0); const scalar_t cwtilde_0 = cwtilde_vars.predictor(node); const scalar_t cwtilde_vel = (cwtilde - cwtilde_0)/m_dt; cwtilde_vars.variable(node) = cwtilde; cwtilde_vars.velocity(node) = cwtilde_vel; // total solid concentration MainVariable& solid_vars = vars->get_solid_concentration(0); const scalar_t cwbar = extr.total_solid_concentration(0); const scalar_t cwbar_0 = solid_vars.predictor(node); const scalar_t cwbar_vel = (cwbar - cwbar_0)/m_dt; solid_vars.variable(node) = cwbar; solid_vars.velocity(node) = cwbar_vel; solid_vars.chemistry_rate(node) = - cwbar_vel; // vapor pressure if (vars->component_has_gas(0)) { MainVariable& pres_vars = vars->get_pressure_main_variables(0); const scalar_t pv = vars->get_vapor_pressure_model()(node, saturation); const scalar_t pv_0 = pres_vars.predictor(node); const scalar_t pv_vel = (pv - pv_0)/m_dt; pres_vars.variable(node) = pv; pres_vars.velocity(node) = pv_vel; const scalar_t transient = ( (pv*porosity*sat_g) - (pv_0*porosity_0*(1.0-saturation_0)) )/ (rt*m_dt); const scalar_t pv_chem_rate = - transient + pres_vars.transport_fluxes(node); pres_vars.chemistry_rate(node) = pv_chem_rate; } } // aqueous components for (index_t component: m_database->range_aqueous_component()) { // total aqueous concentration MainVariable& aqconc = vars->get_aqueous_concentration(component); const scalar_t ctildei = rho_l*extr.total_aqueous_concentration(component); const scalar_t ctildei_0 = aqconc.predictor(node); const scalar_t ctildei_vel = (ctildei - ctildei_0) / m_dt; aqconc.variable(node) = ctildei; aqconc.velocity(node) = ctildei_vel; // total solid concentration MainVariable& solconc = vars->get_solid_concentration(component); const scalar_t cbari = extr.total_solid_concentration(component); const scalar_t cbari_0 = solconc.predictor(node); const scalar_t cbari_vel = (cbari - cbari_0) / m_dt; solconc.variable(node) = cbari; solconc.velocity(node) = cbari_vel; solconc.chemistry_rate(node) = - cbari_vel; // partial pressure if (vars->component_has_gas(component)) { MainVariable& pres_vars = vars->get_pressure_main_variables(component); index_t id_g = vars->get_id_gas(component); const scalar_t pi = extr.fugacity_gas(id_g)*vars->get_total_pressure(); const scalar_t pi_0 = pres_vars.predictor(node); const scalar_t pi_vel = (pi - pi_0) / m_dt; pres_vars.variable(node) = pi; pres_vars.velocity(node) = pi_vel; const scalar_t transient = ( (pi*porosity*sat_g) - (pi_0*porosity_0*(1.0-saturation_0)) ) / (rt * m_dt); const scalar_t pi_chem_rate = - transient + pres_vars.transport_fluxes(node); pres_vars.chemistry_rate(node) = pi_chem_rate; } } } void EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep_one_node( index_t node, UnsaturatedVariables * const vars ) { } EquilibriumStagger::StaggerReturnCode EquilibriumStagger::EquilibriumStaggerImpl::restart_timestep( UnsaturatedVariables* vars ) { int sum_retcode = 0; #ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum_retcode) #endif for (index_t node=0; nodenb_nodes(); ++node) { - if (not m_bcs->is_normal_node(node)) continue; + if (m_bcs->is_gas_node(node)) continue; sum_retcode += compute_one_node(node, vars); } if (sum_retcode != 0) { return StaggerReturnCode::UnknownError; } return StaggerReturnCode::ResidualMinimized; } void EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep( scalar_t dt, UnsaturatedVariables* vars ) { m_dt = dt; // initialize //water { MainVariable& cbar_w = vars->get_solid_concentration(0); cbar_w.predictor = cbar_w.variable; cbar_w.velocity.setZero(); cbar_w.chemistry_rate.setZero(); SecondaryTransientVariable& ctilde_w = vars->get_water_aqueous_concentration(); ctilde_w.predictor = ctilde_w.variable; ctilde_w.velocity.setZero(); if (vars->component_has_gas(0)) { MainVariable& pres_vars = vars->get_pressure_main_variables(0); pres_vars.predictor = pres_vars.variable; pres_vars.velocity.setZero(); pres_vars.chemistry_rate.setZero(); } } // aqueous component for (index_t component: m_database->range_aqueous_component()) { MainVariable& cbar_i = vars->get_solid_concentration(component); cbar_i.predictor = cbar_i.variable; cbar_i.velocity.setZero(); cbar_i.chemistry_rate.setZero(); if (vars->component_has_gas(component)) { MainVariable& pres_vars = vars->get_pressure_main_variables(component); pres_vars.predictor = pres_vars.variable; pres_vars.velocity.setZero(); pres_vars.chemistry_rate.setZero(); } } // porosity { SecondaryTransientVariable& porosity = vars->get_porosity(); porosity.predictor = porosity.variable; porosity.velocity.setZero(); } for (auto node: m_mesh->range_nodes()) { initialize_timestep_one_node(node, vars); } } void EquilibriumStagger::EquilibriumStaggerImpl::initialize( UnsaturatedVariables * const vars ) { // do nothing by default } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/fv_1dof_equation.hpp b/src/reactmicp/systems/unsaturated/fv_1dof_equation.hpp index f3d40a2..bdcdbb2 100644 --- a/src/reactmicp/systems/unsaturated/fv_1dof_equation.hpp +++ b/src/reactmicp/systems/unsaturated/fv_1dof_equation.hpp @@ -1,166 +1,177 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_SYSTEMS_UNSATURATED_FV1DDOFEQUATION_HPP #define SPECMICP_REACTMICP_SYSTEMS_UNSATURATED_FV1DDOFEQUATION_HPP //! \file unsaturated/fv_1dof_equation.hpp //! \brief 1D finite volume equation, with 1 dof per node #include "specmicp_common/types.hpp" #include "dfpm/solver/parabolic_program.hpp" #include "dfpm/meshes/mesh1d.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { +class BoundaryConditions; + //! \brief 1D finite volume equation, with 1DOF per node //! //! This is to be subclassed by the other equations template class FV1DOFEquation: public dfpmsolver::ParabolicProgram> { public: FV1DOFEquation(index_t nb_nodes): m_tot_ndf(nb_nodes) {} Derived* derived() {return static_cast(this);} //! \brief Return the number of equations index_t get_neq() const {return m_neq;} //! \brief Return the number of degrees of freedom per node index_t get_ndf() const {return 1;} //! \brief Return the total number of degrees of freedom index_t get_tot_ndf() const {return m_tot_ndf;} //! \brief Return the id of equation dof index_t id_equation(index_t node) { return derived()->id_equation_impl(node); } mesh::Mesh1D* get_mesh() { return derived()->get_mesh_impl(); } void pre_nodal_residual_hook(index_t node, const Vector& displacement) { return derived()->pre_nodal_residual_hook_impl(node, displacement); } void pre_residual_hook(const Vector& displacement) { return derived()->pre_residual_hook_impl(displacement); } void post_residual_hook(const Vector& displacement) { return derived()->post_residual_hook_impl(displacement); } void residuals_element( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ) { derived()->residuals_element_impl( element, displacement, velocity, element_residual, use_chemistry_rate); } void residuals_element( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual ) { residuals_element( element, displacement, velocity, element_residual, true); } //! \brief Compute the residuals void compute_residuals( const Vector& displacement, const Vector& velocity, Vector& residuals, bool use_chemistry_rate ); //! \brief Compute the residuals void compute_residuals( const Vector& displacement, const Vector& velocity, Vector& residuals ) { compute_residuals(displacement, velocity, residuals, true); } //! \brief Compute the jacobian void compute_jacobian(Vector& displacement, Vector& velocity, Eigen::SparseMatrix& jacobian, scalar_t alphadt ); //! \brief Update the solution void update_solution(const Vector& update, scalar_t lambda, scalar_t alpha_dt, Vector& predictor, Vector& displacement, Vector& velocity); + //! \brief Set the scaling factor of the equation void set_scaling(scalar_t scaling_factor) {m_scaling = scaling_factor;} + //! \brief Return the scaling factor scalar_t get_scaling() {return m_scaling;} + //! \brief Register the number of equations void register_number_equations(scalar_t neq) {m_neq = neq;} + //! \brief Number the equations + void number_equations(const BoundaryConditions * const bcs) { + derived()->number_equations(bcs); + } + + protected: scalar_t m_neq; scalar_t m_tot_ndf; scalar_t m_scaling {1.0}; }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #include "fv_1dof_equation.inl" #endif // SPECMICP_REACTMICP_SYSTEMS_UNSATURATED_FV1DDOFEQUATION_HPP diff --git a/src/reactmicp/systems/unsaturated/pressure_equation.cpp b/src/reactmicp/systems/unsaturated/pressure_equation.cpp index 355cdd6..4a38cd0 100644 --- a/src/reactmicp/systems/unsaturated/pressure_equation.cpp +++ b/src/reactmicp/systems/unsaturated/pressure_equation.cpp @@ -1,270 +1,323 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "pressure_equation.hpp" -#include "transport_constraints.hpp" +#include "boundary_conditions.hpp" #include "variables_box.hpp" #include "dfpm/solver/parabolic_driver.hpp" #include "dfpm/meshes/mesh1d.hpp" #include "specmicp_common/physics/constants.hpp" #include "specmicp_common/physics/maths.hpp" #include "specmicp_common/compat.hpp" #include namespace specmicp { namespace dfpmsolver { // explicit template instanciation template class ParabolicDriver; } //end namespace dfpmsolver namespace reactmicp { namespace systems { namespace unsaturated { struct SPECMICP_DLL_LOCAL PressureEquation::PressureEquationImpl { + uindex_t m_id_component; PressureVariableBox m_vars; std::vector m_ideq; mesh::Mesh1DPtr m_mesh; + std::shared_ptr m_bcs; + mesh::Mesh1D* mesh() {return m_mesh.get();} PressureVariableBox& vars() {return m_vars;} - bool node_has_equation(index_t node) { + bool node_has_equation(uindex_t node) { return m_ideq[node] != no_equation; } - index_t id_equation(index_t node) { + bool is_gas_node(uindex_t node) { + return m_bcs->is_gas_node(node); + } + + index_t id_equation(uindex_t node) { return m_ideq[node]; } - void fix_node(index_t node) { + void fix_node(uindex_t node) { m_ideq[node] = no_equation; } + scalar_t get_bc_gas_flux(uindex_t node) { + return m_bcs->get_flux_gas_dof(node, m_id_component); + } + PressureEquationImpl( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, - PressureVariableBox& the_vars + PressureVariableBox& the_vars, + std::shared_ptr bcs ): + m_id_component(id_component), m_vars(the_vars), m_ideq(the_mesh->nb_nodes(), -5), - m_mesh(the_mesh) + m_mesh(the_mesh), + m_bcs(bcs) {} + + scalar_t get_diff( + index_t node_0, + index_t node_1 + ); + + uindex_t number_equations(); + void compute_transport_rate(scalar_t dt, const Vector& displacement); }; PressureEquation::PressureEquation( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, PressureVariableBox& variables, - const TransportConstraints& constraints + std::shared_ptr bcs ): base(the_mesh->nb_nodes()), - m_impl(make_unique(the_mesh, variables)) + m_impl(make_unique( + id_component, the_mesh, variables, bcs)) { - number_equations(constraints); + number_equations(); } PressureEquation::~PressureEquation() = default; index_t PressureEquation::id_equation_impl(index_t id_dof) { return m_impl->id_equation(id_dof); } mesh::Mesh1D* PressureEquation::get_mesh_impl() { return m_impl->mesh(); } void PressureEquation::pre_nodal_residual_hook_impl( index_t node, const Vector& displacement) {} void PressureEquation::pre_residual_hook_impl(const Vector& displacement) { } void PressureEquation::post_residual_hook_impl(const Vector& displacement) { } +scalar_t PressureEquation::PressureEquationImpl::get_diff( + index_t node_0, + index_t node_1 + ) +{ + scalar_t coeff_diff = 0; + + if (is_gas_node(node_0)) { + coeff_diff = m_vars.resistance_gas_diffusivity(node_0) + * m_vars.relative_gas_diffusivity(node_0); + } else if (is_gas_node(node_1)) { + coeff_diff = m_vars.resistance_gas_diffusivity(node_1) + * m_vars.relative_gas_diffusivity(node_1); + } + else { + const scalar_t coeff_diff_0 = m_vars.resistance_gas_diffusivity(node_0) + * m_vars.relative_gas_diffusivity(node_0); + const scalar_t coeff_diff_1 = m_vars.resistance_gas_diffusivity(node_1) + * m_vars.relative_gas_diffusivity(node_1); + + coeff_diff = average(coeff_diff_0, coeff_diff_1); + } + coeff_diff *= m_vars.binary_diffusion_coefficient; + return coeff_diff; +} + //! \brief Compute the residuals inside 'element' void PressureEquation::residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ) { element_residual.setZero(); mesh::Mesh1D* m_mesh = m_impl->mesh(); PressureVariableBox& vars = m_impl->m_vars; const scalar_t& rt = vars.constants.rt; const scalar_t mass_coeff_0 = m_mesh->get_volume_cell_element(element, 0); const scalar_t mass_coeff_1 = m_mesh->get_volume_cell_element(element, 1); const index_t node_0 = m_mesh->get_node(element, 0); const index_t node_1 = m_mesh->get_node(element, 1); // Diffusion Cw - const scalar_t coeff_diff_0 = vars.resistance_gas_diffusivity(node_0) - * vars.relative_gas_diffusivity(node_0); - const scalar_t coeff_diff_1 = vars.resistance_gas_diffusivity(node_1) - * vars.relative_gas_diffusivity(node_1); - - const scalar_t coeff_diff = vars.binary_diffusion_coefficient * - average(coeff_diff_0, coeff_diff_1); - + const scalar_t coeff_diff = m_impl->get_diff(node_0, node_1); const scalar_t diff_flux = coeff_diff*( displacement(node_1) - displacement(node_0)) / m_mesh->get_dx(element) / rt; // Tot flux - const scalar_t tot_flux = m_mesh->get_face_area(element)*(diff_flux); + scalar_t flux_0 = diff_flux; + scalar_t flux_1 = -diff_flux; + flux_0 += m_impl->get_bc_gas_flux(node_0) / rt; + flux_1 += m_impl->get_bc_gas_flux(node_1) / rt; - const scalar_t flux_0 = tot_flux; - const scalar_t flux_1 = -tot_flux; + const scalar_t section = m_mesh->get_face_area(element); + flux_0 *= section; + flux_1 *= section; // transient if (m_impl->node_has_equation(node_0)) { const scalar_t porosity_0 = vars.porosity(node_0); const scalar_t pressure_0 = displacement(node_0); const scalar_t saturation_0 = 1.0 - vars.liquid_saturation(node_0); const scalar_t transient_0 = mass_coeff_0/rt * ( porosity_0 * saturation_0 * velocity(node_0) + saturation_0 * pressure_0 * vars.porosity.velocity(node_0) - porosity_0 * pressure_0 * vars.liquid_saturation.velocity(node_0) ); auto res = transient_0 - flux_0; if (use_chemistry_rate) { res += mass_coeff_0*vars.partial_pressure.chemistry_rate(node_0); } element_residual(0) = res/get_scaling(); } if (m_impl->node_has_equation(node_1)) { const scalar_t porosity_1 = vars.porosity(node_1); const scalar_t pressure_1 = displacement(node_1); const scalar_t saturation_1 = 1.0 - vars.liquid_saturation(node_1); const scalar_t transient_1 = mass_coeff_1/rt * ( porosity_1 * saturation_1 * velocity(node_1) + saturation_1 * pressure_1 * vars.porosity.velocity(node_1) - porosity_1 * pressure_1 * vars.liquid_saturation.velocity(node_1) ); auto res = transient_1 - flux_1; if (use_chemistry_rate) { res += mass_coeff_1*vars.partial_pressure.chemistry_rate(node_1); } element_residual(1) = res/get_scaling(); } } -void PressureEquation::number_equations(const TransportConstraints& constraints) +uindex_t PressureEquation::PressureEquationImpl::number_equations() { - for (int fixed: constraints.fixed_nodes()) - { - m_impl->fix_node(fixed); - } - index_t neq = 0; - for (index_t node: m_impl->mesh()->range_nodes()) + uindex_t neq = 0; + for (index_t node: m_mesh->range_nodes()) { - if (m_impl->m_ideq[node] == -5) - { - m_impl->m_ideq[node] = neq; + switch (m_bcs->get_bcs_gas_dof(node, m_id_component)){ + case IdBCs::FixedDof: + fix_node(node); + break; + default: + m_ideq[node] = neq; ++neq; } } + return neq; +} + +void PressureEquation::number_equations() +{ + auto neq = m_impl->number_equations(); register_number_equations(neq); } void PressureEquation::compute_transport_rate( scalar_t dt, const Vector& displacement) { m_impl->compute_transport_rate(dt, displacement); } void PressureEquation::PressureEquationImpl::compute_transport_rate( scalar_t dt, const Vector& displacement) { const scalar_t& rt = m_vars.constants.rt; const MainVariable& saturation = m_vars.liquid_saturation; MainVariable& pressure = m_vars.partial_pressure; const SecondaryTransientVariable& porosity = m_vars.porosity; for (index_t node: m_mesh->range_nodes()) { if (! node_has_equation(node)) continue; const scalar_t transient = ( (porosity(node)*(1.0-saturation(node))*displacement(node)) - (porosity.predictor(node)*(1.0-saturation.predictor(node))*pressure.predictor(node)) ) / (rt*dt); const scalar_t chem_rates = ( - pressure.chemistry_rate(node) ); pressure.transport_fluxes(node) = transient - chem_rates; } } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/pressure_equation.hpp b/src/reactmicp/systems/unsaturated/pressure_equation.hpp index da9ebb2..71876d6 100644 --- a/src/reactmicp/systems/unsaturated/pressure_equation.hpp +++ b/src/reactmicp/systems/unsaturated/pressure_equation.hpp @@ -1,110 +1,112 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_UNSATURATED_PRESSUREEQUATION_HPP #define SPECMICP_REACTMICP_UNSATURATED_PRESSUREEQUATION_HPP //! \file unsaturated/pressure_equation.hpp //! \brief The pressure diffusion equation #include "specmicp_common/types.hpp" #include "fv_1dof_equation.hpp" #include "dfpm/meshes/mesh1dfwd.hpp" #include "dfpm/solver/parabolic_program.hpp" #include "variables.hpp" #include namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { -class TransportConstraints; +class BoundaryConditions; //! \brief Pressure diffusion equation - valid for water and aqueous components //! /*! \f[ \frac{1}{R T} \frac{\partial \phi S_g p_i}{\partial t} = \nabla \cdot \left(\frac{D_l}{R T} \nabla p_i \right) - \dot{R}^{g \rightarrow l}_i\f] */ class SPECMICP_DLL_LOCAL PressureEquation: public FV1DOFEquation { using base = FV1DOFEquation; using base::get_scaling; using base::register_number_equations; public: PressureEquation( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, PressureVariableBox& variables, - const TransportConstraints& constraints); + std::shared_ptr bcs + ); ~PressureEquation(); index_t id_equation_impl(index_t id_dof); mesh::Mesh1D* get_mesh_impl(); void pre_nodal_residual_hook_impl(index_t node, const Vector& displacement); void pre_residual_hook_impl(const Vector& displacement); void post_residual_hook_impl(const Vector& displacement); //! \brief Compute the residuals inside 'element' void residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ); void compute_transport_rate(scalar_t dt, const Vector& displacement); private: - void number_equations(const TransportConstraints& constraints); + void number_equations(); struct PressureEquationImpl; std::unique_ptr m_impl; }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_UNSATURATED_PRESSUREEQUATION_HPP diff --git a/src/reactmicp/systems/unsaturated/saturation_equation.cpp b/src/reactmicp/systems/unsaturated/saturation_equation.cpp index 326ca27..0f4921d 100644 --- a/src/reactmicp/systems/unsaturated/saturation_equation.cpp +++ b/src/reactmicp/systems/unsaturated/saturation_equation.cpp @@ -1,357 +1,446 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "saturation_equation.hpp" #include "variables_box.hpp" -#include "transport_constraints.hpp" +#include "boundary_conditions.hpp" #include "dfpm/meshes/mesh1d.hpp" #include "dfpm/solver/parabolic_driver.hpp" #include "specmicp_common/compat.hpp" #include "specmicp_common/physics/maths.hpp" #include #include namespace specmicp { namespace dfpmsolver { // explicit template instanciation template class dfpmsolver::ParabolicDriver; } //end namespace dfpmsolver } //end namespace specmicp namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { -static constexpr index_t no_equation {-1}; -static constexpr index_t no_eq_no_var {-2}; +static constexpr index_t no_equation {-1}; +static constexpr index_t no_eq_no_var {-2}; static constexpr index_t not_initialized {-5}; struct SPECMICP_DLL_LOCAL SaturationEquation::SaturationEquationImpl { mesh::Mesh1DPtr m_mesh; SaturationVariableBox m_vars; std::vector m_ideq; bool m_store_residual_info; + std::shared_ptr m_bcs; SaturationEquationImpl( mesh::Mesh1DPtr the_mesh, - SaturationVariableBox vars): + SaturationVariableBox vars, + std::shared_ptr bcs + ): m_mesh(the_mesh), m_vars(vars), - m_ideq(the_mesh->nb_nodes(), not_initialized) + m_ideq(the_mesh->nb_nodes(), not_initialized), + m_bcs(bcs) {} index_t& id_equation(index_t node) {return m_ideq[node];} bool node_can_flux(index_t node) {return m_ideq[node] != no_eq_no_var;} bool node_has_eq(index_t node) {return m_ideq[node] > no_equation;} void set_store_residual_info() { m_store_residual_info = true; } void reset_store_residual_info() { m_store_residual_info = false; } bool store_residual_info() { return m_store_residual_info; } //! \brief Return a pointer to the mesh mesh::Mesh1D* mesh() {return m_mesh.get();} range_t range_nodes() {return m_mesh->range_nodes();} void set_relative_variables(const Vector& displacement); void set_relative_variables(index_t node, const Vector& displacement); void compute_transport_rate(scalar_t dt, const Vector& displacement); -}; + uindex_t number_equations(); + + scalar_t get_perm( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const; + scalar_t get_diff( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const; + + + scalar_t get_bc_liquid_flux(uindex_t node) { + return m_bcs->get_flux_liquid_dof(node, 0); + } +}; SaturationEquation::~SaturationEquation() = default; -SaturationEquation::SaturationEquation(mesh::Mesh1DPtr the_mesh, - SaturationVariableBox &variables, - const TransportConstraints& constraints +SaturationEquation::SaturationEquation( + uindex_t id_component, + mesh::Mesh1DPtr the_mesh, + SaturationVariableBox& variables, + std::shared_ptr bcs ): base(the_mesh->nb_nodes()), - m_impl(make_unique(the_mesh, variables)) + m_impl(make_unique(the_mesh, variables, bcs)) { - number_equations(constraints); + number_equations(); } -void SaturationEquation::number_equations(const TransportConstraints& constraints) +uindex_t SaturationEquation::SaturationEquationImpl::number_equations() { - - for (index_t node: constraints.fixed_nodes()) - { - m_impl->id_equation(node) = no_equation; - } - for (index_t node: constraints.gas_nodes()) - { - m_impl->id_equation(node) = no_eq_no_var; - } - index_t neq = 0; - for (index_t node: m_impl->range_nodes()) { - if (m_impl->id_equation(node) == not_initialized) - { - m_impl->id_equation(node) = neq; + uindex_t neq = 0; + + for (index_t node: m_mesh->range_nodes()) { + IdBCs bc_l = m_bcs->get_bcs_liquid_dof(node, 0); + + switch (bc_l) { + case IdBCs::FixedDof: + m_ideq[node] = no_equation; + break; + case IdBCs::NoEquation: + m_ideq[node] = no_eq_no_var; + break; + default: + m_ideq[node] = neq; ++neq; } } + return neq; + +} + +void SaturationEquation::number_equations() +{ + + auto neq = m_impl->number_equations(); register_number_equations(neq); } mesh::Mesh1D* SaturationEquation::get_mesh_impl() { return m_impl->mesh(); } index_t SaturationEquation::id_equation_impl(index_t id_dof) { return m_impl->id_equation(id_dof)<0?no_equation:m_impl->id_equation(id_dof); } void SaturationEquation::pre_nodal_residual_hook_impl( index_t node, const Vector& displacement) { m_impl->set_relative_variables(node, displacement); } void SaturationEquation::pre_residual_hook_impl(const Vector& displacement) { m_impl->set_relative_variables(displacement); m_impl->set_store_residual_info(); } void SaturationEquation::post_residual_hook_impl(const Vector& displacement) { m_impl->reset_store_residual_info(); } +scalar_t SaturationEquation::SaturationEquationImpl::get_perm( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const +{ + scalar_t perm; + if (advection_direction == +1) { + perm = m_vars.liquid_permeability(node_0) + * m_vars.relative_liquid_permeability(node_0); + } + else if (advection_direction == -1) { + perm = m_vars.liquid_permeability(node_1) + * m_vars.relative_liquid_permeability(node_1); + } + else { + const scalar_t perm_0 = m_vars.liquid_permeability(node_0) + * m_vars.relative_liquid_permeability(node_0); + const scalar_t perm_1 = m_vars.liquid_permeability(node_1) + * m_vars.relative_liquid_permeability(node_1); + perm = average(perm_0, perm_1); + } + return perm; +} + +scalar_t SaturationEquation::SaturationEquationImpl::get_diff( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const +{ + scalar_t diff; + if (advection_direction == +1) { + diff = m_vars.liquid_diffusivity(node_0) + * m_vars.relative_liquid_diffusivity(node_0); + } + else if (advection_direction == -1) { + diff = m_vars.liquid_diffusivity(node_1) + * m_vars.relative_liquid_diffusivity(node_1); + } + else { + const scalar_t coeff_diff_0 = m_vars.liquid_diffusivity(node_0) + * m_vars.relative_liquid_diffusivity(node_0); + const scalar_t coeff_diff_1 = m_vars.liquid_diffusivity(node_1) + * m_vars.relative_liquid_diffusivity(node_1); + + diff = average( + coeff_diff_0, coeff_diff_1); + } + return diff; +} + // Residuals // ========= void SaturationEquation::residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ) { element_residual.setZero(); mesh::Mesh1D* m_mesh = m_impl->mesh(); SaturationVariableBox& vars = m_impl->m_vars; const scalar_t mass_coeff_0 = m_mesh->get_volume_cell_element(element, 0); const scalar_t mass_coeff_1 = m_mesh->get_volume_cell_element(element, 1); const index_t node_0 = m_mesh->get_node(element, 0); const index_t node_1 = m_mesh->get_node(element, 1); + const scalar_t section = m_mesh->get_face_area(element); scalar_t flux_0 = 0.0; scalar_t flux_1 = 0.0; + index_t advection_direction = 0; if (m_impl->node_can_flux(node_0) and m_impl->node_can_flux(node_1)) { - // Cap pressure gradient - const scalar_t perm_0 = vars.liquid_permeability(node_0) - * vars.relative_liquid_permeability(node_0); - const scalar_t perm_1 = vars.liquid_permeability(node_1) - * vars.relative_liquid_permeability(node_1); - const scalar_t perm = average(perm_0, perm_1); - - const scalar_t aq_coefficient = average( - vars.aqueous_concentration(node_0), - vars.aqueous_concentration(node_1) - ); - const scalar_t cap_pressure_gradient = ( vars.capillary_pressure(node_1) - - vars.capillary_pressure(node_0) - ) / m_mesh->get_dx(element); - const scalar_t advection_flux = (perm/vars.constants.viscosity_liquid_water)*cap_pressure_gradient; - const scalar_t cappres_flux = -aq_coefficient*advection_flux; + // advection -> Pc(S) + scalar_t pc0 = vars.capillary_pressure(node_0); + scalar_t pc1 = vars.capillary_pressure(node_1); + scalar_t aq_coefficient = 0; + scalar_t advection_flux = 0; + scalar_t cappres_flux = 0; + + if (pc1 > pc0) { + advection_direction = +1; + aq_coefficient = vars.aqueous_concentration(node_0); + } + else if (pc0 > pc1) { + advection_direction = -1; + aq_coefficient = vars.aqueous_concentration(node_1); + } - // Diffusion Cw - const scalar_t coeff_diff_0 = vars.liquid_diffusivity(node_0) - * vars.relative_liquid_diffusivity(node_0); - const scalar_t coeff_diff_1 = vars.liquid_diffusivity(node_1) - * vars.relative_liquid_diffusivity(node_1); + if (pc0 != pc1) { + const scalar_t perm = m_impl->get_perm(node_0, node_1, advection_direction); + const scalar_t cap_pressure_gradient = ( + pc1 - pc0 ) / m_mesh->get_dx(element); - const scalar_t coeff_diff = average( - coeff_diff_0, coeff_diff_1); + advection_flux = (perm/vars.constants.viscosity_liquid_water)*cap_pressure_gradient; + cappres_flux = -aq_coefficient*advection_flux; + } + // Diffusion Cw + const scalar_t coeff_diff = m_impl->get_diff(node_0, node_1, advection_direction); const scalar_t aq_flux = coeff_diff*(vars.aqueous_concentration(node_1) - vars.aqueous_concentration(node_0) ) / m_mesh->get_dx(element); // Tot flux - const scalar_t tot_flux = m_mesh->get_face_area(element)*(cappres_flux + aq_flux); + const scalar_t tot_flux = section*(cappres_flux + aq_flux); flux_0 = tot_flux; flux_1 = -tot_flux; + flux_0 += m_impl->get_bc_liquid_flux(node_0) * section; + flux_1 += m_impl->get_bc_liquid_flux(node_1) * section; + // Storage if (m_impl->store_residual_info()) { // advective flux stored by element vars.advection_flux(element) = advection_flux; // fluxes to compute exchange term vars.liquid_saturation.transport_fluxes(node_0) += flux_0; vars.liquid_saturation.transport_fluxes(node_1) += flux_1; } } // transient if (m_impl->node_has_eq(node_0)) { const scalar_t porosity_0 = vars.porosity(node_0); const scalar_t aq_tot_conc_0 = vars.aqueous_concentration(node_0); const scalar_t saturation_0 = displacement(node_0); scalar_t transient_0 = ( porosity_0 * aq_tot_conc_0 * velocity(node_0) + saturation_0 * aq_tot_conc_0 * vars.porosity.velocity(node_0) + porosity_0 * saturation_0 * vars.aqueous_concentration.velocity(node_0) ); auto res = mass_coeff_0*transient_0 - flux_0; if (use_chemistry_rate) { const scalar_t chemistry_0 = vars.liquid_saturation.chemistry_rate(node_0) + vars.solid_concentration.chemistry_rate(node_0) + vars.vapor_pressure.chemistry_rate(node_0) ; res -= mass_coeff_0*chemistry_0; } element_residual(0) = res / get_scaling(); } if (m_impl->node_has_eq(node_1)) { const scalar_t porosity_1 = vars.porosity(node_1); const scalar_t aq_tot_conc_1 = vars.aqueous_concentration(node_1); const scalar_t saturation_1 = displacement(node_1); scalar_t transient_1 = ( porosity_1 * aq_tot_conc_1 * velocity(node_1) + saturation_1 * aq_tot_conc_1 * vars.porosity.velocity(node_1) + porosity_1 * saturation_1 * vars.aqueous_concentration.velocity(node_1) ); auto res = mass_coeff_1*transient_1 - flux_1; if (use_chemistry_rate) { const scalar_t chemistry_1 = vars.liquid_saturation.chemistry_rate(node_1) + vars.solid_concentration.chemistry_rate(node_1) + vars.vapor_pressure.chemistry_rate(node_1) ; res -= mass_coeff_1*chemistry_1; } element_residual(1) = res / get_scaling(); } } void SaturationEquation::set_relative_variables(const Vector& displacement) { return m_impl->set_relative_variables(displacement); } void SaturationEquation::SaturationEquationImpl::set_relative_variables( index_t node, const Vector& displacement ) { if (not node_can_flux(node)) return; const scalar_t saturation = displacement(node); m_vars.relative_liquid_diffusivity(node) = m_vars.relative_liquid_diffusivity_f(node, saturation); m_vars.relative_liquid_permeability(node) = m_vars.relative_liquid_permeability_f(node, saturation); m_vars.capillary_pressure(node) = m_vars.capillary_pressure_f(node, saturation); } void SaturationEquation::SaturationEquationImpl::set_relative_variables(const Vector& displacement) { for (index_t node: m_mesh->range_nodes()) { set_relative_variables(node, displacement); } } void SaturationEquation::compute_transport_rate(scalar_t dt, const Vector& displacement) { m_impl->compute_transport_rate(dt, displacement); } void SaturationEquation::SaturationEquationImpl::compute_transport_rate( scalar_t dt, const Vector& displacement) { MainVariable& saturation = m_vars.liquid_saturation; const MainVariable& solid_conc = m_vars.solid_concentration; const MainVariable& pressure = m_vars.vapor_pressure; const SecondaryTransientVariable& porosity = m_vars.porosity; const SecondaryTransientVariable& aqueous_concentration = m_vars.aqueous_concentration; for (index_t node: m_mesh->range_nodes()) { if (! node_has_eq(node)) continue; const scalar_t transient = ( (porosity(node)*aqueous_concentration(node)*displacement(node)) - (porosity.predictor(node)*aqueous_concentration.predictor(node)*saturation.predictor(node)) ) / dt; const scalar_t chem_rates = ( saturation.chemistry_rate(node) + solid_conc.chemistry_rate(node) + pressure.chemistry_rate(node) ); saturation.transport_fluxes(node) = transient - chem_rates; } } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/saturation_equation.hpp b/src/reactmicp/systems/unsaturated/saturation_equation.hpp index 272bd56..bf6a417 100644 --- a/src/reactmicp/systems/unsaturated/saturation_equation.hpp +++ b/src/reactmicp/systems/unsaturated/saturation_equation.hpp @@ -1,110 +1,113 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_UNSATURATED_SATURATIONEQUATION_HPP #define SPECMICP_REACTMICP_UNSATURATED_SATURATIONEQUATION_HPP //! \file unsaturated/saturation_equation.hpp //! \brief The saturation equation for the unsaturated system #include "specmicp_common/types.hpp" #include "fv_1dof_equation.hpp" #include "dfpm/meshes/mesh1dfwd.hpp" #include "variables.hpp" +#include namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { -class TransportConstraints; +class BoundaryConditions; //! \brief The saturation equation //! //! Solve the transport of liquid water class SPECMICP_DLL_LOCAL SaturationEquation: public FV1DOFEquation { using base = FV1DOFEquation; using base::get_scaling; using base::register_number_equations; public: SaturationEquation( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, SaturationVariableBox& variables, - const TransportConstraints& constraints); + std::shared_ptr bcs + ); ~SaturationEquation(); //! \brief Return the id of equation dof index_t id_equation_impl(index_t id_dof); mesh::Mesh1D* get_mesh_impl(); void pre_nodal_residual_hook_impl(index_t node, const Vector& displacement); void pre_residual_hook_impl(const Vector& displacement); void post_residual_hook_impl(const Vector& displacement); //! \brief Compute the residuals inside 'element' for 'component' void residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ); //! \brief Set the relative variables void set_relative_variables(const Vector& displacement); void compute_transport_rate(scalar_t dt, const Vector& displacement); -private: + void number_equations(); - void number_equations(const TransportConstraints& constraints); +private: struct SaturationEquationImpl; std::unique_ptr m_impl; }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_UNSATURATED_SATURATIONEQUATION_HPP diff --git a/src/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp b/src/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp index 93defee..7098c8b 100644 --- a/src/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp +++ b/src/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp @@ -1,392 +1,539 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "saturation_pressure_equation.hpp" #include "variables_box.hpp" -#include "transport_constraints.hpp" +#include "boundary_conditions.hpp" #include "dfpm/meshes/mesh1d.hpp" #include "dfpm/solver/parabolic_driver.hpp" #include "specmicp_common/compat.hpp" #include "specmicp_common/physics/maths.hpp" #include #include namespace specmicp { namespace dfpmsolver { // explicit template instanciation template class ParabolicDriver; } //end namespace dfpmsolver } //end namespace specmicp namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { -static constexpr index_t no_equation {-1}; -static constexpr index_t no_eq_no_var {-2}; +static constexpr index_t no_equation {-1}; +static constexpr index_t no_eq_no_var {-2}; +static constexpr index_t no_eq_flux_gas {-3}; static constexpr index_t not_initialized {-5}; -struct SPECMICP_DLL_LOCAL SaturationPressureEquation::SaturationPressureEquationImpl +struct SPECMICP_DLL_LOCAL + SaturationPressureEquation::SaturationPressureEquationImpl { mesh::Mesh1DPtr m_mesh; SaturationPressureVariableBox m_vars; std::vector m_ideq; + + std::shared_ptr m_bcs; + bool m_store_residual_info; SaturationPressureEquationImpl( mesh::Mesh1DPtr the_mesh, - SaturationPressureVariableBox vars): + SaturationPressureVariableBox vars, + std::shared_ptr bcs + ): m_mesh(the_mesh), m_vars(vars), - m_ideq(the_mesh->nb_nodes(), not_initialized) + m_ideq(the_mesh->nb_nodes(), not_initialized), + m_bcs(bcs) {} + uindex_t number_equations(); + index_t& id_equation(index_t node) {return m_ideq[node];} - bool node_can_flux(index_t node) {return m_ideq[node] != no_eq_no_var;} - bool node_has_eq(index_t node) {return m_ideq[node] > no_equation;} + bool node_can_flux(index_t node) {return m_ideq[node] != no_eq_no_var;} + bool node_has_eq(index_t node) {return m_ideq[node] > no_equation;} void set_store_residual_info() { m_store_residual_info = true; } void reset_store_residual_info() { m_store_residual_info = false; } bool store_residual_info() { return m_store_residual_info; } //! \brief Return a pointer to the mesh mesh::Mesh1D* mesh() {return m_mesh.get();} range_t range_nodes() {return m_mesh->range_nodes();} void set_relative_variables(const Vector& displacement); void set_relative_variables(index_t node, const Vector& displacement); void compute_transport_rate(scalar_t dt, const Vector& displacement); + + + // diffusion coefficient and all + scalar_t get_perm( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const; + + scalar_t get_diff( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const; + + scalar_t get_diff_gas( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const; + + scalar_t get_bc_liquid_flux(uindex_t node) { + return m_bcs->get_flux_liquid_dof(node, 0); + } + scalar_t get_bc_gas_flux(uindex_t node) { + return m_bcs->get_flux_gas_dof(node, 0); + } }; SaturationPressureEquation::~SaturationPressureEquation() = default; -SaturationPressureEquation::SaturationPressureEquation(mesh::Mesh1DPtr the_mesh, - SaturationPressureVariableBox &variables, - const TransportConstraints& constraints +SaturationPressureEquation::SaturationPressureEquation( + uindex_t id_component, + mesh::Mesh1DPtr the_mesh, + SaturationPressureVariableBox& variables, + std::shared_ptr bcs ): base(the_mesh->nb_nodes()), - m_impl(make_unique(the_mesh, variables)) + m_impl(make_unique(the_mesh, variables, bcs)) { - number_equations(constraints); + number_equations(); } -void SaturationPressureEquation::number_equations(const TransportConstraints& constraints) -{ - for (index_t node: constraints.fixed_nodes()) - { - m_impl->id_equation(node) = no_equation; - } - for (index_t node: constraints.gas_nodes()) - { - m_impl->id_equation(node) = no_eq_no_var; - } - index_t neq = 0; - for (index_t node: m_impl->range_nodes()) { - if (m_impl->id_equation(node) == not_initialized) - { - m_impl->id_equation(node) = neq; +uindex_t +SaturationPressureEquation::SaturationPressureEquationImpl::number_equations() +{ + uindex_t neq = 0; + for (index_t node: m_mesh->range_nodes()) { + IdBCs bc_l = m_bcs->get_bcs_liquid_dof(node, 0); + + switch (bc_l) { + case IdBCs::FixedDof: + m_ideq[node] = no_equation; + break; + case IdBCs::NoEquation: + m_ideq[node] = no_eq_no_var; + break; + default: + m_ideq[node] = neq; ++neq; } } + return neq; +} + +void SaturationPressureEquation::number_equations() +{ + auto neq = m_impl->number_equations(); register_number_equations(neq); } mesh::Mesh1D* SaturationPressureEquation::get_mesh_impl() { return m_impl->mesh(); } index_t SaturationPressureEquation::id_equation_impl(index_t id_dof) { return m_impl->id_equation(id_dof)<0?no_equation:m_impl->id_equation(id_dof); } void SaturationPressureEquation::pre_nodal_residual_hook_impl( index_t node, const Vector& displacement ) { m_impl->set_relative_variables(node, displacement); } void SaturationPressureEquation::pre_residual_hook_impl( const Vector& displacement ) { m_impl->set_relative_variables(displacement); m_impl->set_store_residual_info(); } void SaturationPressureEquation::post_residual_hook_impl( const Vector& displacement ) { m_impl->reset_store_residual_info(); } // Residuals // ========= +scalar_t SaturationPressureEquation::SaturationPressureEquationImpl::get_perm( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const +{ + scalar_t perm = 0.0; + if (advection_direction == +1) { + perm = m_vars.liquid_permeability(node_0) + * m_vars.relative_liquid_permeability(node_0); + } + else if (advection_direction == -1) { + perm = m_vars.liquid_permeability(node_1) + * m_vars.relative_liquid_permeability(node_1); + } + else { + const scalar_t perm_0 = m_vars.liquid_permeability(node_0) + * m_vars.relative_liquid_permeability(node_0); + const scalar_t perm_1 = m_vars.liquid_permeability(node_1) + * m_vars.relative_liquid_permeability(node_1); + perm = average(perm_0, perm_1); + } + return perm; +} + +scalar_t SaturationPressureEquation::SaturationPressureEquationImpl::get_diff( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const +{ + scalar_t diff = 0.0; + if (advection_direction == +1) { + diff = m_vars.liquid_diffusivity(node_0) + * m_vars.relative_liquid_diffusivity(node_0); + } + else if (advection_direction == -1) { + diff = m_vars.liquid_diffusivity(node_1) + * m_vars.relative_liquid_diffusivity(node_1); + } + else { + const scalar_t coeff_diff_0 = m_vars.liquid_diffusivity(node_0) + * m_vars.relative_liquid_diffusivity(node_0); + const scalar_t coeff_diff_1 = m_vars.liquid_diffusivity(node_1) + * m_vars.relative_liquid_diffusivity(node_1); + + diff = average(coeff_diff_0, coeff_diff_1); + } + return diff; +} + +scalar_t SaturationPressureEquation::SaturationPressureEquationImpl::get_diff_gas( + index_t node_0, + index_t node_1, + index_t advection_direction + ) const +{ + scalar_t diff = 0.0; + if (advection_direction == +1) { + diff = m_vars.resistance_gas_diffusivity(node_0) + * m_vars.relative_gas_diffusivity(node_0); + } + else if (advection_direction == -1) { + diff = m_vars.resistance_gas_diffusivity(node_1) + * m_vars.relative_gas_diffusivity(node_1); + } + else { + if (m_bcs->is_gas_node(node_0)) { + diff = m_vars.resistance_gas_diffusivity(node_0) + * m_vars.relative_gas_diffusivity(node_0); + } + else if (m_bcs->is_gas_node(node_1)) { + diff = m_vars.resistance_gas_diffusivity(node_1) + * m_vars.relative_gas_diffusivity(node_1); + } + else { + const scalar_t coeff_diff_0 = m_vars.resistance_gas_diffusivity(node_0) + * m_vars.relative_gas_diffusivity(node_0); + const scalar_t coeff_diff_1 = m_vars.resistance_gas_diffusivity(node_1) + * m_vars.relative_gas_diffusivity(node_1); + + diff = average(coeff_diff_0, coeff_diff_1); + } + } + return diff * m_vars.binary_diffusion_coefficient; +} + void SaturationPressureEquation::residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ) { element_residual.setZero(); mesh::Mesh1D* m_mesh = m_impl->mesh(); SaturationPressureVariableBox& vars = m_impl->m_vars; const scalar_t& rt = vars.constants.rt; const scalar_t mass_coeff_0 = m_mesh->get_volume_cell_element(element, 0); const scalar_t mass_coeff_1 = m_mesh->get_volume_cell_element(element, 1); + const scalar_t section = m_mesh->get_face_area(element); + const index_t node_0 = m_mesh->get_node(element, 0); const index_t node_1 = m_mesh->get_node(element, 1); scalar_t flux_0 = 0.0; scalar_t flux_1 = 0.0; scalar_t tot_flux = 0.0; scalar_t advection_flux = 0.0; + index_t advection_direction = 0; + if (m_impl->node_can_flux(node_0) and m_impl->node_can_flux(node_1)) { - // Cap pressure gradient - const scalar_t perm_0 = vars.liquid_permeability(node_0) - * vars.relative_liquid_permeability(node_0); - const scalar_t perm_1 = vars.liquid_permeability(node_1) - * vars.relative_liquid_permeability(node_1); - const scalar_t perm = average(perm_0, perm_1); + // advection -> Pc(S) + scalar_t pc0 = vars.capillary_pressure(node_0); + scalar_t pc1 = vars.capillary_pressure(node_1); + + if (pc1 > pc0) advection_direction = +1; + else if (pc0 > pc1) advection_direction = -1; + + const scalar_t perm = m_impl->get_perm( + node_0, node_1, advection_direction); const scalar_t aq_coefficient = average( vars.aqueous_concentration(node_0), vars.aqueous_concentration(node_1) ); - const scalar_t cap_pressure_gradient = ( vars.capillary_pressure(node_1) - - vars.capillary_pressure(node_0) + const scalar_t cap_pressure_gradient = ( pc1 - pc0 ) / m_mesh->get_dx(element); - advection_flux = (perm/vars.constants.viscosity_liquid_water)*cap_pressure_gradient; + advection_flux = ( perm / vars.constants.viscosity_liquid_water + ) * cap_pressure_gradient; const scalar_t cappres_flux = -aq_coefficient*advection_flux; - // Diffusion Cw - const scalar_t coeff_diff_0 = vars.liquid_diffusivity(node_0) - * vars.relative_liquid_diffusivity(node_0); - const scalar_t coeff_diff_1 = vars.liquid_diffusivity(node_1) - * vars.relative_liquid_diffusivity(node_1); + // diffusion Cw + const scalar_t coeff_diff = m_impl->get_diff(node_0, node_1, advection_direction); - const scalar_t coeff_diff = average( - coeff_diff_0, coeff_diff_1); - - const scalar_t aq_flux = coeff_diff*(vars.aqueous_concentration(node_1) - - vars.aqueous_concentration(node_0) - ) / m_mesh->get_dx(element); + const scalar_t aq_flux = coeff_diff*( + vars.aqueous_concentration(node_1) + - vars.aqueous_concentration(node_0) + ) / m_mesh->get_dx(element); tot_flux = (cappres_flux + aq_flux); } // Diffusion pressure pv - const scalar_t coeff_diff_gas_0 = vars.resistance_gas_diffusivity(node_0) - * vars.relative_gas_diffusivity(node_0); - const scalar_t coeff_diff_gas_1 = vars.resistance_gas_diffusivity(node_1) - * vars.relative_gas_diffusivity(node_1); - - const scalar_t coeff_diff_gas = vars.binary_diffusion_coefficient * - average(coeff_diff_gas_0, coeff_diff_gas_1); + const scalar_t coeff_diff_gas = m_impl->get_diff_gas( + node_0, node_1, advection_direction); const scalar_t diff_flux_gas = coeff_diff_gas*( - vars.partial_pressure(node_1) - vars.partial_pressure(node_0)) + vars.partial_pressure(node_1) + - vars.partial_pressure(node_0) + ) / m_mesh->get_dx(element) / rt; // Tot flux tot_flux += diff_flux_gas; tot_flux *= m_mesh->get_face_area(element); - flux_0 = tot_flux; - flux_1 = -tot_flux; + flux_0 = tot_flux; + flux_1 = -tot_flux; + + // liquid and gas fluxes (0 if no fluxes) + flux_0 += m_impl->get_bc_liquid_flux(node_0); + flux_0 += m_impl->get_bc_gas_flux( node_0) / rt; + flux_1 += m_impl->get_bc_liquid_flux(node_1); + flux_1 += m_impl->get_bc_gas_flux( node_1) / rt; + + //if (m_impl->m_bcs->is_gas_node(node_0)) { + // flux_1 += 2.2*( + // vars.partial_pressure(node_1) + // - vars.partial_pressure(node_0) + // )/rt; + //} + + flux_0 *= section; + flux_1 *= section; // Storage if (m_impl->store_residual_info()) { // advective flux stored by element vars.advection_flux(element) = advection_flux; // fluxes to compute exchange term vars.liquid_saturation.transport_fluxes(node_0) += flux_0; vars.liquid_saturation.transport_fluxes(node_1) += flux_1; } // transient if (m_impl->node_has_eq(node_0)) { const scalar_t porosity_0 = vars.porosity(node_0); const scalar_t aq_tot_conc_0 = vars.aqueous_concentration(node_0); const scalar_t saturation_0 = displacement(node_0); scalar_t transient_0 = ( porosity_0 * aq_tot_conc_0 * velocity(node_0) + saturation_0 * aq_tot_conc_0 * vars.porosity.velocity(node_0) + porosity_0 * saturation_0 * vars.aqueous_concentration.velocity(node_0) ); auto res = mass_coeff_0*transient_0 - flux_0; if (use_chemistry_rate) { const scalar_t chemistry_0 = vars.liquid_saturation.chemistry_rate(node_0) + vars.solid_concentration.chemistry_rate(node_0) + vars.partial_pressure.chemistry_rate(node_0) ; res -= mass_coeff_0*chemistry_0; } element_residual(0) = res / get_scaling(); } - if (m_impl->node_has_eq(node_1)) { const scalar_t porosity_1 = vars.porosity(node_1); const scalar_t aq_tot_conc_1 = vars.aqueous_concentration(node_1); const scalar_t saturation_1 = displacement(node_1); scalar_t transient_1 = ( porosity_1 * aq_tot_conc_1 * velocity(node_1) + saturation_1 * aq_tot_conc_1 * vars.porosity.velocity(node_1) + porosity_1 * saturation_1 * vars.aqueous_concentration.velocity(node_1) ); auto res = mass_coeff_1*transient_1 - flux_1; if (use_chemistry_rate) { const scalar_t chemistry_1 = vars.liquid_saturation.chemistry_rate(node_1) + vars.solid_concentration.chemistry_rate(node_1) + vars.partial_pressure.chemistry_rate(node_1) ; res -= mass_coeff_1*chemistry_1; } element_residual(1) = res / get_scaling(); } } void SaturationPressureEquation::set_relative_variables(const Vector& displacement) { return m_impl->set_relative_variables(displacement); } void SaturationPressureEquation::SaturationPressureEquationImpl::set_relative_variables( index_t node, const Vector& displacement ) { if (not node_can_flux(node)) return; const scalar_t saturation = displacement(node); m_vars.relative_liquid_diffusivity(node) = m_vars.relative_liquid_diffusivity_f(node, saturation); m_vars.relative_liquid_permeability(node) = m_vars.relative_liquid_permeability_f(node, saturation); m_vars.capillary_pressure(node) = m_vars.capillary_pressure_f(node, saturation); //m_vars.partial_pressure(node) = m_vars.partial_pressure_f(node, saturation); } void SaturationPressureEquation::SaturationPressureEquationImpl::set_relative_variables( const Vector& displacement ) { for (index_t node: m_mesh->range_nodes()) { set_relative_variables(node, displacement); } } void SaturationPressureEquation::compute_transport_rate( scalar_t dt, const Vector& displacement ) { m_impl->compute_transport_rate(dt, displacement); } void SaturationPressureEquation::SaturationPressureEquationImpl::compute_transport_rate( scalar_t dt, const Vector& displacement) { MainVariable& saturation = m_vars.liquid_saturation; const MainVariable& solid_conc = m_vars.solid_concentration; const MainVariable& pressure = m_vars.partial_pressure; const SecondaryTransientVariable& porosity = m_vars.porosity; const SecondaryTransientVariable& aqueous_concentration = m_vars.aqueous_concentration; for (index_t node: m_mesh->range_nodes()) { if (! node_has_eq(node)) continue; const scalar_t transient = ( (porosity(node)*aqueous_concentration(node)*displacement(node)) - (porosity.predictor(node)*aqueous_concentration.predictor(node)*saturation.predictor(node)) ) / dt; const scalar_t chem_rates = ( saturation.chemistry_rate(node) + solid_conc.chemistry_rate(node) + pressure.chemistry_rate(node) ); saturation.transport_fluxes(node) = transient - chem_rates; } } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/src/reactmicp/systems/unsaturated/saturation_pressure_equation.hpp b/src/reactmicp/systems/unsaturated/saturation_pressure_equation.hpp index 3637491..867340c 100644 --- a/src/reactmicp/systems/unsaturated/saturation_pressure_equation.hpp +++ b/src/reactmicp/systems/unsaturated/saturation_pressure_equation.hpp @@ -1,110 +1,115 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_UNSATURATED_SATURATIONPRESSUREEQUATION_HPP #define SPECMICP_REACTMICP_UNSATURATED_SATURATIONPRESSUREEQUATION_HPP //! \file unsaturated/saturation_pressure_equation.hpp //! \brief The saturation and pressure equation for the unsaturated system #include "specmicp_common/types.hpp" #include "fv_1dof_equation.hpp" #include "dfpm/meshes/mesh1dfwd.hpp" #include "variables.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace unsaturated { -class TransportConstraints; +class SaturationPressureVariableBox; +class BoundaryConditions; //! \brief The saturation equation //! //! Solve the transport of liquid water class SPECMICP_DLL_LOCAL SaturationPressureEquation: public FV1DOFEquation { using base = FV1DOFEquation; using base::get_scaling; using base::register_number_equations; public: SaturationPressureEquation( + uindex_t id_component, mesh::Mesh1DPtr the_mesh, SaturationPressureVariableBox& variables, - const TransportConstraints& constraints); + std::shared_ptr bcs + ); ~SaturationPressureEquation(); //! \brief Return the id of equation dof index_t id_equation_impl(index_t id_dof); mesh::Mesh1D* get_mesh_impl(); void pre_nodal_residual_hook_impl(index_t node, const Vector& displacement); void pre_residual_hook_impl(const Vector& displacement); void post_residual_hook_impl(const Vector& displacement); //! \brief Compute the residuals inside 'element' for 'component' void residuals_element_impl( index_t element, const Vector& displacement, const Vector& velocity, Eigen::Vector2d& element_residual, bool use_chemistry_rate ); //! \brief Set the relative variables void set_relative_variables(const Vector& displacement); + //! \brief Compute the transport rates void compute_transport_rate(scalar_t dt, const Vector& displacement); -private: + //! \brief Number the equations + void number_equations(); - void number_equations(const TransportConstraints& constraints); +private: struct SaturationPressureEquationImpl; std::unique_ptr m_impl; }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_UNSATURATED_SATURATIONPRESSUREEQUATION_HPP diff --git a/src/reactmicp/systems/unsaturated/transport_constraints.hpp b/src/reactmicp/systems/unsaturated/transport_constraints.hpp deleted file mode 100644 index 786dfc8..0000000 --- a/src/reactmicp/systems/unsaturated/transport_constraints.hpp +++ /dev/null @@ -1,98 +0,0 @@ -/* ============================================================================= - - Copyright (c) 2014 - 2016 - F. Georget Princeton University - All rights reserved. - - -Redistribution and use in source and binary forms, with or without modification, -are permitted provided that the following conditions are met: - - 1. Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimer. - - 2. Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - 3. Neither the name of the copyright holder nor the names of its - contributors may be used to endorse or promote products derived from this - software without specific prior written permission. - -* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR -ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON -ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * - -============================================================================= */ - -#ifndef SPECMICP_REACTMICP_TRANSPORTCONSTRAINTS_HPP -#define SPECMICP_REACTMICP_TRANSPORTCONSTRAINTS_HPP - -//! \file unsaturated/transport_constraints.hpp -//! \brief Constraints for the transport equations - -#include "specmicp_common/types.hpp" -#include -#include - -namespace specmicp { -namespace reactmicp { -namespace systems { -namespace unsaturated { - -//! \brief Constraints for the transport equations -class SPECMICP_DLL_PUBLIC TransportConstraints -{ -public: - TransportConstraints() {} - - //! \brief Set the fixed nodes - void set_fixed_nodes(const std::vector& fixed_nodes) { - m_fixed_nodes = fixed_nodes; - } - - //! \brief Set the gas nodes - void set_gas_nodes(const std::vector& gas_nodes) { - m_gas_nodes = gas_nodes; - } - - //! \brief Add a fixed node - void add_fixed_node(index_t node) { - m_fixed_nodes.push_back(node); - } - - //! \brief Add a gas node - //! - //! A gas node is a node with S_l = 0 - void add_gas_node(index_t node) { - m_gas_nodes.push_back(node); - } - - //! \brief Return the list of fixed nodes - const std::vector& fixed_nodes() const { - return m_fixed_nodes; - } - - //! \brief Return the list of gas nodes - const std::vector& gas_nodes() const { - return m_gas_nodes; - } - -private: - std::vector m_fixed_nodes; - std::vector m_gas_nodes; -}; - -} //end namespace unsaturated -} //end namespace systems -} //end namespace reactmicp -} //end namespace specmicp - -#endif // SPECMICP_REACTMICP_TRANSPORTCONSTRAINTS_HPP diff --git a/src/reactmicp/systems/unsaturated/transport_stagger.cpp b/src/reactmicp/systems/unsaturated/transport_stagger.cpp index 9039d66..eac647e 100644 --- a/src/reactmicp/systems/unsaturated/transport_stagger.cpp +++ b/src/reactmicp/systems/unsaturated/transport_stagger.cpp @@ -1,1096 +1,1092 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "transport_stagger.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "variables.hpp" #include "variables_box.hpp" -#include "transport_constraints.hpp" +#include "boundary_conditions.hpp" #include "saturation_equation.hpp" #include "saturation_pressure_equation.hpp" #include "pressure_equation.hpp" #include "aqueous_equation.hpp" #include "aqueous_pressure_equation.hpp" #include "specmicp_database/database_holder.hpp" #include "dfpm/solver/parabolic_driver.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/config.h" +#include #include namespace specmicp { namespace dfpmsolver { extern template class ParabolicDriver; extern template class ParabolicDriver; extern template class ParabolicDriver; extern template class ParabolicDriver; extern template class ParabolicDriver; } //end namespace dfpmsolver namespace reactmicp { namespace systems { namespace unsaturated { //! \brief Type of the equation //! \internal enum class EquationType { Saturation, LiquidAqueous, Pressure }; //! \brief Format type to a string (for message error) //! \internal static std::string to_string(EquationType eq_type); //! \brief Format DFPM solver retcode to a string static std::string to_string(dfpmsolver::ParabolicDriverReturnCode retcode); // forward declaration of wrapper classes // They are defined at the bootom of this file class TaskBase; class SaturationTask; using VariablesBase = solver::VariablesBase; // =================================== // // // // Declaration of the implementation // // // // =================================== // //! \brief Implementation class for the Unsaturated transport stagger //! \internal //! //! This class does all the work //! It was designed to be OpenMP compatible class UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl { public: // Constructor // ----------- UnsaturatedTransportStaggerImpl( UnsaturatedVariablesPtr variables, std::shared_ptr boundary_conditions, bool merge_saturation_pressure, bool merge_aqueous_pressure ); // Main functions // -------------- //! \brief Return the residual scalar_t get_residual(UnsaturatedVariables * const vars); //! \brief Return the first residual (start of timestep) scalar_t get_residual_0(UnsaturatedVariables * const vars); //! \brief Return the update (norm of velocity vector) scalar_t get_update(UnsaturatedVariables * const vars); //! \brief Initialize timestep void initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars); //! \brief Restart the timestep solver::StaggerReturnCode restart_timestep(UnsaturatedVariables* vars); // Timestep management // ------------------- //! \brief Save the timestep void save_dt(scalar_t dt) {m_dt = dt;} //! \brief Return the timestep scalar_t get_dt() {return m_dt;} // Task index // ---------- index_t& id_aqueous_task(index_t component) { return m_aq_equation_task[static_cast(component)]; } index_t id_aqueous_task(index_t component) const { return m_aq_equation_task[static_cast(component)]; } index_t& id_gas_task(index_t component) { return m_aq_equation_task[static_cast(component)]; } - index_t id_gas_tasj(index_t component) const { + index_t id_gas_task(index_t component) const { return m_aq_equation_task[static_cast(component)]; } //! \brief Return the options of the saturation equation dfpmsolver::ParabolicDriverOptions* get_saturation_options(); //! \brief Return the options of aqueous equation of component dfpmsolver::ParabolicDriverOptions* get_aqueous_options(index_t component); //! \brief Return the options of gas equation of component dfpmsolver::ParabolicDriverOptions* get_gas_options(index_t component); private: scalar_t m_norm_0 {1.0}; // timestep management scalar_t m_dt {-1.0}; //! \brief The saturation equations std::unique_ptr m_saturation_equation; // Equation and solver std::vector> m_equation_list; std::vector m_aq_equation_task; std::vector m_gas_equation_task; std::vector m_merged_gas; }; // Main functions // =============== // call of the implementation UnsaturatedTransportStagger::UnsaturatedTransportStagger( std::shared_ptr variables, std::shared_ptr boundary_conditions, bool merge_saturation_pressure, bool merge_aqueous_pressure ): m_impl(utils::make_pimpl( variables, boundary_conditions, merge_saturation_pressure, merge_aqueous_pressure )) { } UnsaturatedTransportStagger::~UnsaturatedTransportStagger() = default; static UnsaturatedVariables * const get_var(VariablesBase * const var) { return static_cast(var); } void UnsaturatedTransportStagger::initialize_timestep( scalar_t dt, VariablesBase* var ) { m_impl->initialize_timestep(dt, get_var(var)); } solver::StaggerReturnCode UnsaturatedTransportStagger::restart_timestep(VariablesBase * const var) { return m_impl->restart_timestep(get_var(var)); } scalar_t UnsaturatedTransportStagger::get_residual(VariablesBase * const var) { return m_impl->get_residual(get_var(var)); } scalar_t UnsaturatedTransportStagger::get_residual_0(VariablesBase * const var) { return m_impl->get_residual_0(get_var(var)); } scalar_t UnsaturatedTransportStagger::get_update(VariablesBase * const var) { return m_impl->get_update(get_var(var)); } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::get_saturation_options() { return m_impl->get_saturation_options(); } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::get_aqueous_options(index_t component) { return m_impl->get_aqueous_options(component); } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::get_gas_options(index_t component) { return m_impl->get_gas_options(component); } // =============================== // // =============================== // // // // Implementation details // // ---------------------- // // // // =============================== // // =============================== // // 2 main sections // - Solver wrappers : wrapper around 1 couple equation/solver // - UnsaturatedTransportStaggerImpl : call the wrappers // =============================== // // // // Solver wrappers // // // // =============================== // // This wrappers are used to abstract the residual computation // and the methods to solve every equations //! \brief Base class for a task //! //! A task is how we solve governing equations, //! and obtain residuals and update //! //! \internal class SPECMICP_DLL_LOCAL TaskBase { public: TaskBase(index_t component, EquationType eq_type): m_type(eq_type), m_component(component) {} virtual ~TaskBase() {} //! \brief Compute the squared residuals virtual scalar_t compute_squared_residual( UnsaturatedVariables * const vars ) = 0; //! \brief Compute the squared residuals at the beginning of the timestep virtual scalar_t compute_squared_residual_0( UnsaturatedVariables * const vars ) = 0; //! \brief Compute the squared update of the variables virtual scalar_t compute_squared_update( UnsaturatedVariables * const vars ) = 0; //! \brief Initialize the timestep virtual void initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars ) = 0; //! \brief Solve the governing equation virtual dfpmsolver::ParabolicDriverReturnCode restart_timestep( UnsaturatedVariables * const vars ) = 0; //! \brief Return the component index (in the db) index_t component() {return m_component;} //! \brief The equation type EquationType equation_type() {return m_type;} private: EquationType m_type; index_t m_component; }; //! \brief Base class for a equation solver //! //! \internal template class SPECMICP_DLL_LOCAL EquationTask: public TaskBase { public: using EqT = typename traits::EqT; using SolverT = typename traits::SolverT; EquationTask( index_t component, mesh::Mesh1DPtr the_mesh, typename traits::VariableBoxT& variables, - const TransportConstraints& constraints + std::shared_ptr bcs ): TaskBase(component, traits::equation_type), - m_equation(the_mesh, variables, constraints), + m_equation(component, the_mesh, variables, bcs), m_solver(m_equation) {} EquationTask( index_t component, mesh::Mesh1DPtr the_mesh, typename traits::VariableBoxT& variables, - const TransportConstraints& constraints, + std::shared_ptr bcs, scalar_t scaling ): TaskBase(component, traits::equation_type), - m_equation(the_mesh, variables, constraints), + m_equation(component, the_mesh, variables, bcs), m_solver(m_equation) { m_equation.set_scaling(scaling); } Derived* derived() {return static_cast(this);} // This function must be implemented by subclass MainVariable& get_var(UnsaturatedVariables * const vars) { return derived()->get_var(vars); } scalar_t compute_squared_residual( UnsaturatedVariables * const vars ) override { Vector residuals; const MainVariable& main = get_var(vars); m_equation.compute_residuals(main.variable, main.velocity, residuals, true); return residuals.squaredNorm()/m_norm_square_0; } scalar_t compute_squared_residual_0( UnsaturatedVariables * const vars ) override { Vector residuals; const MainVariable& main = get_var(vars); m_equation.compute_residuals(main.variable, main.velocity, residuals, false); m_norm_square_0 = residuals.squaredNorm(); if (m_norm_square_0 < 1e-8) m_norm_square_0 = 1.0; return 1.0; } scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { const MainVariable& main = get_var(vars); return main.velocity.squaredNorm(); } void initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars ) override { m_dt = dt; MainVariable& main = get_var(vars); main.predictor = main.variable; main.transport_fluxes.setZero(); main.velocity.setZero(); m_solver.initialize_timestep(dt, get_var(vars).variable); } dfpmsolver::ParabolicDriverReturnCode restart_timestep( UnsaturatedVariables * const vars ) override { MainVariable& main = get_var(vars); m_solver.velocity() = main.velocity; auto retcode = m_solver.restart_timestep(main.variable); if (retcode > dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet) { main.velocity = m_solver.velocity(); } return retcode; } dfpmsolver::ParabolicDriverOptions& get_options() { return m_solver.get_options(); } const dfpmsolver::ParabolicDriverOptions& get_options() const { return m_solver.get_options(); } scalar_t get_dt() {return m_dt;} protected: scalar_t m_norm_square_0 {-1}; scalar_t m_dt {-1}; EqT m_equation; SolverT m_solver; }; //! \brief Traits struct for the SaturationTask class //! //! \internal struct SPECMICP_DLL_LOCAL SaturationTaskTraits { using EqT = SaturationEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = SaturationVariableBox; static constexpr EquationType equation_type {EquationType::Saturation}; }; //! \brief Wrapper for the saturation equation solver //! //! \internal class SPECMICP_DLL_LOCAL SaturationTask: public EquationTask { using base = EquationTask; public: SaturationTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints + std::shared_ptr bcs ): EquationTask( - component, the_mesh, variables, constraints) + component, the_mesh, variables, bcs) {} SaturationTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints, + std::shared_ptr bcs, scalar_t scaling ): EquationTask( - component, the_mesh, variables, constraints, scaling) + component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_liquid_saturation(); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the SaturationPressureTask class //! //! \internal struct SPECMICP_DLL_LOCAL SaturationPressureTaskTraits { using EqT = SaturationPressureEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = SaturationPressureVariableBox; static constexpr EquationType equation_type {EquationType::Saturation}; }; //! \brief Wrapper for the saturation equation solver //! //! \internal class SPECMICP_DLL_LOCAL SaturationPressureTask: public EquationTask { using base = EquationTask; public: SaturationPressureTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationPressureTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints + std::shared_ptr bcs ): EquationTask( - component, the_mesh, variables, constraints) + component, the_mesh, variables, bcs) {} SaturationPressureTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationPressureTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints, + std::shared_ptr bcs, scalar_t scaling ): EquationTask( - component, the_mesh, variables, constraints, scaling) + component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_liquid_saturation(); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the LiquidAqueousTask class //! //! \internal struct SPECMICP_DLL_LOCAL LiquidAqueousTaskTraits { using EqT = AqueousTransportEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = LiquidAqueousComponentVariableBox; static constexpr EquationType equation_type {EquationType::LiquidAqueous}; }; //! \brief Wrapper for the liquid transport of aqueous component equation //! //! \internal class SPECMICP_DLL_LOCAL LiquidAqueousTask: public EquationTask { using base = EquationTask; public: LiquidAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidAqueousTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints + std::shared_ptr bcs ): EquationTask( - component, the_mesh, variables, constraints) + component, the_mesh, variables, bcs) {} LiquidAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidAqueousTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints, + std::shared_ptr bcs, scalar_t scaling ): EquationTask( - component, the_mesh, variables, constraints, scaling) + component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_aqueous_concentration(component()); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the LiquidGasAqueousTask class //! //! \internal struct SPECMICP_DLL_LOCAL LiquidGasAqueousTaskTraits { using EqT = AqueousGasTransportEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = LiquidGasAqueousVariableBox; static constexpr EquationType equation_type {EquationType::LiquidAqueous}; }; //! \brief Wrapper for the liquid transport of aqueous component equation //! //! \internal class SPECMICP_DLL_LOCAL LiquidGasAqueousTask: public EquationTask { using base = EquationTask; public: LiquidGasAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidGasAqueousTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints + std::shared_ptr bcs ): - base(component, the_mesh, variables, constraints) + base(component, the_mesh, variables, bcs) {} LiquidGasAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidGasAqueousTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints, + std::shared_ptr bcs, scalar_t scaling ): - base(component, the_mesh, variables, constraints, scaling) + base(component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_aqueous_concentration(component()); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the Pressure Task traits //! //! \internal struct SPECMICP_DLL_LOCAL PressureTaskTraits { using EqT = PressureEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = PressureVariableBox; static constexpr EquationType equation_type {EquationType::Pressure}; }; //! \brief Wrapper for the pressure equation solver //! //! \internal class SPECMICP_DLL_LOCAL PressureTask: public EquationTask { using base = EquationTask; public: PressureTask( index_t component, mesh::Mesh1DPtr the_mesh, PressureTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints + std::shared_ptr bcs ): EquationTask( - component, the_mesh, variables, constraints) + component, the_mesh, variables, bcs) {} PressureTask( index_t component, mesh::Mesh1DPtr the_mesh, PressureTaskTraits::VariableBoxT&& variables, - const TransportConstraints& constraints, + std::shared_ptr bcs, scalar_t scaling ): EquationTask( - component, the_mesh, variables, constraints, scaling) + component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_pressure_main_variables(component()); } }; // ================================== // // // // UnsaturatedTransportStaggerImpl // // // // ================================== // // constructor // =========== UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::UnsaturatedTransportStaggerImpl( UnsaturatedVariablesPtr variables, - std::shared_ptr boundary_conditions, + std::shared_ptr bcs, bool merge_saturation_pressure, bool merge_aqueous_pressure ): m_merged_gas(variables->get_database()->nb_component(), false) { database::RawDatabasePtr raw_db = variables->get_database(); mesh::Mesh1DPtr the_mesh = variables->get_mesh(); m_aq_equation_task = std::vector( raw_db->nb_component(), no_equation); m_gas_equation_task = std::vector( raw_db->nb_component(), no_equation); - - TransportConstraints constraints; - constraints.set_fixed_nodes(boundary_conditions->get_fixed_nodes()); - constraints.set_gas_nodes(boundary_conditions->get_gas_nodes()); - dfpmsolver::ParabolicDriverOptions* opts = nullptr; // Saturation equations // ==================== // // There is 2 choices for the saturation equation // Either SaturationPressureEquation or SaturationPressure if (merge_saturation_pressure and variables->component_has_gas(0)) { m_saturation_equation = make_unique( 0, the_mesh, variables->get_saturation_pressure_variables(), - constraints, + bcs, variables->get_aqueous_scaling(0) ); m_merged_gas[0] = true; opts = &(static_cast( m_saturation_equation.get())->get_options()); } else { m_saturation_equation = make_unique( 0, the_mesh, variables->get_saturation_variables(), - constraints, + bcs, variables->get_aqueous_scaling(0) ); opts = &(static_cast( m_saturation_equation.get())->get_options()); } opts->residuals_tolerance = 1e-8; opts->absolute_tolerance = 1e-12; opts->step_tolerance = 1e-12; opts->threshold_stationary_point = 1e-12; opts->maximum_iterations = 500; opts->sparse_solver = sparse_solvers::SparseSolver::SparseLU; opts->quasi_newton = 3; opts->maximum_step_length = 0.1; opts->max_iterations_at_max_length = 5000; const index_t size = raw_db->nb_aqueous_components() + variables->nb_gas(); m_equation_list.reserve(size); // Liquid aqueous diffusion-advection // ================================== for (index_t id: raw_db->range_aqueous_component()) { if (merge_aqueous_pressure and variables->component_has_gas(id)) { m_equation_list.emplace_back( make_unique( id, the_mesh, variables->get_liquid_gas_aqueous_variables(id), - constraints, + bcs, variables->get_aqueous_scaling(id)) ); m_merged_gas[id] = true; } else { m_equation_list.emplace_back( make_unique( id, the_mesh, variables->get_liquid_aqueous_component_variables(id), - constraints, + bcs, variables->get_aqueous_scaling(id)) ); } id_aqueous_task(id) = m_equation_list.size()-1; dfpmsolver::ParabolicDriverOptions& opts = *get_aqueous_options(id); opts.maximum_iterations = 500; opts.residuals_tolerance = 1e-8; opts.absolute_tolerance = 1e-12; opts.step_tolerance = 1e-12; opts.threshold_stationary_point = 1e-12; opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU; opts.maximum_step_length= 1e6; } // Water partial pressure // ====================== // // Depending on the saturation equation chosen we may or may not include // this equations if (variables->component_has_gas(0) and (not m_merged_gas[0])) { m_equation_list.emplace_back( make_unique(0, the_mesh, variables->get_pressure_variables(0), - constraints, + bcs, variables->get_gaseous_scaling(0) )); id_gas_task(0) = m_equation_list.size()-1; dfpmsolver::ParabolicDriverOptions& opts = *get_gas_options(0); opts.maximum_iterations = 500; opts.residuals_tolerance = 1e-8; opts.absolute_tolerance = 1e-12; opts.step_tolerance = 1e-12; opts.threshold_stationary_point = 1e-12; opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU; opts.maximum_step_length= 1e6; } // Gaseous diffusion for (index_t id: raw_db->range_aqueous_component()) { if (variables->component_has_gas(id) and (not m_merged_gas[id])) { m_equation_list.emplace_back( make_unique(id, the_mesh, variables->get_pressure_variables(id), - constraints, + bcs, variables->get_gaseous_scaling(id) )); id_gas_task(id) = m_equation_list.size()-1; dfpmsolver::ParabolicDriverOptions& opts = *get_gas_options(id); opts.maximum_iterations = 500; opts.residuals_tolerance = 1e-8; opts.absolute_tolerance = 1e-12; opts.step_tolerance = 1e-12; opts.threshold_stationary_point = 1e-12; opts.sparse_solver = sparse_solvers::SparseSolver::SparseLU; opts.maximum_step_length= 1e6; } } } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_saturation_options() { dfpmsolver::ParabolicDriverOptions* opts = nullptr; if (m_merged_gas[0]) { opts = &static_cast( m_saturation_equation.get())->get_options(); } else { opts = &static_cast( m_saturation_equation.get())->get_options(); } return opts; } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_aqueous_options( index_t component) { dfpmsolver::ParabolicDriverOptions* opts = nullptr; auto id = id_aqueous_task(component); if (id != no_equation) { if (m_merged_gas[component]) { opts = &static_cast( m_equation_list[id].get())->get_options(); } else { opts = &static_cast( m_equation_list[id].get())->get_options(); } } return opts; } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_gas_options( index_t component) { dfpmsolver::ParabolicDriverOptions* opts = nullptr; auto id = id_gas_task(component); if (id != no_equation) { opts = &static_cast( m_equation_list[id].get())->get_options(); } return opts; } // Residuals // ========= scalar_t UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual( UnsaturatedVariables * const vars ) { auto sum = m_saturation_equation->compute_squared_residual(vars); #ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum) #endif for (std::size_t ideq=0; ideqcompute_squared_residual(vars); } auto norm = std::sqrt(sum); //std::cout << "; " << norm << std::endl; return norm; } scalar_t UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual_0( UnsaturatedVariables * const vars ) { return m_norm_0; } scalar_t UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_update( UnsaturatedVariables * const vars ) { scalar_t sum = m_saturation_equation->compute_squared_update(vars); #ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum) #endif for (std::size_t ideq=0; ideqcompute_squared_update(vars); } return std::sqrt(sum); } // Solving the equations // ===================== void UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars ) { save_dt(dt); vars->get_advection_flux().set_zero(); m_saturation_equation->initialize_timestep(dt, vars); scalar_t sum = m_saturation_equation->compute_squared_residual_0(vars); #ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum) #endif for (std::size_t ideq=0; ideqinitialize_timestep(dt, vars); sum += task->compute_squared_residual_0(vars); } m_norm_0 = 1.0;//std::sqrt(sum); } solver::StaggerReturnCode UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep( UnsaturatedVariables * const vars ) { dfpmsolver::ParabolicDriverReturnCode retcode; bool flag_fail = false; bool flag_error_minimized = false; // Saturation equation retcode = m_saturation_equation->restart_timestep(vars); if (dfpmsolver::has_failed(retcode)) { WARNING << "Failed to solve saturation equation, return code :" << to_string(retcode); flag_fail = true; } if (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorMinimized) flag_error_minimized = true; // Other equation if (not flag_fail) { vars->set_relative_variables(); #ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(or: flag_fail) \ reduction(and: flag_error_minimized) #endif for (std::size_t ideq=0; ideqrestart_timestep(vars); if (dfpmsolver::has_failed(retcode)) { WARNING << "Equation of type '" << to_string(task->equation_type()) << "' for component " << task->component() << " has failed with return code : " << to_string(retcode); flag_fail = true; } flag_error_minimized = flag_error_minimized and (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorMinimized); } } // Return code solver::StaggerReturnCode return_code = solver::StaggerReturnCode::ResidualMinimized; if (flag_error_minimized) { return_code = solver::StaggerReturnCode::ErrorMinimized; } if (flag_fail) return_code = solver::StaggerReturnCode::UnknownError; return return_code; } // ================================ // // // // Helper functions // // // // ================================ // static std::string to_string(EquationType eq_type) { std::string str; switch (eq_type) { case EquationType::LiquidAqueous: str = "Liquid advection-diffusion"; break; case EquationType::Pressure: str = "Gaseous diffusion"; break; case EquationType::Saturation: str = "Saturation equation"; break; default: str = "Unknown equation"; break; } return str; } std::string to_string(dfpmsolver::ParabolicDriverReturnCode retcode) { using RetCode = dfpmsolver::ParabolicDriverReturnCode; std::string str; switch (retcode) { case RetCode::ResidualMinimized: str = "ResidualMinimized"; break; case RetCode::ErrorMinimized: str = "ErrorMinimized"; break; case RetCode::NotConvergedYet: str = "NotConvergedYet"; break; case RetCode::StationaryPoint: str = "StationaryPoint"; break; case RetCode::ErrorLinearSystem: str = "ErrorLinearSystem"; break; case RetCode::MaxStepTakenTooManyTimes: str = "MaxStepTakenTooManyTimes"; break; case RetCode::MaxIterations: str = "MaxIterations"; break; case RetCode::LinesearchFailed: str = "LinesearchFailed"; default: str = "Unknown return code"; break; } return str; } } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp diff --git a/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp b/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp index b757408..56d35cf 100644 --- a/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp +++ b/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp @@ -1,183 +1,183 @@ #include "catch.hpp" -#include "reactmicp/systems/unsaturated/transport_constraints.hpp" +#include "reactmicp/systems/unsaturated/boundary_conditions.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 "specmicp_database/database.hpp" #include "dfpm/meshes/generic_mesh1d.hpp" #include "dfpm/solver/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", "[transport],[aqueous]") { // begin initialisation index_t nb_nodes = 10; scalar_t dx = 0.1; scalar_t cross_section = 1.0; 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(D); vars_interface.set_relative_liquid_diffusivity_model(one); vars_interface.set_advection_flux(v); - TransportConstraints constraints; - constraints.add_fixed_node(0); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs->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); + AqueousTransportEquation equation(2, the_mesh, aq_vars, bcs); 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*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*D-v)); CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1*D)); 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*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*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); + AqueousTransportEquation equation(2, the_mesh, aq_vars, bcs); dfpmsolver::ParabolicDriver solver(equation); 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 << "\n ----- " << std::endl; Vector save = aq_concentration.variable; aq_concentration.set_constant(1.0); //aq_concentration(6) = 2.0; aq_concentration(7) = 2.0; aq_concentration.velocity.setZero(); aq_concentration.transport_fluxes.setZero(); std::cout << aq_concentration.variable << "\n ----- " << std::endl; vars_interface.set_advection_flux(-v); - constraints = TransportConstraints(); - constraints.add_fixed_node(9); + auto bcs2 = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs2->add_fixed_node(9); - AqueousTransportEquation equation2(the_mesh, aq_vars, constraints); + AqueousTransportEquation equation2(2, the_mesh, aq_vars, bcs2); dfpmsolver::ParabolicDriver solver2(equation2); solver2.get_options().sparse_solver = sparse_solvers::SparseSolver::SparseLU; 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 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", "HCO3[-]"}); raw_data = thedatabase.get_database(); raw_data->freeze_db(); } return raw_data; } TEST_CASE("Aqueous pressure equation", "[transport],[aqueous],[gas]") { // begin initialisation index_t nb_nodes = 10; scalar_t dx = 0.1; scalar_t cross_section = 1.0; 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(); index_t id_co2 = raw_data->get_id_component("HCO3[-]"); VariablesInterface vars_interface(the_mesh, raw_data, {0, id_co2}); Vector aq_conc(nb_nodes); aq_conc.setConstant(1.0); aq_conc(2) = 2.0; vars_interface.set_aqueous_concentration(id_co2, aq_conc); vars_interface.set_liquid_saturation(1.0); vars_interface.set_porosity(0.5); vars_interface.set_liquid_diffusivity(D); vars_interface.set_relative_liquid_diffusivity_model(one); vars_interface.set_advection_flux(v); - TransportConstraints constraints; - constraints.add_fixed_node(0); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs->add_fixed_node(0); // end initialisation SECTION("Simple Residuals") { REQUIRE(id_co2 != no_equation); UnsaturatedVariablesPtr vars = vars_interface.get_variables(); LiquidGasAqueousVariableBox aq_vars = vars->get_liquid_gas_aqueous_variables(id_co2); MainVariable& aq_concentration = aq_vars.aqueous_concentration; - AqueousGasTransportEquation equation(the_mesh, aq_vars, constraints); + AqueousGasTransportEquation equation(id_co2, the_mesh, aq_vars, bcs); 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*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*D-v)); CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1*D)); 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*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*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 (and gas) component transport \n --------------\n"; UnsaturatedVariablesPtr vars = vars_interface.get_variables(); LiquidGasAqueousVariableBox aq_vars = vars->get_liquid_gas_aqueous_variables(id_co2); MainVariable& aq_concentration = aq_vars.aqueous_concentration; - AqueousGasTransportEquation equation(the_mesh, aq_vars, constraints); + AqueousGasTransportEquation equation(id_co2, the_mesh, aq_vars, bcs); dfpmsolver::ParabolicDriver solver(equation); 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 << "\n ----- " << std::endl; Vector save = aq_concentration.variable; aq_concentration.set_constant(1.0); //aq_concentration(6) = 2.0; aq_concentration(7) = 2.0; aq_concentration.velocity.setZero(); aq_concentration.transport_fluxes.setZero(); std::cout << aq_concentration.variable << "\n ----- " << std::endl; vars_interface.set_advection_flux(-v); - constraints = TransportConstraints(); - constraints.add_fixed_node(9); + auto bcs2 = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs2->add_fixed_node(9); - AqueousGasTransportEquation equation2(the_mesh, aq_vars, constraints); + AqueousGasTransportEquation equation2(id_co2, the_mesh, aq_vars, bcs2); dfpmsolver::ParabolicDriver solver2(equation2); solver2.get_options().sparse_solver = sparse_solvers::SparseSolver::SparseLU; 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; iadd_gas_node(0); vars_interface.set_relative_gas_diffusivity_model(one); vars_interface.set_resistance_gas_diffusivity(1.0); vars_interface.set_binary_gas_diffusivity(id_co2, 0.01); Vector partial_pressure_co2(nb_nodes); partial_pressure_co2.setConstant(10); partial_pressure_co2(0) = 100; vars_interface.set_partial_pressure(id_co2, partial_pressure_co2); aq_conc.setConstant(1.0); aq_conc(0) = 1.0; vars_interface.set_aqueous_concentration(id_co2, aq_conc); UnsaturatedVariablesPtr vars = vars_interface.get_variables(); LiquidGasAqueousVariableBox aq_vars = vars->get_liquid_gas_aqueous_variables(id_co2); MainVariable& aq_concentration = aq_vars.aqueous_concentration; - AqueousGasTransportEquation equation(the_mesh, aq_vars, constraints); + AqueousGasTransportEquation equation(id_co2, the_mesh, aq_vars, bcs); dfpmsolver::ParabolicDriver solver(equation); 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 << "\n ----- " << std::endl; Vector save = aq_concentration.variable; for (index_t i=0; i= 1.0); } } } diff --git a/tests/reactmicp/systems/unsaturated/boundary_conditions.cpp b/tests/reactmicp/systems/unsaturated/boundary_conditions.cpp index f6e7083..18158de 100644 --- a/tests/reactmicp/systems/unsaturated/boundary_conditions.cpp +++ b/tests/reactmicp/systems/unsaturated/boundary_conditions.cpp @@ -1,38 +1,52 @@ #include #include "reactmicp/systems/unsaturated/boundary_conditions.hpp" #include "specmicp/adimensional/adimensional_system_structs.hpp" using namespace specmicp; using namespace specmicp::reactmicp::systems::unsaturated; -TEST_CASE("Boundary conditions", "[bcs],[constraitns]") +TEST_CASE("Boundary conditions", "[bcs],[constraints]") { SECTION("Boundary conditions") { - auto bcs = BoundaryConditions::make(5); + auto bcs = BoundaryConditions::make(5, 5); + bcs->add_gas_node(0); bcs->add_fixed_node(4); + bcs->set_fixed_flux_liquid_dof(1, 0, 3.0); + bcs->set_fixed_flux_gas_dof( 1, 2, 4.0); CHECK(bcs->is_gas_node(0) == true); - CHECK(bcs->is_fixed_node(0) == false); - CHECK(bcs->is_normal_node(0) == false); - CHECK(bcs->is_normal_node(1) == true); - CHECK(bcs->is_normal_node(2) == true); - CHECK(bcs->is_normal_node(3) == true); + CHECK(bcs->get_bcs_liquid_dof(1, 0) == IdBCs::FixedFlux); + CHECK(bcs->get_bcs_gas_dof( 1, 0) == IdBCs::NormalNode); + CHECK(bcs->get_bcs_liquid_dof(2, 0) == IdBCs::NormalNode); + CHECK(bcs->get_bcs_gas_dof( 2, 0) == IdBCs::NormalNode); + CHECK(bcs->get_bcs_liquid_dof(3, 0) == IdBCs::NormalNode); + CHECK(bcs->get_bcs_gas_dof( 3, 0) == IdBCs::NormalNode); + CHECK(bcs->get_bcs_liquid_dof(4, 0) == IdBCs::FixedDof); + CHECK(bcs->get_bcs_gas_dof( 4, 0) == IdBCs::FixedDof); + CHECK(bcs->get_bcs_liquid_dof(1, 1) == IdBCs::NormalNode); + CHECK(bcs->get_bcs_liquid_dof(2, 1) == IdBCs::NormalNode); + CHECK(bcs->get_bcs_gas_dof(4, 0) == IdBCs::FixedDof); + CHECK(bcs->get_bcs_gas_dof(4, 2) == IdBCs::FixedDof); + CHECK(bcs->get_bcs_gas_dof(4, 3) == IdBCs::FixedDof); + CHECK(bcs->get_bcs_gas_dof(1, 2) == IdBCs::FixedFlux); + + + CHECK(bcs->get_flux_liquid_dof(1, 0) == 3.0); + CHECK(bcs->get_flux_gas_dof( 1, 0) == 0.0); + CHECK(bcs->get_flux_gas_dof( 1, 2) == 4.0); - CHECK(bcs->is_fixed_node(4) == true); - CHECK(bcs->is_gas_node(4) == false); - CHECK(bcs->is_normal_node(4) == false); AdimensionalSystemConstraints& fixed_compo = bcs->fork_constraint(4, "fixed_compo"); fixed_compo.set_inert_volume_fraction(0.2); AdimensionalSystemConstraints& not_fixed = bcs->get_constraint(0); CHECK(not_fixed.inert_volume_fraction == 0.0); CHECK(bcs->get_constraint("fixed_compo").inert_volume_fraction == 0.2); } } diff --git a/tests/reactmicp/systems/unsaturated/configuration.cpp b/tests/reactmicp/systems/unsaturated/configuration.cpp index 0775c8b..876cf49 100644 --- a/tests/reactmicp/systems/unsaturated/configuration.cpp +++ b/tests/reactmicp/systems/unsaturated/configuration.cpp @@ -1,124 +1,126 @@ #include #include "reactmicp/io/configuration_unsaturated.hpp" #include "reactmicp/systems/unsaturated/boundary_conditions.hpp" #include "specmicp/adimensional/adimensional_system_solver_structs.hpp" #include "reactmicp/systems/unsaturated/equilibrium_stagger.hpp" #include "specmicp_common/io/safe_config.hpp" #include "specmicp_database/database.hpp" const char * conf_bcs = R"plop( boundary_conditions: gas_nodes: 0 fixed_nodes: 9 constraints: default: charge_keeper: "H[+]" others: - nodes: 9 name: fixed_node fork_from: default fixed_fugacity: - label_component: "CO2" label_gas: "CO2(g)" amount: 4e-4 )plop"; const char * conf_options = R"plop( equilibrium_options: default: residual_tolerance: 1e-10 step_tolerance: 1e-10 cutoff_total_concentration: 1e-14 others: - name: totally_reacted fork_from: default residual_tolerance: 1e-8 step_tolerance: 1e-8 - name: difficult node: 9 residual_tolerance: 1e-12 step_tolerance: 1e-12 )plop"; using namespace specmicp; - +using namespace specmicp::reactmicp::systems::unsaturated; 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", "HCO3[-]"} ); std::map swap = {{"HCO3[-]", "CO2"},}; thedatabase.swap_components(swap); raw_data = thedatabase.get_database(); raw_data->freeze_db(); } return raw_data; } TEST_CASE("Boundary conditions configuration", "[io],[configuration],[bcs]") { SECTION("Configuration") { auto bcs_yaml = io::YAMLConfigFile::load_from_string(conf_bcs); auto raw_data = get_database(); auto bcs = io::configure_unsaturated_boundary_conditions( 10, raw_data.get(), bcs_yaml.get_section("boundary_conditions") ); CHECK(bcs->is_gas_node(0) == true); - CHECK(bcs->is_fixed_node(9) == true); + CHECK(bcs->get_bcs_liquid_dof(9, 0) == IdBCs::FixedDof); + CHECK(bcs->get_bcs_gas_dof( 9, 0) == IdBCs::FixedDof); + CHECK(bcs->get_bcs_liquid_dof(9, 2) == IdBCs::FixedDof); const auto& default_const = bcs->get_constraint("default"); CHECK(default_const.charge_keeper == raw_data->get_id_component("H[+]")); REQUIRE(bcs->has_constraint("fixed_node")); const auto& other_const = bcs->get_constraint("fixed_node"); CHECK(other_const.charge_keeper == raw_data->get_id_component("H[+]")); const auto& from_node = bcs->get_constraint(9); CHECK(from_node.fixed_fugacity_cs.size() == 1); CHECK(from_node.fixed_fugacity_cs[0].id_gas == other_const.fixed_fugacity_cs[0].id_gas ); } } TEST_CASE("Equilibrium options", "[io],[configuration],[options]") { SECTION("Configuration") { auto opts_yaml = io::YAMLConfigFile::load_from_string(conf_options); units::UnitsSet units_set; units_set.length = units::LengthUnit::centimeter; auto opts = io::configure_unsaturated_equilibrium_options( 10, units_set, opts_yaml.get_section("equilibrium_options") ); const AdimensionalSystemSolverOptions& default_opts = opts->get("default"); CHECK(default_opts.solver_options.fvectol == 1e-10); CHECK(default_opts.solver_options.steptol == 1e-10); CHECK(default_opts.system_options.cutoff_total_concentration == 1e-14); const AdimensionalSystemSolverOptions& reacted_opts = opts->get("totally_reacted"); CHECK(reacted_opts.solver_options.fvectol == 1e-8); CHECK(reacted_opts.solver_options.steptol == 1e-8); CHECK(reacted_opts.system_options.cutoff_total_concentration == 1e-14); const AdimensionalSystemSolverOptions& diff_opts = opts->get("difficult"); CHECK(diff_opts.solver_options.fvectol == 1e-12); CHECK(diff_opts.solver_options.steptol == 1e-12); CHECK(diff_opts.system_options.cutoff_total_concentration == AdimensionalSystemOptions().cutoff_total_concentration); } } diff --git a/tests/reactmicp/systems/unsaturated/equilibrium_stagger.cpp b/tests/reactmicp/systems/unsaturated/equilibrium_stagger.cpp index d59351f..9b18dfc 100644 --- a/tests/reactmicp/systems/unsaturated/equilibrium_stagger.cpp +++ b/tests/reactmicp/systems/unsaturated/equilibrium_stagger.cpp @@ -1,192 +1,192 @@ #include "catch.hpp" #include "specmicp_database/database.hpp" #include "dfpm/meshes/generic_mesh1d.hpp" #include "reactmicp/systems/unsaturated/equilibrium_stagger.hpp" #include "reactmicp/systems/unsaturated/variables_interface.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" #include "reactmicp/systems/unsaturated/boundary_conditions.hpp" #include "specmicp/adimensional/adimensional_system_solver.hpp" #include "specmicp/adimensional/adimensional_system_solver_structs.hpp" #include "specmicp/adimensional/adimensional_system_solution_extractor.hpp" #include "specmicp/problem_solver/formulation.hpp" #include "specmicp/problem_solver/dissolver.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/cached_vector.hpp" #include #include using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::unsaturated; 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", "HCO3[-]"} ); std::map swap = {{"HCO3[-]", "CO2"},}; thedatabase.swap_components(swap); raw_data = thedatabase.get_database(); raw_data->freeze_db(); } return raw_data; } static AdimensionalSystemSolution get_init_solution(database::RawDatabasePtr the_database) { Formulation formulation; formulation.set_mass_solution(500); formulation.add_mineral("C2S", 1000.0); formulation.add_mineral("C3S", 4000.0); formulation.add_aqueous_species("CO2", 20); Dissolver dissolver(the_database); Vector tot_conc = dissolver.dissolve(formulation, false); AdimensionalSystemConstraints constraints; constraints.total_concentrations = tot_conc; AdimensionalSystemSolver adim_solver(the_database, constraints); adim_solver.get_options().solver_options.set_tolerance(1e-14); Vector x; adim_solver.initialise_variables(x, 0.8, -3); micpsolver::MiCPPerformance perf = adim_solver.solve(x); specmicp_assert(perf.return_code > micpsolver::MiCPSolverReturnCode::Success); return adim_solver.get_raw_solution(x); } TEST_CASE("EquilibriumConstraints", "[Equilibrium],[Constraints]") { init_logger(&std::cout, logger::Warning); SECTION("Initialisation") { - auto bcs = BoundaryConditions::make(10); + auto bcs = BoundaryConditions::make(10, 5); bcs->add_fixed_node(0); bcs->add_fixed_node(9); AdimensionalSystemConstraints& def_constraint = bcs->get_constraint("default"); def_constraint.set_charge_keeper(2); AdimensionalSystemConstraints& fixed_constraint = bcs->fork_constraint(0, "fixed_compo"); fixed_constraint.set_charge_keeper(3); bcs->set_constraint(9, "fixed_compo"); CHECK(bcs->get_constraint(0).charge_keeper == 3); CHECK(bcs->get_constraint(1).charge_keeper == 2); CHECK(bcs->get_constraint(2).charge_keeper == 2); CHECK(bcs->get_constraint(3).charge_keeper == 2); CHECK(bcs->get_constraint(4).charge_keeper == 2); CHECK(bcs->get_constraint(5).charge_keeper == 2); CHECK(bcs->get_constraint(6).charge_keeper == 2); CHECK(bcs->get_constraint(7).charge_keeper == 2); CHECK(bcs->get_constraint(8).charge_keeper == 2); CHECK(bcs->get_constraint(9).charge_keeper == 3); } } TEST_CASE("Equilibrium stagger", "[Equilibrium],[Stagger]") { init_logger(&std::cout, logger::Warning); index_t nb_nodes {10}; scalar_t dx {1.0}; scalar_t cross_section {5.0}; mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d({dx, nb_nodes, cross_section}); database::RawDatabasePtr raw_data = get_database(); - auto bcs = BoundaryConditions::make(the_mesh->nb_nodes()); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); auto options = EquilibriumOptionsVector::make(the_mesh->nb_nodes()); options->get("default").solver_options.set_tolerance(1e-14); index_t id_co2 = raw_data->get_id_component_from_element("C"); VariablesInterface interface(the_mesh, raw_data, {id_co2}); units::UnitsSet units_set; AdimensionalSystemSolution solution = get_init_solution(raw_data); AdimensionalSystemSolutionExtractor extr(solution, raw_data, units_set); scalar_t rho_l = extr.density_water(); interface.set_porosity(extr.porosity()); scalar_t init_saturation = extr.saturation_water(); interface.set_liquid_saturation(init_saturation); scalar_t init_water_aqconc = rho_l*extr.total_aqueous_concentration(0); interface.set_water_aqueous_concentration(init_water_aqconc); scalar_t init_water_solidconc = extr.total_solid_concentration(0); interface.set_solid_concentration(0, init_water_solidconc); for (index_t comp: raw_data->range_aqueous_component()) { interface.set_aqueous_concentration(comp, rho_l*extr.total_aqueous_concentration(comp)); interface.set_solid_concentration(comp, extr.total_solid_concentration(comp)); } // pv co2 UnsaturatedVariablesPtr vars = interface.get_variables(); scalar_t pv_co2 = extr.fugacity_gas(vars->get_id_gas(id_co2) )*vars->get_total_pressure(); std::cout << "pv co2 : " << pv_co2 << " tot pressure : " << vars->get_total_pressure() << std::endl; interface.set_partial_pressure(id_co2, pv_co2); scalar_t init_co2_aqconc = rho_l*extr.total_aqueous_concentration(id_co2); scalar_t init_co2_solidconc = extr.total_solid_concentration(id_co2); for (index_t node: the_mesh->range_nodes()) { vars->set_adim_solution(node, solution); } SECTION("Initialisation") { EquilibriumStagger stagger(vars, bcs, options); // Check that nothing has changed after we solved an already solved problem stagger.initialize(vars.get()); stagger.initialize_timestep(1.0, vars.get()); solver::StaggerReturnCode retcode = stagger.restart_timestep(vars.get()); CHECK(retcode > solver::StaggerReturnCode::NotConvergedYet); for (index_t node: the_mesh->range_nodes()) { CHECK(interface.get_liquid_saturation()(node) == init_saturation); CHECK(interface.get_water_aqueous_concentration()(node)/init_water_aqconc == Approx(1.0).epsilon(1e-14)); CHECK(interface.get_solid_concentration(0)(node) / init_water_solidconc == Approx(1.0).epsilon(1e-14)); CHECK(interface.get_aqueous_concentration(id_co2)(node)/init_co2_aqconc == Approx(1.0).epsilon(1e-14)); CHECK(interface.get_solid_concentration(id_co2)(node) / init_co2_solidconc == Approx(1.0).epsilon(1e-14)); CHECK(interface.get_partial_pressure(id_co2)(node) == Approx(pv_co2).epsilon(1e-12)); } } } diff --git a/tests/reactmicp/systems/unsaturated/pressure_equation.cpp b/tests/reactmicp/systems/unsaturated/pressure_equation.cpp index ce95e19..0844f6a 100644 --- a/tests/reactmicp/systems/unsaturated/pressure_equation.cpp +++ b/tests/reactmicp/systems/unsaturated/pressure_equation.cpp @@ -1,173 +1,173 @@ #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/transport_constraints.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 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); - - TransportConstraints constraints; - constraints.add_fixed_node(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(the_mesh, pressure_vars, constraints); + 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 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(the_mesh, pressure_vars, constraints); + PressureEquation equation(0, the_mesh, pressure_vars, bcs); dfpmsolver::ParabolicDriver 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(the_mesh, pressure_vars, constraints); + PressureEquation equation(0, the_mesh, pressure_vars, bcs); dfpmsolver::ParabolicDriver 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; } } diff --git a/tests/reactmicp/systems/unsaturated/saturation_equation.cpp b/tests/reactmicp/systems/unsaturated/saturation_equation.cpp index 7ab0c6f..2093a5c 100644 --- a/tests/reactmicp/systems/unsaturated/saturation_equation.cpp +++ b/tests/reactmicp/systems/unsaturated/saturation_equation.cpp @@ -1,247 +1,247 @@ #include "catch.hpp" -#include "reactmicp/systems/unsaturated/transport_constraints.hpp" +#include "reactmicp/systems/unsaturated/boundary_conditions.hpp" #include "reactmicp/systems/unsaturated/saturation_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 using namespace specmicp; using namespace specmicp::mesh; using namespace specmicp::reactmicp::systems::unsaturated; 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 one(index_t _, scalar_t __) {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("Saturation 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); SECTION("Simple Residuals") { // 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(); SaturationVariableBox sat_vars = vars->get_saturation_variables(); MainVariable& saturation = sat_vars.liquid_saturation; - TransportConstraints constraints; - constraints.add_fixed_node(0); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs->add_fixed_node(0); - SaturationEquation equation(the_mesh, sat_vars, constraints); + SaturationEquation 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 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") { // 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(); SaturationVariableBox sat_vars = vars->get_saturation_variables(); MainVariable& saturation = sat_vars.liquid_saturation; - TransportConstraints constraints; - constraints.add_fixed_node(0); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs->add_fixed_node(0); - SaturationEquation equation(the_mesh, sat_vars, constraints); + SaturationEquation equation(0, the_mesh, sat_vars, bcs); 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()) { 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); } } std::cout << saturation.variable << "\n ---------" << std::endl; std::cout << sat_vars.advection_flux.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()) { 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); 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)); } std::cout << saturation.variable << std::endl; } SECTION("Solving - layered") { std::cout << "-------- \n Layered system \n ---------" << std::endl; 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_2); UnsaturatedVariablesPtr vars = vars_interface.get_variables(); SaturationVariableBox sat_vars = vars->get_saturation_variables(); MainVariable& saturation = sat_vars.liquid_saturation; - TransportConstraints constraints; - constraints.add_fixed_node(0); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs->add_fixed_node(0); - SaturationEquation equation(the_mesh, sat_vars, constraints); + SaturationEquation equation(0, the_mesh, sat_vars, bcs); dfpmsolver::ParabolicDriver solver(equation); dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(0.01, saturation.variable); REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); std::cout << saturation.variable << "\n ---------" << std::endl; std::cout << sat_vars.advection_flux.variable << "\n ---------" << std::endl; retcode = solver.solve_timestep(0.01, saturation.variable); REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); std::cout << saturation.variable << "\n ---------" << std::endl; for (int i=0; i<500; ++i) { retcode = solver.solve_timestep(0.01, saturation.variable); REQUIRE((int) retcode > (int) dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet); } REQUIRE(saturation.variable(3) == Approx(saturation.variable(4)).epsilon(1e-7)); REQUIRE(saturation.variable(4) != saturation.variable(5)); REQUIRE(saturation.variable(5) == Approx(saturation.variable(6)).epsilon(1e-7)); std::cout << saturation.variable << std::endl; } } diff --git a/tests/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp b/tests/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp index b5665cb..2c86719 100644 --- a/tests/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp +++ b/tests/reactmicp/systems/unsaturated/saturation_pressure_equation.cpp @@ -1,317 +1,317 @@ #include "catch.hpp" -#include "reactmicp/systems/unsaturated/transport_constraints.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 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); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs->add_fixed_node(0); - SaturationPressureEquation equation(the_mesh, sat_vars, constraints); + 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 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); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs->add_fixed_node(0); - SaturationPressureEquation equation(the_mesh, sat_vars, constraints); + SaturationPressureEquation equation(0, the_mesh, sat_vars, bcs); 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) { 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); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); + bcs->add_fixed_node(0); + bcs->add_gas_node(0); - SaturationPressureEquation equation(the_mesh, sat_vars, constraints); + SaturationPressureEquation equation(0, the_mesh, sat_vars, bcs); 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()) { 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; } } diff --git a/tests/reactmicp/systems/unsaturated/transport_stagger.cpp b/tests/reactmicp/systems/unsaturated/transport_stagger.cpp index cd13bd1..4f75c58 100644 --- a/tests/reactmicp/systems/unsaturated/transport_stagger.cpp +++ b/tests/reactmicp/systems/unsaturated/transport_stagger.cpp @@ -1,169 +1,168 @@ #include "catch.hpp" -//#include "reactmicp/systems/unsaturated/transport_constraints.hpp" #include "reactmicp/systems/unsaturated/boundary_conditions.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" #include "reactmicp/systems/unsaturated/variables_box.hpp" #include "reactmicp/systems/unsaturated/transport_stagger.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 "reactmicp/io/reactive_transport.hpp" #include "specmicp_common/log.hpp" #include using namespace specmicp; using namespace specmicp::mesh; using namespace specmicp::reactmicp::systems::unsaturated; static scalar_t identity(index_t _, scalar_t saturation) { return saturation; } static scalar_t one_minus(index_t _, scalar_t saturation) { return 1.0 - saturation; } static scalar_t decreasing_linear(index_t _, scalar_t saturation) { return 10*(1-saturation); } 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", "HCO3[-]"}); std::map swap = {{"HCO3[-]", "CO2"},}; thedatabase.swap_components(swap); raw_data = thedatabase.get_database(); raw_data->freeze_db(); } return raw_data; } TEST_CASE("Transport stagger", "[transport][stagger]") { init_logger(&std::cout, logger::Warning); SECTION("Initialization") { std::cout << "\n ------------- \n Transport stagger \n ------------- \n"; 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(); index_t id_co2 = raw_data->get_id_component_from_element("C"); VariablesInterface vars_interface( the_mesh, raw_data, {0, id_co2} ); Vector sat_init(nb_nodes); sat_init.setConstant(0.71); sat_init(0) = 0.0; sat_init(1) = 0.6; vars_interface.set_liquid_saturation(sat_init); Vector pres_0_init(nb_nodes); pres_0_init.setConstant(1.0); pres_0_init(0) = 0.5; vars_interface.set_partial_pressure(0, pres_0_init); Vector pres_5_init(nb_nodes); pres_5_init.setConstant(1.0); pres_5_init(0) = 0.5; vars_interface.set_partial_pressure(id_co2, pres_5_init); Vector aq_concentration(nb_nodes); aq_concentration.setConstant(1.0); aq_concentration(1) = 0.5; aq_concentration(0) = 0.0; for (index_t aqcomp: raw_data->range_aqueous_component()) { vars_interface.set_aqueous_concentration(aqcomp, aq_concentration); } vars_interface.set_porosity(0.3); vars_interface.set_water_aqueous_concentration(1.0); vars_interface.set_liquid_permeability(1e-4); vars_interface.set_liquid_diffusivity(1.0); vars_interface.set_resistance_gas_diffusivity(1.0); vars_interface.set_binary_gas_diffusivity(0, 1.0); vars_interface.set_capillary_pressure_model(decreasing_linear); vars_interface.set_relative_liquid_diffusivity_model(identity); vars_interface.set_relative_liquid_permeability_model(identity); vars_interface.set_relative_gas_diffusivity_model(one_minus); - auto bcs = BoundaryConditions::make(the_mesh->nb_nodes()); + auto bcs = BoundaryConditions::make(the_mesh->nb_nodes(), raw_data->nb_component()); bcs->add_fixed_node(0); bcs->add_gas_node(0); UnsaturatedVariablesPtr vars = vars_interface.get_variables(); UnsaturatedTransportStagger stagger(vars, bcs); scalar_t res_0 = stagger.get_residual_0(vars.get()); std::cout << "res_0 : " << res_0 << std::endl; stagger.initialize_timestep(0.01, vars.get()); auto retcode = stagger.restart_timestep(vars.get()); REQUIRE(static_cast(retcode) > 0); std::cout << "Retcode : " << (int) retcode << std::endl; scalar_t res = stagger.get_residual(vars.get()); std::cout << "res : " << res << std::endl; //REQUIRE(res < 1e-6); scalar_t sat1 = vars_interface.get_liquid_saturation().norm(); scalar_t presw1 = vars_interface.get_partial_pressure(0).norm(); scalar_t presc1 = vars_interface.get_partial_pressure(id_co2).norm(); scalar_t aqc1 = vars_interface.get_aqueous_concentration(2).norm(); retcode = stagger.restart_timestep(vars.get()); std::cout << "Retcode : " << (int) retcode << std::endl; res = stagger.get_residual(vars.get()); std::cout << "res : " << res << std::endl; scalar_t sat2 = vars_interface.get_liquid_saturation().norm(); scalar_t presw2 = vars_interface.get_partial_pressure(0).norm(); scalar_t presc2 = vars_interface.get_partial_pressure(id_co2).norm(); scalar_t aqc2 = vars_interface.get_aqueous_concentration(2).norm(); REQUIRE(sat1 == Approx(sat2)); REQUIRE(presw1 == Approx(presw2)); REQUIRE(presc1 == Approx(presc2)); REQUIRE(aqc1 == Approx(aqc2)); for (int i=0; i<100; ++i) { stagger.initialize_timestep(0.01, vars.get()); retcode = stagger.restart_timestep(vars.get()); REQUIRE(static_cast(retcode) > 0); } std::cout << vars->get_liquid_saturation().variable << "\n - \n"; std::cout << vars->get_aqueous_concentration(2).variable << "\n - \n"; std::cout << vars->get_aqueous_concentration(3).variable << "\n - \n"; std::cout << vars->get_pressure_main_variables(0).variable << "\n - \n"; std::cout << vars->get_pressure_main_variables(id_co2).variable << "\n -" << std::endl; } }