diff --git a/CMakeLists.txt b/CMakeLists.txt index 422a163..21fc043 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,335 +1,336 @@ project(specmicp) cmake_minimum_required(VERSION 2.8.8) # Options : # - SPECMICP_USE_OPENMP : bool : use openmp to parallelize code # # ####################### External Package ##################################### set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/") find_package(Boost REQUIRED) include_directories(${Boost_INCLUDE_DIR}) find_package(Eigen3 REQUIRED) # This module comes from the Eigen3 Package include_directories(${EIGEN3_INCLUDE_DIR}) if(EXISTS "${EIGEN3_INCLUDE_DIR}/unsupported/Eigen/src/IterativeSolvers/GMRES.h") add_definitions(-DEIGEN_UNSUPPORTED_FOUND) INCLUDE_DIRECTORIES("${EIGEN3_INCLUDE_DIR}/unsupported/") endif() ######################## Compilation flags ################################# if(UNIX) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -std=c++11 -fuse-ld=gold") SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -pedantic") SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O2 -DNDEBUG") if(SPECMICP_USE_OPENMP) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp") add_definitions(-DSPECMICP_USE_OPENMP) endif() message(STATUS "c++ flags ${CMAKE_CXX_FLAGS}") else() message(WARNING "not tested !") endif() ################ Directories ##################### # 3rd party set(JSONCPP_DIR ${PROJECT_SOURCE_DIR}/3rdparty/jsoncpp/) add_custom_target(3rdparty_header SOURCES ${JSONCPP_DIR}/json/json.h ${JSONCPP_DIR}/json/json-forwards.h ) include_directories(${PROJECT_SOURCE_DIR}/3rdparty/jsoncpp/) set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) # include include_directories(${PROJECT_SOURCE_DIR}) add_custom_target(common SOURCES ${PROJECT_SOURCE_DIR}/common.hpp ${PROJECT_SOURCE_DIR}/database.hpp ) ################ Modules ########################## # -------- MiCPSolver ---------------- set(MICPSOLVER_DIR ${PROJECT_SOURCE_DIR}/micpsolver) add_custom_target(specmisolver_inc SOURCES ${MICPSOLVER_DIR}/micpsolver_base.hpp ${MICPSOLVER_DIR}/micpsolver_base.inl ${MICPSOLVER_DIR}/micpsolver.hpp ${MICPSOLVER_DIR}/micpsolver.inl ${MICPSOLVER_DIR}/micpsolver_min.hpp ${MICPSOLVER_DIR}/micpsolver_min.inl ${MICPSOLVER_DIR}/micpsolverold.hpp ${MICPSOLVER_DIR}/micpsolverold.inl ${MICPSOLVER_DIR}/micpsolver_structs.hpp ${MICPSOLVER_DIR}/micpprog.hpp ${MICPSOLVER_DIR}/ncp_function.hpp ${MICPSOLVER_DIR}/estimate_cond_number.hpp ) # ------- SpecMiCP --------------- set(SPECMICP_DIR ${PROJECT_SOURCE_DIR}/specmicp) set(KINETICS_DIR ${SPECMICP_DIR}/kinetics) set(ADIMENSIONAL_DIR ${SPECMICP_DIR}/adimensional) set(PROBLEM_SOLVER_DIR ${SPECMICP_DIR}/problem_solver) add_custom_target(specmic_inc SOURCES #${SPECMICP_DIR}/.hpp ${SPECMICP_DIR}/simulation_box.hpp ${SPECMICP_DIR}/reduced_system_structs.hpp ${SPECMICP_DIR}/reduced_system_solver_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solution.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solver_structs.hpp ${PROBLEM_SOLVER_DIR}/formulation.hpp ${KINETICS_DIR}/kinetic_model.hpp ${KINETICS_DIR}/kinetic_variables.hpp ${KINETICS_DIR}/kinetic_system_solver_structs.hpp ) set(SPECMICP_LIBFILE ${SPECMICP_DIR}/equilibrium_data.cpp ${SPECMICP_DIR}/reduced_system.cpp ${SPECMICP_DIR}/reduced_system_solver.cpp ${SPECMICP_DIR}/reaction_path.cpp #${SPECMICP_DIR}/speciation_solver.cpp #${SPECMICP_DIR}/simple_composition.cpp ${ADIMENSIONAL_DIR}/adimensional_system.cpp ${ADIMENSIONAL_DIR}/adimensional_system_solver.cpp ${ADIMENSIONAL_DIR}/adimensional_system_solution_extractor.cpp ${ADIMENSIONAL_DIR}/equilibrium_curve.cpp ${PROBLEM_SOLVER_DIR}/dissolver.cpp # kinetic ${KINETICS_DIR}/kinetic_system.cpp ${KINETICS_DIR}/kinetic_system_solver.cpp ) add_library(specmicp ${SPECMICP_LIBFILE}) # -------- Database ---------------- set(DATABASE_DIR ${PROJECT_SOURCE_DIR}/database) set(DATABASE_LIBFILE ${DATABASE_DIR}/data_container.cpp ${DATABASE_DIR}/database.cpp ${DATABASE_DIR}/reader.cpp ${DATABASE_DIR}/selector.cpp ${DATABASE_DIR}/mineral_selector.cpp ${DATABASE_DIR}/canonicalizer.cpp ${DATABASE_DIR}/switch_basis.cpp ${JSONCPP_DIR}/jsoncpp.cpp ) add_custom_target(data_incl SOURCES ${DATABASE_DIR}/common_def.hpp ${DATABASE_DIR}/module.hpp ) add_library(specmicp_database ${DATABASE_LIBFILE}) # --------- ODE int ------------------ set(ODEINT_DIR ${PROJECT_SOURCE_DIR}/odeint) add_custom_target(odeint_incl SOURCES ${ODEINT_DIR}/odeint_types.hpp ${ODEINT_DIR}/butcher_tableau.hpp ${ODEINT_DIR}/runge_kutta_step.hpp ${ODEINT_DIR}/runge_kutta_step.inl ${ODEINT_DIR}/runge_kutta_step_structs.hpp ) # --------- DFPMSolver ----------------- set(DFPMSOLVER_DIR ${PROJECT_SOURCE_DIR}/dfpmsolver) add_custom_target(dfpmsolver_incl SOURCES ${DFPMSOLVER_DIR}/driver.hpp ${DFPMSOLVER_DIR}/driver.inl ${DFPMSOLVER_DIR}/driver_structs.hpp ${DFPMSOLVER_DIR}/dfpm_program.hpp # # Parabolic solver ${DFPMSOLVER_DIR}/parabolic_driver.hpp ${DFPMSOLVER_DIR}/parabolic_driver.inl ${DFPMSOLVER_DIR}/parabolic_program.hpp ${DFPMSOLVER_DIR}/parabolic_structs.hpp ) # -------- DFPM ------------------------ set (DFPM_DIR ${PROJECT_SOURCE_DIR}/dfpm) add_custom_target(dfpm_incl SOURCES ${DFPM_DIR}/types.hpp ${DFPM_DIR}/1dtransport/diffusion_parameters.hpp ${DFPM_DIR}/mesh.hpp ${DFPM_DIR}/meshes/mesh1dfwd.hpp ${DFPM_DIR}/meshes/mesh1d.hpp ${DFPM_DIR}/meshes/uniform_mesh1d.hpp ${DFPM_DIR}/meshes/axisymmetric_uniform_mesh1d.hpp ) set(DFPMLIB ${DFPM_DIR}/1dtransport/diffusion.cpp ) add_library(dfpm ${DFPMLIB}) # --------- ReactMiCP ----------------- set(REACTMICP_DIR ${PROJECT_SOURCE_DIR}/reactmicp) set(REACTMICP_SOLVER_DIR ${REACTMICP_DIR}/solver/) add_custom_target(reactmicp_incl SOURCES ${REACTMICP_DIR}/common.hpp # ${REACTMICP_DIR}/micpsolvers/parabolicprogram.hpp ${REACTMICP_DIR}/micpsolvers/micpsolver_structs.hpp ${REACTMICP_DIR}/micpsolvers/micpsolver.hpp ${REACTMICP_DIR}/micpsolvers/micpsolver.inl # ${REACTMICP_DIR}/systems/secondary_variables/secondary.hpp ${REACTMICP_DIR}/systems/secondary_variables/secondary.inl ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_solver_structs.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters_base.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_parameters.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/boundary_conditions.hpp ${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 ) set(REACTMICP_LIBFILES #${REACTMICP_DIR}/systems/diffusion/diffusion.cpp #${REACTMICP_DIR}/systems/diffusion/diffusion_numbering.cpp #${REACTMICP_DIR}/systems/diffusion/diffusion_secondary.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_solver.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_neutrality_solver.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_variables.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/variables.cpp ${REACTMICP_DIR}/equilibrium_curve/chemistry.cpp ${REACTMICP_DIR}/equilibrium_curve/eqcurve_extractor.cpp ${REACTMICP_DIR}/equilibrium_curve/eqcurve_coupler.cpp ${REACTMICP_DIR}/equilibrium_curve/eqcurve_solid_transport.cpp # reactive transport solver ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.cpp ) add_library(reactmicp ${REACTMICP_LIBFILES}) # ------- Utils ------------------ set(UTILS_DIR ${PROJECT_SOURCE_DIR}/utils) add_custom_target(utils_inc SOURCES ${UTILS_DIR}/log.hpp ${UTILS_DIR}/sparse_solver.hpp ${UTILS_DIR}/sparse_solver_structs.hpp ${UTILS_DIR}/options_handler.hpp ${UTILS_DIR}/perfs_handler.hpp ${UTILS_DIR}/moving_average.hpp ) # ------- Physics ---------------- set(PHYSICS_DIR ${PROJECT_SOURCE_DIR}/physics) add_custom_target(physics_inc SOURCES ${PHYSICS_DIR}/constants.hpp ${PHYSICS_DIR}/laws.hpp ${PHYSICS_DIR}/units.hpp ) # ------- Data ------------------- set(DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data) set(DATABASE_LIST ${DATA_DIR}/specmicp_database.js ${DATA_DIR}/cemdata_specmicp.js ) add_custom_target(data SOURCES ${DATABASE_LIST}) file(INSTALL ${DATABASE_LIST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data/) # ------ Documentation ----------- set( DOCS_LIST ${CMAKE_CURRENT_SOURCE_DIR}/README.md ${CMAKE_CURRENT_SOURCE_DIR}/INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/TODO ${CMAKE_CURRENT_SOURCE_DIR}/COPYING ) add_custom_target(docs SOURCES ${DOCS_LIST}) # add a target to generate API documentation with Doxygen find_package(Doxygen) if(DOXYGEN_FOUND) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/doc/Doxyfile.in ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile @ONLY) add_custom_target(doc ${DOXYGEN_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMENT "Generating API documentation with Doxygen" VERBATIM ) add_custom_target(citations SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib) file(INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/doc/) endif(DOXYGEN_FOUND) # --------- Python interface ---------- add_subdirectory( cython ) # --------- Test ---------------------- add_subdirectory( tests ) # --------- Examples ------------------ add_subdirectory( examples ) diff --git a/src/reactmicp/solver/reactive_transport_solver.cpp b/src/reactmicp/solver/reactive_transport_solver.cpp index ab898aa..8aea09d 100644 --- a/src/reactmicp/solver/reactive_transport_solver.cpp +++ b/src/reactmicp/solver/reactive_transport_solver.cpp @@ -1,15 +1,142 @@ #include "reactive_transport_solver.hpp" -#include "staggers_base/transport_stagger_base.hpp" -#include "staggers_base/chemistry_stagger_base.hpp" -#include "staggers_base/upscaling_stagger_base.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()) + { + return retcode; + } + // else + while (retcode == ReactiveTransportReturnCode::NotConvergedYet) + { + std::cout << "iter : " << residuals.nb_iterations << std::endl; + 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; + } + 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) + { + return ReactiveTransportReturnCode::TransportFailure; + } + StaggerReturnCode chemistry_ret_code = m_chemistry_stagger->restart_timestep(variables); + 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) + { + 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/reactive_transport_solver.hpp b/src/reactmicp/solver/reactive_transport_solver.hpp index dc5797a..df82396 100644 --- a/src/reactmicp/solver/reactive_transport_solver.hpp +++ b/src/reactmicp/solver/reactive_transport_solver.hpp @@ -1,53 +1,77 @@ #ifndef SPECMICP_REACTMICP_SOLVER_REACTIVETRANSPORTSOLVER_HPP #define SPECMICP_REACTMICP_SOLVER_REACTIVETRANSPORTSOLVER_HPP //! \file reactive_transport_solver.hpp The reactive transport solver #include "common.hpp" #include #include "reactive_transport_solver_structs.hpp" #include "utils/options_handler.hpp" #include "utils/perfs_handler.hpp" namespace specmicp { namespace reactmicp { namespace solver { // forward declarations class VariablesBase; using VariablesBasePtr = std::shared_ptr; -class TransportStagger; -class ChemistryStagger; -class UpscalingStagger; -using TransportStaggerPtr = std::shared_ptr; -using ChemistryStaggerPtr = std::shared_ptr; -using UpscalingStaggerPtr = std::shared_ptr; +class TransportStaggerBase; +class ChemistryStaggerBase; +class UpscalingStaggerBase; +using TransportStaggerPtr = std::shared_ptr; +using ChemistryStaggerPtr = std::shared_ptr; +using UpscalingStaggerPtr = std::shared_ptr; + +namespace internal { +struct ReactiveTransportResiduals; +} // end namespace internal //! \brief The reactive transport solver class ReactiveTransportSolver: public OptionsHandler, public PerformanceHandler { public: - //! \brief Solve a timestep - ReactiveTransportReturnCode solve_timestep(scalar_t timestep, - VariablesBasePtr variables); + ReactiveTransportSolver( + TransportStaggerPtr transport_stagger, + ChemistryStaggerPtr chemistry_stagger, + UpscalingStaggerPtr upscaling_stagger + ): + m_transport_stagger(transport_stagger), + m_chemistry_stagger(chemistry_stagger), + m_upscaling_stagger(upscaling_stagger) + {} + //! \brief Solve a timestep + ReactiveTransportReturnCode solve_timestep( + scalar_t timestep, + VariablesBasePtr variables + ); +private: + //! \brief One iteration inside the timestep + ReactiveTransportReturnCode one_iteration(VariablesBasePtr variables, + internal::ReactiveTransportResiduals& residuals + ); //! \brief Check the convergence - ReactiveTransportReturnCode check_convergence(VariablesBasePtr variables); + ReactiveTransportReturnCode check_convergence( + VariablesBasePtr variables, + const internal::ReactiveTransportResiduals& residuals + ); private: TransportStaggerPtr m_transport_stagger; //!< The transport stagger ChemistryStaggerPtr m_chemistry_stagger; //!< The chemistry stagger UpscalingStaggerPtr m_upscaling_stagger; //!< The upscaling stagger + }; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_REACTIVETRANSPORTSOLVER_HPP diff --git a/src/reactmicp/solver/reactive_transport_solver_structs.hpp b/src/reactmicp/solver/reactive_transport_solver_structs.hpp index 932ed5b..c960590 100644 --- a/src/reactmicp/solver/reactive_transport_solver_structs.hpp +++ b/src/reactmicp/solver/reactive_transport_solver_structs.hpp @@ -1,42 +1,55 @@ #ifndef SPECMICP_REACTMICP_SOLVER_REACTIVETRANSPORTSOLVERSTRUCTS_HPP #define SPECMICP_REACTMICP_SOLVER_REACTIVETRANSPORTSOLVERSTRUCTS_HPP #include "common.hpp" namespace specmicp { namespace reactmicp { namespace solver { enum class ReactiveTransportReturnCode { - UpscalingFailure, - ChemistryFailure, - TransportFailure, - MaximumIterationsReached, - StationaryPoint, - NotConvergedYet, - ResidualMinimized, - ErrorMinimized + UpscalingFailure = -13, + ChemistryFailure = -12, + TransportFailure = -11, + StaggerFailure = -10, + MaximumIterationsReached = - 2, + StationaryPoint = - 1, + NotConvergedYet = 0, + ResidualMinimized = 1, + ErrorMinimized = 2 }; struct ReactiveTransportOptions { scalar_t residuals_tolerance; scalar_t absolute_residuals_tolerance; scalar_t step_tolerance; scalar_t good_enough_tolerance; - index_t maximum_iterations; + index_t maximum_iterations; + bool implicit_upscaling; + + bool is_snia() {return maximum_iterations <= 1;} + + ReactiveTransportOptions(): + residuals_tolerance(1e-4), + absolute_residuals_tolerance(1e-16), + step_tolerance(1e-10), + good_enough_tolerance(1e-2), + maximum_iterations(100), + implicit_upscaling(false) + {} }; struct ReactiveTransportPerformance { index_t nb_iterations; ReactiveTransportReturnCode return_code; }; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_REACTIVETRANSPORTSOLVERSTRUCTS_HPP diff --git a/src/reactmicp/solver/staggers_base/chemistry_stagger_base.hpp b/src/reactmicp/solver/staggers_base/chemistry_stagger_base.hpp index f45b91a..217ef3f 100644 --- a/src/reactmicp/solver/staggers_base/chemistry_stagger_base.hpp +++ b/src/reactmicp/solver/staggers_base/chemistry_stagger_base.hpp @@ -1,35 +1,37 @@ #ifndef SPECMICP_REACTMICP_SOLVER_CHEMISTRYSTAGGERBASE_HPP #define SPECMICP_REACTMICP_SOLVER_CHEMISTRYSTAGGERBASE_HPP //! \file chemistry_stagger_base.hpp The base class for a chemistry 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 ChemistryStaggerBase { +public: + virtual ~ChemistryStaggerBase() {} //! \brief Initialize the stagger at the beginning of the computation virtual void initialize(VariablesBasePtr var) {} //! \brief Initialize the stagger at the beginning of an iteration virtual void initialize_timestep(scalar_t dt, VariablesBasePtr var) = 0; //! \brief Solve the equation for the timestep virtual StaggerReturnCode restart_timestep(VariablesBasePtr var) = 0; }; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_CHEMISTRYSTAGGERBASE_HPP diff --git a/src/reactmicp/solver/staggers_base/decl.inl b/src/reactmicp/solver/staggers_base/decl.inl index fb03d87..338db31 100644 --- a/src/reactmicp/solver/staggers_base/decl.inl +++ b/src/reactmicp/solver/staggers_base/decl.inl @@ -1,23 +1,23 @@ #ifndef SPECMICP_REACTMICP_SOLVER_DECL_INL #define SPECMICP_REACTMICP_SOLVER_DECL_INL //! \file decl.inl common declaration needed for the *stagger_base headers #include "common.hpp" #include namespace specmicp { namespace reactmicp { namespace solver { // forward declaration class VariablesBase; using VariablesBasePtr = std::shared_ptr; -class StaggerReturnCode; +enum class StaggerReturnCode; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_DECL_INL diff --git a/src/reactmicp/solver/staggers_base/staggers_base.hpp b/src/reactmicp/solver/staggers_base/staggers_base.hpp new file mode 100644 index 0000000..4c81c25 --- /dev/null +++ b/src/reactmicp/solver/staggers_base/staggers_base.hpp @@ -0,0 +1,7 @@ +//! \file staggers_base.hpp Include headers of the staggers base directory + +#include "variables_base.hpp" +#include "transport_stagger_base.hpp" +#include "chemistry_stagger_base.hpp" +#include "upscaling_stagger_base.hpp" +#include "stagger_structs.hpp" diff --git a/src/reactmicp/solver/staggers_base/transport_stagger_base.hpp b/src/reactmicp/solver/staggers_base/transport_stagger_base.hpp index cf4136c..3148075 100644 --- a/src/reactmicp/solver/staggers_base/transport_stagger_base.hpp +++ b/src/reactmicp/solver/staggers_base/transport_stagger_base.hpp @@ -1,37 +1,43 @@ #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 virtual void initialize(VariablesBasePtr var) {} //! \brief Initialize the stagger at the beginning of an iteration virtual void initialize_timestep(scalar_t dt, VariablesBasePtr var) = 0; //! \brief Solve the equation for the timestep virtual StaggerReturnCode restart_timestep(VariablesBasePtr var) = 0; //! \brief Compute the residuals norm - virtual scalar_t residual(VariablesBasePtr var) = 0; + virtual scalar_t get_residual(VariablesBasePtr var) = 0; + //! \brief Compute the residuals norm + virtual scalar_t get_residual_0(VariablesBasePtr var) = 0; + + //! \brief Obtain the update + virtual scalar_t get_update(VariablesBasePtr var) = 0; }; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_TRANSPORTSTAGGERBASE_HPP diff --git a/src/reactmicp/solver/staggers_base/upscaling_stagger_base.hpp b/src/reactmicp/solver/staggers_base/upscaling_stagger_base.hpp index 0cf43ea..5a504c1 100644 --- a/src/reactmicp/solver/staggers_base/upscaling_stagger_base.hpp +++ b/src/reactmicp/solver/staggers_base/upscaling_stagger_base.hpp @@ -1,34 +1,36 @@ #ifndef SPECMICP_REACTMICP_SOLVER_UPSCALINGSTAGGERBASE_HPP #define SPECMICP_REACTMICP_SOLVER_UPSCALINGSTAGGERBASE_HPP //! \file upscaling_stagger_base.hpp The base class for the upscaling 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 an upscaling stagger //! class UpscalingStaggerBase { +public: + virtual ~UpscalingStaggerBase() {} //! \brief Initialize the stagger at the beginning of the computation virtual void initialize(VariablesBasePtr var) {} //! \brief Initialize the stagger at the beginning of an iteration virtual void initialize_timestep(scalar_t dt, VariablesBasePtr var) = 0; //! \brief Solve the equation for the timestep virtual StaggerReturnCode restart_timestep(VariablesBasePtr var) = 0; }; } // end namespace solver } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SOLVER_UPSCALINGSTAGGERBASE_HPP