diff --git a/src/reactmicp/solver/staggers_base/parabolic_driver_help.hpp b/src/reactmicp/solver/staggers_base/parabolic_driver_help.hpp new file mode 100644 index 0000000..00de421 --- /dev/null +++ b/src/reactmicp/solver/staggers_base/parabolic_driver_help.hpp @@ -0,0 +1,87 @@ +/* ============================================================================= + + Copyright (c) 2014-2017 F. Georget Princeton University + Copyright (c) 2017 F. Georget EPFL + 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 REACTMICP_SOLVER_PARABOLICDRIVERHELP_HPP +#define REACTMICP_SOLVER_PARABOLICDRIVERHELP_HPP + +#include "stagger_structs.hpp" +#include "dfpm/solver/parabolic_structs.hpp" + +namespace specmicp { +namespace reactmicp { +namespace solver { + +//! \brief This function transform a parabolic driver Return code to the reactmicp stagger return code +inline StaggerReturnCode parabolic2StaggerReturnCode( + dfpmsolver::ParabolicDriverReturnCode retcode + ) +{ + using ParRC = dfpmsolver::ParabolicDriverReturnCode; + using StaRC = solver::StaggerReturnCode; + + StaRC ret; + + switch (retcode) + { + case ParRC::ResidualMinimized: + ret = StaRC::ResidualMinimized; + break; + case ParRC::ErrorMinimized: + ret = StaRC::ErrorMinimized; + break; + default: // error return code + switch(retcode) + { + case ParRC::MaxIterations: + ret = StaRC::MaximumIterationsReached; + break; + case ParRC::StationaryPoint: + ret = StaRC::StationaryPoint; + break; + case ParRC::NotConvergedYet: + ret = StaRC::NotConvergedYet; + break; + default: + ret = StaRC::UnknownError; + } + } + + return ret; +} + +} // namespace solver +} // namespace reactmicp +} // namespace specmicp + +#endif // REACTMICP_SOLVER_PARABOLICDRIVERHELP_HPP diff --git a/src/reactmicp/systems/saturated_react/transport_stagger.cpp b/src/reactmicp/systems/saturated_react/transport_stagger.cpp index 4aa3439..ea3fd5f 100644 --- a/src/reactmicp/systems/saturated_react/transport_stagger.cpp +++ b/src/reactmicp/systems/saturated_react/transport_stagger.cpp @@ -1,262 +1,223 @@ /* ============================================================================= 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 "variables.hpp" #include "reactmicp/solver/staggers_base/stagger_structs.hpp" #include "transport_program.hpp" #include "dfpm/solver/parabolic_driver.hpp" +#include "reactmicp/solver/staggers_base/parabolic_driver_help.hpp" #include "specmicp_common/compat.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace satdiff { // Declaration of the implementation of the stagger struct SaturatedTransportStagger::SaturatedTransportStaggerImpl { scalar_t m_dt {-1}; scalar_t m_residual_0 {-1}; SaturatedDiffusion m_program; dfpmsolver::ParabolicDriver m_solver; SaturatedTransportStaggerImpl(SaturatedVariables* variables, std::vector& list_fixed_nodes): m_program(variables, list_fixed_nodes), m_solver(m_program) {} SaturatedTransportStaggerImpl( SaturatedVariables* variables, std::vector& list_fixed_nodes, std::map& list_slave_nodes, std::vector& list_immobile_components): m_program(variables, list_fixed_nodes, list_slave_nodes, list_immobile_components), m_solver(m_program) {} void initialize_timestep( scalar_t dt, SaturatedVariables * const variables ); solver::StaggerReturnCode restart_timestep( SaturatedVariables * const variables ); scalar_t get_residual( SaturatedVariables * const var ); scalar_t get_residual_0() {return m_residual_0;} dfpmsolver::ParabolicDriverOptions& get_options() { return m_solver.get_options(); } - - //! \brief Translate a ParabolicDriverReturnCode to a StaggerReturnCode - static solver::StaggerReturnCode parabolic2StaggerReturnCode( - dfpmsolver::ParabolicDriverReturnCode retcode); }; SaturatedTransportStagger::SaturatedTransportStagger( SaturatedVariablesPtr variables, std::vector list_fixed_nodes ): m_impl(make_unique( variables.get(), list_fixed_nodes) ) { } SaturatedTransportStagger::SaturatedTransportStagger( SaturatedVariablesPtr variables, std::vector list_fixed_nodes, std::map list_slave_nodes, std::vector list_immobile_components): m_impl(make_unique( variables.get(), list_fixed_nodes, list_slave_nodes, list_immobile_components) ) { } SaturatedTransportStagger::~SaturatedTransportStagger() = default; using VarConstPtr = SaturatedVariables * const; //!< Const pointer to the variables //! \brief Initialize the stagger at the beginning of an iteration void SaturatedTransportStagger::initialize_timestep( scalar_t dt, VariablesBase * const var ) { VarConstPtr true_var = static_cast(var); m_impl->initialize_timestep(dt, true_var); } //! \brief Solve the equation for the timestep solver::StaggerReturnCode SaturatedTransportStagger::restart_timestep(VariablesBase* var) { SaturatedVariables* true_var = static_cast(var); return m_impl->restart_timestep(true_var); } scalar_t SaturatedTransportStagger::get_update(VariablesBase * const var) { return static_cast(var)->velocity().norm(); } //! \brief Compute the residuals norm scalar_t SaturatedTransportStagger::get_residual(VariablesBase * const var) { VarConstPtr true_var = static_cast(var); return m_impl->get_residual(true_var); } scalar_t SaturatedTransportStagger::get_residual_0(VariablesBase * const _) { return m_impl->get_residual_0(); } dfpmsolver::ParabolicDriverOptions& SaturatedTransportStagger::options_solver() { return m_impl->get_options(); } // ################ // // // // Implementation // // // // ############### // void SaturatedTransportStagger::SaturatedTransportStaggerImpl::initialize_timestep( scalar_t dt, SaturatedVariables * const var ) { m_dt = dt; // Set the predictor var->predictor() = var->displacement(); // Reset the values var->velocity().setZero(); var->transport_rate().setZero(); m_program.set_variables(var); // Initialize the DFPM solver m_solver.initialize_timestep(dt, var->displacement()); // Compute the initial residual without external rate : Eigen::VectorXd residuals; m_program.compute_residuals(var->displacement(), var->velocity(), residuals, m_program.NO_CHEMISTRY_RATE); m_residual_0 = residuals.norm(); } //! \brief Solve the equation for the timestep solver::StaggerReturnCode SaturatedTransportStagger::SaturatedTransportStaggerImpl::restart_timestep( SaturatedVariables * const var ) { m_solver.velocity() = var->velocity(); dfpmsolver::ParabolicDriverReturnCode retcode = m_solver.restart_timestep(var->displacement()); // copy variables if successful if (retcode > dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet) { var->velocity() = m_solver.velocity(); } - return parabolic2StaggerReturnCode(retcode); + return solver::parabolic2StaggerReturnCode(retcode); } scalar_t SaturatedTransportStagger::SaturatedTransportStaggerImpl::get_residual( SaturatedVariables * const var ) { Eigen::VectorXd residuals; m_program.compute_residuals(var->displacement(), var->velocity(), residuals); return residuals.norm(); } - -//! \brief Translate a ParabolicDriverReturnCode to a StaggerReturnCode -solver::StaggerReturnCode -SaturatedTransportStagger::SaturatedTransportStaggerImpl::parabolic2StaggerReturnCode( - dfpmsolver::ParabolicDriverReturnCode retcode - ) -{ - using ParRC = dfpmsolver::ParabolicDriverReturnCode; - using StaRC = solver::StaggerReturnCode; - - switch (retcode) - { - case ParRC::ResidualMinimized: - return StaRC::ResidualMinimized; - break; - case ParRC::ErrorMinimized: - return StaRC::ErrorMinimized; - break; - default: // error - switch(retcode) - { - case ParRC::MaxIterations: - return StaRC::MaximumIterationsReached; - break; - case ParRC::StationaryPoint: - return StaRC::StationaryPoint; - break; - case ParRC::NotConvergedYet: - return StaRC::NotConvergedYet; - break; - default: - return StaRC::UnknownError; - } - } -} - } // end namespace satdiff } // end namespace systems } // end namespace reactmicp } // end namespace specmicp