diff --git a/src/reactmicp/solver/reactive_transport_solver.cpp b/src/reactmicp/solver/reactive_transport_solver.cpp index eccfb00..2854517 100644 --- a/src/reactmicp/solver/reactive_transport_solver.cpp +++ b/src/reactmicp/solver/reactive_transport_solver.cpp @@ -1,287 +1,290 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "reactive_transport_solver.hpp" #include "staggers_base/staggers_base.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/timer.hpp" namespace specmicp { namespace reactmicp { namespace solver { namespace internal { // This contains internal information for the reactive transport solver // It is used to check the convergence struct ReactiveTransportResiduals { index_t nb_iterations; scalar_t transport_residual_0; scalar_t transport_residuals; scalar_t update; ReactiveTransportResiduals(): nb_iterations(0), transport_residual_0(-1), transport_residuals(-1), update(-1) {} }; } // end namespace internal // ReactiveTransportSolver:: // Solve a timestep // // Attention : This function uses goto to handle return code and performance ReactiveTransportReturnCode ReactiveTransportSolver::solve_timestep( scalar_t timestep, VariablesBasePtr variables ) { return solve_timestep(timestep, variables.get()); } ReactiveTransportReturnCode ReactiveTransportSolver::solve_timestep( scalar_t timestep, VariablesBase* variables ) { // start the timer Timer tot_timer; tot_timer.start(); // initialization internal::ReactiveTransportResiduals residuals; reset_perfs(); get_perfs().timestep = timestep; // copy of variables // ? m_transport_stagger->initialize_timestep(timestep, variables); m_chemistry_stagger->initialize_timestep(timestep, variables); m_upscaling_stagger->initialize_timestep(timestep, variables); // residuals.transport_residual_0 = m_transport_stagger->get_residual_0(variables); ReactiveTransportReturnCode retcode = one_iteration(variables, residuals); // check for failure if (retcode < ReactiveTransportReturnCode::StaggerFailure) { ERROR << "Failed to solve the iteration, return code : " << (int) retcode; goto set_return; } retcode = check_convergence(variables, residuals, retcode); ++residuals.nb_iterations; // if sequential non-iterative algorithm if (get_options().is_snia()) { goto end; } // else while (retcode == ReactiveTransportReturnCode::NotConvergedYet) { retcode = one_iteration(variables, residuals); // check for failure if (retcode < ReactiveTransportReturnCode::StaggerFailure) { ERROR << "Failed to solve the iteration, return code : " << (int) retcode; goto set_return; } retcode = check_convergence(variables, residuals, retcode); ++residuals.nb_iterations; } // wrapup end: // upscaling, if needed if (not get_options().implicit_upscaling) { Timer timer; timer.start(); StaggerReturnCode upscaling_ret_code = m_upscaling_stagger->restart_timestep(variables); if (upscaling_ret_code <= StaggerReturnCode::NotConvergedYet) { retcode = ReactiveTransportReturnCode::UpscalingFailure; goto set_return; } else if (upscaling_ret_code == StaggerReturnCode::UserTermination) retcode = ReactiveTransportReturnCode::UserTermination; timer.stop(); m_timer.upscaling_time += timer.elapsed_time(); } // record performance and return code set_return: - tot_timer.stop(); + if (retcode <= ReactiveTransportReturnCode::NotConvergedYet) { + m_transport_stagger->print_debug_information(variables); + } get_perfs().total_time = tot_timer.elapsed_time(); get_perfs().nb_iterations = residuals.nb_iterations; get_perfs().return_code = retcode; get_perfs().residuals = residuals.transport_residuals/residuals.transport_residual_0; get_perfs().update = residuals.update; + tot_timer.stop(); return retcode; } ReactiveTransportReturnCode ReactiveTransportSolver::one_iteration( VariablesBasePtr variables, internal::ReactiveTransportResiduals& residuals ) { return one_iteration(variables.get(), residuals); } ReactiveTransportReturnCode ReactiveTransportSolver::one_iteration( VariablesBase* variables, internal::ReactiveTransportResiduals& residuals ) { Timer timer; bool bypass = false; bool user_termination = false; // Transport // --------- timer.start(); StaggerReturnCode transport_ret_code = m_transport_stagger->restart_timestep(variables); if (transport_ret_code <= StaggerReturnCode::NotConvergedYet) { return ReactiveTransportReturnCode::TransportFailure; } else if (transport_ret_code == StaggerReturnCode::ErrorMinimized) { bypass = true; } timer.stop(); const scalar_t ttime = timer.elapsed_time(); m_timer.transport_time += ttime; get_perfs().transport_time += ttime; // Chemistry // --------- timer.start(); StaggerReturnCode chemistry_ret_code = m_chemistry_stagger->restart_timestep(variables); if (chemistry_ret_code <= StaggerReturnCode::NotConvergedYet) { return ReactiveTransportReturnCode::ChemistryFailure; } timer.stop(); const scalar_t ctime = timer.elapsed_time(); m_timer.chemistry_time += ctime; get_perfs().chemistry_time += ctime; // Upscaling // --------- if (get_options().implicit_upscaling) { timer.start(); StaggerReturnCode upscaling_ret_code = m_upscaling_stagger->restart_timestep(variables); if (upscaling_ret_code <= StaggerReturnCode::NotConvergedYet) { return ReactiveTransportReturnCode::UpscalingFailure; } else if (upscaling_ret_code == StaggerReturnCode::UserTermination) { user_termination = true; } timer.stop(); m_timer.upscaling_time += timer.elapsed_time(); } // Final residuals // --------------- residuals.transport_residuals = m_transport_stagger->get_residual(variables); residuals.update = m_transport_stagger->get_update(variables); if (user_termination) return ReactiveTransportReturnCode::UserTermination; else if (bypass) return ReactiveTransportReturnCode::TransportBypass; else return ReactiveTransportReturnCode::NotConvergedYet; } // Check the convergence ReactiveTransportReturnCode ReactiveTransportSolver::check_convergence( VariablesBase* _, const internal::ReactiveTransportResiduals& residuals, ReactiveTransportReturnCode iteration_return_code ) { const scalar_t relative_residual = residuals.transport_residuals/residuals.transport_residual_0; // Residual if (relative_residual < get_options().residuals_tolerance or residuals.transport_residuals < get_options().absolute_residuals_tolerance) { return ReactiveTransportReturnCode::ResidualMinimized; } // Step else if (residuals.update < get_options().step_tolerance) { if (relative_residual < get_options().good_enough_tolerance) { return ReactiveTransportReturnCode::ErrorMinimized; } else return ReactiveTransportReturnCode::StationaryPoint; } else if (iteration_return_code == ReactiveTransportReturnCode::TransportBypass) { return ReactiveTransportReturnCode::TransportBypass; } else if (iteration_return_code == ReactiveTransportReturnCode::UserTermination) { return ReactiveTransportReturnCode::UserTermination; } // Number of iterations else if (residuals.nb_iterations >= get_options().maximum_iterations) { if (relative_residual < get_options().good_enough_tolerance) { return ReactiveTransportReturnCode::GoodEnough; } return ReactiveTransportReturnCode::MaximumIterationsReached; } return ReactiveTransportReturnCode::NotConvergedYet; } } // end namespace solver } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/solver/staggers_base/transport_stagger_base.hpp b/src/reactmicp/solver/staggers_base/transport_stagger_base.hpp index b3d3b1b..a9d13b6 100644 --- a/src/reactmicp/solver/staggers_base/transport_stagger_base.hpp +++ b/src/reactmicp/solver/staggers_base/transport_stagger_base.hpp @@ -1,112 +1,117 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_SOLVER_TRANSPORTSTAGGERBASE_HPP #define SPECMICP_REACTMICP_SOLVER_TRANSPORTSTAGGERBASE_HPP //! \file transport_stagger_base.hpp The base class for the transport stagger // following file include main data types and forward declaration needed for a stagger #include "decl.inl" namespace specmicp { namespace reactmicp { namespace solver { //! \brief The base class for a transport stagger //! class TransportStaggerBase { public: virtual ~TransportStaggerBase() {} //! \brief Initialize the stagger at the beginning of the computation //! //! \param var raw pointer to the variables virtual void initialize(VariablesBase * const var) {} //! \brief Initialize the stagger at the beginning of the computation //! //! \param var shared_ptr to the variables void initialize(std::shared_ptr var) { return initialize(var.get()); } //! \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 raw pointer to the variables virtual void initialize_timestep( scalar_t dt, VariablesBase * const var ) = 0; //! \brief Solve the equation for the timetep //! //! \param var raw pointer to the variables virtual StaggerReturnCode restart_timestep( VariablesBase * const var ) = 0; //! \brief Compute the residuals norm //! //! \param var raw pointer to the variables virtual scalar_t get_residual( VariablesBase * const var ) = 0; //! \brief Compute the residuals norm //! //! \param var raw pointer to the variables virtual scalar_t get_residual_0( VariablesBase * const var ) = 0; //! \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 virtual scalar_t get_update( VariablesBase * const var ) = 0; + + //! \brief Print debug information when timestep has failed + virtual void print_debug_information( + VariablesBase * const var + ) {} }; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_TRANSPORTSTAGGERBASE_HPP diff --git a/src/reactmicp/systems/unsaturated/transport_stagger.cpp b/src/reactmicp/systems/unsaturated/transport_stagger.cpp index 1ae4959..0fa681c 100644 --- a/src/reactmicp/systems/unsaturated/transport_stagger.cpp +++ b/src/reactmicp/systems/unsaturated/transport_stagger.cpp @@ -1,1081 +1,1109 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "transport_stagger.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "variables.hpp" #include "variables_box.hpp" #include "boundary_conditions.hpp" #include "saturation_equation.hpp" #include "saturation_pressure_equation.hpp" #include "pressure_equation.hpp" #include "aqueous_equation.hpp" #include "aqueous_pressure_equation.hpp" #include "specmicp_database/database_holder.hpp" #include "dfpm/solver/parabolic_driver.hpp" #include "specmicp_common/log.hpp" #include "specmicp_common/config.h" #include #include namespace specmicp { namespace dfpmsolver { extern template class ParabolicDriver; extern template class ParabolicDriver; extern template class ParabolicDriver; extern template class ParabolicDriver; extern template class ParabolicDriver; } //end namespace dfpmsolver namespace reactmicp { namespace systems { namespace unsaturated { //! \brief Type of the equation //! \internal enum class EquationType { Saturation, LiquidAqueous, Pressure }; //! \brief Format type to a string (for message error) //! \internal static std::string to_string(EquationType eq_type); //! \brief Format DFPM solver retcode to a string static std::string to_string(dfpmsolver::ParabolicDriverReturnCode retcode); // forward declaration of wrapper classes // They are defined at the bootom of this file class TaskBase; class SaturationTask; using VariablesBase = solver::VariablesBase; // =================================== // // // // Declaration of the implementation // // // // =================================== // //! \brief Implementation class for the Unsaturated transport stagger //! \internal //! //! This class does all the work //! It was designed to be OpenMP compatible class UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl { public: // Constructor // ----------- UnsaturatedTransportStaggerImpl( UnsaturatedVariablesPtr variables, std::shared_ptr boundary_conditions, bool merge_saturation_pressure, bool merge_aqueous_pressure ); // Main functions // -------------- //! \brief Return the residual scalar_t get_residual(UnsaturatedVariables * const vars); //! \brief Return the first residual (start of timestep) scalar_t get_residual_0(UnsaturatedVariables * const vars); //! \brief Return the update (norm of velocity vector) scalar_t get_update(UnsaturatedVariables * const vars); //! \brief Initialize timestep void initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars); //! \brief Restart the timestep solver::StaggerReturnCode restart_timestep(UnsaturatedVariables* vars); // Timestep management // ------------------- //! \brief Save the timestep void save_dt(scalar_t dt) {m_dt = dt;} //! \brief Return the timestep scalar_t get_dt() {return m_dt;} // Task index // ---------- index_t& id_aqueous_task(index_t component) { return m_aq_equation_task[static_cast(component)]; } index_t id_aqueous_task(index_t component) const { return m_aq_equation_task[static_cast(component)]; } index_t& id_gas_task(index_t component) { return m_aq_equation_task[static_cast(component)]; } index_t id_gas_task(index_t component) const { return m_aq_equation_task[static_cast(component)]; } //! \brief Return the options of the saturation equation dfpmsolver::ParabolicDriverOptions* get_saturation_options(); //! \brief Return the options of aqueous equation of component dfpmsolver::ParabolicDriverOptions* get_aqueous_options(index_t component); //! \brief Return the options of gas equation of component dfpmsolver::ParabolicDriverOptions* get_gas_options(index_t component); + + void print_debug_information( + UnsaturatedVariables * const var + ); private: scalar_t m_norm_0 {1.0}; // timestep management scalar_t m_dt {-1.0}; //! \brief The saturation equations std::unique_ptr m_saturation_equation; // Equation and solver std::vector> m_equation_list; std::vector m_aq_equation_task; std::vector m_gas_equation_task; std::vector m_merged_gas; }; // Main functions // =============== // call of the implementation UnsaturatedTransportStagger::UnsaturatedTransportStagger( std::shared_ptr variables, std::shared_ptr boundary_conditions, bool merge_saturation_pressure, bool merge_aqueous_pressure ): m_impl(utils::make_pimpl( variables, boundary_conditions, merge_saturation_pressure, merge_aqueous_pressure )) { } UnsaturatedTransportStagger::~UnsaturatedTransportStagger() = default; static UnsaturatedVariables * const get_var(VariablesBase * const var) { return static_cast(var); } void UnsaturatedTransportStagger::initialize_timestep( scalar_t dt, VariablesBase* var ) { m_impl->initialize_timestep(dt, get_var(var)); } solver::StaggerReturnCode UnsaturatedTransportStagger::restart_timestep(VariablesBase * const var) { return m_impl->restart_timestep(get_var(var)); } scalar_t UnsaturatedTransportStagger::get_residual(VariablesBase * const var) { return m_impl->get_residual(get_var(var)); } scalar_t UnsaturatedTransportStagger::get_residual_0(VariablesBase * const var) { return m_impl->get_residual_0(get_var(var)); } scalar_t UnsaturatedTransportStagger::get_update(VariablesBase * const var) { return m_impl->get_update(get_var(var)); } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::get_saturation_options() { return m_impl->get_saturation_options(); } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::get_aqueous_options(index_t component) { return m_impl->get_aqueous_options(component); } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::get_gas_options(index_t component) { return m_impl->get_gas_options(component); } +void UnsaturatedTransportStagger::print_debug_information( + VariablesBase * const var + ) +{ + m_impl->print_debug_information(get_var(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) {} virtual ~TaskBase() {} //! \brief Compute the squared residuals virtual scalar_t compute_squared_residual( UnsaturatedVariables * const vars ) = 0; //! \brief Compute the squared residuals at the beginning of the timestep virtual scalar_t compute_squared_residual_0( UnsaturatedVariables * const vars ) = 0; //! \brief Compute the squared update of the variables virtual scalar_t compute_squared_update( UnsaturatedVariables * const vars ) = 0; //! \brief Initialize the timestep virtual void initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars ) = 0; //! \brief Solve the governing equation virtual dfpmsolver::ParabolicDriverReturnCode restart_timestep( UnsaturatedVariables * const vars ) = 0; //! \brief Return the component index (in the db) index_t component() {return m_component;} //! \brief The equation type EquationType equation_type() {return m_type;} private: EquationType m_type; index_t m_component; }; //! \brief Base class for a equation solver //! //! \internal template class SPECMICP_DLL_LOCAL EquationTask: public TaskBase { public: using EqT = typename traits::EqT; using SolverT = typename traits::SolverT; EquationTask( index_t component, mesh::Mesh1DPtr the_mesh, typename traits::VariableBoxT& variables, std::shared_ptr bcs ): TaskBase(component, traits::equation_type), m_equation(component, the_mesh, variables, bcs), m_solver(m_equation) {} EquationTask( index_t component, mesh::Mesh1DPtr the_mesh, typename traits::VariableBoxT& variables, std::shared_ptr bcs, scalar_t scaling ): TaskBase(component, traits::equation_type), m_equation(component, the_mesh, variables, bcs), m_solver(m_equation) { m_equation.set_scaling(scaling); } Derived* derived() {return static_cast(this);} // This function must be implemented by subclass MainVariable& get_var(UnsaturatedVariables * const vars) { return derived()->get_var(vars); } scalar_t compute_squared_residual( UnsaturatedVariables * const vars ) override { Vector residuals; const MainVariable& main = get_var(vars); m_equation.compute_residuals(main.variable, main.velocity, residuals, true); return residuals.squaredNorm()/m_norm_square_0; } scalar_t compute_squared_residual_0( UnsaturatedVariables * const vars ) override { Vector residuals; const MainVariable& main = get_var(vars); m_equation.compute_residuals(main.variable, main.velocity, residuals, false); m_norm_square_0 = residuals.squaredNorm(); if (m_norm_square_0 < 1e-12) m_norm_square_0 = 1.0; return 1.0; } scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { const MainVariable& main = get_var(vars); return main.velocity.squaredNorm(); } void initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars ) override { m_dt = dt; MainVariable& main = get_var(vars); main.predictor = main.variable; main.transport_fluxes.setZero(); main.velocity.setZero(); m_solver.initialize_timestep(dt, get_var(vars).variable); } dfpmsolver::ParabolicDriverReturnCode restart_timestep( UnsaturatedVariables * const vars ) override { MainVariable& main = get_var(vars); m_solver.velocity() = main.velocity; auto retcode = m_solver.restart_timestep(main.variable); if (retcode > dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet) { main.velocity = m_solver.velocity(); } return retcode; } dfpmsolver::ParabolicDriverOptions& get_options() { return m_solver.get_options(); } const dfpmsolver::ParabolicDriverOptions& get_options() const { return m_solver.get_options(); } scalar_t get_dt() {return m_dt;} protected: scalar_t m_norm_square_0 {-1}; scalar_t m_dt {-1}; EqT m_equation; SolverT m_solver; }; //! \brief Traits struct for the SaturationTask class //! //! \internal struct SPECMICP_DLL_LOCAL SaturationTaskTraits { using EqT = SaturationEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = SaturationVariableBox; static constexpr EquationType equation_type {EquationType::Saturation}; }; //! \brief Wrapper for the saturation equation solver //! //! \internal class SPECMICP_DLL_LOCAL SaturationTask: public EquationTask { using base = EquationTask; public: SaturationTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs ): EquationTask( component, the_mesh, variables, bcs) {} SaturationTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs, scalar_t scaling ): EquationTask( component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_liquid_saturation(); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the SaturationPressureTask class //! //! \internal struct SPECMICP_DLL_LOCAL SaturationPressureTaskTraits { using EqT = SaturationPressureEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = SaturationPressureVariableBox; static constexpr EquationType equation_type {EquationType::Saturation}; }; //! \brief Wrapper for the saturation equation solver //! //! \internal class SPECMICP_DLL_LOCAL SaturationPressureTask: public EquationTask { using base = EquationTask; public: SaturationPressureTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationPressureTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs ): EquationTask( component, the_mesh, variables, bcs) {} SaturationPressureTask( index_t component, mesh::Mesh1DPtr the_mesh, SaturationPressureTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs, scalar_t scaling ): EquationTask( component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_liquid_saturation(); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the LiquidAqueousTask class //! //! \internal struct SPECMICP_DLL_LOCAL LiquidAqueousTaskTraits { using EqT = AqueousTransportEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = LiquidAqueousComponentVariableBox; static constexpr EquationType equation_type {EquationType::LiquidAqueous}; }; //! \brief Wrapper for the liquid transport of aqueous component equation //! //! \internal class SPECMICP_DLL_LOCAL LiquidAqueousTask: public EquationTask { using base = EquationTask; public: LiquidAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidAqueousTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs ): EquationTask( component, the_mesh, variables, bcs) {} LiquidAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidAqueousTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs, scalar_t scaling ): EquationTask( component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_aqueous_concentration(component()); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the LiquidGasAqueousTask class //! //! \internal struct SPECMICP_DLL_LOCAL LiquidGasAqueousTaskTraits { using EqT = AqueousGasTransportEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = LiquidGasAqueousVariableBox; static constexpr EquationType equation_type {EquationType::LiquidAqueous}; }; //! \brief Wrapper for the liquid transport of aqueous component equation //! //! \internal class SPECMICP_DLL_LOCAL LiquidGasAqueousTask: public EquationTask { using base = EquationTask; public: LiquidGasAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidGasAqueousTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs ): base(component, the_mesh, variables, bcs) {} LiquidGasAqueousTask( index_t component, mesh::Mesh1DPtr the_mesh, LiquidGasAqueousTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs, scalar_t scaling ): base(component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_aqueous_concentration(component()); } // Also take into account the solid total concentration scalar_t compute_squared_update( UnsaturatedVariables * const vars ) override { scalar_t solid_update = vars->get_solid_concentration(component()).velocity.squaredNorm(); return solid_update + base::compute_squared_update(vars); } }; //! \brief Traits struct for the Pressure Task traits //! //! \internal struct SPECMICP_DLL_LOCAL PressureTaskTraits { using EqT = PressureEquation; using SolverT = dfpmsolver::ParabolicDriver; using VariableBoxT = PressureVariableBox; static constexpr EquationType equation_type {EquationType::Pressure}; }; //! \brief Wrapper for the pressure equation solver //! //! \internal class SPECMICP_DLL_LOCAL PressureTask: public EquationTask { using base = EquationTask; public: PressureTask( index_t component, mesh::Mesh1DPtr the_mesh, PressureTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs ): EquationTask( component, the_mesh, variables, bcs) {} PressureTask( index_t component, mesh::Mesh1DPtr the_mesh, PressureTaskTraits::VariableBoxT&& variables, std::shared_ptr bcs, scalar_t scaling ): EquationTask( component, the_mesh, variables, bcs, scaling) {} MainVariable& get_var(UnsaturatedVariables * const vars) { return vars->get_pressure_main_variables(component()); } }; // ================================== // // // // UnsaturatedTransportStaggerImpl // // // // ================================== // // constructor // =========== UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::UnsaturatedTransportStaggerImpl( UnsaturatedVariablesPtr variables, std::shared_ptr bcs, bool merge_saturation_pressure, bool merge_aqueous_pressure ): m_merged_gas(variables->get_database()->nb_component(), false) { database::RawDatabasePtr raw_db = variables->get_database(); mesh::Mesh1DPtr the_mesh = variables->get_mesh(); m_aq_equation_task = std::vector( raw_db->nb_component(), no_equation); m_gas_equation_task = std::vector( raw_db->nb_component(), no_equation); // Saturation equations // ==================== // // There is 2 choices for the saturation equation // Either SaturationPressureEquation or SaturationPressure if (merge_saturation_pressure and variables->component_has_gas(0)) { m_saturation_equation = make_unique( 0, the_mesh, variables->get_saturation_pressure_variables(), bcs, variables->get_aqueous_scaling(0) ); m_merged_gas[0] = true; } else { m_saturation_equation = make_unique( 0, the_mesh, variables->get_saturation_variables(), bcs, variables->get_aqueous_scaling(0) ); } const index_t size = raw_db->nb_aqueous_components() + variables->nb_gas(); m_equation_list.reserve(size); // Liquid aqueous diffusion-advection // ================================== for (index_t id: raw_db->range_aqueous_component()) { if (merge_aqueous_pressure and variables->component_has_gas(id)) { m_equation_list.emplace_back( make_unique( id, the_mesh, variables->get_liquid_gas_aqueous_variables(id), bcs, variables->get_aqueous_scaling(id)) ); m_merged_gas[id] = true; } else { m_equation_list.emplace_back( make_unique( id, the_mesh, variables->get_liquid_aqueous_component_variables(id), bcs, variables->get_aqueous_scaling(id)) ); } id_aqueous_task(id) = m_equation_list.size()-1; } // Water partial pressure // ====================== // // Depending on the saturation equation chosen we may or may not include // this equations if (variables->component_has_gas(0) and (not m_merged_gas[0])) { m_equation_list.emplace_back( make_unique(0, the_mesh, variables->get_pressure_variables(0), bcs, variables->get_gaseous_scaling(0) )); id_gas_task(0) = m_equation_list.size()-1; } // Gaseous diffusion for (index_t id: raw_db->range_aqueous_component()) { if (variables->component_has_gas(id) and (not m_merged_gas[id])) { m_equation_list.emplace_back( make_unique(id, the_mesh, variables->get_pressure_variables(id), bcs, variables->get_gaseous_scaling(id) )); id_gas_task(id) = m_equation_list.size()-1; } } } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_saturation_options() { dfpmsolver::ParabolicDriverOptions* opts = nullptr; if (m_merged_gas[0]) { opts = &static_cast( m_saturation_equation.get())->get_options(); } else { opts = &static_cast( m_saturation_equation.get())->get_options(); } return opts; } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_aqueous_options( index_t component) { dfpmsolver::ParabolicDriverOptions* opts = nullptr; auto id = id_aqueous_task(component); if (id != no_equation) { if (m_merged_gas[component]) { opts = &static_cast( m_equation_list[id].get())->get_options(); } else { opts = &static_cast( m_equation_list[id].get())->get_options(); } } return opts; } dfpmsolver::ParabolicDriverOptions* UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_gas_options( index_t component) { dfpmsolver::ParabolicDriverOptions* opts = nullptr; auto id = id_gas_task(component); if (id != no_equation) { opts = &static_cast( m_equation_list[id].get())->get_options(); } return opts; } // Residuals // ========= scalar_t UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual( UnsaturatedVariables * const vars ) { scalar_t sum = 0.0; #ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum) #endif for (std::size_t ideq=0; ideqcompute_squared_residual(vars); } sum += m_saturation_equation->compute_squared_residual(vars); auto norm = std::sqrt(sum); //std::cout << "; " << norm << std::endl; return norm; } scalar_t UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_residual_0( UnsaturatedVariables * const vars ) { return m_norm_0; } scalar_t UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::get_update( UnsaturatedVariables * const vars ) { scalar_t sum = 0.0; #ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum) #endif for (std::size_t ideq=0; ideqcompute_squared_update(vars); } sum += m_saturation_equation->compute_squared_update(vars); return std::sqrt(sum); } // Solving the equations // ===================== void UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::initialize_timestep( scalar_t dt, UnsaturatedVariables * const vars ) { // auto& aqwtilde = vars->get_water_aqueous_concentration(); // auto& saturation = vars->get_liquid_saturation(); // vars->set_aqueous_scaling(0, saturation.variable.maxCoeff()); // for (auto aqc: vars->get_database()->range_aqueous_component()) { // auto& aqtilde = vars->get_aqueous_concentration(aqc); // vars->set_aqueous_scaling(0, aqtilde.variable.maxCoeff()); // } save_dt(dt); vars->get_advection_flux().set_zero(); vars->get_water_aqueous_concentration().velocity.setZero(); vars->get_water_aqueous_concentration().predictor = vars->get_water_aqueous_concentration().variable; vars->get_porosity().velocity.setZero(); vars->get_porosity().predictor = vars->get_porosity().variable; vars->set_relative_variables(); m_saturation_equation->initialize_timestep(dt, vars); scalar_t sum = 0.0; #ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(+: sum) #endif for (std::size_t ideq=0; ideqinitialize_timestep(dt, vars); sum += task->compute_squared_residual_0(vars); } sum += m_saturation_equation->compute_squared_residual_0(vars); m_norm_0 = 1.0;//std::sqrt(sum); } solver::StaggerReturnCode UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::restart_timestep( UnsaturatedVariables * const vars ) { dfpmsolver::ParabolicDriverReturnCode retcode; //auto& satu = vars->get_liquid_saturation(); //for (index_t id=0; id 0.95) { // satu.variable(id) = 0.95; // WARNING << "Reset saturation at node : " << id; // } //} bool flag_fail = false; bool flag_error_minimized = false; // Saturation equation retcode = m_saturation_equation->restart_timestep(vars); if (dfpmsolver::has_failed(retcode)) { WARNING << "Failed to solve saturation equation, return code :" << to_string(retcode); flag_fail = true; } if (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorMinimized) flag_error_minimized = true; // Other equation if (not flag_fail) { vars->set_relative_variables(); #ifdef SPECMICP_HAVE_OPENMP #pragma omp parallel for \ reduction(or: flag_fail) \ reduction(and: flag_error_minimized) #endif for (std::size_t ideq=0; ideqrestart_timestep(vars); if (dfpmsolver::has_failed(retcode)) { WARNING << "Equation of type '" << to_string(task->equation_type()) << "' for component " << task->component() << " has failed with return code : " << to_string(retcode); flag_fail = true; } flag_error_minimized = flag_error_minimized and (retcode == dfpmsolver::ParabolicDriverReturnCode::ErrorMinimized); } } // Return code solver::StaggerReturnCode return_code = solver::StaggerReturnCode::ResidualMinimized; if (flag_error_minimized) { return_code = solver::StaggerReturnCode::ErrorMinimized; } if (flag_fail) return_code = solver::StaggerReturnCode::UnknownError; return return_code; } + +void UnsaturatedTransportStagger::UnsaturatedTransportStaggerImpl::print_debug_information( + UnsaturatedVariables * const vars + ) +{ + ERROR << "Current residuals per equation : \n" + << "Saturation equation : " << std::sqrt(m_saturation_equation->compute_squared_residual(vars)); + + for (std::size_t ideq=0; ideqequation_type()) << " - " << task->component() + << " : " + << std::sqrt(task->compute_squared_residual(vars)); + } +} + // ================================ // // // // 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 index 82fb459..c5e09d8 100644 --- a/src/reactmicp/systems/unsaturated/transport_stagger.hpp +++ b/src/reactmicp/systems/unsaturated/transport_stagger.hpp @@ -1,152 +1,157 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_REACTMICP_UNSATURATED_TRANSPORTSTAGGER_HPP #define SPECMICP_REACTMICP_UNSATURATED_TRANSPORTSTAGGER_HPP //! \file unsaturated/transport_stagger.hpp //! \brief The unsaturated transport stagger #include "reactmicp/solver/staggers_base/transport_stagger_base.hpp" #include "specmicp_common/pimpl_ptr.hpp" #include "boundary_conditions.hpp" namespace specmicp { namespace dfpmsolver { struct ParabolicDriverOptions; } //end namespace dfpmsolver namespace reactmicp { namespace systems { namespace unsaturated { class UnsaturatedVariables; class TransportConstraints; class SPECMICP_DLL_PUBLIC UnsaturatedTransportStagger: public solver::TransportStaggerBase { using VariablesBase = solver::VariablesBase; using StaggerReturnCode = solver::StaggerReturnCode; public: UnsaturatedTransportStagger( std::shared_ptr variables, std::shared_ptr constraints, bool merge_saturation_pressure=false, bool merge_aqueous_pressure=false ); //! \brief Return a shared pointer to the transport stagger //! //! \param variables shared_pointer to the unsaturated variables //! \param constraints The transport constraints //! \param merge_saturation_pressure true if the saturation //! and water vapor pressure equations must be merged static std::shared_ptr make( std::shared_ptr variables, std::shared_ptr boundary_conditions, bool merge_saturation_pressure = false, bool merge_aqueous_pressure = false ) { return std::make_shared( variables, boundary_conditions, merge_saturation_pressure, merge_aqueous_pressure ); } ~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, VariablesBase * const var) override; //! \brief Solve the equation for the timetep //! //! \param var shared_ptr to the variables StaggerReturnCode restart_timestep(VariablesBase * const var) override; //! \brief Compute the residuals norm //! //! \param var shared_ptr to the variables scalar_t get_residual(VariablesBase * const var) override; //! \brief Compute the residuals norm //! //! \param var shared_ptr to the variables scalar_t get_residual_0(VariablesBase * const 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(VariablesBase * const var) override; + //! \brief Print debug information + virtual void print_debug_information( + VariablesBase * const var + ) override; + //! Return options of saturation equation //! //! \return the options or nullptr if the equation does not exist dfpmsolver::ParabolicDriverOptions* get_saturation_options(); //! return options of aqueous transport equation //! //! \return the options or nullptr if the equation does not exist dfpmsolver::ParabolicDriverOptions* get_aqueous_options(index_t component); //! return options of gas transport equation //! //! \return the options or nullptr if the equation does not exist dfpmsolver::ParabolicDriverOptions* get_gas_options(index_t component); private: // The implementation details class UnsaturatedTransportStaggerImpl; utils::pimpl_ptr m_impl; }; } //end namespace unsaturated } //end namespace systems } //end namespace reactmicp } //end namespace specmicp #endif // SPECMICP_REACTMICP_UNSATURATED_TRANSPORTSTAGGER_HPP