diff --git a/CMakeLists.txt b/CMakeLists.txt
index a8a1789..3298d26 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,245 +1,232 @@
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_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src)
# include
include_directories(${EIGEN3_INCLUDE_DIR})
include_directories(${PROJECT_SOURCE_DIR})
-add_custom_target(common
+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)
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}/mesh.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/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/transport_program.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
)
# ------- 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 )
diff --git a/src/reactmicp/mesh.hpp b/src/reactmicp/mesh.hpp
new file mode 100644
index 0000000..a361e0f
--- /dev/null
+++ b/src/reactmicp/mesh.hpp
@@ -0,0 +1,12 @@
+/*-------------------------------------------------------
+
+ - Module : reactmicp
+ - File : mesh.hpp
+ - Author : Fabien Georget
+
+ Copyright (c) 2014, Fabien Georget, Princeton University
+
+---------------------------------------------------------*/
+
+#include "meshes/uniform_mesh1d.hpp"
+#include "meshes/axisymmetric_uniform_mesh1d.hpp"
diff --git a/src/reactmicp/meshes/axisymmetric_uniform_mesh1d.hpp b/src/reactmicp/meshes/axisymmetric_uniform_mesh1d.hpp
new file mode 100644
index 0000000..74b5ad6
--- /dev/null
+++ b/src/reactmicp/meshes/axisymmetric_uniform_mesh1d.hpp
@@ -0,0 +1,72 @@
+/*-------------------------------------------------------
+
+ - Module : reactmicp/meshes
+ - File : axisymmetric_uniform_mesh1d.hpp
+ - Author : Fabien Georget
+
+ Copyright (c) 2014, Fabien Georget, Princeton University
+
+---------------------------------------------------------*/
+
+#ifndef SPECMICP_REACTMICP_MESHES_AXISYMMETRICMESH1D_HPP
+#define SPECMICP_REACTMICP_MESHES_AXISYMMETRICMESH1D_HPP
+
+#include "mesh1d.hpp"
+
+// \file axisymmetric_uniform_mesh1d.hpp Uniform axisymmetric 1D mesh
+
+namespace specmicp {
+namespace reactmicp {
+namespace mesh {
+
+//! \brief A uniform 1D mesh
+class AxisymmetricUniformMesh1D: public Mesh1D
+{
+public:
+ AxisymmetricUniformMesh1D(index_t nb_nodes, scalar_t radius, scalar_t height):
+ Mesh1D(nb_nodes),
+ m_radius(radius),
+ m_dx(radius/(nb_nodes-1)),
+ m_height(height)
+ {}
+
+ scalar_t get_radius_node(index_t node) {
+ return m_radius-m_dx*node;
+ }
+ //! Return the radius of a face of an element
+ scalar_t get_radius_face(index_t element) {
+ return get_radius_node(get_node(element, 1))+m_dx/2;
+ }
+
+ scalar_t get_dx(index_t _) {return m_dx;}
+ scalar_t get_face_area(index_t element) {
+ return 2*M_PI*get_radius_face(element)*m_height;}
+ scalar_t get_volume_element(index_t element) {return M_PI*m_height*(
+ std::pow(get_radius_node(get_node(element, 0)),2)
+ -std::pow(get_radius_node(get_node(element, 1)),2));}
+ scalar_t get_volume_cell(index_t node) {
+ if (node ==0)
+ return M_PI*m_height*(std::pow(m_radius,2)-std::pow(get_radius_face(0),2));
+ else if (node == nb_nodes()-1)
+ return M_PI*std::pow(m_dx/2.0,2)*m_height;
+ else
+ return M_PI*m_height*(std::pow(get_radius_face(node-1),2) - std::pow(get_radius_face(node),2));
+ }
+
+private:
+ scalar_t m_radius;
+ scalar_t m_dx;
+ scalar_t m_height;
+};
+
+//! \brief Factory method to build a pointer to a uniform 1D mesh
+inline Mesh1DPtr axisymmetric_uniform_mesh1d(index_t nb_nodes, scalar_t dx, scalar_t cross_section)
+{
+ return std::static_pointer_cast<Mesh1D>(std::make_shared<AxisymmetricUniformMesh1D>(nb_nodes, dx, cross_section));
+}
+
+} // end namespace mesh
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_MESHES_AXISYMMETRICMESH1D_HPP
diff --git a/src/reactmicp/meshes/mesh1d.hpp b/src/reactmicp/meshes/mesh1d.hpp
index f639804..fbbf074 100644
--- a/src/reactmicp/meshes/mesh1d.hpp
+++ b/src/reactmicp/meshes/mesh1d.hpp
@@ -1,116 +1,95 @@
/*-------------------------------------------------------
- Module : reactmicp/meshes
- File : mesh1d.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_REACTMICP_MESH_MESH1D_HPP
#define SPECMICP_REACTMICP_MESH_MESH1D_HPP
#include "common.hpp"
#include <memory>
namespace specmicp {
namespace reactmicp {
namespace mesh {
+//! AbstractBase class for a 1D mesh
class Mesh1D
{
public:
const index_t nen = 2; //!< Number of elemental nodes
Mesh1D(index_t nb_nodes): m_nb_nodes(nb_nodes) {}
virtual ~Mesh1D() {}
//! \brief Return the number of elements
index_t nb_elements() {return m_nb_nodes-1;}
//! \brief Return the number of nodes
index_t nb_nodes() {return m_nb_nodes;}
//! \brief Return the node from local nod enumber and element
index_t get_node(index_t element, index_t index)
{
specmicp_assert(index < nen);
return element+index;
}
//! \brief Return the length of an element
virtual scalar_t get_dx(index_t element) = 0;
//! \brief Return the area of the face at the middle of an element
virtual scalar_t get_face_area(index_t element) = 0;
//! \brief Return the volume of an element
virtual scalar_t get_volume_element(index_t element) = 0;
//! \brief Return the volume of a cell (element of the dual mesh)
virtual scalar_t get_volume_cell(index_t node) = 0;
//! \brief Range over the elements
range_t range_elements() {return boost::irange((index_t) 0, nb_elements());}
//! \brief Range over the nodes
range_t range_nodes() {return boost::irange((index_t) 0, nb_nodes());}
//! \brief Range over the elemental nodes
range_t range_nen() {return boost::irange((index_t) 0, nen);}
private:
index_t m_nb_nodes; //!< Number of elements
};
-//! a one dimensional mesh
-class UniformMesh1D: public Mesh1D
-{
-public:
- UniformMesh1D(index_t nb_nodes, scalar_t dx, scalar_t cross_section):
- Mesh1D(nb_nodes), m_dx(dx), m_crosssection(cross_section)
- {}
-
- scalar_t get_dx(index_t _) {return m_dx;}
- scalar_t get_face_area(index_t _) { return m_crosssection;}
- scalar_t get_volume_element(index_t _) {return m_dx*m_crosssection;}
- scalar_t get_volume_cell(index_t node) {
- if (node ==0 or node == nb_nodes()-1) return m_dx*m_crosssection/2;
- else return m_dx*m_crosssection;
- }
-
-private:
- scalar_t m_dx;
- scalar_t m_crosssection;
-};
-
+//! \brief type of a pointer to a mesh
+//!
+//! This is a smart pointer allowing the mesh to be shared freely
using Mesh1DPtr = std::shared_ptr<Mesh1D>;
-inline Mesh1DPtr uniform_mesh1d(index_t nb_nodes, scalar_t dx, scalar_t cross_section)
-{
- return std::static_pointer_cast<Mesh1D>(std::make_shared<UniformMesh1D>(nb_nodes, dx, cross_section));
-}
} // end namespace mesh
} // end namespace reactmicp
} // end namespace specmicp
#endif // SPECMICP_REACTMICP_MESH_MESH1D_HPP
diff --git a/src/reactmicp/meshes/uniform_mesh1d.hpp b/src/reactmicp/meshes/uniform_mesh1d.hpp
new file mode 100644
index 0000000..60746aa
--- /dev/null
+++ b/src/reactmicp/meshes/uniform_mesh1d.hpp
@@ -0,0 +1,53 @@
+/*-------------------------------------------------------
+
+ - Module : reactmicp/meshes
+ - File : uniform_mesh1d
+ - Author : Fabien Georget
+
+ Copyright (c) 2014, Fabien Georget, Princeton University
+
+---------------------------------------------------------*/
+
+#ifndef SPECMICP_REACTMICP_MESH_UNIFORMMESH1D_HPP
+#define SPECMICP_REACTMICP_MESH_UNIFORMMESH1D_HPP
+
+#include "mesh1d.hpp"
+
+// \file uniform_mesh1d.hpp Uniform 1D mesh
+
+namespace specmicp {
+namespace reactmicp {
+namespace mesh {
+
+//! \brief A uniform 1D mesh
+class UniformMesh1D: public Mesh1D
+{
+public:
+ UniformMesh1D(index_t nb_nodes, scalar_t dx, scalar_t cross_section):
+ Mesh1D(nb_nodes), m_dx(dx), m_crosssection(cross_section)
+ {}
+
+ scalar_t get_dx(index_t _) {return m_dx;}
+ scalar_t get_face_area(index_t _) { return m_crosssection;}
+ scalar_t get_volume_element(index_t _) {return m_dx*m_crosssection;}
+ scalar_t get_volume_cell(index_t node) {
+ if (node ==0 or node == nb_nodes()-1) return m_dx*m_crosssection/2;
+ else return m_dx*m_crosssection;
+ }
+
+private:
+ scalar_t m_dx;
+ scalar_t m_crosssection;
+};
+
+//! \brief Factory method to build a pointer to a uniform 1D mesh
+inline Mesh1DPtr uniform_mesh1d(index_t nb_nodes, scalar_t dx, scalar_t cross_section)
+{
+ return std::static_pointer_cast<Mesh1D>(std::make_shared<UniformMesh1D>(nb_nodes, dx, cross_section));
+}
+
+} // end namespace mesh
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_MESH_UNIFORMMESH1D_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 5361f66..df8bf86 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,136 +1,151 @@
##################### External packages #########################
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
)
################## Configuration ################################
set(PROJECT_TEST_DIR ${CMAKE_CURRENT_SOURCE_DIR})
ExternalProject_Get_Property(catch source_dir)
set(CATCH_INCLUDE_DIR ${source_dir}/include CACHE INTERNAL "Path to include folder for Catch")
include_directories(${CATCH_INCLUDE_DIR})
include_directories(${PROJECT_TEST_DIR})
################## Tests ########################################
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
#===========
+# Database
+# --------
+
add_custom_target(catch_header SOURCES 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)
+# Mesh
+# ----
+
+add_executable(test_meshes
+ ${PROJECT_TEST_DIR}/reactmip/meshes/test_meshes.cpp
+ ${PROJECT_TEST_DIR}/reactmip/meshes/test_uniform_mesh_1d.cpp
+ ${PROJECT_TEST_DIR}/reactmip/meshes/test_axisymmetric_uniform_mesh_1d.cpp
+)
+
+# Saturated Diffusion
+# -------------------
+
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/reactive_transport.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/tests/reactmip/meshes/test_axisymmetric_uniform_mesh_1d.cpp b/tests/reactmip/meshes/test_axisymmetric_uniform_mesh_1d.cpp
new file mode 100644
index 0000000..f86d518
--- /dev/null
+++ b/tests/reactmip/meshes/test_axisymmetric_uniform_mesh_1d.cpp
@@ -0,0 +1,42 @@
+#include "catch.hpp"
+
+#include "reactmicp/meshes/axisymmetric_uniform_mesh1d.hpp"
+
+using namespace specmicp::reactmicp::mesh;
+
+TEST_CASE("Axisymmetric uniform mesh 1d", "[Reactmicp, mesh, axisymmetric, uniform, 1D]") {
+
+ SECTION("test uniform mesh 1d") {
+
+ double nb_node = 6;
+ double radius = 5;
+ double dx = 1;
+ double height = 2;
+
+ AxisymmetricUniformMesh1D mesh(nb_node, radius, height);
+
+ REQUIRE(mesh.nb_nodes() == nb_node);
+ REQUIRE(mesh.nb_elements() == nb_node - 1);
+
+ REQUIRE(mesh.get_dx(0) == dx);
+ REQUIRE(mesh.get_radius_face(0) == 4.5);
+ REQUIRE(mesh.get_radius_face(1) == 3.5);
+ REQUIRE(std::abs(mesh.get_face_area(0) - 2*M_PI*4.5*height) < 1e-10);
+ REQUIRE(std::abs(mesh.get_volume_element(0) - M_PI*height*(5*5-4*4)) < 1e-10);
+ REQUIRE(std::abs(mesh.get_volume_cell(0) - M_PI*height*(5*5-4.5*4.5)) < 1e-10);
+ REQUIRE(std::abs(mesh.get_volume_element(1) - M_PI*height*(4*4-3*3)) < 1e-10);
+ REQUIRE(std::abs(mesh.get_volume_cell(1) - M_PI*height*(4.5*4.5-3.5*3.5)) < 1e-10);
+
+ int check_nb_node = 0;
+ for (int node: mesh.range_nodes()) {++check_nb_node;}
+ REQUIRE(check_nb_node == nb_node);
+
+ int check_nb_element = 0;
+ for (int element: mesh.range_elements()) {++check_nb_element;}
+ REQUIRE(check_nb_element == mesh.nb_elements());
+
+ int check_nen = 0;
+ for (int node: mesh.range_nen()) {++check_nen;}
+ REQUIRE(check_nen == 2);
+ }
+}
diff --git a/tests/reactmip/meshes/test_meshes.cpp b/tests/reactmip/meshes/test_meshes.cpp
new file mode 100644
index 0000000..178916e
--- /dev/null
+++ b/tests/reactmip/meshes/test_meshes.cpp
@@ -0,0 +1,3 @@
+#define CATCH_CONFIG_MAIN
+#include "catch.hpp"
+
diff --git a/tests/reactmip/meshes/test_uniform_mesh_1d.cpp b/tests/reactmip/meshes/test_uniform_mesh_1d.cpp
new file mode 100644
index 0000000..2530a00
--- /dev/null
+++ b/tests/reactmip/meshes/test_uniform_mesh_1d.cpp
@@ -0,0 +1,36 @@
+#include "catch.hpp"
+
+#include "reactmicp/meshes/uniform_mesh1d.hpp"
+
+using namespace specmicp::reactmicp::mesh;
+
+TEST_CASE("Uniform mesh 1d", "[Reactmicp, mesh, uniform, 1D]") {
+
+ SECTION("test uniform mesh 1d") {
+
+ double nb_node = 5;
+ double dx = 1;
+ double cross = 2;
+
+ UniformMesh1D mesh(nb_node, dx, cross);
+
+ REQUIRE(mesh.nb_nodes() == nb_node);
+ REQUIRE(mesh.nb_elements() == nb_node - 1);
+
+ REQUIRE(mesh.get_dx(0) == dx);
+ REQUIRE(mesh.get_face_area(2) == cross);
+ REQUIRE(mesh.get_volume_element(2) == dx*cross);
+
+ int check_nb_node = 0;
+ for (int node: mesh.range_nodes()) {++check_nb_node;}
+ REQUIRE(check_nb_node == nb_node);
+
+ int check_nb_element = 0;
+ for (int element: mesh.range_elements()) {++check_nb_element;}
+ REQUIRE(check_nb_element == mesh.nb_elements());
+
+ int check_nen = 0;
+ for (int node: mesh.range_nen()) {++check_nen;}
+ REQUIRE(check_nen == 2);
+ }
+}
diff --git a/tests/reactmip/systems/saturated_diffusion/reactive_transport.cpp b/tests/reactmip/systems/saturated_diffusion/reactive_transport.cpp
index cf5caf2..355d3af 100644
--- a/tests/reactmip/systems/saturated_diffusion/reactive_transport.cpp
+++ b/tests/reactmip/systems/saturated_diffusion/reactive_transport.cpp
@@ -1,95 +1,96 @@
#include "catch.hpp"
#include <iostream>
#include "utils.hpp"
+#include "reactmicp/meshes/uniform_mesh1d.hpp"
#include "reactmicp/systems/saturated_diffusion/boundary_conditions.hpp"
#include "reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp"
#include <fstream>
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems;
using namespace specmicp::reactmicp::systems::siasaturated;
TEST_CASE("reactive transport solver", "[reactive transport, speciation, dfpm, solver]") {
specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
specmicp::logger::ErrFile::stream() = &std::cerr;
int nb_nodes = 10;
std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 0.1, 2.5);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-4, 0.2);
EquilibriumState initial_state = sample_carbo_composition(database, 0.02);
// scaling {
double volume = initial_state.volume_minerals();
volume /= (1 - parameters->porosity(1) );
double scale = themesh->get_volume_cell(1)*1e-6 / volume;
initial_state.scale_condensed(scale);
// }
SIABoundaryConditions bcs(nb_nodes);
bcs.list_initial_states.push_back(sample_carbo_composition(database));
bcs.list_initial_states.push_back(blank_composition(database));
bcs.bs_types[0] = BoundaryConditionType::FixedComposition;
bcs.initial_states[0] = 1;
SECTION("Initialization") {
SIASaturatedReactiveTransportSolver solver(themesh, database, parameters);
solver.apply_boundary_conditions(bcs);
SIASaturatedVariables& variables = solver.get_variables();
REQUIRE(std::abs(variables.mineral_amount(1, 1)
- bcs.list_initial_states[0].moles_mineral(1)) < 1e-10 );
REQUIRE(std::abs(variables.component_concentration(1, 1)
- bcs.list_initial_states[0].molality_component(1)) < 1e-10 );
REQUIRE(std::abs(variables.mineral_amount(0, 1)
- bcs.list_initial_states[1].moles_mineral(1)) < 1e-10 );
REQUIRE(std::abs(variables.component_concentration(0, 1)
- bcs.list_initial_states[1].molality_component(1)) < 1e-10 );
}
SECTION("Solving") {
std::ofstream output;
output.open("out.dat");
SIASaturatedReactiveTransportSolver solver(themesh, database, parameters);
solver.apply_boundary_conditions(bcs);
for (auto it=database->labels_basis.begin(); it!=database->labels_basis.end(); ++it)
{
std::cout << *it << " ";
}
std::cout << std::endl;
std::cout << database->nu_mineral << std::endl;
solver.use_sia(50);
double dt = 100;
double total=0 ;
for (int k=0; k<10000; ++k)
{
std::cout << " iterations : " << k << std::endl;
SIASaturatedReactiveTransportSolverReturnCode retcode = solver.solve_timestep(dt);
REQUIRE(retcode == SIASaturatedReactiveTransportSolverReturnCode::Success);
SIASaturatedVariables& variables = solver.get_variables();
output << total << " " << total/(3600.0*24.0);
for(int node: themesh->range_nodes())
{
output << " " << variables.nodal_component_update_total_amount(node, 3);
}
output << std::endl;
total+=dt;
}
output.close();
}
}
diff --git a/tests/reactmip/systems/saturated_diffusion/solver.cpp b/tests/reactmip/systems/saturated_diffusion/solver.cpp
index b32bcb2..2bb898b 100644
--- a/tests/reactmip/systems/saturated_diffusion/solver.cpp
+++ b/tests/reactmip/systems/saturated_diffusion/solver.cpp
@@ -1,90 +1,91 @@
#include "catch.hpp"
#include "utils.hpp"
#include "dfpmsolver/parabolic_driver.hpp"
+#include "reactmicp/meshes/uniform_mesh1d.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_program.hpp"
#include "reactmicp/systems/saturated_diffusion/variables.hpp"
#include "specmicp/reduced_system_solver.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") {
std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
int nb_nodes = 5;
EquilibriumState initial_state = sample_carbo_composition(database);
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 1e-3, 0.001);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
SIASaturatedVariables variables(database, themesh, parameters);
SaturatedDiffusionProgram the_program(themesh, database, parameters);
Eigen::VectorXd& concentrations = variables.total_mobile_concentrations();
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->nb_nodes(); ++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->get_volume_cell(node));
}
}
retcode = solver.restart_timestep(concentrations);
REQUIRE(retcode == dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
}
SECTION("micpsolver") {
std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
database::Database database_handler(database);
int nb_nodes = 5;
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 1e-3, 0.001);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-8, 0.2);
SIASaturatedVariables variables(database, themesh, parameters);
EquilibriumState initial_state = sample_carbo_composition(database);
variables.update_composition(1, initial_state);
Eigen::VectorXd totconc;
variables.nodal_update_total_amount(1, totconc);
totconc(database_handler.component_label_to_id("HCO3[-]")) += 0.05;
totconc(database_handler.component_label_to_id("HO[-]")) -= 0.05;
ReducedSystemSolverOptions options;
options.solver_options.max_factorization_step = 1;
ReducedSystemSolver solver(database, totconc, options);
Eigen::VectorXd spec_var(variables.speciation_variables().col(1));
specmicp::micpsolver::MiCPPerformance perf = solver.solve(spec_var);
REQUIRE((int) perf.return_code == (int) specmicp::micpsolver::MiCPSolverReturnCode::ResidualMinimized);
EquilibriumState final_compo = solver.get_solution(spec_var);
variables.update_composition(1, final_compo);
}
}
diff --git a/tests/reactmip/systems/saturated_diffusion/transport_program.cpp b/tests/reactmip/systems/saturated_diffusion/transport_program.cpp
index 6672847..b7697dc 100644
--- a/tests/reactmip/systems/saturated_diffusion/transport_program.cpp
+++ b/tests/reactmip/systems/saturated_diffusion/transport_program.cpp
@@ -1,26 +1,27 @@
#include "catch.hpp"
#include "utils.hpp"
+#include "reactmicp/meshes/uniform_mesh1d.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);
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 1e-3, 0.001);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-6, 0.2);
SECTION("Initialization") {
SaturatedDiffusionProgram(themesh, database, parameters);
}
}
diff --git a/tests/reactmip/systems/saturated_diffusion/variables.cpp b/tests/reactmip/systems/saturated_diffusion/variables.cpp
index cc1cd65..1d2845d 100644
--- a/tests/reactmip/systems/saturated_diffusion/variables.cpp
+++ b/tests/reactmip/systems/saturated_diffusion/variables.cpp
@@ -1,64 +1,64 @@
#include "catch.hpp"
#include "utils.hpp"
#include "reactmicp/systems/saturated_diffusion/variables.hpp"
-#include "reactmicp/meshes/mesh1d.hpp"
+#include "reactmicp/meshes/uniform_mesh1d.hpp"
using namespace specmicp;
using namespace specmicp::reactmicp;
using namespace specmicp::reactmicp::systems::siasaturated;
TEST_CASE("SIA Variables", "[Variables, Secondary variables, SIA]") {
std::shared_ptr<database::DataContainer> database = get_test_carbo_database();
int nb_nodes = 5;
mesh::Mesh1DPtr themesh = mesh::uniform_mesh1d(nb_nodes, 1e-3, 0.001);
auto parameters = std::make_shared<SaturatedDiffusionTransportParameters>(nb_nodes, 1e-6, 0.2);
SECTION("Initialization") {
SIASaturatedVariables variables(database, themesh, parameters);
REQUIRE_NOTHROW( variables.component_concentration(0, 2) );
REQUIRE_NOTHROW( variables.component_concentration(2, 2) );
REQUIRE_NOTHROW( variables.component_concentration(4, 4) );
REQUIRE_NOTHROW( variables.component_concentration(3, 1) );
REQUIRE_NOTHROW( variables.log_component_concentration(3, 3) );
REQUIRE_NOTHROW( variables.log_gamma_component(4, 2) );
}
SECTION("Ionic strength computation") {
SIASaturatedVariables variables(database, themesh, parameters);
database::Database database_handler(database);
variables.log_component_concentration(1, database_handler.component_label_to_id("Ca[2+]")) = -3;
variables.log_component_concentration(1, database_handler.component_label_to_id("HO[-]")) = -3;
variables.log_component_concentration(1, database_handler.component_label_to_id("HCO3[-]")) = -3;
variables.log_component_concentration(1, database_handler.component_label_to_id("SiO(OH)3[-]")) = -3;
variables.nodal_ionic_strength(1);
REQUIRE( variables.ionic_strength(1) == (4*1e-3 +3*1e-3)*0.5 );
variables.secondary_concentration(1, database_handler.aqueous_label_to_id("CO3[2-]")) = 1e-2;
variables.nodal_ionic_strength(1);
REQUIRE( variables.ionic_strength(1) == (4*1e-3 +3*1e-3 + 4*1e-2)*0.5 );
variables.secondary_concentration(1, database_handler.aqueous_label_to_id("CO2")) = 1e-1;
variables.nodal_ionic_strength(1);
REQUIRE( variables.ionic_strength(1) == (4*1e-3 +3*1e-3 + 4*1e-2)*0.5 );
}
SECTION("Update") {
SIASaturatedVariables variables(database, themesh, parameters);
EquilibriumState composition = sample_carbo_composition(database);
variables.update_composition(1, composition);
REQUIRE(variables.component_concentration(1, 1) == composition.molality_component(1));
REQUIRE(variables.total_mobile_concentration(1, 1) == composition.total_aqueous_concentration_component(1));
}
}