diff --git a/src/reactmicp/CMakeLists.txt b/src/reactmicp/CMakeLists.txt index 17e5be1..d1fce14 100644 --- a/src/reactmicp/CMakeLists.txt +++ b/src/reactmicp/CMakeLists.txt @@ -1,163 +1,164 @@ # 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 + ${REACTMICP_SYSTEMS_DIR}/unsaturated/transport_stagger.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/transport_stagger.cpp b/src/reactmicp/systems/unsaturated/transport_stagger.cpp new file mode 100644 index 0000000..d7a08ad --- /dev/null +++ b/src/reactmicp/systems/unsaturated/transport_stagger.cpp @@ -0,0 +1,653 @@ +/*------------------------------------------------------------------------------- + +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 "transport_stagger.hpp" +#include "../../solver/staggers_base/stagger_structs.hpp" +#include "variables.hpp" +#include "variables_box.hpp" + +#include "saturation_equation.hpp" +#include "pressure_equation.hpp" +#include "aqueous_equation.hpp" + +#include "../../../database/database_holder.hpp" + +#include "../../../dfpmsolver/parabolic_driver.hpp" + +#include "../../../utils/log.hpp" + +namespace specmicp { +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 VariablesBasePtr = solver::VariablesBasePtr; + + +// =================================== // +// // +// 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, + const TransportConstraints& constraints + ); + + // Main functions + // -------------- + + //! \brief Return the residual + scalar_t get_residual(UnsaturatedVariables* vars); + //! \brief Return the first residual (start of timestep) + scalar_t get_residual_0(UnsaturatedVariables* vars); + //! \brief Return the update (norm of velocity vector) + scalar_t get_update(UnsaturatedVariables* vars); + + //! \brief Initialize timestep + void initialize_timestep( + scalar_t dt, + UnsaturatedVariables* 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;} + +private: + + // 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; + + +}; + +// Main functions +// =============== + +// call of the implementation + + +UnsaturatedTransportStagger::UnsaturatedTransportStagger( + std::shared_ptr variables, + const TransportConstraints& constraints + ): + m_impl(make_unique(variables, constraints)) +{ +} + +UnsaturatedTransportStagger::~UnsaturatedTransportStagger() = default; + + +static UnsaturatedVariables* get_var(VariablesBasePtr var) { + return static_cast(var.get()); +} + +void +UnsaturatedTransportStagger::initialize_timestep( + scalar_t dt, + VariablesBasePtr var + ) +{ + UnsaturatedVariables* true_var = get_var(var); + m_impl->initialize_timestep(dt, true_var); + +} + +solver::StaggerReturnCode +UnsaturatedTransportStagger::restart_timestep(VariablesBasePtr var) +{ + UnsaturatedVariables* true_var = get_var(var); + return m_impl->restart_timestep(true_var); +} + +scalar_t +UnsaturatedTransportStagger::get_residual(VariablesBasePtr var) +{ + UnsaturatedVariables* true_var = get_var(var); + return m_impl->get_residual(true_var); +} + +scalar_t +UnsaturatedTransportStagger::get_residual_0(VariablesBasePtr var) +{ + UnsaturatedVariables* true_var = get_var(var); + return m_impl->get_residual_0(true_var); +} + +scalar_t +UnsaturatedTransportStagger::get_update(VariablesBasePtr var) +{ + UnsaturatedVariables* true_var = get_var(var); + return m_impl->get_update(true_var); +} + +// =============================== // +// =============================== // +// // +// 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) + {} + + + //! \brief Compute the squared residuals + virtual scalar_t compute_squared_residual(UnsaturatedVariables* vars) = 0; + //! \brief Compute the squared residuals at the beginning of the timestep + virtual scalar_t compute_squared_residual_0(UnsaturatedVariables* vars) = 0; + //! \brief Compute the squared update of the variables + virtual scalar_t compute_squared_update(UnsaturatedVariables* vars) = 0; + + //! \brief Initialize the timestep + virtual void initialize_timestep(scalar_t dt, UnsaturatedVariables* vars) = 0; + //! \brief Solve the governing equation + virtual dfpmsolver::ParabolicDriverReturnCode restart_timestep(UnsaturatedVariables* 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 + ): + TaskBase(component, traits::equation_type), + m_equation(the_mesh, variables, constraints), + m_solver(m_equation) + {} + + Derived* derived() {return static_cast(this);} + + MainVariable& get_var(UnsaturatedVariables* vars) { + return derived()->get_var(vars); + } + + scalar_t compute_squared_residual(UnsaturatedVariables* vars) override { + Vector residuals; + const MainVariable& main = get_var(vars); + m_equation.compute_residuals(main.variable, main.velocity, residuals, true); + return residuals.squaredNorm(); + } + scalar_t compute_squared_residual_0(UnsaturatedVariables* vars) override { + Vector residuals; + const MainVariable& main = get_var(vars); + m_equation.compute_residuals(main.variable, main.velocity, residuals, false); + return residuals.squaredNorm(); + } + scalar_t compute_squared_update(UnsaturatedVariables *vars) override { + const MainVariable& main = get_var(vars); + return main.velocity.squaredNorm(); + } + + void initialize_timestep(scalar_t dt, UnsaturatedVariables* vars) override { + m_solver.initialize_timestep(dt, get_var(vars).variable); + } + + dfpmsolver::ParabolicDriverReturnCode restart_timestep(UnsaturatedVariables *vars) override { + MainVariable& main = get_var(vars); + m_solver.velocity() = main.velocity; + dfpmsolver::ParabolicDriverReturnCode 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();} + +private: + 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 + ): + EquationTask(component, the_mesh, variables, constraints) + {} + + MainVariable& get_var(UnsaturatedVariables* vars) { + return vars->get_liquid_saturation(); + } + + // Also take into account the solid total concentration + scalar_t compute_squared_update(UnsaturatedVariables *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 + ): + EquationTask(component, the_mesh, variables, constraints) + {} + + MainVariable& get_var(UnsaturatedVariables *vars) { + return vars->get_aqueous_concentration(component()); + } + + // Also take into account the solid total concentration + scalar_t compute_squared_update(UnsaturatedVariables *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 +{ +public: + PressureTask( + index_t component, + mesh::Mesh1DPtr the_mesh, + PressureTaskTraits::VariableBoxT&& variables, + const TransportConstraints& constraints + ): + EquationTask(component, the_mesh, variables, constraints) + {} + + MainVariable& get_var(UnsaturatedVariables* vars) { + return vars->get_pressure_main_variables(component()); + } +}; + +// ================================== // +// // +// UnsaturatedTransportStaggerImpl // +// // +// ================================== // + +// constructor +// =========== + +UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::UnsaturatedTransportStaggerImpl( + UnsaturatedVariablesPtr variables, + const TransportConstraints& constraints + ) +{ + database::RawDatabasePtr raw_db = variables->get_database(); + mesh::Mesh1DPtr the_mesh = variables->get_mesh(); + + // Saturation equations + m_saturation_equation = make_unique(0, the_mesh, + variables->get_saturation_variables(), + constraints); + + 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()) + { + m_equation_list.emplace_back( + make_unique(id, the_mesh, + variables->get_liquid_aqueous_component_variables(id), + constraints) + ); + } + // Gaseous diffusion + for (index_t id: raw_db->range_component()) + { + if (variables->component_has_gas(id)) { + m_equation_list.emplace_back( + make_unique(id, the_mesh, + variables->get_pressure_variables(id), + constraints + ) + ); + + } + } +} + +// Residuals +// ========= + +scalar_t +UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual( + UnsaturatedVariables* vars + ) +{ + scalar_t sum = m_saturation_equation->compute_squared_residual(vars); + for (auto& task: m_equation_list) { + sum += task->compute_squared_residual(vars); + } + return std::sqrt(sum); +} + +scalar_t +UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual_0( + UnsaturatedVariables* vars + ) +{ + scalar_t sum = m_saturation_equation->compute_squared_residual_0(vars); + for (auto& task: m_equation_list) { + sum += task->compute_squared_residual_0(vars); + } + return std::sqrt(sum); +} + +scalar_t +UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_update( + UnsaturatedVariables* vars + ) +{ + scalar_t sum = m_saturation_equation->compute_squared_update(vars); + for (auto& task: m_equation_list) { + sum += task->compute_squared_update(vars); + } + return std::sqrt(sum); +} + +// Solving the equations +// ===================== + +void +UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::initialize_timestep( + scalar_t dt, + UnsaturatedVariables* vars + ) +{ + save_dt(dt); + + vars->get_advection_flux().set_zero(); + + m_saturation_equation->initialize_timestep(dt, vars); + for (auto& task: m_equation_list) { + task->initialize_timestep(dt, vars); + } + +} + +solver::StaggerReturnCode +UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep( + UnsaturatedVariables* 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 or flag_error_minimized)) + { + vars->set_relative_variables(); + for (auto& task: m_equation_list) { + retcode = task->restart_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; + } + if (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorMinimized) + flag_error_minimized = true; + } + } + + // 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/src/reactmicp/systems/unsaturated/transport_stagger.hpp b/src/reactmicp/systems/unsaturated/transport_stagger.hpp new file mode 100644 index 0000000..27c1b18 --- /dev/null +++ b/src/reactmicp/systems/unsaturated/transport_stagger.hpp @@ -0,0 +1,100 @@ +/*------------------------------------------------------------------------------- + +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_TRANSPORTSTAGGER_HPP +#define SPECMICP_REACTMICP_UNSATURATED_TRANSPORTSTAGGER_HPP + +#include "../../solver/staggers_base/transport_stagger_base.hpp" + + +namespace specmicp { +namespace reactmicp { +namespace systems { +namespace unsaturated { + +class UnsaturatedVariables; +class TransportConstraints; + +class UnsaturatedTransportStagger: public solver::TransportStaggerBase +{ + using VariablesBasePtr = solver::VariablesBasePtr; + using StaggerReturnCode = solver::StaggerReturnCode; + +public: + UnsaturatedTransportStagger( + std::shared_ptr variables, + const TransportConstraints& constraints + ); + + ~UnsaturatedTransportStagger(); + + //! \brief Initialize the stagger at the beginning of an iteration + //! + //! This is where the first residual may be computed, the predictor saved, ... + //! \param dt the duration of the timestep + //! \param var shared_ptr to the variables + void initialize_timestep(scalar_t dt, VariablesBasePtr var) override; + + //! \brief Solve the equation for the timetep + //! + //! \param var shared_ptr to the variables + StaggerReturnCode restart_timestep(VariablesBasePtr var) override; + + //! \brief Compute the residuals norm + //! + //! \param var shared_ptr to the variables + scalar_t get_residual(VariablesBasePtr var) override; + //! \brief Compute the residuals norm + //! + //! \param var shared_ptr to the variables + scalar_t get_residual_0(VariablesBasePtr var) override; + + //! \brief Obtain the norm of the step size + //! + //! This is used to check if the algorithm has reach a stationary points. + //! It should look like : return main_variables.norm() + //! + //! \param var shared_ptr to the variables + scalar_t get_update(VariablesBasePtr var) override; + +private: + // The implementation details + class UnsaturatedTransportStaggerImpl; + std::unique_ptr m_impl; +}; + +} //end namespace unsaturated +} //end namespace systems +} //end namespace reactmicp +} //end namespace specmicp + +#endif // SPECMICP_REACTMICP_UNSATURATED_TRANSPORTSTAGGER_HPP diff --git a/tests/reactmicp/CMakeLists.txt b/tests/reactmicp/CMakeLists.txt index dc8f7da..f628e43 100644 --- a/tests/reactmicp/CMakeLists.txt +++ b/tests/reactmicp/CMakeLists.txt @@ -1,63 +1,64 @@ # 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 + systems/unsaturated/transport_stagger.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/transport_stagger.cpp b/tests/reactmicp/systems/unsaturated/transport_stagger.cpp new file mode 100644 index 0000000..6f3240a --- /dev/null +++ b/tests/reactmicp/systems/unsaturated/transport_stagger.cpp @@ -0,0 +1,155 @@ +#include "catch.hpp" + +#include "reactmicp/systems/unsaturated/transport_constraints.hpp" +#include "reactmicp/systems/unsaturated/variables.hpp" +#include "reactmicp/systems/unsaturated/variables_box.hpp" +#include "reactmicp/systems/unsaturated/transport_stagger.hpp" + +#include "database/database.hpp" +#include "dfpm/meshes/generic_mesh1d.hpp" + +#include "dfpmsolver/parabolic_driver.hpp" + +#include "reactmicp/io/reactive_transport.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]") { + 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(); + std::vector has_gas = {true, false, false, false, false, true}; + std::shared_ptr vars = std::make_shared( + the_mesh, raw_data, has_gas); + + vars->set_capillary_pressure_model(decreasing_linear); + vars->set_relative_liquid_diffusivity_model(identity); + vars->set_relative_liquid_permeability_model(identity); + vars->set_relative_gas_diffusivity_model(one_minus); + + SaturationVariableBox sat_vars = vars->get_saturation_variables(); + MainVariable& saturation = sat_vars.liquid_saturation; + + saturation.variable.setConstant(0.61); + saturation.variable(0) = 0.0; + saturation.variable(1) = 0.6; + + saturation.predictor = saturation.variable; + + vars->get_porosity().set_constant(0.3); + vars->get_water_aqueous_concentration().set_constant(1.0); + vars->get_liquid_permeability().set_constant(1e-4); + vars->get_liquid_diffusivity().set_constant(1); + vars->get_gas_diffusivity().set_constant(1); + + MainVariable& pressure_water = vars->get_pressure_main_variables(0); + pressure_water.set_constant(1.0); + pressure_water.variable(0) = 0.5; + pressure_water.predictor = pressure_water.variable; + + MainVariable& pressure_co2 = vars->get_pressure_main_variables(5); + pressure_co2.set_constant(0.0); + pressure_co2.variable(0) = 1.0; + pressure_co2.predictor = pressure_co2.variable; + + for (index_t aqcomp: raw_data->range_aqueous_component()) + { + MainVariable& aq_concentration = vars->get_aqueous_concentration(aqcomp); + + aq_concentration.set_constant(1.0); + aq_concentration(1) = 0.5; + aq_concentration(0) = 0.0; + + aq_concentration.predictor = aq_concentration.variable; + } + + TransportConstraints constraints; + constraints.add_fixed_node(0); + constraints.add_gas_node(0); + + + UnsaturatedTransportStagger stagger(vars, constraints); + + scalar_t res_0 = stagger.get_residual_0(vars); + + std::cout << "res_0 : " << res_0 << std::endl; + + stagger.initialize_timestep(0.01, vars); + auto retcode = stagger.restart_timestep(vars); + REQUIRE(static_cast(retcode) > 0); + std::cout << "Retcode : " << (int) retcode << std::endl; + + scalar_t res = stagger.get_residual(vars); + std::cout << "res : " << res << std::endl; + //REQUIRE(res < 1e-6); + + scalar_t sat1 = saturation.variable.norm(); + scalar_t presw1 = pressure_water.variable.norm(); + scalar_t presc1 = pressure_co2.variable.norm(); + scalar_t aqc1 = vars->get_aqueous_concentration(2).variable.norm(); + + + retcode = stagger.restart_timestep(vars); + std::cout << "Retcode : " << (int) retcode << std::endl; + + res = stagger.get_residual(vars); + std::cout << "res : " << res << std::endl; + + scalar_t sat2 = saturation.variable.norm(); + scalar_t presw2 = pressure_water.variable.norm(); + scalar_t presc2 = pressure_co2.variable.norm(); + scalar_t aqc2 = vars->get_aqueous_concentration(2).variable.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); + retcode = stagger.restart_timestep(vars); + REQUIRE(static_cast(retcode) > 0); + } + + std::cout << vars->get_liquid_saturation().variable << "\n -" << std::endl; + std::cout << vars->get_aqueous_concentration(2).variable << "\n -" << std::endl; + std::cout << vars->get_aqueous_concentration(3).variable << "\n -" << std::endl; + std::cout << vars->get_pressure_main_variables(0).variable << "\n -" << std::endl; + std::cout << vars->get_pressure_main_variables(5).variable << "\n -" << std::endl; + + } +}