diff --git a/CMakeLists.txt b/CMakeLists.txt
index de20380..00aed33 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,242 +1,252 @@
project(specmicp)
cmake_minimum_required(VERSION 2.8.8)
####################### External Package #####################################
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/")
find_package(Eigen3 REQUIRED) # This module comes from the Eigen3 Package
######################## 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")
message(STATUS "c++ flags ${CMAKE_CXX_FLAGS}")
else()
message(WARNING "not tested !")
endif()
################ Directories #####################
set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src)
# include
include_directories(${EIGEN3_INCLUDE_DIR})
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)
add_custom_target(specmic_inc SOURCES
#${SPECMICP_DIR}/.hpp
${SPECMICP_DIR}/simulation_box.hpp
+ ${SPECMICP_DIR}/reduced_system_solver_structs.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
+
+ # 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
)
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
)
# --------- ReactMiCP -----------------
set(REACTMICP_DIR ${PROJECT_SOURCE_DIR}/reactmicp)
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}/mesh.hpp
${REACTMICP_DIR}/meshes/mesh1dfwd.hpp
${REACTMICP_DIR}/meshes/mesh1d.hpp
${REACTMICP_DIR}/meshes/uniform_mesh1d.hpp
${REACTMICP_DIR}/meshes/axisymmetric_uniform_mesh1d.hpp
${REACTMICP_DIR}/systems/secondary_variables/secondary.hpp
${REACTMICP_DIR}/systems/secondary_variables/secondary.inl
${REACTMICP_DIR}/systems/saturated_diffusion/options.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
)
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
)
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
)
# ------- 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/odeint/runge_kutta_step.hpp b/src/odeint/runge_kutta_step.hpp
index ff49f64..2eb5cf4 100644
--- a/src/odeint/runge_kutta_step.hpp
+++ b/src/odeint/runge_kutta_step.hpp
@@ -1,178 +1,186 @@
/*-------------------------------------------------------
- Module : odeint
- File : runge_kutta_step
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, 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:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* 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.
* Neither the name of the Princeton University 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 OWNER 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_ODEINT_RUNGEKUTTASTEP_HPP
#define SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP
//! \file runge_kutta_step algorithm for embedded runge kutta
#include "odeint_types.hpp"
#include "butcher_tableau.hpp"
+#include "runge_kutta_step_structs.hpp"
namespace specmicp {
namespace odeint {
-struct EmbeddedRungeKuttaStepOptions
-{
- scalar_t safety;
- scalar_t p_grow;
- scalar_t p_shrink;
- scalar_t errcon;
- scalar_t max_try;
- scalar_t eps;
-
- EmbeddedRungeKuttaStepOptions():
- safety(0.9), p_grow(-0.2), p_shrink(-0.25), errcon(1.89e-4), max_try(20), eps(1e-5)
- {}
-};
-
struct StepLength
{
static constexpr scalar_t not_available = -1.0;
scalar_t test;
scalar_t did;
scalar_t next;
StepLength(scalar_t h_test):
test(h_test), did(not_available), next(not_available)
{}
void get_next() {test=next; did=not_available; next=not_available;}
bool is_next_available() {return next != not_available;}
};
//! \brief Base class for computing a step of the Range Kutta algorithm
template <int dim, class butcher_tableau_t>
class EmbeddedRungeKuttaStep
{
public:
using vec_t = vector_t<dim>;
- EmbeddedRungeKuttaStep(const rhs_f<dim>& rhs,
- const butcher_tableau_t& butcher):
+ EmbeddedRungeKuttaStep(
+ const rhs_f<dim>& rhs,
+ const butcher_tableau_t& butcher):
+ m_rhs(rhs),
+ m_butcher(butcher)
+ {}
+
+ EmbeddedRungeKuttaStep(
+ const rhs_f<dim>& rhs,
+ const butcher_tableau_t& butcher,
+ EmbeddedRungeKuttaStepOptions& options
+ ):
+ m_options(options),
m_rhs(rhs),
m_butcher(butcher)
{}
//! \brief compute a step
//!
//! \param y initial value
//! \param dydx initial value of the right hand site
//! \param x value of the integration variable
//! \param timestep the integration step
//! \param y_out value at the end of the step (y+timestep)
//! \param y_err estimation of the truncation error
void run_step(
const vec_t& y,
const vec_t& dydx,
scalar_t x,
scalar_t timestep,
vec_t& y_out,
vec_t& y_err);
protected:
//! Return the value of the right hand side at x for value y
void get_rhs(
scalar_t x,
const vector_t<dim>& y,
vec_t& dydx
) {return m_rhs(x, y, dydx);}
//! Return the Butcher tableau
butcher_tableau_t& butcher() {return m_butcher;}
EmbeddedRungeKuttaStepOptions& get_options() {return m_options;}
const EmbeddedRungeKuttaStepOptions& get_options() const {return m_options;}
private:
EmbeddedRungeKuttaStepOptions m_options;
rhs_f<dim> m_rhs;
butcher_tableau_t m_butcher;
};
//! Embedded Runge Kutta of fifth order
template <int dim, class butcher_tableau_t>
class EmbeddedRungeKuttaStep45: EmbeddedRungeKuttaStep<dim, butcher_tableau_t>
{
public:
using base = EmbeddedRungeKuttaStep<dim, butcher_tableau_t>;
using vec_t = typename base::vec_t;
using base::butcher;
using base::get_rhs;
using base::get_options;
//check that template argument are consistent
static_assert(butcher_tableau_t::RK_order == 5, "The butcher tableau must correspond to a Runge-Kutta method of order 5");
static_assert(butcher_tableau_t::RK_evaluations == 6 || butcher_tableau_t::RK_evaluations == 7, "The number of functions evaluations must be 6 or 7 for an embedded Runge-Kutta method");
EmbeddedRungeKuttaStep45(const rhs_f<dim>& rhs,
const butcher_tableau_t& butcher):
base(rhs, butcher)
{
}
+
+ EmbeddedRungeKuttaStep45(
+ const rhs_f<dim>& rhs,
+ const butcher_tableau_t& butcher,
+ EmbeddedRungeKuttaStepOptions options
+ ):
+ base(rhs, butcher, options)
+ {
+ }
+
void run_step(const vec_t& y, const vec_t& dydx, scalar_t x, scalar_t timestep, vec_t& y_out, vec_t& y_err);
void rk_step(vector_t<dim> &y,
const vector_t<dim>& dydx,
double& x,
StepLength& stepl);
};
//! typedef to create a Dormant-Prince method
template <int dim>
using DormandPrinceStep = EmbeddedRungeKuttaStep45<dim, ButcherTableauDormandPrince_t>;
template <int dim>
DormandPrinceStep<dim> get_dormand_prince_step(const rhs_f<dim>& rhs) {return DormandPrinceStep<dim>(rhs, butcher_dormand_prince45);}
//! typedef to create a Cash-Karp method
template <int dim>
using CashKarpStep = EmbeddedRungeKuttaStep45<dim, ButcherTableauCashKarp_t>;
template <int dim>
CashKarpStep<dim> get_cash_karp_step(const rhs_f<dim>& rhs) {return CashKarpStep<dim>(rhs, butcher_cash_karp45);}
template <int dim, class T>
T get_rk_step(const rhs_f<dim>& rhs)
{
if (std::is_same<T, CashKarpStep<dim>>::value) return get_cash_karp_step(rhs);
else if (std::is_same<T, DormandPrinceStep<dim>>::value) return get_dormand_prince_step(rhs);
}
} // end namespace odeint
} // end namespace specmicp
#include "runge_kutta_step.inl"
#endif // SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP
diff --git a/src/odeint/runge_kutta_step_structs.hpp b/src/odeint/runge_kutta_step_structs.hpp
new file mode 100644
index 0000000..b2710da
--- /dev/null
+++ b/src/odeint/runge_kutta_step_structs.hpp
@@ -0,0 +1,26 @@
+#ifndef SPECMICP_ODEINT_RUNGEKUTTASTEPSTRUCTS_HPP
+#define SPECMICP_ODEINT_RUNGEKUTTASTEPSTRUCTS_HPP
+
+namespace specmicp {
+namespace odeint {
+
+struct EmbeddedRungeKuttaStepOptions
+{
+ scalar_t safety;
+ scalar_t p_grow;
+ scalar_t p_shrink;
+ scalar_t errcon;
+ scalar_t max_try;
+ scalar_t eps;
+
+ EmbeddedRungeKuttaStepOptions():
+ safety(0.9), p_grow(-0.2), p_shrink(-0.25), errcon(1.89e-4), max_try(20), eps(1e-5)
+ {}
+};
+
+
+} // end namespace odeint
+} // end namespace specmicp
+
+
+#endif// SPECMICP_ODEINT_RUNGEKUTTASTEPSTRUCTS_HPP
diff --git a/src/specmicp/kinetics/kinetic_model.hpp b/src/specmicp/kinetics/kinetic_model.hpp
new file mode 100644
index 0000000..ea5a7cb
--- /dev/null
+++ b/src/specmicp/kinetics/kinetic_model.hpp
@@ -0,0 +1,47 @@
+#ifndef SPECMICP_SPECMICP_KINETICS_KINETICMODEL_HPP
+#define SPECMICP_SPECMICP_KINETICS_KINETICMODEL_HPP
+
+#include "common.hpp"
+#include "kinetic_variables.hpp"
+
+namespace specmicp {
+namespace kinetics {
+
+// \brief Base class for a kinetic model
+class KineticModel
+{
+public:
+ KineticModel(index_t nb_equation):
+ m_neq(nb_equation)
+ {}
+
+ virtual ~KineticModel() {}
+
+ //! \brief Return the number of kinetic minerals included in the problem
+ index_t get_neq() {return m_list_species.size();}
+
+ //! \brief Compute the kinetic rates and store them in dydt
+ virtual void compute_rate(
+ scalar_t t,
+ const Vector& y,
+ KineticVariables& variables,
+ Vector& dydt
+ ) = 0;
+
+
+ //! \brief Index of species corresponding to equation 'id_equation'
+ index_t index_species(index_t id_equation) {return m_list_species[id_equation];}
+
+ //! \brief Iterator over the species vector
+ std::vector<index_t>::iterator species_begin() {return m_list_species.begin();}
+ std::vector<index_t>::iterator species_end() {return m_list_species.end();}
+
+private:
+ std::vector<index_t> list_species;
+};
+
+} // end namespace kinetics
+} // end namespace specmicp
+
+
+#endif //SPECMICP_SPECMICP_KINETICS_KINETICMODEL_HPP
diff --git a/src/specmicp/kinetics/kinetic_system.cpp b/src/specmicp/kinetics/kinetic_system.cpp
new file mode 100644
index 0000000..b04c5c5
--- /dev/null
+++ b/src/specmicp/kinetics/kinetic_system.cpp
@@ -0,0 +1,59 @@
+#include "kinetic_system.hpp"
+
+#include "specmicp/reduced_system_solver.hpp"
+#include "kinetic_model.hpp"
+
+namespace specmicp {
+namespace kinetics {
+
+void KineticSystem::update_total_concentrations(scalar_t dt, const Vector& y)
+{
+ Vector update_minerals = (y - m_variables.moles_minerals())/dt;
+ RawDatabasePtr database = m_variables.get_database();
+ m_variables.total_concentrations() = m_variables.total_concentrations_initial();
+ for (auto it=m_model->species_begin(); it!=m_model->species_end(); ++it)
+ {
+ for (index_t component: database->range_component())
+ {
+ m_variables.total_concentration(component) +=
+ database->nu_mineral_kinetic(*it, component)*update_minerals(it);
+ }
+ }
+}
+
+void KineticSystem::compute_rates(scalar_t t, const Vector& y, Vector& dydt)
+{
+ m_model->compute_rate(t, y, m_variables, dydt);
+}
+
+
+void KineticSystem::compute_equilibrium(ReducedSystemSolverOptions options)
+{
+ ReducedSystemSolver solver;
+ Vector variables;
+ specmicp::micpsolver::MiCPPerformance perf;
+ if (m_variables.equilibrium_state().is_valid())
+ {
+ solver = ReducedSystemSolver(m_variables.get_database(),
+ m_variables.total_concentrations(),
+ m_variables.equilibrium_state(),
+ options);
+ variables = m_variables.equilibrium_state().main_variables();
+ perf = solver.solve(variables);
+ }
+ else
+ {
+ solver = ReducedSystemSolver(m_variables.get_database(),
+ m_variables.total_concentrations(),
+ options);
+ perf = solver.solve(variables, true);
+ }
+ if ((int) perf.return_code < 0)
+ {
+ std::cout << "Failed to solve the system ! Err code " << (int) perf.return_code << std::endl;
+ }
+ m_variables.update_equilibrium_state(solver.get_solution(variables));
+}
+
+} // end namespace kinetics
+} // end namespace specmicp
diff --git a/src/specmicp/kinetics/kinetic_system.hpp b/src/specmicp/kinetics/kinetic_system.hpp
new file mode 100644
index 0000000..e18c60d
--- /dev/null
+++ b/src/specmicp/kinetics/kinetic_system.hpp
@@ -0,0 +1,48 @@
+#ifndef SPECMICP_SPECMICP_KINETICS_KINETICSYSTEM_HPP
+#define SPECMICP_SPECMICP_KINETICS_KINETICSYSTEM_HPP
+
+#include <memory>
+#include "common.hpp"
+#include "database.hpp"
+#include "kinetic_variables.hpp"
+#include "specmicp/reduced_system_solver_structs.hpp"
+
+namespace specmicp {
+
+namespace kinetics {
+
+class KineticModel;
+
+//! \brief Kinetics System,
+//!
+//! Wrapper for a kinetic model
+class KineticSystem
+{
+public:
+
+ KineticSystem(std::shared_ptr<KineticModel> model
+ ):
+ m_model(model)
+ {}
+ //! \brief Compute the kinetics rates to be solved
+ //!
+ //! Use the kinetic model provided by the user
+ void compute_rates(scalar_t x, const Vector& y, Vector& dydx);
+
+ //! \brief Compute the equilibrium state of the solution
+ void compute_equilibrium(ReducedSystemSolverOptions options=ReducedSystemSolverOptions());
+
+ //! Update the total concentrations
+ void update_total_concentrations(scalar_t dt, const Vector& y);
+
+ KineticVariables& variables() {return m_variables;}
+private:
+ KineticVariables m_variables;
+ std::shared_ptr<KineticModel> m_model;
+};
+
+} // end namespace kinetics
+} // end namespace specmicp
+
+
+#endif //SPECMICP_SPECMICP_KINETICS_KINETICSYSTEM_HPP
diff --git a/src/specmicp/kinetics/kinetic_system_solver.cpp b/src/specmicp/kinetics/kinetic_system_solver.cpp
new file mode 100644
index 0000000..e69de29
diff --git a/src/specmicp/kinetics/kinetic_system_solver.hpp b/src/specmicp/kinetics/kinetic_system_solver.hpp
new file mode 100644
index 0000000..beb6af9
--- /dev/null
+++ b/src/specmicp/kinetics/kinetic_system_solver.hpp
@@ -0,0 +1,4 @@
+#ifndef SPECMICP_SPECMICP_KINETICS_KINETICSYSTEMSOLVER_HPP
+#define SPECMICP_SPECMICP_KINETICS_KINETICSYSTEMSOLVER_HPP
+
+#endif //SPECMICP_SPECMICP_KINETICS_KINETICSYSTEMSOLVER_HPP
diff --git a/src/specmicp/kinetics/kinetic_system_solver_structs.hpp b/src/specmicp/kinetics/kinetic_system_solver_structs.hpp
new file mode 100644
index 0000000..1b72cf5
--- /dev/null
+++ b/src/specmicp/kinetics/kinetic_system_solver_structs.hpp
@@ -0,0 +1,20 @@
+#ifndef SPECMICP_SPECMICP_KINETICS_KINETICSYSTEMSOLVERSTRUCTS
+#define SPECMICP_SPECMICP_KINETICS_KINETICSYSTEMSOLVERSTRUCTS
+
+#include "specmicp/reduced_system_solver_structs.hpp"
+#include "odeint/runge_kutta_step_structs.hpp"
+
+namespace specmicp {
+namespace kinetics {
+
+struct KineticSystemSolverOptions
+{
+ odeint::EmbeddedRungeKuttaStepOptions rk_options;
+ ReducedSystemSolverOptions speciation_options;
+};
+
+
+} // end namespace kinetics
+} // end namespace specmicp
+
+#endif // SPECMICP_SPECMICP_KINETICS_KINETICSYSTEMSOLVERSTRUCTS
diff --git a/src/specmicp/kinetics/kinetic_variables.hpp b/src/specmicp/kinetics/kinetic_variables.hpp
new file mode 100644
index 0000000..fc4cf78
--- /dev/null
+++ b/src/specmicp/kinetics/kinetic_variables.hpp
@@ -0,0 +1,75 @@
+#ifndef SPECMICP_SPECMICP_KINETICS_KINETICVARIABLES_HPP
+#define SPECMICP_SPECMICP_KINETICS_KINETICVARIABLES_HPP
+
+#include "common.hpp"
+#include "specmicp/equilibrium_data.hpp"
+
+namespace specmicp {
+namespace kinetics {
+
+//! \brief Variables used in a kinetic computation
+class KineticVariables
+{
+public:
+ KineticVariables(
+ const Vector& total_concentrations,
+ const Vector& mineral_moles):
+ m_mineral_kinetics(mineral_moles),
+ m_total_concentrations_initial(total_concentrations),
+ m_total_concentrations(total_concentrations)
+ {}
+
+ KineticVariables(
+ const Vector& total_concentrations,
+ const Vector& mineral_moles,
+ const EquilibriumState& equilibrium_solution
+ ):
+ KineticVariables(total_concentrations, mineral_moles),
+ m_equilibrium(equilibrium_solution)
+ {}
+
+ //! \brief Return the number of species considered in the kinetics computation
+ index_t get_neq() {return m_mineral_kinetics.rows();}
+ //! \brief Return the raw database
+ RawDatabasePtr get_database() {return m_equilibrium.get_database();}
+
+
+ //! \brief Return mole number of mineral "mineral kinetic"
+ scalar_t moles_mineral(index_t mineral_kinetic) const {return m_mineral_kinetics(mineral_kinetic);}
+ //! \brief Return a const reference to the vector of moles number
+ const Vector& moles_minerals() const {return m_mineral_kinetics;}
+ //! \brief Return a reference to the vector of moles number
+ Vector& moles_minerals() {return m_mineral_kinetics;}
+
+ //! \brief Return the total concentration (at equilibrium) of 'component'
+ scalar_t total_concentration(index_t component) const {return m_total_concentrations(component);}
+ //! \brief Return a const reference to the total concentration (at equilibrium) vector
+ const Vector& total_concentrations() const {return m_total_concentrations;}
+ //! \brief Return a reference to the total concentration at equilibrium vector
+ Vector& total_concentrations() {return m_total_concentrations; }
+
+ //! \brief Return the total concentration (at equilibrium) of 'component'
+ scalar_t total_concentration_initial(index_t component) const {
+ return m_total_concentrations_initial(component);}
+ //! \brief Return a const reference to the total concentration (at equilibrium) vector
+ const Vector& total_concentrations_initial() const {return m_total_concentrations_initial;}
+ //! \brief Return a reference to the total concentration at equilibrium vector
+ Vector& total_concentrations_initial() {return m_total_concentrations_initial; }
+
+ //! \brief Const Reference to the equilibrium state
+ const EquilibriumState& equilibrium_state() const {return m_equilibrium;}
+ //! \brief Reference to the equilibrium State
+ EquilibriumState& equilibrium_state() {return m_equilibrium;}
+ //! \brief Update the equilibrium state
+ void update_equilibrium_state(const EquilibriumState& other) {m_equilibrium = other;}
+private:
+ Vector m_mineral_kinetics; //!< Number of moles of minerals at equilibrium
+ Vector m_total_concentrations_initial; //!< Initial total concentrations at equilibrium
+ Vector m_total_concentrations; //!< total concentration at equilibrium
+ EquilibriumState m_equilibrium; //!< Current equilibrium state
+};
+
+} // end namespace kinetics
+} // end namespace specmicp
+
+#endif // SPECMICP_SPECMICP_KINETICS_KINETICVARIABLES_HPP
diff --git a/src/specmicp/reduced_system_solver.hpp b/src/specmicp/reduced_system_solver.hpp
index b9d019f..0c5f6cb 100644
--- a/src/specmicp/reduced_system_solver.hpp
+++ b/src/specmicp/reduced_system_solver.hpp
@@ -1,145 +1,125 @@
/*-------------------------------------------------------
- Module : specmicp
- File : reduced_system_solver.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, 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:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* 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.
* Neither the name of the Princeton University 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 OWNER 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_SPECMICP_REDUCEDSYSTEMSOLVER_HPP
#define SPECMICP_SPECMICP_REDUCEDSYSTEMSOLVER_HPP
#include "database.hpp"
#include "reduced_system.hpp"
-#include "micpsolver/micpsolver_structs.hpp"
+#include "reduced_system_solver_structs.hpp"
//! \file reduced_system_solver.hpp Solve a reduced system
namespace specmicp {
// forward declaration
class EquilibriumState;
-//! \brief Options for the Equilibrium solver
-//!
-//! Most of the options are contained in the MiCP solver options
-struct ReducedSystemSolverOptions
-{
- micpsolver::NCPfunction ncp_function; //!< the NCP function to use to solve the problem
- micpsolver::MiCPSolverOptions solver_options; //!< Options of the MiCP solver
- bool conservation_water; //!< Solve the conservation of liquid water
- index_t charge_keeper; //!< Species that explicitely conserve the charge
- bool allow_restart; //!< allow the restarting if the problem failed the first part
-
- ReducedSystemSolverOptions():
- ncp_function(micpsolver::NCPfunction::penalizedFB),
- solver_options(),
- conservation_water(true),
- charge_keeper(no_species),
- allow_restart(true)
- {}
-};
-
//! \brief Solver a reduced system
//!
//! Take care of possibly non-existing component in the system
class ReducedSystemSolver
{
public:
//! Default constructor - you need to use another method to initialize it correctly
ReducedSystemSolver() {}
ReducedSystemSolver(RawDatabasePtr data,
const Vector& totconc,
const SimulationBox& simulbox,
ReducedSystemSolverOptions options=ReducedSystemSolverOptions()
);
ReducedSystemSolver(RawDatabasePtr data,
const Vector& totconc,
ReducedSystemSolverOptions options=ReducedSystemSolverOptions()
);
ReducedSystemSolver(RawDatabasePtr data,
const Vector& totconc,
const EquilibriumState& previous_solution,
ReducedSystemSolverOptions options=ReducedSystemSolverOptions()
);
ReducedSystemSolver(RawDatabasePtr data,
const Vector& totconc,
const SimulationBox& simulbox,
const EquilibriumState& previous_solution,
ReducedSystemSolverOptions options=ReducedSystemSolverOptions()
);
//! \brief solve the problem using initial guess x
//!
//! \param[in,out] x in -> initial guess, out -> solution
//! \param init if true, the algorithm guess a starting point
micpsolver::MiCPPerformance solve(Vector& x, bool init=false);
//! \brief Return the system used for the computation
std::shared_ptr<ReducedSystem> get_system() {return m_system;}
//! \brief Return the solution in a manageable form
EquilibriumState get_solution(const Eigen::VectorXd& x);
//! \brief Return a const reference to the options
const ReducedSystemSolverOptions& get_options() const {return m_options;}
//! \brief Return a reference to the options
ReducedSystemSolverOptions& get_options() {return m_options;}
private:
//! \brief set up the true variable vector
void set_true_variable_vector(const Vector& x);
//! \brief set up the true solution vector
//!
//! add zero components
void set_return_vector(Vector& x);
//! \brief solve the problem
micpsolver::MiCPPerformance solve();
ReducedSystemSolverOptions m_options; //! The options
RawDatabasePtr m_data; //! The raw database
std::shared_ptr<ReducedSystem> m_system; //! The system to solve
Vector m_var; //! Copy of the solution vector
};
//! \brief Solve a reduced system, function provided for convenience
inline micpsolver::MiCPPerformance solve_reduced_system(std::shared_ptr<database::DataContainer> data,
Vector& totconc,
Vector& x)
{
return ReducedSystemSolver(data, totconc).solve(x);
}
} // end namespace specmicp
#endif // SPECMICP_SPECMICP_REDUCEDSYSTEMSOLVER_HPP
diff --git a/src/specmicp/reduced_system_solver_structs.hpp b/src/specmicp/reduced_system_solver_structs.hpp
new file mode 100644
index 0000000..b251052
--- /dev/null
+++ b/src/specmicp/reduced_system_solver_structs.hpp
@@ -0,0 +1,32 @@
+#ifndef SPECMICP_SPECMICP_REDUCEDSYSTEMSOLVERSTRUCTS_HPP
+#define SPECMICP_SPECMICP_REDUCEDSYSTEMSOLVERSTRUCTS_HPP
+
+#include "common.hpp"
+#include "micpsolver/micpsolver_structs.hpp"
+
+namespace specmicp {
+
+//! \brief Options for the Equilibrium solver
+//!
+//! Most of the options are contained in the MiCP solver options
+struct ReducedSystemSolverOptions
+{
+ micpsolver::NCPfunction ncp_function; //!< the NCP function to use to solve the problem
+ micpsolver::MiCPSolverOptions solver_options; //!< Options of the MiCP solver
+ bool conservation_water; //!< Solve the conservation of liquid water
+ index_t charge_keeper; //!< Species that explicitely conserve the charge
+ bool allow_restart; //!< allow the restarting if the problem failed the first part
+
+ ReducedSystemSolverOptions():
+ ncp_function(micpsolver::NCPfunction::penalizedFB),
+ solver_options(),
+ conservation_water(true),
+ charge_keeper(no_species),
+ allow_restart(true)
+ {}
+};
+
+
+} // end namespace specmicp
+
+#endif //SPECMICP_SPECMICP_REDUCEDSYSTEMSOLVERSTRUCTS_HPP