diff --git a/src/reactmicp/solver/reactive_transport_solver.cpp b/src/reactmicp/solver/reactive_transport_solver.cpp index 83e81d9..fd02aeb 100644 --- a/src/reactmicp/solver/reactive_transport_solver.cpp +++ b/src/reactmicp/solver/reactive_transport_solver.cpp @@ -1,151 +1,151 @@ #include "reactive_transport_solver.hpp" #include "staggers_base/staggers_base.hpp" #include "utils/log.hpp" #include namespace specmicp { namespace reactmicp { namespace solver { namespace internal { 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 ReactiveTransportReturnCode ReactiveTransportSolver::solve_timestep( scalar_t timestep, VariablesBasePtr variables ) { internal::ReactiveTransportResiduals residuals; // initialization // 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; return retcode; } retcode = check_convergence(variables, residuals); ++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; return retcode; } retcode = check_convergence(variables, residuals); ++residuals.nb_iterations; } // wrapup end: // upscaling, if needed if (not get_options().implicit_upscaling) m_upscaling_stagger->restart_timestep(variables); // record performance get_perfs().nb_iterations = residuals.nb_iterations; get_perfs().return_code = retcode; get_perfs().residuals = residuals.transport_residuals/residuals.transport_residual_0; return retcode; } ReactiveTransportReturnCode ReactiveTransportSolver::one_iteration( VariablesBasePtr variables, internal::ReactiveTransportResiduals& residuals ) { StaggerReturnCode transport_ret_code = m_transport_stagger->restart_timestep(variables); - if (transport_ret_code <= StaggerReturnCode::NotconvergedYet) + if (transport_ret_code <= StaggerReturnCode::NotConvergedYet) { return ReactiveTransportReturnCode::TransportFailure; } StaggerReturnCode chemistry_ret_code = m_chemistry_stagger->restart_timestep(variables); - if (chemistry_ret_code <= StaggerReturnCode::NotconvergedYet) + if (chemistry_ret_code <= StaggerReturnCode::NotConvergedYet) { return ReactiveTransportReturnCode::ChemistryFailure; } if (get_options().implicit_upscaling) { StaggerReturnCode upscaling_ret_code = m_upscaling_stagger->restart_timestep(variables); - if (upscaling_ret_code <= StaggerReturnCode::NotconvergedYet) + if (upscaling_ret_code <= StaggerReturnCode::NotConvergedYet) { return ReactiveTransportReturnCode::UpscalingFailure; } } residuals.transport_residuals = m_transport_stagger->get_residual(variables); residuals.update = m_transport_stagger->get_update(variables); } // Check the convergence ReactiveTransportReturnCode ReactiveTransportSolver::check_convergence( VariablesBasePtr _, const internal::ReactiveTransportResiduals& residuals ) { if (residuals.transport_residuals/residuals.transport_residual_0 < get_options().residuals_tolerance or residuals.transport_residuals < get_options().absolute_residuals_tolerance) { return ReactiveTransportReturnCode::ResidualMinimized; } if (residuals.update < get_options().step_tolerance) { if (residuals.transport_residuals/residuals.transport_residual_0 < get_options().good_enough_tolerance) { return ReactiveTransportReturnCode::ErrorMinimized; } else return ReactiveTransportReturnCode::StationaryPoint; } if (residuals.nb_iterations >= get_options().maximum_iterations) { return ReactiveTransportReturnCode::MaximumIterationsReached; } return ReactiveTransportReturnCode::NotConvergedYet; } } // end namespace solver } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/solver/staggers_base/stagger_structs.hpp b/src/reactmicp/solver/staggers_base/stagger_structs.hpp index ef681f2..6e4b9d3 100644 --- a/src/reactmicp/solver/staggers_base/stagger_structs.hpp +++ b/src/reactmicp/solver/staggers_base/stagger_structs.hpp @@ -1,29 +1,29 @@ #ifndef SPECMICP_REACTMICP_SOLVER_STAGGERSTRUCTS_HPP #define SPECMICP_REACTMICP_SOLVER_STAGGERSTRUCTS_HPP //! \file stagger_structs.hpp common structs for the staggers namespace specmicp { namespace reactmicp { namespace solver { //! \brief Return code error used by a stagger //! //! Success of the stagger should be tested as : //! return code > NotConvergedYet enum class StaggerReturnCode { LolThatsNotSupposedToHappen, //!< Mainly for debugging purposes UnknownError, //!< Generic error, when used, an entry in the log is necessary MaximumIterationsReached, //!< Maximum number of iterations reached in the code StationaryPoint, //!< Stagger is stuck in a Stationary point - NotconvergedYet, //!< The stagger has not converged yet + NotConvergedYet, //!< The stagger has not converged yet ResidualMinimized, //!< The residuals are minimized ErrorMinimized //!< The error is minimized }; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_STAGGERSTRUCTS_HPP