diff --git a/CMakeLists.txt b/CMakeLists.txt
index eaf41d2..984d80d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,330 +1,352 @@
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
find_package(CxxTest) # Necessary if you want the tests
# catch for unit testing
include(ExternalProject)
find_package(Git REQUIRED)
ExternalProject_Add(
catch
PREFIX ${CMAKE_BINARY_DIR}/catch
GIT_REPOSITORY https://github.com/philsquared/Catch.git
TIMEOUT 10
UPDATE_COMMAND ${GIT_EXECUTABLE} pull
CONFIGURE_COMMAND ""
BUILD_COMMAND ""
INSTALL_COMMAND ""
LOG_DOWNLOAD ON
)
######################## Compilation flags #################################
if(UNIX)
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -std=c++11")
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_TEST_DIR ${PROJECT_SOURCE_DIR}/tests)
ExternalProject_Get_Property(catch source_dir)
set(CATCH_INCLUDE_DIR ${source_dir}/include CACHE INTERNAL "Path to include folder for Catch")
set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src)
# include
include_directories(${EIGEN3_INCLUDE_DIR})
include_directories(${PROJECT_SOURCE_DIR})
include_directories(${CATCH_INCLUDE_DIR})
include_directories(${PROJECT_TEST_DIR})
################ 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)
add_custom_target(specmic_inc SOURCES
#${SPECMICP_DIR}/.hpp
${SPECMICP_DIR}/simulation_box.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
)
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
)
+# --------- 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
+ #
+ ${DFPMSOLVER_DIR}/sparse_solver.hpp
+ ${DFPMSOLVER_DIR}/sparse_solver_structs.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}/meshes/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/variables.hpp
+ ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters.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/equation_numbering.cpp
${REACTMICP_DIR}/systems/saturated_diffusion/transport_program.cpp
${REACTMICP_DIR}/systems/saturated_diffusion/speciation_program.cpp
)
add_library(reactmicp ${REACTMICP_LIBFILES})
# ------- Utils ------------------
set(UTILS_DIR ${PROJECT_SOURCE_DIR}/utils)
add_custom_target(utils_inc SOURCES
${UTILS_DIR}/log.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 ----------------------
enable_testing(true)
# CXX Test
# ========
if(CXXTEST_FOUND)
set(CXXTEST_USE_PYTHON TRUE)
include_directories(${CXXTEST_INCLUDE_DIR})
enable_testing()
# MiCP Solver
# ------------
CXXTEST_ADD_TEST(test_cond_number
test_cond_number.cpp
${PROJECT_TEST_DIR}/micpsolver/test_cond_number.hpp)
CXXTEST_ADD_TEST(test_ncp_funtion
test_ncp_function.cpp
${PROJECT_TEST_DIR}/micpsolver/test_ncp_function.hpp)
CXXTEST_ADD_TEST(test_micpsolver
test_micpsolver.cpp
${PROJECT_TEST_DIR}/micpsolver/test_micpsolver.hpp)
# Database
# --------
CXXTEST_ADD_TEST(test_data_reader
test_data_reader.cpp
${PROJECT_TEST_DIR}/database/test_data_reader.hpp)
target_link_libraries(test_data_reader specmicp_database jsoncpp)
CXXTEST_ADD_TEST(test_data_selector
test_data_selector.cpp
${PROJECT_TEST_DIR}/database/test_data_selector.hpp)
target_link_libraries(test_data_selector specmicp_database jsoncpp)
# Odeint
# ------
CXXTEST_ADD_TEST(test_butchertableau
test_butchertableau.cpp
${PROJECT_TEST_DIR}/odeint/test_butchertableau.hpp)
CXXTEST_ADD_TEST(test_embeddedrungekutta
test_embeddedrungekutta.cpp
${PROJECT_TEST_DIR}/odeint/test_embeddedrungekutta.hpp)
endif()
# Test 'Perso'
# ============
add_executable(thermocarbo ${PROJECT_TEST_DIR}/specmicp/thermocarbo.cpp)
target_link_libraries(thermocarbo specmicp specmicp_database jsoncpp)
add_executable(carboalu ${PROJECT_TEST_DIR}/specmicp/carboalu.cpp)
target_link_libraries(carboalu specmicp specmicp_database jsoncpp)
add_executable(alkaliactivated ${PROJECT_TEST_DIR}/specmicp/alkaliactivated.cpp)
target_link_libraries(alkaliactivated specmicp specmicp_database jsoncpp)
add_executable(rate_nicoleau ${PROJECT_TEST_DIR}/specmicp/rate_nicoleau.cpp)
target_link_libraries(rate_nicoleau specmicp specmicp_database jsoncpp)
add_executable(kinetics ${PROJECT_TEST_DIR}/specmicp/kinetics.cpp)
target_link_libraries(kinetics specmicp specmicp_database jsoncpp)
set(DATA_TEST ${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3sm.csv
${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st1.csv
${PROJECT_TEST_DIR}/specmicp/data_test/data_nicoleau_c3st2.csv
)
add_custom_target(data_nicoleau SOURCES ${DATA_TEST})
file(INSTALL ${DATA_TEST}
DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data_test/)
# Catch test
#===========
add_custom_target(catch_header SOURCES
${PROJECT_TEST_DIR}/test_database.hpp
)
add_executable(test_reactmicp_diffusion
${PROJECT_TEST_DIR}/reactmip/systems/diffusion/test_reactmicp_diffusion.cpp
${PROJECT_TEST_DIR}/reactmip/systems/diffusion/diffusion.cpp
${PROJECT_TEST_DIR}/reactmip/systems/diffusion/numbering.cpp
${PROJECT_TEST_DIR}/reactmip/systems/diffusion/secondary.cpp
${PROJECT_TEST_DIR}/reactmip/systems/diffusion/solving.cpp
${PROJECT_TEST_DIR}/reactmip/systems/diffusion/utils.cpp
)
target_link_libraries(test_reactmicp_diffusion reactmicp specmicp specmicp_database jsoncpp)
add_executable(test_reactmicp_saturated_diffusion
${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/test_reactmicp_saturated_diffusion.cpp
${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/variables.cpp
+ ${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/solver.cpp
+ ${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/transport_program.cpp
${PROJECT_TEST_DIR}/reactmip/systems/saturated_diffusion/utils.cpp
)
target_link_libraries(test_reactmicp_saturated_diffusion reactmicp specmicp specmicp_database jsoncpp)
diff --git a/src/dfpmsolver/dfpm_program.hpp b/src/dfpmsolver/dfpm_program.hpp
new file mode 100644
index 0000000..09d9436
--- /dev/null
+++ b/src/dfpmsolver/dfpm_program.hpp
@@ -0,0 +1,27 @@
+#ifndef SPECMICP_DFPMSOLVER_DFPMPROGRAM_HPP
+#define SPECMICP_DFPMSOLVER_DFPMPROGRAM_HPP
+
+namespace specmicp {
+namespace dfpmsolver {
+
+//! \brief Base class for a program
+template <class Derived>
+class DFPMProgram
+{
+public:
+ //! \brief Return a pointer to the true object
+ Derived* derived() {static_cast<Derived*>(this);}
+
+ //! \brief Return the number of equations
+ int get_neq() {return derived()->get_neq();}
+
+ //! \brief Return the number of degrees of freedom per node
+ int get_ndf() {return derived()->get_ndf();}
+
+
+};
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
+
+#endif // SPECMICP_DFPMSOLVER_DFPMPROGRAM_HPP
diff --git a/src/dfpmsolver/driver.hpp b/src/dfpmsolver/driver.hpp
new file mode 100644
index 0000000..52d0a30
--- /dev/null
+++ b/src/dfpmsolver/driver.hpp
@@ -0,0 +1,76 @@
+#ifndef SPECMICP_DFPMSOLVER_DRIVER_HPP
+#define SPECMICP_DFPMSOLVER_DRIVER_HPP
+
+namespace specmicp {
+namespace dfpmsolver {
+
+//! \brief Base class for a driver
+//!
+//! Options should be a subclass of DriverOptions
+//! Performance should be a subclass of DriverPerformance
+//!
+template <class Program, class Options, class Performance>
+class Driver
+{
+public:
+ Driver(Program& the_program):
+ m_program(the_program),
+ m_options()
+ {}
+
+ Driver(Program& the_program,
+ Options some_options):
+ m_program(the_program),
+ m_options(some_options)
+ {}
+
+ // basic program management
+ // ------------------------
+
+ //! Return the program
+ Program& program() {return m_program;}
+
+ //! Return the number of equations
+ int get_neq() {return m_program.get_neq();}
+
+ // common process
+ // --------------
+ //! \brief rescale the update if needed
+ //!
+ //! Return the step length
+ double is_step_too_long(Eigen::VectorXd& update);
+
+
+ // options
+ // -------
+
+ //! \brief Return a read/write reference to the options
+ Options& get_options() {return m_options;}
+ //! \brief Return a read-only reference to the options
+ const Options& get_options() const {return m_options;}
+
+ // performance
+ // -----------
+
+ //! \brief Return a const reference to the performance
+ const Performance& get_perfs() const {return m_performance;}
+
+protected:
+ //! \brief Read/write access to the performance
+ Performance& get_perfs() {return m_performance;}
+
+private:
+ Program& m_program; //!< The program
+ Performance m_performance; //!< The performance
+ Options m_options; //!< The options
+};
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
+
+// implementation
+// ==============
+
+#include "driver.inl"
+
+#endif // SPECMICP_DFPMSOLVER_DRIVER_HPP
diff --git a/src/dfpmsolver/driver.inl b/src/dfpmsolver/driver.inl
new file mode 100644
index 0000000..778178e
--- /dev/null
+++ b/src/dfpmsolver/driver.inl
@@ -0,0 +1,22 @@
+#include "driver.hpp" // syntaxic coloration only
+
+namespace specmicp {
+namespace dfpmsolver {
+
+
+template <class Program, class Options, class Performance>
+double Driver<Program, Options, Performance>::is_step_too_long(Eigen::VectorXd& update)
+{
+ double steplength = update.norm();
+ if (steplength > get_options().maximum_step_length)
+ {
+ get_perfs().maximum_step_taken = true;
+ update = get_options().maximum_step_length / steplength * update;
+ steplength = get_options().maximum_step_length;
+ }
+ return steplength;
+}
+
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
diff --git a/src/dfpmsolver/driver_structs.hpp b/src/dfpmsolver/driver_structs.hpp
new file mode 100644
index 0000000..6fb407e
--- /dev/null
+++ b/src/dfpmsolver/driver_structs.hpp
@@ -0,0 +1,49 @@
+#ifndef SPECMICP_DFPMSOLVER_DRIVERSTRUCTS_HPP
+#define SPECMICP_DFPMSOLVER_DRIVERSTRUCTS_HPP
+
+#include "sparse_solver_structs.hpp"
+
+namespace specmicp {
+namespace dfpmsolver {
+
+
+//! \brief Options of a driver
+//!
+struct DriverOptions
+{
+ double residuals_tolerance; //!< Tolerance for the residual
+ double step_tolerance; //!< Tolerance for the minimum step length
+ double threshold_stationary_point; //!< if ||R||>threshold, the point is classified as stationary
+ int maximum_iterations; //!< Maximum iterations allowed
+ double maximum_step_length; //!< Maximum step length allowed
+ int max_iterations_at_max_length; //!< Maximum number of iterations at maximum step length
+ double coeff_accept_newton_step; //!< Accept Newton step if enough progress is made
+ SparseSolver sparse_solver; //!< The sparse solver to use
+
+ DriverOptions():
+ residuals_tolerance(1e-6),
+ step_tolerance(1e-8),
+ threshold_stationary_point(1e-4),
+ maximum_iterations(200),
+ maximum_step_length(1e3),
+ max_iterations_at_max_length(50),
+ coeff_accept_newton_step(0.9),
+ sparse_solver(SparseSolver::SparseQR)
+ {}
+};
+
+//! \brief Performance of a driver
+struct DriverPerformance
+{
+ int nb_call_residuals;
+ int nb_call_jacobian;
+ int nb_iterations;
+ int nb_consecutive_max_step_taken;
+ int nb_max_step_taken;
+ bool maximum_step_taken;
+};
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
+
+#endif // SPECMICP_DFPMSOLVER_DRIVERSTRUCTS_HPP
diff --git a/src/dfpmsolver/parabolic_driver.hpp b/src/dfpmsolver/parabolic_driver.hpp
new file mode 100644
index 0000000..9ad8335
--- /dev/null
+++ b/src/dfpmsolver/parabolic_driver.hpp
@@ -0,0 +1,100 @@
+#ifndef SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
+#define SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+#include "driver.hpp"
+#include "parabolic_structs.hpp"
+
+namespace specmicp {
+namespace dfpmsolver {
+
+//! \brief The parabolic driver
+//!
+//! Parabolic driver for finite element programs
+template <class Program>
+class ParabolicDriver: public Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>
+{
+public:
+ using base=Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>;
+ using base::program;
+ using base::get_neq;
+ using base::get_options;
+ using base::get_perfs;
+
+
+ ParabolicDriver(Program& the_program):
+ Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>(the_program)
+ {}
+
+ //! \brief Solve a timestep of length dt
+ ParabolicDriverReturnCode solve_timestep(double dt, Eigen::VectorXd& displacement);
+ //! \brief Restart the current timestep
+ ParabolicDriverReturnCode restart_timestep(Eigen::VectorXd& displacement);
+ //! \brief Check if the solution has been found
+ ParabolicDriverReturnCode check_convergence(double norm_update);
+
+ //! \brief Compute the residuals
+ void compute_residuals(const Eigen::VectorXd& displacement,
+ const Eigen::VectorXd& velocity,
+ Eigen::VectorXd& residual)
+ {
+ program().compute_residuals(displacement, velocity, residual);
+ get_perfs().nb_call_residuals += 1;
+ }
+ //! \brief Compute the jacobian
+ void compute_jacobian(Eigen::VectorXd& displacement,
+ Eigen::VectorXd& velocity,
+ Eigen::SparseMatrix<double>& jacobian
+ )
+ {
+ program().compute_jacobian(displacement, velocity, jacobian, get_options().alpha*m_current_dt);
+ get_perfs().nb_call_jacobian += 1;
+ }
+ //! \brief Return the norm of the current residuals
+ double norm_residuals() {return m_residuals.norm();}
+
+ Eigen::VectorXd& get_velocity() {return m_velocity;}
+
+private:
+ void update_variable(
+ const Eigen::VectorXd& update,
+ double lambda,
+ Eigen::VectorXd& displacement
+ );
+ //! \brief Set the predictor
+ void set_predictor(Eigen::VectorXd& displacement);
+ //! \brief solve the linear system
+ int solve_linear_system(Eigen::VectorXd& update);
+ //! \brief perform the linesearch
+ int linesearch(Eigen::VectorXd& update,
+ Eigen::VectorXd& displacements);
+
+ double compute_residuals_linesearch(Eigen::VectorXd& update,
+ double lambda,
+ Eigen::VectorXd& displacement
+ );
+
+ //! \brief Initialize the computation
+ void initialize_timestep(double dt, Eigen::VectorXd& displacement);
+
+ double m_norm_0;
+ double m_current_dt;
+ Eigen::VectorXd m_gradient;
+ Eigen::VectorXd m_velocity;
+ Eigen::VectorXd m_residuals;
+ Eigen::VectorXd m_predictor;
+ Eigen::SparseMatrix<double> m_jacobian;
+
+};
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
+
+
+// implementation
+// ==============
+#include "parabolic_driver.inl"
+
+#endif // SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
diff --git a/src/dfpmsolver/parabolic_driver.inl b/src/dfpmsolver/parabolic_driver.inl
new file mode 100644
index 0000000..db4c060
--- /dev/null
+++ b/src/dfpmsolver/parabolic_driver.inl
@@ -0,0 +1,245 @@
+#include "parabolic_driver.hpp" // for syntaxic coloration only
+#include "utils/log.hpp"
+#include "sparse_solver.hpp"
+
+namespace specmicp {
+namespace dfpmsolver {
+
+template <class Program>
+ParabolicDriverReturnCode ParabolicDriver<Program>::solve_timestep(double dt, Eigen::VectorXd& displacement)
+{
+ initialize_timestep(dt, displacement);
+ compute_residuals(displacement, m_velocity, m_residuals);
+ m_norm_0 = norm_residuals();
+ return restart_timestep(displacement);
+}
+
+template <class Program>
+void ParabolicDriver<Program>::initialize_timestep(double dt, Eigen::VectorXd& displacement)
+{
+ m_residuals = Eigen::VectorXd::Zero(get_neq());
+ m_current_dt = dt;
+ set_predictor(displacement);
+ m_velocity.setZero();
+ program().apply_bc(dt, displacement, m_velocity);
+}
+
+template <class Program>
+void ParabolicDriver<Program>::set_predictor(Eigen::VectorXd& displacement)
+{
+ if (get_options().alpha < 1)
+ m_predictor = displacement + (1-get_options().alpha)*m_current_dt*m_velocity;
+ else
+ m_predictor = displacement;
+}
+
+template <class Program>
+ParabolicDriverReturnCode ParabolicDriver<Program>::restart_timestep(Eigen::VectorXd& displacement)
+{
+ ParabolicDriverReturnCode return_code = ParabolicDriverReturnCode::NotConvergedYet;
+ while (return_code == ParabolicDriverReturnCode::NotConvergedYet)
+ {
+ Eigen::VectorXd update(get_neq());
+ compute_residuals(displacement, m_velocity, m_residuals);
+ return_code = check_convergence(-1.0);
+ if (return_code != ParabolicDriverReturnCode::NotConvergedYet) break;
+ compute_jacobian(displacement, m_velocity, m_jacobian);
+ m_gradient = m_jacobian.transpose()*m_residuals;
+ int retcode = solve_linear_system(update);
+ if (retcode <1)
+ {
+ ERROR << "Error when solving linear system : " << retcode << std::endl;
+ return_code = ParabolicDriverReturnCode::ErrorLinearSystem;
+ }
+ linesearch(update, displacement);
+ }
+ return return_code;
+
+}
+
+template <class Program>
+ParabolicDriverReturnCode ParabolicDriver<Program>::check_convergence(double norm_update)
+{
+ ParabolicDriverReturnCode termcode = ParabolicDriverReturnCode::NotConvergedYet;
+ const int nb_iterations = get_perfs().nb_iterations;
+ const double norm_residuals = m_residuals.norm();
+ WARNING << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0;
+ if (norm_residuals/m_norm_0 < get_options().residuals_tolerance)
+ {
+ termcode = ParabolicDriverReturnCode::ResidualMinimized;
+ }
+ else if (nb_iterations > 0 and norm_update >0 and norm_update < get_options().step_tolerance)
+ {
+ if (norm_residuals > get_options().threshold_stationary_point)
+ {
+ ERROR << "Stationary point detected !";
+ termcode = ParabolicDriverReturnCode::StationaryPoint;
+ }
+ WARNING << "Error is minimized - may indicate a stationnary point";
+ termcode = ParabolicDriverReturnCode::ErrorMinimized;
+
+ }
+ else if (nb_iterations > get_options().maximum_iterations)
+ {
+ ERROR << "Maximum number of iteration reached (" << get_options().maximum_iterations << ")";
+ termcode = ParabolicDriverReturnCode::MaxIterations;
+ }
+ else if (get_perfs().maximum_step_taken)
+ {
+ get_perfs().nb_consecutive_max_step_taken += 1;
+ get_perfs().nb_max_step_taken += 1;
+ if (get_perfs().nb_consecutive_max_step_taken == get_options().max_iterations_at_max_length) {
+ ERROR << "Divergence detected - Maximum step length taken two many times";
+ termcode = ParabolicDriverReturnCode::MaxStepTakenTooManyTimes;
+ }
+ }
+ else
+ {
+ get_perfs().nb_consecutive_max_step_taken = 0;
+ }
+ return termcode;
+}
+
+template <class Program>
+int ParabolicDriver<Program>::solve_linear_system(Eigen::VectorXd& update)
+{
+ return solve_sparse_linear_system(m_residuals, m_jacobian, update, get_options().sparse_solver);
+}
+
+template <class Program>
+double ParabolicDriver<Program>::compute_residuals_linesearch(
+ Eigen::VectorXd& update,
+ double lambda,
+ Eigen::VectorXd& displacement
+ )
+{
+ Eigen::VectorXd velocity(m_velocity);
+ Eigen::VectorXd residual = Eigen::VectorXd::Zero(get_neq());
+ program().update_solution(update, lambda, get_options().alpha*m_current_dt, m_predictor, displacement, velocity);
+ compute_residuals(displacement, velocity, residual);
+ return 0.5*residual.squaredNorm();
+}
+
+template <class Program>
+void ParabolicDriver<Program>::update_variable(
+ const Eigen::VectorXd& update,
+ double lambda,
+ Eigen::VectorXd& displacement
+ )
+{
+ program().update_solution(update, lambda, get_options().alpha*m_current_dt,
+ m_predictor, displacement, m_velocity);
+}
+
+template <class Program>
+int ParabolicDriver<Program>::linesearch(
+ Eigen::VectorXd& update,
+ Eigen::VectorXd& displacements)
+{
+ // References
+ // ----------
+ // - Algo A6.3.1 : Dennis and Schnabel (1983)
+ // - Nocedal & Wrigth (2006)
+ DEBUG << "Linesearch";
+ Eigen::VectorXd xp(displacements);
+ double fcp;
+ get_perfs().maximum_step_taken = false;
+ int retcode = 2;
+ const double alpha = 1e-4;
+ double newtlen = base::is_step_too_long(update);
+ double init_slope = m_gradient.dot(update);
+ double rellength = update.cwiseAbs().maxCoeff();
+ const double minlambda = get_options().step_tolerance / rellength;
+ double lambda = 1.0;
+ double lambda_prev = lambda;
+
+ double merit_value = 0.5*m_residuals.squaredNorm();
+ // new residual
+ fcp = compute_residuals_linesearch(update, lambda, xp);
+ // Skip linesearch if enough progress is done
+ // ------------------------------------------
+ if (fcp < get_options().coeff_accept_newton_step *merit_value)
+ {
+ WARNING << "Skip linesearch ";
+ update_variable(update, lambda, displacements);
+ return 0;
+ }
+ // The linesearch
+ // --------------
+ double fc = merit_value;
+ double fcp_prev;
+ int cnt = 0;
+ do
+ {
+ WARNING << "cnt : " <<cnt << " - lambda : " << lambda << " # " << fc << " # " << fcp << " # " << init_slope;
+ if (fcp <= fc + alpha*lambda*init_slope) //pg760 Fachinei2003
+ {
+ retcode = 0;
+ if (lambda ==1 and (newtlen > 0.99 * get_options().maximum_step_length)) {
+ get_perfs().maximum_step_taken = true;
+ }
+ break;
+ }
+ else if (lambda < minlambda)
+ {
+ retcode = 1;
+ break;
+ }
+ else
+ {
+ //WARNING << "denom : " << fcp - fc -init_slope;
+ // Select a new step length
+ // - - - - - - - - - - - -
+ double lambdatmp;
+ if (cnt == 0) { // only a quadratic at the first
+ lambdatmp = - init_slope / (2*(fcp - fc -init_slope));
+ }
+ else
+ {
+ const double factor = 1 /(lambda - lambda_prev);
+ const double x1 = fcp - fc - lambda*init_slope;
+ const double x2 = fcp_prev - fc - lambda_prev*init_slope;
+
+ const double a = factor * ( x1/(lambda*lambda) - x2/(lambda_prev*lambda_prev));
+ const double b = factor * ( -x1*lambda_prev/(lambda*lambda) + x2*lambda/(lambda_prev*lambda_prev));
+
+ if (a == 0)
+ { // cubic interpolation is in fact a quadratic interpolation
+ lambdatmp = - init_slope/(2*b);
+ }
+ else
+ {
+ const double disc = b*b-3*a*init_slope;
+ lambdatmp = (-b+std::sqrt(disc))/(3*a);
+ }
+ if (lambdatmp > 0.5*lambda ) lambdatmp = 0.5*lambda;
+ }
+ //WARNING << "lambdatmp : " << lambdatmp;
+ lambda_prev = lambda;
+ fcp_prev = fcp;
+ if (lambdatmp < 0.1*lambda) {
+ lambda = 0.1 * lambda;
+ } else {
+ lambda = lambdatmp;
+ }
+ }
+ if (not std::isfinite(lambda))
+ {
+ ERROR << "Lambda is non finite - we stop";
+ throw std::runtime_error("Lambda is nonfinite : "+std::to_string(lambda));
+ }
+ fcp = compute_residuals_linesearch(update, lambda, xp);
+ //xp.velocity = x.velocity;
+ ++cnt;
+ } while(retcode == 2 and cnt < 100);
+ //WARNING << "Lambda : " << lambda << " - iterations : " << cnt;
+ if (cnt == 100)
+ {
+ ERROR << "Too much linesearch iterations ! We stop";
+ }
+ update_variable(update, lambda, displacements);
+ return retcode;
+}
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
diff --git a/src/dfpmsolver/parabolic_program.hpp b/src/dfpmsolver/parabolic_program.hpp
new file mode 100644
index 0000000..5a4d783
--- /dev/null
+++ b/src/dfpmsolver/parabolic_program.hpp
@@ -0,0 +1,63 @@
+#ifndef SPECMICP_DFPMSOLVER_PARABOLICPROGRAM_HPP
+#define SPECMICP_DFPMSOLVER_PARABOLICPROGRAM_HPP
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include "dfpm_program.hpp"
+
+namespace specmicp {
+namespace dfpmsolver {
+
+//! Concept class for a program
+template <class Derived>
+class ParabolicProgram: DFPMProgram<ParabolicProgram<Derived>>
+{
+public:
+ Derived* derived() {static_cast<Derived*>(this);}
+
+ //! \brief Return the number of equations
+ int get_neq() {return derived()->get_neq();}
+ //! \brief Return the number of degrees of freedom per node
+ int get_ndf() {return derived()->get_ndf();}
+
+ //! \brief Compute the residuals
+ void compute_residuals(const Eigen::VectorXd& displacement,
+ const Eigen::VectorXd& velocity,
+ Eigen::VectorXd& residual
+ )
+ {
+ derived()->compute_residuals(displacement, velocity, residual);
+ }
+ //! \brief Compute the jacobian
+ void compute_jacobian(Eigen::VectorXd& displacement,
+ Eigen::VectorXd& velocity,
+ Eigen::SparseMatrix<double>& jacobian,
+ double alphadt
+ )
+ {
+ derived()->compute_jacobian(displacement, velocity, jacobian);
+ }
+ //! \brief Update the solutions
+ void update_solution(const Eigen::VectorXd& update,
+ double lambda,
+ double alpha_dt,
+ Eigen::VectorXd& predictor,
+ Eigen::VectorXd& displacement,
+ Eigen::VectorXd& velocity)
+ {
+ derived()->update_solution(update, lambda, alpha_dt, predictor, displacement, velocity);
+ }
+ //! \brief Apply boundary conditions to the velocity vector
+ //!
+ //! by default do nothing.
+ void apply_bc(double dt,
+ const Eigen::VectorXd& displacement,
+ Eigen::VectorXd velocity)
+ {}
+
+};
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
+
+#endif // SPECMICP_DFPMSOLVER_PARABOLICPROGRAM_HPP
diff --git a/src/dfpmsolver/parabolic_structs.hpp b/src/dfpmsolver/parabolic_structs.hpp
new file mode 100644
index 0000000..19f1085
--- /dev/null
+++ b/src/dfpmsolver/parabolic_structs.hpp
@@ -0,0 +1,42 @@
+#ifndef SPECMICP_DFPMSOLVER_PARABOLCSTRUCTS_HPP
+#define SPECMICP_DFPMSOLVER_PARABOLICSTRUCTS_HPP
+
+#include "driver_structs.hpp"
+
+namespace specmicp {
+namespace dfpmsolver {
+
+
+//! \brief Options of a parabolic driver
+//!
+struct ParabolicDriverOptions: public DriverOptions
+{
+ double alpha; //!< Implicit/Cranck-Nicholson parameter
+
+ ParabolicDriverOptions():
+ DriverOptions(),
+ alpha(1.0)
+ {}
+};
+
+//! \brief Performance of the parabolic driver
+struct ParabolicDriverPerformance: public DriverPerformance
+{
+};
+
+//! \brief Return codes
+enum class ParabolicDriverReturnCode
+{
+ MaxIterations = -4, //!< Maximum number of iterations reached
+ MaxStepTakenTooManyTimes = -3, //!< Maximum step taken too many times (divergence)
+ ErrorLinearSystem = -2, //!< Error when solving the linear system
+ StationaryPoint = -1, //!< Stationnary points are detected
+ NotConvergedYet = 0, //!< Problem is not converged
+ ResidualMinimized = 1, //!< The residual is minimized (Success)
+ ErrorMinimized = 2 //!< Error is minimized (may be good)
+};
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
+
+#endif // SPECMICP_DFPMSOLVER_PARABOLICSTRUCTS_HPP
diff --git a/src/dfpmsolver/sparse_solver.hpp b/src/dfpmsolver/sparse_solver.hpp
new file mode 100644
index 0000000..e6d8a1d
--- /dev/null
+++ b/src/dfpmsolver/sparse_solver.hpp
@@ -0,0 +1,97 @@
+#ifndef SPECMICP_DFPMSOLVER_SPARSESOLVER_HPP
+#define SPECMICP_DFPMSOLVER_SPARSESOLVER_HPP
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+#include "sparse_solver_structs.hpp"
+
+namespace specmicp {
+namespace dfpmsolver {
+
+// header
+// ======
+
+//! \brief Solve a sparse linear system
+template <class Vector, class Matrix>
+int solve_sparse_linear_system(Vector& residual, Matrix& jacobian, Vector& update, SparseSolver solver);
+
+
+//! \brief Solve a sparse linear system using the QR decomposition
+template <class Vector, class Matrix>
+int solve_linear_system_sparseqr(Vector& residual, Matrix& jacobian, Vector& update);
+
+
+//! \brief Solve a sparse linear system using the LU decomposition
+template <class Vector, class Matrix>
+int solve_linear_system_sparselu(Vector& residual, Matrix& jacobian, Vector& update);
+
+
+//! \brief Solve a sparse linear system using the BiCGSTAB algorithm
+template <class Vector, class Matrix>
+int solve_linear_system_bicgstab(Vector& residual, Matrix& jacobian, Vector& update);
+
+// Implementation
+// ==============
+
+template <class Vector, class Matrix>
+int solve_sparse_linear_system(Vector& residual, Matrix& jacobian, Vector& update, SparseSolver solver)
+{
+ if (solver == SparseSolver::SparseQR) {
+ return solve_linear_system_sparseqr(residual, jacobian, update);
+ } else if (solver == SparseSolver::SparseLU) {
+ return solve_linear_system_sparselu(residual, jacobian, update);
+ }
+ else {
+ return solve_linear_system_bicgstab(residual, jacobian, update);
+ }
+}
+
+template <class Vector, class Matrix>
+int solve_linear_system_sparseqr(Vector& residual, Matrix& jacobian, Vector& update)
+{
+ Eigen::SparseQR<Matrix, Eigen::COLAMDOrdering<int>> solver(jacobian);
+ if (solver.info() != Eigen::Success)
+ {
+ return -2;
+ }
+ update = solver.solve(-residual);
+ if (solver.info() != Eigen::Success)
+ {
+ return -1;
+ }
+ return 1;
+}
+
+template <class Vector, class Matrix>
+int solve_linear_system_sparselu(Vector& residual, Matrix& jacobian, Vector& update)
+{
+ Eigen::SparseLU<Matrix> solver(jacobian);
+ if (solver.info() != Eigen::Success)
+ {
+ return -2;
+ }
+ update = solver.solve(-residual);
+ if (solver.info() != Eigen::Success)
+ {
+ return -1;
+ }
+ return 1;
+}
+
+template <class Vector, class Matrix>
+int solve_linear_system_bicgstab(Vector& residual, Matrix& jacobian, Vector& update)
+{
+ Eigen::BiCGSTAB<Matrix, Eigen::IncompleteLUT<typename Matrix::Scalar>> solver(jacobian);
+ update = solver.solve(-residual);
+ if (solver.info() != Eigen::Success)
+ {
+ return -1;
+ }
+ return 1;
+}
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
+
+#endif //SPECMICP_DFPMSOLVER_SPARSESOLVER_HPP
diff --git a/src/dfpmsolver/sparse_solver_structs.hpp b/src/dfpmsolver/sparse_solver_structs.hpp
new file mode 100644
index 0000000..50cec47
--- /dev/null
+++ b/src/dfpmsolver/sparse_solver_structs.hpp
@@ -0,0 +1,20 @@
+#ifndef SPECMICP_DFPMSOLVER_SPARSESOLVERSTRUCT_HPP
+#define SPECMICP_DFPMSOLVER_SPARSESOLVERSTRUCT_HPP
+
+
+namespace specmicp {
+namespace dfpmsolver {
+
+//! \brief Available sparse solver
+enum class SparseSolver
+{
+ SparseLU, //!< Eigen Sparse LU solver
+ SparseQR, //!< Eigen Sparse QR solver
+ BiCGSTAB //!< Eigne BiCGSTAB solver
+};
+
+
+} // end namespace dfpmsolver
+} // end namespace specmicp
+
+#endif //SPECMICP_DFPMSOLVER_SPARSESOLVERSTRUCT_HPP
diff --git a/src/reactmicp/systems/saturated_diffusion/transport_parameters.hpp b/src/reactmicp/systems/saturated_diffusion/transport_parameters.hpp
new file mode 100644
index 0000000..0a61019
--- /dev/null
+++ b/src/reactmicp/systems/saturated_diffusion/transport_parameters.hpp
@@ -0,0 +1,38 @@
+#ifndef SPECMICP_REACTMICP_SATURATEDDIFFUSION_TRANSPORTPARAMETERS_HPP
+#define SPECMICP_REACTMICP_SATURATEDDIFFUSION_TRANSPORTPARAMETERS_HPP
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace siasaturated {
+
+
+//! \brief Parameters for the diffusion system
+class SaturatedDiffusionTransportParameters
+{
+public:
+
+ SaturatedDiffusionTransportParameters(double diffusion_coefficient, double porosity):
+ m_diff(diffusion_coefficient), m_poro(porosity)
+ {}
+
+ //! \brief Density of water (kg/m^3)
+ double density_water() {return 1e3;}
+
+ //! \brief Diffusion coefficient (element by element) (m^2/s)
+ double diffusion_coefficient(int element) {return m_diff;}
+
+ //! \brief Porosity (at a node (function of composition))
+ double porosity(int node) {return m_poro;}
+
+private:
+ double m_diff;
+ double m_poro;
+};
+
+} // end namespace siasaturated
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_SATURATEDDIFFUSION_TRANSPORTPARAMETERS_HPP
diff --git a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp
index 4810681..2f6689d 100644
--- a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp
+++ b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp
@@ -1,7 +1,205 @@
#include "transport_program.hpp"
+#include <iostream>
+
+#define EPS_J 1e-8
namespace specmicp {
namespace reactmicp {
+namespace systems {
+namespace siasaturated {
+
+void SaturatedDiffusionProgram::number_equations(
+ std::vector<int> list_fixed_nodes
+ )
+{
+ m_ideq = Eigen::VectorXd::Zero(get_ndf()*m_mesh->nnodes());
+ for (auto it=list_fixed_nodes.begin(); it<list_fixed_nodes.end(); ++it)
+ {
+ for (int component: m_data->range_aqueous_component())
+ {
+ m_ideq(get_dof_component(*it, component)) = no_equation;
+ }
+ }
+
+ int neq = 0;
+ for (int dof=0; dof<m_ideq.rows(); ++dof)
+ {
+ if (m_ideq(dof) != no_equation)
+ {
+ m_ideq(dof) = neq;
+ ++neq;
+ }
+ }
+ m_neq = neq;
+}
+
+void SaturatedDiffusionProgram::compute_residuals(
+ const Eigen::VectorXd& displacement,
+ const Eigen::VectorXd& velocity,
+ Eigen::VectorXd& residual
+ )
+{
+ residual = Eigen::VectorXd::Zero(get_neq());
+ for (ind_t element: m_mesh->range_elements())
+ {
+ element_residual_transport(element, displacement, velocity, residual);
+ }
+}
+
+
+void SaturatedDiffusionProgram::element_residual_transport(
+ ind_t element,
+ const Eigen::VectorXd& displacement,
+ const Eigen::VectorXd& velocity,
+ Eigen::VectorXd& residual)
+{
+ Eigen::VectorXd element_residual(2);
+ for (ind_t component: m_data->range_aqueous_component())
+ {
+ element_residual_transport_component(element, component, displacement, velocity, element_residual);
+
+ for (int en: m_mesh->range_nen())
+ {
+ const ind_t node = m_mesh->get_node(element, en);
+ const ind_t id = m_ideq(get_dof_component(node, component));
+ if (id != no_equation) {residual(id) += element_residual(en);}
+ }
+ }
+}
+
+void SaturatedDiffusionProgram::element_residual_transport_component(
+ ind_t element,
+ int component,
+ const Eigen::VectorXd& displacement,
+ const Eigen::VectorXd& velocity,
+ Vector& element_residual)
+{
+
+ Eigen::Matrix2d mass, jacob;
+ Eigen::Vector2d evelocity, edisplacement;
+ double mass_coeff = -(m_param->density_water() *
+ m_param->porosity(0) *
+ m_mesh->crosssection()*m_mesh->get_dx(element)/2);
+ mass << 1, 0, 0, 1;
+ mass *= mass_coeff;
+ double flux_coeff = -( m_mesh->crosssection() / m_mesh->get_dx(element)
+ * m_param->porosity(0)
+ * m_param->diffusion_coefficient(element)
+ * m_param->density_water()
+ );
+ jacob << 1, -1, -1, 1;
+ jacob *= flux_coeff;
+ evelocity << velocity(get_dof_component(m_mesh->get_node(element, 0), component)),
+ velocity(get_dof_component(m_mesh->get_node(element, 1), component));
+
+ edisplacement << displacement(get_dof_component(m_mesh->get_node(element, 0), component)),
+ displacement(get_dof_component(m_mesh->get_node(element, 1), component));
+
+ element_residual = mass*evelocity+jacob*edisplacement;
+
+ for (int en: m_mesh->range_nen())
+ {
+ element_residual(en) += (m_mesh->crosssection()*m_mesh->get_dx(element)/2
+ *external_flux(m_mesh->get_node(element, en), component));
+ }
+}
+
+
+void SaturatedDiffusionProgram::compute_jacobian(
+ Eigen::VectorXd& displacement,
+ Eigen::VectorXd& velocity,
+ Eigen::SparseMatrix<double>& jacobian,
+ double alphadt
+ )
+{
+
+ list_triplet_t jacob;
+ const ind_t ncomp = m_data->nb_component-1;
+ const ind_t estimation = m_mesh->nnodes()*(ncomp*m_mesh->nen);
+ jacob.reserve(estimation);
+ for (ind_t element: m_mesh->range_elements())
+ {
+ element_jacobian_transport(element, displacement, velocity, jacob, alphadt);
+ }
+ jacobian = SparseMatrix(get_neq(), get_neq());
+ jacobian.setFromTriplets(jacob.begin(), jacob.end());
+}
+
+// Transport
+// =========
+
+void SaturatedDiffusionProgram::element_jacobian_transport(
+ ind_t element,
+ Eigen::VectorXd& displacement,
+ Eigen::VectorXd& velocity,
+ list_triplet_t& jacobian,
+ double alphadt)
+{
+ for (int component: m_data->range_aqueous_component())
+ {
+ Eigen::VectorXd element_residual_orig(Eigen::VectorXd::Zero(2));
+ element_residual_transport_component(element, component, displacement, velocity, element_residual_orig);
+
+ for (int en: m_mesh->range_nen())
+ {
+ Eigen::VectorXd element_residual(Eigen::VectorXd::Zero(2));
+
+ const ind_t node = m_mesh->get_node(element, en);
+ const int dof = get_dof_component(node, component);
+ const int idc = m_ideq(dof);
+
+ if (idc == no_equation) continue;
+
+ const double tmp_v = velocity(dof);
+ const double tmp_d = displacement(dof);
+
+ double h = EPS_J*std::abs(tmp_v);
+ if (h == 0) h = EPS_J;
+ velocity(dof) = tmp_v + h;
+ h = velocity(dof) - tmp_v;
+
+ displacement(dof) = tmp_d + alphadt*h;
+
+ element_residual_transport_component(element, component, displacement, velocity, element_residual);
+ velocity(dof) = tmp_v;
+ displacement(dof) = tmp_d;
+
+ for (int enr: m_mesh->range_nen())
+ {
+ const ind_t noder = m_mesh->get_node(element, enr);
+ const ind_t idr = m_ideq(get_dof_component(noder, component));
+
+ if (idr == no_equation) continue;
+ jacobian.push_back(triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h));
+ }
+ }
+
+ }
+}
+
+void SaturatedDiffusionProgram::update_solution(
+ const Eigen::VectorXd& update,
+ double lambda,
+ double alpha_dt,
+ Eigen::VectorXd& predictor,
+ Eigen::VectorXd& displacement,
+ Eigen::VectorXd& velocity
+ )
+{
+ for (int node: m_mesh->range_nodes())
+ {
+ for (int component: m_data->range_aqueous_component())
+ {
+ const int dof = get_dof_component(node, component);
+ const int id = m_ideq(dof);
+ if (id == no_equation) continue;
+ velocity(dof) += update(id);
+ }
+ }
+ displacement = predictor + alpha_dt*velocity;
+}
+} // end namespace siasaturated
+} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
diff --git a/src/reactmicp/systems/saturated_diffusion/transport_program.hpp b/src/reactmicp/systems/saturated_diffusion/transport_program.hpp
index 8e2e74c..cd2744c 100644
--- a/src/reactmicp/systems/saturated_diffusion/transport_program.hpp
+++ b/src/reactmicp/systems/saturated_diffusion/transport_program.hpp
@@ -1,10 +1,137 @@
#ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP
#define SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP
+#include <memory>
+
+#include "database/data_container.hpp"
+#include "dfpmsolver/parabolic_program.hpp"
+#include "transport_parameters.hpp"
+#include "reactmicp/meshes/mesh1d.hpp"
+
namespace specmicp {
namespace reactmicp {
+namespace systems {
+namespace siasaturated {
+
+class SaturatedDiffusionProgram: dfpmsolver::ParabolicProgram<SaturatedDiffusionProgram>
+{
+public:
+ SaturatedDiffusionProgram(
+ std::shared_ptr<mesh::Mesh1D> the_mesh,
+ std::shared_ptr<database::DataContainer> the_data,
+ std::shared_ptr<SaturatedDiffusionTransportParameters> the_param
+ ):
+ m_mesh(the_mesh),
+ m_data(the_data),
+ m_param(the_param),
+ m_external_flux(Eigen::VectorXd::Zero(the_mesh->nnodes()*(the_data->nb_component-1)))
+ {}
+
+ const int no_equation = -1;
+ //! \brief Return the number of equations
+ int get_neq() const {return m_neq;}
+ //! \brief Return the number of degrees of freedom per node
+ int get_ndf() const {return m_data->nb_component-1;}
+
+ // number the equations
+ void number_equations(std::vector<int> list_fixed_nodes);
+
+ //! \brief Compute the residuals
+ void compute_residuals(const Eigen::VectorXd& displacement,
+ const Eigen::VectorXd& velocity,
+ Eigen::VectorXd& residual
+ );
+
+ //! \brief Compute the jacobian
+ void compute_jacobian(Eigen::VectorXd& displacement,
+ Eigen::VectorXd& velocity,
+ Eigen::SparseMatrix<double> &jacobian,
+ double alphadt);
+
+ //! \brief Update the solutions
+ void update_solution(const Eigen::VectorXd& update,
+ double lambda,
+ double alpha_dt,
+ Eigen::VectorXd& predictor,
+ Eigen::VectorXd& displacement,
+ Eigen::VectorXd& velocity);
+
+ //! \brief Apply boundary conditions to the velocity vector
+ //!
+ //! by default do nothing.
+ void apply_bc(double dt,
+ const Eigen::VectorXd& displacement,
+ Eigen::VectorXd velocity)
+ {}
+
+ //! Return the degree of freedom number of 'component' at 'node'
+ int get_dof_component(int node, int component) const {
+ assert(component >=1 and component < m_data->nb_component);
+ assert(node >= 0 and node < m_mesh->nnodes());
+ return (component - 1)+get_ndf()*node;
+ }
+ //! Return the degree of freedom number of 'component' at 'node'
+ int get_dof_mineral(int node, int mineral) const {
+ assert(mineral >= 0 and mineral < m_data->nb_mineral);
+ assert(node >= 0 and node < m_mesh->nnodes());
+ return mineral+m_data->nb_mineral*node;
+ }
+ //! Return the total mobile concentration
+ double& get_total_mobile_concentration(int node, int component, Eigen::VectorXd& concentrations) const
+ {
+ return concentrations.coeffRef(get_dof_component(node, component));
+ }
+ //! Return the total mobile concentration
+ double get_total_mobile_concentration(int node, int component, const Eigen::VectorXd& concentrations) const
+ {
+ return concentrations.coeff(get_dof_component(node, component));
+ }
+ //! Return the external flux
+ Eigen::VectorXd& external_flux() {return m_external_flux;}
+ //! Return the external flux for 'component' at 'node'
+ double& external_flux(int node, int component) {return m_external_flux(get_dof_component(node, component));}
+
+private:
+
+ void element_residual_transport(
+ ind_t element,
+ const Eigen::VectorXd& displacement,
+ const Eigen::VectorXd& velocity,
+ Eigen::VectorXd& residual);
+
+ void element_residual_transport_component(
+ ind_t element,
+ int component,
+ const Eigen::VectorXd& displacement,
+ const Eigen::VectorXd& velocity,
+ Vector& element_residual);
+
+ double nodal_mineral_transient_term_transport(
+ ind_t node,
+ ind_t component);
+
+
+ void element_jacobian_transport(
+ ind_t element,
+ Eigen::VectorXd& displacement,
+ Eigen::VectorXd& velocity,
+ list_triplet_t& jacobian,
+ double alphadt);
+
+ int m_neq;
+ Eigen::VectorXd m_ideq;
+
+ std::shared_ptr<mesh::Mesh1D> m_mesh;
+ std::shared_ptr<database::DataContainer> m_data;
+ std::shared_ptr<SaturatedDiffusionTransportParameters> m_param;
+
+ Eigen::VectorXd m_external_flux;
+
+};
+} // end namespace siasaturated
+} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
#endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP
diff --git a/tests/reactmip/systems/saturated_diffusion/solver.cpp b/tests/reactmip/systems/saturated_diffusion/solver.cpp
new file mode 100644
index 0000000..5a33c39
--- /dev/null
+++ b/tests/reactmip/systems/saturated_diffusion/solver.cpp
@@ -0,0 +1,51 @@
+#include "catch.hpp"
+
+#include "utils.hpp"
+#include "dfpmsolver/parabolic_driver.hpp"
+#include "reactmicp/systems/saturated_diffusion/transport_program.hpp"
+
+using namespace specmicp;
+using namespace specmicp::reactmicp;
+using namespace specmicp::reactmicp::systems;
+using namespace specmicp::reactmicp::systems::siasaturated;
+
+
+TEST_CASE("dfpm solver", "[dfpm, solver]") {
+ SECTION("transport program") {
+ //specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
+ //specmicp::logger::ErrFile::stream() = &std::cerr;
+
+ std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
+ int nb_nodes = 5;
+
+ EquilibriumState initial_state = sample_carbo_composition(database);
+
+ std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(nb_nodes-1, 1e-3, 0.001);
+ auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(1e-8, 0.2);
+
+ SaturatedDiffusionProgram the_program(themesh, database, parameters);
+
+ Eigen::VectorXd concentrations(themesh->nnodes()*the_program.get_ndf());
+ concentrations.setOnes();
+ concentrations.block(0, 0, the_program.get_ndf(), 1).setZero();
+ the_program.number_equations({0,});
+
+ dfpmsolver::ParabolicDriver<SaturatedDiffusionProgram> solver(the_program);
+ solver.get_velocity() = Eigen::VectorXd::Zero(concentrations.rows());
+ double dt = 100;
+ dfpmsolver::ParabolicDriverReturnCode retcode = solver.solve_timestep(dt, concentrations);
+ REQUIRE(retcode == dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
+
+ for (int node=1; node<themesh->nnodes(); ++node)
+ {
+ for (int component: database->range_aqueous_component())
+ {
+ the_program.external_flux(node, component) =
+ parameters->porosity(0)*(1.0 - concentrations(the_program.get_dof_component(node, component)))
+ / (parameters->density_water()*dt*themesh->cell_volume(node));
+ }
+ }
+ retcode = solver.restart_timestep(concentrations);
+ REQUIRE(retcode == dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
+ }
+}
diff --git a/tests/reactmip/systems/saturated_diffusion/transport_program.cpp b/tests/reactmip/systems/saturated_diffusion/transport_program.cpp
new file mode 100644
index 0000000..86e540b
--- /dev/null
+++ b/tests/reactmip/systems/saturated_diffusion/transport_program.cpp
@@ -0,0 +1,26 @@
+#include "catch.hpp"
+
+#include "utils.hpp"
+#include "reactmicp/systems/saturated_diffusion/transport_program.hpp"
+
+using namespace specmicp;
+using namespace specmicp::reactmicp;
+using namespace specmicp::reactmicp::systems;
+using namespace specmicp::reactmicp::systems::siasaturated;
+
+TEST_CASE("Transport program", "[Transport, Diffusion, Finite Element, Solver]") {
+
+ std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
+ int nb_nodes = 5;
+
+ EquilibriumState initial_state = sample_carbo_composition(database);
+
+ std::shared_ptr<mesh::Mesh1D> themesh = std::make_shared<mesh::Mesh1D>(nb_nodes-1, 1e-3, 0.001);
+ auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(1e-6, 0.2);
+
+ SECTION("Initialization") {
+ SaturatedDiffusionProgram(themesh, database, parameters);
+ }
+
+
+}