diff --git a/src/reactmicp/CMakeLists.txt b/src/reactmicp/CMakeLists.txt index 56e4a01..b309c6f 100644 --- a/src/reactmicp/CMakeLists.txt +++ b/src/reactmicp/CMakeLists.txt @@ -1,161 +1,162 @@ # ReactMiCP # ============= set(REACTMICP_SOLVER_DIR solver) set(REACTMICP_SYSTEMS_DIR systems) set(REACTMICP_EQUILIBRIUMCURVE equilibrium_curve) add_custom_target(reactmicp_incl SOURCES ${REACTMICP_SOLVER_DIR}/reactive_transport_solver_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/variables_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/transport_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/chemistry_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/upscaling_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/stagger_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/decl.inl ${REACTMICP_SOLVER_DIR}/staggers_base/staggers_base.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/variablesfwd.hpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables_sub.hpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables_box.hpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/transport_constraints.hpp ) set(REACTMICP_LIBFILES ${REACTMICP_EQUILIBRIUMCURVE}/chemistry.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_extractor.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_coupler.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_solid_transport.cpp # Flexible reactive transport solver # ---------------------------------- ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.cpp ${REACTMICP_SOLVER_DIR}/timestepper.cpp ${REACTMICP_SOLVER_DIR}/runner.cpp # ${REACTMICP_SYSTEMS_DIR}/saturated_react/variables.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_program.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/diffusive_upscaling_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/variables.cpp ${REACTMICP_SYSTEMS_DIR}/unsaturated/saturation_equation.cpp + ${REACTMICP_SYSTEMS_DIR}/unsaturated/pressure_equation.cpp io/reactive_transport.cpp io/saturated_react.cpp io/configuration.cpp ) set_visibility_flag(${REACTMICP_LIBFILES}) set_pgo_flag(${REACTMICP_LIBFILES}) if(OPENMP_FOUND) if(SPECMICP_USE_OPENMP) set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.cpp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" ) set_source_files_properties(${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.cpp PROPERTIES COMPILE_DEFINITIONS "SPECMICP_USE_OPENMP" ) endif() endif() add_library(objreactmicp OBJECT ${REACTMICP_LIBFILES}) set_property(TARGET objreactmicp PROPERTY POSITION_INDEPENDENT_CODE 1) add_library(reactmicp SHARED $) target_link_libraries(reactmicp dfpm) install(TARGETS reactmicp LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR} ) # include files # ------------- set(REACTMICP_STAGGERBASE_INCLUDE_LIST ${REACTMICP_SOLVER_DIR}/staggers_base/variables_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/transport_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/chemistry_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/upscaling_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/stagger_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/decl.inl ${REACTMICP_SOLVER_DIR}/staggers_base/staggers_base.hpp ) set(REACTMICP_SOLVER_INCLUDE_LIST ${REACTMICP_SOLVER_DIR}/reactive_transport_solver_structs.hpp ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.hpp ${REACTMICP_SOLVER_DIR}/timestepper.hpp ${REACTMICP_SOLVER_DIR}/runner.hpp ) set(REACTMICP_SATURATEDREACT_INCLUDE_LIST ${REACTMICP_SYSTEMS_DIR}/saturated_react/variablesfwd.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/variables.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_program.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/kinetic_stagger.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/diffusive_upscaling_stagger.hpp ) set(REACTMICP_EQUILIBRIUMCURVE_INCLUDE_LIST ${REACTMICP_EQUILIBRIUMCURVE}/chemistry.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_extractor.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_coupler.cpp ${REACTMICP_EQUILIBRIUMCURVE}/eqcurve_solid_transport.cpp ) set(REACTMICP_IO_INCLUDE_LIST io/reactive_transport.hpp io/saturated_react.hpp io/configuration.hpp ) install(FILES ${REACTMICP_STAGGERBASE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/solver/staggers_base ) install(FILES ${REACTMICP_SOLVER_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/solver ) install(FILES ${REACTMICP_SATURATEDREACT_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/systems/saturated_react ) install(FILES ${REACTMICP_EQUILIBRIUMCURVE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/equilibrium_curve ) install(FILES ${REACTMICP_IO_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/io ) # static libraries # ---------------- if(SPECMICP_BUILD_STATIC) add_library(reactmicp_static STATIC $) install(TARGETS reactmicp_static ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR} ) else() add_library(reactmicp_static EXCLUDE_FROM_ALL STATIC $) endif() set_target_properties(reactmicp_static PROPERTIES OUTPUT_NAME reactmicp) target_link_libraries(reactmicp_static dfpm_static) diff --git a/src/reactmicp/systems/unsaturated/pressure_equation.cpp b/src/reactmicp/systems/unsaturated/pressure_equation.cpp new file mode 100644 index 0000000..7bce335 --- /dev/null +++ b/src/reactmicp/systems/unsaturated/pressure_equation.cpp @@ -0,0 +1,328 @@ +/*------------------------------------------------------------------------------- + +Copyright (c) 2015 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 "../../../dfpm/meshes/mesh1d.hpp" +#include "../../../utils/compat.hpp" +#include "variables_box.hpp" +#include "../../../physics/constants.hpp" + + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace unsaturated { + +struct SPECMICP_DLL_LOCAL PressureEquation::PressureEquationImpl +{ + bool m_store_residual_info; + mesh::Mesh1DPtr m_mesh; + PressureVariableBox m_vars; + std::vector m_ideq; + + bool store_residual_info() { + return m_store_residual_info; + } + + void set_store_residual_info() { + m_store_residual_info = true; + } + void reset_store_residual_info() { + m_store_residual_info = false; + } + + mesh::Mesh1D* mesh() {return m_mesh.get();} + PressureVariableBox& vars() {return m_vars;} + + bool node_has_equation(index_t node) { + return m_ideq[node] != no_equation; + } + + index_t id_equation(index_t node) { + return m_ideq[node]; + } + + void fix_node(index_t node) { + m_ideq[node] = no_equation; + } + + PressureEquationImpl( + mesh::Mesh1DPtr the_mesh, + PressureVariableBox& the_vars + ): + m_mesh(the_mesh), + m_vars(the_vars), + m_ideq(the_mesh->nb_nodes(), -5) + {} + + +}; + +PressureEquation::PressureEquation( + mesh::Mesh1DPtr the_mesh, + PressureVariableBox& variables, + const TransportConstraints& constraints + ): + m_tot_ndf(the_mesh->nb_nodes()), + m_impl(make_unique(the_mesh, variables)) +{ + number_equations(constraints); +} + +PressureEquation::~PressureEquation() = default; + +index_t PressureEquation::id_equation(index_t id_dof) +{ + return m_impl->id_equation(id_dof); +} + +//! \brief Compute the residuals inside 'element' +void PressureEquation::residuals_element( + index_t element, + const Vector& displacement, + const Vector& velocity, + Eigen::Vector2d& element_residual, + bool use_chemistry_rate + ) +{ + static constexpr scalar_t rt = constants::gas_constant*(273.16+25); + + element_residual.setZero(); + + mesh::Mesh1D* m_mesh = m_impl->mesh(); + PressureVariableBox& 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); + + // Diffusion Cw + const scalar_t coeff_diff_0 = vars.gas_diffusivity(node_0) + * vars.relative_gas_diffusivity(node_0); + const scalar_t coeff_diff_1 = vars.gas_diffusivity(node_1) + * vars.relative_gas_diffusivity(node_1); + + const scalar_t coeff_diff = 1.0/(0.5/coeff_diff_0 + 0.5/coeff_diff_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); + + const scalar_t flux_0 = tot_flux; + const scalar_t flux_1 = -tot_flux; + + // Storage + if (m_impl->store_residual_info()) + { + // fluxes to compute exchange term + vars.partial_pressure.transport_fluxes(node_0) += flux_0; + vars.partial_pressure.transport_fluxes(node_1) += flux_1; + } + + + // 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); + + scalar_t transient_0 = 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); + transient_0 /= rt; + + if (use_chemistry_rate) + { + const scalar_t chemistry_0 = - vars.partial_pressure.chemistry_rate(node_0); + transient_0 -= chemistry_0; + } + element_residual(0) = mass_coeff_0*transient_0 - flux_0; + + } + + 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); + + scalar_t transient_1 = 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); + transient_1 /= rt; + + if (use_chemistry_rate) + { + const scalar_t chemistry_1 = - vars.partial_pressure.chemistry_rate(node_1); + transient_1 -= chemistry_1; + } + element_residual(1) = mass_coeff_1*transient_1 - flux_1; + } +} + +//! \brief Compute the residuals +void PressureEquation::compute_residuals( + const Vector& displacement, + const Vector& velocity, + Vector& residuals, + bool use_chemistry_rate + ) +{ + mesh::Mesh1D* m_mesh = m_impl->mesh(); + residuals.setZero(get_neq()); + m_impl->set_store_residual_info(); + Eigen::Vector2d element_residual; + for(index_t element: m_mesh->range_elements()) + { + residuals_element(element, displacement, velocity, element_residual, use_chemistry_rate); + + const index_t node_0 = m_mesh->get_node(element, 0); + if (m_impl->node_has_equation(node_0)) + { + residuals(m_impl->id_equation(node_0)) += element_residual(0); + } + const index_t node_1 = m_mesh->get_node(element, 1); + if (m_impl->node_has_equation(node_1)) + { + residuals(m_impl->id_equation(node_1)) += element_residual(1); + } + } + m_impl->reset_store_residual_info(); +} + +//! \brief Compute the jacobian +void PressureEquation::compute_jacobian(Vector& displacement, + Vector& velocity, + Eigen::SparseMatrix& jacobian, + scalar_t alphadt + ) +{ + mesh::Mesh1D* m_mesh = m_impl->mesh(); + dfpm::list_triplet_t jacob; + const index_t estimation = 3*get_neq(); + jacob.reserve(estimation); + // assume relative variables are set + for (index_t element: m_mesh->range_elements()) + { + Eigen::Vector2d element_residual_orig; + residuals_element(element, displacement, velocity, element_residual_orig, false); + + for (index_t enodec=0; enodec<2; ++enodec) + { + const index_t nodec = m_mesh->get_node(element, enodec); + if (not m_impl->node_has_equation(nodec)) continue; + + const scalar_t tmp_d = displacement(nodec); + const scalar_t tmp_v = velocity(nodec); + + scalar_t h = eps_jacobian*std::abs(tmp_v); + if (h<1e-4*eps_jacobian) h = eps_jacobian; + velocity(nodec) = tmp_v + h; + h = velocity(nodec) - tmp_v; + + displacement(nodec) = tmp_d + alphadt*h; + + Eigen::Vector2d element_residual; + residuals_element(element, displacement, velocity, element_residual, false); + + displacement(nodec) = tmp_d; + velocity(nodec) = tmp_v; + + for (index_t enoder=0; enoder<2; ++enoder) + { + const index_t noder = m_mesh->get_node(element, enoder); + if (not m_impl->node_has_equation(noder)) continue; + + jacob.push_back(dfpm::triplet_t( + m_impl->id_equation(noder), + m_impl->id_equation(nodec), + (element_residual(enoder) - element_residual_orig(enoder))/h + )); + } + } + } + jacobian = Eigen::SparseMatrix(get_neq(), get_neq()); + jacobian.setFromTriplets(jacob.begin(), jacob.end()); +} + +//! \brief Update the solution +void PressureEquation::update_solution(const Vector& update, + scalar_t lambda, + scalar_t alpha_dt, + Vector& predictor, + Vector& displacement, + Vector& velocity) +{ + for (index_t node: m_impl->mesh()->range_nodes()) + { + if (m_impl->node_has_equation(node)) + { + velocity(node) += lambda*update(m_impl->id_equation(node)); + } + displacement = predictor + alpha_dt*velocity; + } +} + + + +void PressureEquation::number_equations(const TransportConstraints& constraints) +{ + for (int fixed: constraints.fixed_nodes()) + { + m_impl->fix_node(fixed); + } + index_t neq = 0; + for (index_t node: m_impl->mesh()->range_nodes()) + { + if (m_impl->m_ideq[node] == -5) + { + m_impl->m_ideq[node] = neq; + ++neq; + } + } + m_neq = neq; +} + +} //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 new file mode 100644 index 0000000..42bfbc3 --- /dev/null +++ b/src/reactmicp/systems/unsaturated/pressure_equation.hpp @@ -0,0 +1,141 @@ +/*------------------------------------------------------------------------------- + +Copyright (c) 2015 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 + +#include "../../../types.hpp" +#include "../../../dfpmsolver/parabolic_program.hpp" + +#include "../../../dfpm/meshes/mesh1dfwd.hpp" + +#include "variables.hpp" +#include + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace unsaturated { + +class TransportConstraints; + +//! \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 dfpmsolver::ParabolicProgram +{ +public: + PressureEquation( + mesh::Mesh1DPtr the_mesh, + PressureVariableBox& variables, + const TransportConstraints& constraints); + + ~PressureEquation(); + + //! \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 m_ndf;} + //! \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 id_dof); + + //! \brief Compute the residuals inside 'element' + void residuals_element( + index_t element, + const Vector& displacement, + const Vector& velocity, + Eigen::Vector2d& element_residual, + bool use_chemistry_rate + ); + //! \brief Compute the residuals inside 'element' + 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); + + +private: + void number_equations(const TransportConstraints& constraints); + + index_t m_neq; + index_t m_ndf {1}; + index_t m_tot_ndf; + + 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/tests/reactmicp/CMakeLists.txt b/tests/reactmicp/CMakeLists.txt index 5b521d5..1289b29 100644 --- a/tests/reactmicp/CMakeLists.txt +++ b/tests/reactmicp/CMakeLists.txt @@ -1,61 +1,62 @@ # Mesh # ---- set(test_meshes_files meshes/test_meshes.cpp meshes/test_uniform_mesh_1d.cpp meshes/test_generic1d.cpp meshes/test_axisymmetric_uniform_mesh_1d.cpp meshes/test_configuration.cpp ) add_catch_test(NAME meshes SOURCES "${test_meshes_files}" LINK_LIBRARIES dfpm_static specmicp_common_static ${YAML_LIBRARIES}) FILE(COPY meshes/unif_mesh.yaml meshes/ramp_mesh.yaml meshes/unif_axisym.yaml DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) # Reactmicp : Reactive transport solver # ------------------------------------- set(test_reactmicp_files solver/test_reactive_transport_solver.cpp solver/test_coupling.cpp ) add_catch_test(NAME reactmicp SOURCES ${test_reactmicp_files} LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}) # Saturated diffusion using new reactive transport solver # ------------------------------------------------------- set(test_reactmicp_saturated_files systems/saturated_react/test_reactmicp_saturated_react.cpp systems/saturated_react/speciation_system.cpp systems/saturated_react/variables.cpp systems/saturated_react/transport_program.cpp systems/saturated_react/equilibrium_stagger.cpp ) add_catch_test(NAME reactmicp_saturated_react SOURCES ${test_reactmicp_saturated_files} LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}) # Unsaturated reactive transport system # ------------------------------------- set(test_reactmicp_unsaturated_files systems/unsaturated/test_reactmicp_unsaturated.cpp systems/unsaturated/variables.cpp systems/unsaturated/saturation_equation.cpp + systems/unsaturated/pressure_equation.cpp ) set(TEST_CEMDATA_PATH \"../../data/cemdata.yaml\") set_source_files_properties(${test_reactmicp_unsaturated_files} PROPERTIES COMPILE_DEFINITIONS "TEST_CEMDATA_PATH=${TEST_CEMDATA_PATH}" ) add_catch_test(NAME reactmicp_unsaturated SOURCES ${test_reactmicp_unsaturated_files} LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}) diff --git a/tests/reactmicp/systems/unsaturated/pressure_equation.cpp b/tests/reactmicp/systems/unsaturated/pressure_equation.cpp new file mode 100644 index 0000000..6715c6e --- /dev/null +++ b/tests/reactmicp/systems/unsaturated/pressure_equation.cpp @@ -0,0 +1,140 @@ +#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/transport_constraints.hpp" + +#include "database/database.hpp" +#include "dfpm/meshes/generic_mesh1d.hpp" + +#include "dfpmsolver/parabolic_driver.hpp" + +#include "physics/constants.hpp" + +#include + +using namespace specmicp; +using namespace specmicp::mesh; +using namespace specmicp::reactmicp::systems::unsaturated; + +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") { + + SECTION("Simple Residuals") { + 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(); + UnsaturatedVariables vars(the_mesh, raw_data, {true, false, false, false, false}); + + + PressureVariableBox pressure_vars = vars.get_vapor_pressure_variables(); + + MainVariable& pressure = pressure_vars.partial_pressure; + + vars.get_liquid_saturation().set_constant(0.5); + pressure.set_constant(0.5); + pressure.variable(2) = 0.8; + + vars.get_porosity().set_constant(0.5); + vars.get_relative_gas_diffusivity().set_constant(1.0); + vars.get_gas_diffusivity().set_constant(1.0); + + TransportConstraints constraints; + constraints.add_fixed_node(0); + + PressureEquation equation(the_mesh, pressure_vars, constraints); + + 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"; + 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(); + UnsaturatedVariables vars(the_mesh, raw_data, {true, false, false, false, false}); + + + PressureVariableBox pressure_vars = vars.get_vapor_pressure_variables(); + + MainVariable& pressure = pressure_vars.partial_pressure; + + vars.get_liquid_saturation().set_constant(0.5); + pressure.set_constant(0.5); + pressure.variable(2) = 0.8; + + vars.get_porosity().set_constant(0.5); + vars.get_relative_gas_diffusivity().set_constant(1.0); + vars.get_gas_diffusivity().set_constant(0.01); + + TransportConstraints constraints; + constraints.add_fixed_node(0); + + PressureEquation equation(the_mesh, pressure_vars, constraints); + + + 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; + + } +}