diff --git a/src/reactmicp/CMakeLists.txt b/src/reactmicp/CMakeLists.txt index b309c6f..17e5be1 100644 --- a/src/reactmicp/CMakeLists.txt +++ b/src/reactmicp/CMakeLists.txt @@ -1,162 +1,163 @@ # 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 + ${REACTMICP_SYSTEMS_DIR}/unsaturated/aqueous_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/aqueous_equation.cpp b/src/reactmicp/systems/unsaturated/aqueous_equation.cpp new file mode 100644 index 0000000..9033564 --- /dev/null +++ b/src/reactmicp/systems/unsaturated/aqueous_equation.cpp @@ -0,0 +1,356 @@ +/*------------------------------------------------------------------------------- + +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 "aqueous_equation.hpp" +#include "variables_box.hpp" +#include "transport_constraints.hpp" + +#include "../../../dfpm/meshes/mesh1d.hpp" +#include "../../../utils/compat.hpp" +#include "../../../physics/constants.hpp" + +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 not_initialized {-5}; + +struct SPECMICP_DLL_LOCAL AqueousTransportEquation::AqueousTransportEquationImpl +{ + bool m_store_residual_info; + mesh::Mesh1DPtr m_mesh; + LiquidAqueousComponentVariableBox 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();} + 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( + mesh::Mesh1DPtr the_mesh, + LiquidAqueousComponentVariableBox the_vars + ): + m_mesh(the_mesh), + m_vars(the_vars), + m_ideq(the_mesh->nb_nodes(), not_initialized) + {} + + +}; + +AqueousTransportEquation::AqueousTransportEquation( + mesh::Mesh1DPtr the_mesh, + LiquidAqueousComponentVariableBox& variables, + const TransportConstraints& constraints): + m_tot_ndf(the_mesh->nb_nodes()), + m_impl(make_unique(the_mesh, variables)) +{ + number_equations(constraints); +} + +AqueousTransportEquation::~AqueousTransportEquation() = default; + +index_t AqueousTransportEquation::id_equation(index_t id_dof) +{ + return m_impl->id_equation(id_dof); +} + +//! \brief Compute the residuals inside 'element' +void AqueousTransportEquation::residuals_element( + 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)) + { + // 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 = 1.0/(0.5/coeff_diff_0 + 0.5/coeff_diff_1); + + 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) + { + flux_0 = -vars.advection_flux(element)*diff_aq; + } + else + { + flux_1 = -vars.advection_flux(element)*diff_aq; + } + flux_0 += diff_flux; + flux_0 *= m_mesh->get_face_area(element); + flux_1 -= diff_flux; + flux_1 *= m_mesh->get_face_area(element); + + + // Storage + if (m_impl->store_residual_info()) + { + // fluxes to compute exchange term + vars.aqueous_concentration.transport_fluxes(node_0) += flux_0; + vars.aqueous_concentration.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 aq_tot_conc_0 = displacement(node_0); + const scalar_t saturation_0 = vars.saturation(node_0); + + 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); + + 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); + 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 aq_tot_conc_1 = displacement(node_1); + const scalar_t saturation_1 = vars.saturation(node_1); + + 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); + + 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); + transient_1 -= chemistry_1; + } + element_residual(1) = mass_coeff_1*transient_1 - flux_1; + } +} + +//! \brief Compute the residuals +void AqueousTransportEquation::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 AqueousTransportEquation::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 AqueousTransportEquation::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 AqueousTransportEquation::number_equations(const TransportConstraints& constraints) +{ + 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; + ++neq; + } + } + m_neq = neq; +} + +} //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 new file mode 100644 index 0000000..18e1b50 --- /dev/null +++ b/src/reactmicp/systems/unsaturated/aqueous_equation.hpp @@ -0,0 +1,144 @@ +/*------------------------------------------------------------------------------- + +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_AQUEOUSEQUATION_HPP +#define SPECMICP_REACTMICP_UNSATURATED_AQUEOUSEQUATION_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 Transport for liquid transport of aqueous component +class SPECMICP_DLL_LOCAL AqueousTransportEquation: + public dfpmsolver::ParabolicProgram +{ + static constexpr bool chemistry_rate = true; + static constexpr bool no_chemistry_rate = false; + +public: + AqueousTransportEquation( + mesh::Mesh1DPtr the_mesh, + LiquidAqueousComponentVariableBox& variables, + const TransportConstraints& constraints); + + ~AqueousTransportEquation(); + + //! \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 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/tests/reactmicp/CMakeLists.txt b/tests/reactmicp/CMakeLists.txt index 1289b29..dc8f7da 100644 --- a/tests/reactmicp/CMakeLists.txt +++ b/tests/reactmicp/CMakeLists.txt @@ -1,62 +1,63 @@ # 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 + systems/unsaturated/aqueous_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/aqueous_equation.cpp b/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp new file mode 100644 index 0000000..22ff9a0 --- /dev/null +++ b/tests/reactmicp/systems/unsaturated/aqueous_equation.cpp @@ -0,0 +1,182 @@ +#include "catch.hpp" + +#include "reactmicp/systems/unsaturated/transport_constraints.hpp" +#include "reactmicp/systems/unsaturated/saturation_equation.hpp" +#include "reactmicp/systems/unsaturated/aqueous_equation.hpp" +#include "reactmicp/systems/unsaturated/variables.hpp" +#include "reactmicp/systems/unsaturated/variables_box.hpp" + +#include "database/database.hpp" +#include "dfpm/meshes/generic_mesh1d.hpp" + +#include "dfpmsolver/parabolic_driver.hpp" + +#include + +using namespace specmicp; +using namespace specmicp::mesh; +using namespace specmicp::reactmicp::systems::unsaturated; + +static 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") { + 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}); + + + LiquidAqueousComponentVariableBox aq_vars = vars.get_liquid_aqueous_component_variables(2); + MainVariable& aq_concentration = aq_vars.aqueous_concentration; + + aq_concentration.set_constant(1.0); + aq_concentration(2) = 2.0; + + vars.get_porosity().set_constant(0.5); + vars.get_liquid_diffusivity().set_constant(1.0); + vars.get_relative_liquid_diffusivity().set_constant(1.0); + + vars.get_advection_flux().set_constant(0.3); + + TransportConstraints constraints; + constraints.add_fixed_node(0); + + AqueousTransportEquation equation(the_mesh, aq_vars, constraints); + + REQUIRE(equation.get_neq() == 9); + + Vector residuals; + equation.compute_residuals(aq_concentration.variable, aq_concentration.velocity, residuals); + + REQUIRE(residuals.rows() == equation.get_neq()); + + CHECK(residuals(4) == Approx(0.0).epsilon(1e-10)); + CHECK(residuals(0) == Approx(-(2.0-1.0)/0.1).epsilon(1e-10)); + CHECK(residuals(1) == Approx(+2*(2.0-1.0)/0.1-0.3*(1.0-2.0)).epsilon(1e-10)); + CHECK(residuals(2) == Approx(-(2.0-1.0)/0.1-0.3*(2.0-1.0)).epsilon(1e-10)); + + CHECK(residuals(0) == - aq_concentration.transport_fluxes(1)); + CHECK(aq_concentration.transport_fluxes(0) == 0.0); + + Eigen::SparseMatrix jacobian; + equation.compute_jacobian(aq_concentration.variable, aq_concentration.velocity, jacobian, 1.0); + + CHECK(jacobian.rows() == equation.get_neq()); + CHECK(jacobian.cols() == equation.get_neq()); + + CHECK(jacobian.coeff(1, 8) == Approx(0.0).epsilon(1e-10)); + CHECK(jacobian.coeff(1, 0) == Approx(-1/0.1-0.3)); + CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1)); + + vars.get_advection_flux().set_constant(-0.3); + + aq_concentration.transport_fluxes.setZero(); + equation.compute_residuals(aq_concentration.variable, aq_concentration.velocity, residuals); + + REQUIRE(residuals.rows() == equation.get_neq()); + + CHECK(residuals(4) == Approx(0.0).epsilon(1e-10)); + CHECK(residuals(0) == Approx(-(2.0-1.0)/0.1-0.3*(2.0-1.0)).epsilon(1e-10)); + CHECK(residuals(1) == Approx(+2*(2.0-1.0)/0.1+0.3*(2.0-1.0)).epsilon(1e-10)); + CHECK(residuals(2) == Approx(-(2.0-1.0)/0.1).epsilon(1e-10)); + + CHECK(residuals(0) == - aq_concentration.transport_fluxes(1)); + CHECK(aq_concentration.transport_fluxes(0) == 0.0); + + equation.compute_jacobian(aq_concentration.variable, aq_concentration.velocity, jacobian, 1.0); + + CHECK(jacobian.rows() == equation.get_neq()); + CHECK(jacobian.cols() == equation.get_neq()); + + CHECK(jacobian.coeff(1, 8) == Approx(0.0).epsilon(1e-10)); + CHECK(jacobian.coeff(1, 0) == Approx(-1/0.1)); + CHECK(jacobian.coeff(0, 1) == Approx(-1/0.1-0.3)); + + } + + SECTION("Solving") { + std::cout << "---------------- \n Aqueous component transport \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}); + + + LiquidAqueousComponentVariableBox aq_vars = vars.get_liquid_aqueous_component_variables(2); + MainVariable& aq_concentration = aq_vars.aqueous_concentration; + + aq_concentration.set_constant(1.0); + aq_concentration(2) = 2.0; + aq_concentration(3) = 2.0; + + vars.get_porosity().set_constant(0.5); + vars.get_liquid_saturation().set_constant(0.5); + vars.get_liquid_diffusivity().set_constant(0.01); + vars.get_relative_liquid_diffusivity().set_constant(1.0); + + vars.get_advection_flux().set_constant(0.3); + + TransportConstraints constraints; + constraints.add_fixed_node(0); + + AqueousTransportEquation equation(the_mesh, aq_vars, constraints); + + dfpmsolver::ParabolicDriver solver(equation); + solver.get_options().sparse_solver = sparse_solvers::SparseSolver::GMRES; + + dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(0.01, aq_concentration.variable); + REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); + + std::cout << aq_concentration.variable << 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(); + + vars.get_advection_flux().set_constant(-0.3); + + constraints = TransportConstraints(); + constraints.add_fixed_node(9); + + AqueousTransportEquation equation2(the_mesh, aq_vars, constraints); + + dfpmsolver::ParabolicDriver solver2(equation2); + + retcode = solver2.solve_timestep(0.01, aq_concentration.variable); + REQUIRE((int) retcode == (int) dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized); + + + std::cout << aq_concentration.variable << std::endl; + + for (int i=0; i