diff --git a/CMakeLists.txt b/CMakeLists.txt index 8042f81..9fe70f2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,235 +1,236 @@ 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") 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) 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/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/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 ) # --------- Examples ------------------ add_subdirectory( examples ) diff --git a/src/reactmicp/meshes/mesh1dfwd.hpp b/src/reactmicp/meshes/mesh1dfwd.hpp new file mode 100644 index 0000000..e3cb365 --- /dev/null +++ b/src/reactmicp/meshes/mesh1dfwd.hpp @@ -0,0 +1,30 @@ +/*------------------------------------------------------- + + - Module : reactmicp/meshes + - File : mesh1dfwd.hpp + - Author : Fabien Georget + + Copyright (c) 2014, Fabien Georget, Princeton University + +---------------------------------------------------------*/ + +#ifndef SPECMICP_REACTMICP_MESH_MESH1DFWD +#define SPECMICP_REACTMICP_MESH_MESH1DFWD + +#include + +// define the forward declaration of a 1D mesh + +namespace specmicp { +namespace reactmicp { + +namespace mesh { + +class Mesh1D; +using Mesh1DPtr = std::shared_ptr; + +} // end namespace mesh +} // end namespace reactmicp +} // end namespace specmicp + +#endif // SPECMICP_REACTMICP_MESH_MESH1DFWD diff --git a/src/reactmicp/systems/saturated_diffusion/options.hpp b/src/reactmicp/systems/saturated_diffusion/options.hpp index 47f8c38..a78eeaa 100644 --- a/src/reactmicp/systems/saturated_diffusion/options.hpp +++ b/src/reactmicp/systems/saturated_diffusion/options.hpp @@ -1,65 +1,69 @@ #ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_OPTIONS_HPP #define SPECMICP_REACTMICP_SATURATED_DIFFUSION_OPTIONS_HPP #include "common.hpp" #include "dfpmsolver/parabolic_structs.hpp" +#include "micpsolver/micpsolver_structs.hpp" +#include "specmicp/reduced_system_solver.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace siasaturated { //! \brief Return codes for the SIASaturatedReactiveTransportSolver enum class SIASaturatedReactiveTransportSolverReturnCode { LolThatsNotSuposedToHappen = -30, //!< Well, have fun debugging that... // Transport error codes ErrorInTransport = -20, //!< Error in the transport solver MaxIterationsInTransport = -21, //!< Maximum iterations reached in the transport solver StationaryPointsInTransport = -22, //!< Stationnary points in the transport solver // Speciation error codes ErrorInSpeciation = -10, //!< Error in the speciation solver MaxIterationsInSpeciation = -11, //!< Maximum iterations reached in the speciation solver StationaryPointsInSpeciation = -12, //!< Stationnary points in the speciation solver // Coupling algorithm MaxIterations = -1, //!< Maximum iterations of the coupling solver StationaryPoints = -2, //!< Stationary points reached in the coupling solver // The good ones NotConvergedYet = 0, //!< Problem is not converged yet Success = 1 //!< The problem is succesfully solved }; //! \brief Options for the SIASaturatedReactiveTransportSolver struct SIASaturatedReactiveTransportSolverOptions { int max_iterations; //!< Maximum number of fixed-point iterations scalar_t residual_tolerance; //!< The tolerance for fixed-point convergence scalar_t good_enough_transport_residual_tolerance; //!< If transport failed, this value is used to check if it is 'good' enough scalar_t threshold_update_transport; //!< Threshold to accept a step if the update is too small scalar_t threshold_update_disable_speciation; //!< Threshold to disable next computation of speciation at a node dfpmsolver::ParabolicDriverOptions transport_options; //!< Options to use for the transport solver + specmicp::ReducedSystemSolverOptions speciation_options; //!< Options for the speciation solver + SIASaturatedReactiveTransportSolverOptions(): max_iterations(50), residual_tolerance(1e-3), good_enough_transport_residual_tolerance(1e-2), threshold_update_transport(1e-10), threshold_update_disable_speciation(1e-10) {} //! \brief Use a sequential non-iterative algorithm void use_snia() {max_iterations=1;} //! \brief Use a sequential iterative algorithm void use_sia(int nb_iter=50) {max_iterations=nb_iter;} //! \brief Set the tolerance void set_tolerance(scalar_t tol) {residual_tolerance = tol;} }; } // end namespace siasaturated } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_OPTIONS_HPP diff --git a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp index b9bc227..e5825d0 100644 --- a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp +++ b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp @@ -1,262 +1,277 @@ //#define EIGEN_DEFAULT_IO_FORMAT IOFormat(8) #include "reactive_transport_solver.hpp" #include "specmicp/reduced_system_solver.hpp" +#include "boundary_conditions.hpp" +#include "reactmicp/meshes/mesh1d.hpp" #ifdef SPECMICP_USE_OPENMP #include #endif // SPECMICP_USE_OPENMP #include namespace specmicp { namespace reactmicp { namespace internals { //! \brief Store residual information struct SIAResiduals { index_t nb_iterations; scalar_t transport_residual_0; scalar_t transport_residual; scalar_t transport_update; Vector speciation_residual; Vector speciation_update; SIAResiduals(index_t nb_node): nb_iterations(0), speciation_residual(Vector::Zero(nb_node)), speciation_update(Vector::Zero(nb_node)) {} }; } // end namespace internals +SIASaturatedReactiveTransportSolver::SIASaturatedReactiveTransportSolver( + std::shared_ptr the_mesh, + RawDatabasePtr the_database, + std::shared_ptr transport_parameters +): + m_mesh(the_mesh), + m_database(the_database), + m_variables(the_database, the_mesh, transport_parameters), + m_transport_program(the_mesh, the_database, transport_parameters), + m_transport_solver(m_transport_program), + m_bcs(the_mesh->nb_nodes()) +{} + void SIASaturatedReactiveTransportSolver::apply_boundary_conditions(SIABoundaryConditions bcs) { specmicp_assert(bcs.bs_types.size() == (unsigned int) m_mesh->nnodes()); specmicp_assert(bcs.initial_states.size() == (unsigned int) m_mesh->nnodes()); std::vector list_fixed_nodes; for (index_t node: m_mesh->range_nodes()) { switch (bcs.bs_types[node]) { case BoundaryConditionType::FixedComposition: m_bcs(node) = 1; list_fixed_nodes.push_back(node); break; default: m_bcs(node) = 0; break; } m_variables.update_composition(node, bcs.list_initial_states[bcs.initial_states[node]]); m_transport_program.number_equations(list_fixed_nodes); } } scalar_t SIASaturatedReactiveTransportSolver::residuals_transport() { Eigen::VectorXd residuals; m_transport_solver.compute_residuals(m_variables.total_mobile_concentrations(), residuals); return residuals.norm(); } SIASaturatedReactiveTransportSolverReturnCode convert_dfpm_retcode(dfpmsolver::ParabolicDriverReturnCode retcode) { switch (retcode) { case dfpmsolver::ParabolicDriverReturnCode::MaxIterations: return SIASaturatedReactiveTransportSolverReturnCode::MaxIterationsInTransport; case dfpmsolver::ParabolicDriverReturnCode::StationaryPoint: return SIASaturatedReactiveTransportSolverReturnCode::StationaryPointsInTransport; default: return SIASaturatedReactiveTransportSolverReturnCode::ErrorInTransport; } } SIASaturatedReactiveTransportSolverReturnCode SIASaturatedReactiveTransportSolver::solve_timestep(scalar_t dt) { retcode_t return_code = retcode_t::NotConvergedYet; internals::SIAResiduals residuals_storage(m_mesh->nb_nodes()); m_flag_compute_speciation = Eigen::VectorXi::Ones(m_mesh->nb_nodes()); // initialization transport program and transport solver Eigen::VectorXd transport_variable(m_variables.total_mobile_concentrations()); m_transport_solver.get_options() = get_options().transport_options; m_transport_solver.get_velocity() = Eigen::VectorXd::Zero(transport_variable.rows()); m_transport_program.external_flux().setZero(); // residuals_storage.transport_residual_0 = residuals_transport(); m_minerals_start_iteration = m_variables.speciation_variables().block( m_database->nb_component, 0, m_database->nb_mineral, m_mesh->nb_nodes()); // solve transport // --------------- dfpmsolver::ParabolicDriverReturnCode retcode = m_transport_solver.solve_timestep(dt, transport_variable); if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized) { ERROR << "Failed to solve 1st transport problem; return code " << (int) retcode; return convert_dfpm_retcode(retcode); } residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual; residuals_storage.transport_update = m_transport_solver.get_perfs().current_update; // update m_variables.total_mobile_concentrations() = transport_variable; // Speciation // ---------- restart_speciation_iterations(dt, residuals_storage); // Convergence // ----------- return_code = check_convergence(residuals_storage); residuals_storage.nb_iterations +=1; // If SNIA, we stop if (get_options().max_iterations == 1) return retcode_t::Success; // SIA algorithm // ------------- while (return_code == retcode_t::NotConvergedYet) { //std::cout << m_flag_compute_speciation.sum() << std::endl; transport_variable = m_variables.total_mobile_concentrations(); retcode = m_transport_solver.restart_timestep(transport_variable); if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized) { if (retcode == dfpmsolver::ParabolicDriverReturnCode::MaxIterations and m_transport_solver.get_perfs().current_residual/residuals_storage.transport_residual_0 < get_options().good_enough_transport_residual_tolerance) { restart_speciation_iterations(dt, residuals_storage); return_code = retcode_t::Success; break; } ERROR << "Failed to solve transport problem; return code " << (int) retcode; return convert_dfpm_retcode(retcode); } residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual; residuals_storage.transport_update = m_transport_solver.get_perfs().current_update; m_variables.total_mobile_concentrations() = transport_variable; residuals_storage.nb_iterations +=1; restart_speciation_iterations(dt, residuals_storage); return_code = check_convergence(residuals_storage); } std::cout << "End fixed point : return code " << (int) return_code << " - nb iter : " <nnodes(); ++node) #else for (index_t node: m_mesh->range_nodes()) #endif { if (m_bcs(node) == 0 and m_flag_compute_speciation(node)) { solve_speciation_node(dt, node, residuals_storage); } } #ifdef _OPENMP } #endif } void SIASaturatedReactiveTransportSolver::solve_speciation_node( scalar_t dt, index_t node, internals::SIAResiduals& residuals_storage ) { Eigen::VectorXd totconc; m_variables.nodal_update_total_amount(node, totconc); - ReducedSystemSolver solver(m_database, totconc); - solver.get_options().solver_options.max_iter = 50; - solver.get_options().solver_options.factor_gradient_search_direction = 400; - solver.get_options().solver_options.use_scaling = false; + ReducedSystemSolver solver(m_database, totconc, get_options().speciation_options); + //solver.get_options().solver_options.max_iter = 50; + //solver.get_options().solver_options.factor_gradient_search_direction = 400; + //solver.get_options().solver_options.use_scaling = false; Eigen::VectorXd variables = m_variables.speciation_variables(node); micpsolver::MiCPPerformance perf = solver.solve(variables); if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized) { ERROR << "Failed to solve speciation solver for node " << node << "; return code = " << (int) perf.return_code; } residuals_storage.speciation_residual(node) = perf.current_residual; residuals_storage.speciation_update(node) = perf.current_update; if (perf.current_update < get_options().threshold_update_disable_speciation) m_flag_compute_speciation(node) = 0; EquilibriumState solution = solver.get_solution(variables); m_variables.update_composition(node, solution); // mineral const scalar_t denom = 1e-3*m_mesh->get_volume_cell(node)*dt; for (index_t component: m_database->range_aqueous_component()) { scalar_t sum = 0.0; for (index_t mineral: m_database->range_mineral()) { if (m_database->nu_mineral(mineral, component) == 0.0) continue; sum -= m_database->nu_mineral(mineral, component)*( m_variables.mineral_amount(node, mineral) -m_minerals_start_iteration(mineral, node) ); } m_transport_program.external_flux(node, component) = sum/denom; } } SIASaturatedReactiveTransportSolverReturnCode SIASaturatedReactiveTransportSolver::check_convergence( internals::SIAResiduals& residuals_storage ) { // compute residuals of transport const scalar_t norm_transport = residuals_transport(); retcode_t return_code = retcode_t::NotConvergedYet; if (norm_transport/residuals_storage.transport_residual_0 < get_options().residual_tolerance) { return_code = retcode_t::Success; } else if (residuals_storage.transport_update < get_options().threshold_update_transport) { return_code = retcode_t::Success; } else if (residuals_storage.nb_iterations >= get_options().max_iterations) { WARNING << "Maximum number of iterations (" << get_options().max_iterations << ") reached during fixed-point iterations"; return_code = retcode_t::MaxIterations; } return return_code; } } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp index 4701145..15cce2e 100644 --- a/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp +++ b/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.hpp @@ -1,91 +1,99 @@ #ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_REACTIVETRANSPORTSOLVER_HPP #define SPECMICP_REACTMICP_SATURATED_DIFFUSION_REACTIVETRANSPORTSOLVER_HPP #include "common.hpp" #include "options.hpp" #include "variables.hpp" #include "transport_program.hpp" #include "dfpmsolver/parabolic_driver.hpp" -#include "boundary_conditions.hpp" namespace specmicp { namespace reactmicp { + using namespace systems::siasaturated; +// Forward declarations +// -------------------- +namespace systems { +namespace siasaturated { + +class SIABoundaryConditions; + +} // end namespace siasaturated +} // end namespace systems + namespace internals { - struct SIAResiduals; // forward declaration + struct SIAResiduals; } // end namespace internals +// end forward declarations + + +// Class declaration +// ================= //! \brief SIA solver for saturated reactive transport problem //! class SIASaturatedReactiveTransportSolver { public: using retcode_t = SIASaturatedReactiveTransportSolverReturnCode; SIASaturatedReactiveTransportSolver( std::shared_ptr the_mesh, RawDatabasePtr the_database, std::shared_ptr transport_parameters - ): - m_mesh(the_mesh), - m_database(the_database), - m_variables(the_database, the_mesh, transport_parameters), - m_transport_program(the_mesh, the_database, transport_parameters), - m_transport_solver(m_transport_program), - m_bcs(the_mesh->nb_nodes()) - {} + ); void apply_boundary_conditions(SIABoundaryConditions bcs); SIASaturatedVariables& get_variables() {return m_variables;} //! \brief Solve a timestep SIASaturatedReactiveTransportSolverReturnCode solve_timestep(scalar_t dt); //! \brief Return a read-write reference to the options SIASaturatedReactiveTransportSolverOptions& get_options() {return m_options;} //! \brief Return a read-only reference to the options const SIASaturatedReactiveTransportSolverOptions& get_options() const {return m_options;} //! \brief Use the sequential iterative algorithm void use_sia(int nb_iterations) {get_options().use_sia(nb_iterations);} //! \brief Use the sequential non-iterative algorithm void use_snia() {get_options().use_snia();} private: //! \brief Restart the transport problem void restart_transport_iterations(internals::SIAResiduals& residuals_storage); //! \brief Restart the speciation problems void restart_speciation_iterations(scalar_t dt, internals::SIAResiduals& residuals_storage); //! \brief Solve the speciation problem at one node void solve_speciation_node(scalar_t dt, index_t node, internals::SIAResiduals& residuals_storage); retcode_t check_convergence(internals::SIAResiduals& residuals_storage); scalar_t residuals_transport(); SIASaturatedReactiveTransportSolverOptions m_options; mesh::Mesh1DPtr m_mesh; database::RawDatabasePtr m_database; SIASaturatedVariables m_variables; Eigen::MatrixXd m_minerals_start_iteration; // transport program SaturatedDiffusionProgram m_transport_program; dfpmsolver::ParabolicDriver m_transport_solver; // Eigen::VectorXi m_bcs; Eigen::VectorXi m_flag_compute_speciation; }; } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_REACTIVETRANSPORTSOLVER_HPP diff --git a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp index f33e4a5..f02212a 100644 --- a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp +++ b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp @@ -1,203 +1,218 @@ #include "transport_program.hpp" +#include "transport_parameters.hpp" +#include "reactmicp/meshes/mesh1d.hpp" + #include #define EPS_J 1e-8 namespace specmicp { namespace reactmicp { namespace systems { namespace siasaturated { +SaturatedDiffusionProgram::SaturatedDiffusionProgram( + std::shared_ptr the_mesh, + std::shared_ptr the_data, + std::shared_ptr the_param + ): + m_mesh(the_mesh), + m_data(the_data), + m_param(the_param), + m_external_flux(Vector::Zero(the_mesh->nb_nodes()*(the_data->nb_component-1))) +{} + + void SaturatedDiffusionProgram::number_equations(std::vector list_fixed_nodes ) { m_ideq = Eigen::VectorXd::Zero(get_ndf()*m_mesh->nb_nodes()); for (auto it=list_fixed_nodes.begin(); itrange_aqueous_component()) { m_ideq(get_dof_component(*it, component)) = no_equation; } } index_t neq = 0; for (index_t dof=0; dofrange_elements()) { element_residual_transport(element, displacement, velocity, residual); } } void SaturatedDiffusionProgram::element_residual_transport(index_t element, const Eigen::VectorXd& displacement, const Eigen::VectorXd& velocity, Eigen::VectorXd& residual) { Eigen::VectorXd element_residual(2); for (index_t component: m_data->range_aqueous_component()) { element_residual_transport_component(element, component, displacement, velocity, element_residual); for (index_t en: m_mesh->range_nen()) { const index_t node = m_mesh->get_node(element, en); const index_t id = m_ideq(get_dof_component(node, component)); if (id != no_equation) {residual(id) += element_residual(en);} } } } void SaturatedDiffusionProgram::element_residual_transport_component(index_t element, index_t component, const Vector &displacement, const Vector &velocity, Vector& element_residual) { Eigen::Matrix jacob; Eigen::Matrix evelocity, edisplacement; scalar_t mass_coeff = -(m_param->density_water() * m_mesh->get_volume_element(element)/2); const scalar_t porosity = (m_param->porosity(m_mesh->get_node(element, 0)) + m_param->porosity(m_mesh->get_node(element, 1)) )/2; const scalar_t inv_diff_coeff = (0.5/m_param->diffusion_coefficient(m_mesh->get_node(element, 0)) + 0.5/m_param->diffusion_coefficient(m_mesh->get_node(element, 1))); scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element) * porosity * m_param->density_water() * 1.0/inv_diff_coeff ); jacob << 1, -1, -1, 1; jacob *= flux_coeff; evelocity << velocity(get_dof_component(m_mesh->get_node(element, 0), component))* mass_coeff*m_param->porosity(m_mesh->get_node(element, 0)), velocity(get_dof_component(m_mesh->get_node(element, 1), component))* mass_coeff*m_param->porosity(m_mesh->get_node(element, 0)); 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 = evelocity+jacob*edisplacement; for (index_t en: m_mesh->range_nen()) { element_residual(en) += (m_mesh->get_volume_element(element)/2 *external_flux(m_mesh->get_node(element, en), component)); } } void SaturatedDiffusionProgram::compute_jacobian( Vector& displacement, Vector& velocity, Eigen::SparseMatrix& jacobian, scalar_t alphadt ) { list_triplet_t jacob; const index_t ncomp = m_data->nb_component-1; const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen); jacob.reserve(estimation); for (index_t element: m_mesh->range_elements()) { element_jacobian_transport(element, displacement, velocity, jacob, alphadt); } jacobian = Eigen::SparseMatrix(get_neq(), get_neq()); jacobian.setFromTriplets(jacob.begin(), jacob.end()); } // Transport // ========= void SaturatedDiffusionProgram::element_jacobian_transport( index_t element, Vector& displacement, Vector& velocity, list_triplet_t& jacobian, scalar_t alphadt) { for (index_t 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 (index_t en: m_mesh->range_nen()) { Eigen::VectorXd element_residual(Eigen::VectorXd::Zero(2)); const index_t node = m_mesh->get_node(element, en); const index_t dof = get_dof_component(node, component); const index_t idc = m_ideq(dof); if (idc == no_equation) continue; const scalar_t tmp_v = velocity(dof); const scalar_t tmp_d = displacement(dof); scalar_t 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 (index_t enr: m_mesh->range_nen()) { const index_t noder = m_mesh->get_node(element, enr); const index_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 Vector &update, scalar_t lambda, scalar_t alpha_dt, Vector &predictor, Vector &displacement, Vector &velocity ) { for (index_t node: m_mesh->range_nodes()) { for (index_t component: m_data->range_aqueous_component()) { const index_t dof = get_dof_component(node, component); const index_t id = m_ideq(dof); if (id == no_equation) continue; velocity(dof) += lambda*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 2aa3615..bcdc008 100644 --- a/src/reactmicp/systems/saturated_diffusion/transport_program.hpp +++ b/src/reactmicp/systems/saturated_diffusion/transport_program.hpp @@ -1,135 +1,133 @@ #ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP #define SPECMICP_REACTMICP_SATURATED_DIFFUSION_TRANSPORTPROGRAM_HPP #include #include "database/data_container.hpp" #include "dfpmsolver/parabolic_program.hpp" -#include "transport_parameters.hpp" -#include "reactmicp/meshes/mesh1d.hpp" +#include "reactmicp/meshes/mesh1dfwd.hpp" + + namespace specmicp { namespace reactmicp { namespace systems { namespace siasaturated { +class SaturatedDiffusionTransportParameters; + //! \brief Transport program for saturated diffusion class SaturatedDiffusionProgram: dfpmsolver::ParabolicProgram { public: using triplet_t = Eigen::Triplet; //!< Triplet type for building the transport matrix using list_triplet_t = std::vector; //!< List of triplet for building the sparse matrix SaturatedDiffusionProgram( std::shared_ptr the_mesh, std::shared_ptr the_data, std::shared_ptr the_param - ): - m_mesh(the_mesh), - m_data(the_data), - m_param(the_param), - m_external_flux(Vector::Zero(the_mesh->nb_nodes()*(the_data->nb_component-1))) - {} + ); //! \brief Return the number of equations index_t get_neq() const {return m_neq;} //! \brief Return the number of degrees of freedom per node index_t get_ndf() const {return m_data->nb_component-1;} // number the equations void number_equations(std::vector list_fixed_nodes); //! \brief Compute the residuals void compute_residuals(const Vector &displacement, const Vector &velocity, Vector &residual ); //! \brief Compute the jacobian void compute_jacobian(Vector& displacement, Vector& velocity, Eigen::SparseMatrix &jacobian, scalar_t alphadt); //! \brief Update the solutions void update_solution(const Vector& update, scalar_t lambda, scalar_t alpha_dt, Vector& predictor, Vector& displacement, Vector& velocity); //! \brief Apply boundary conditions to the velocity vector //! //! by default do nothing. void apply_bc(scalar_t dt, const Vector& displacement, Vector velocity) {} //! Return the degree of freedom number of 'component' at 'node' index_t get_dof_component(index_t node, index_t component) const { specmicp_assert(component >=1 and component < m_data->nb_component); specmicp_assert(node >= 0 and node < m_mesh->nnodes()); return (component - 1)+get_ndf()*node; } //! Return the degree of freedom number of 'component' at 'node' index_t get_dof_mineral(index_t node, index_t mineral) const { specmicp_assert(mineral >= 0 and mineral < m_data->nb_mineral); specmicp_assert(node >= 0 and node < m_mesh->nnodes()); return mineral+m_data->nb_mineral*node; } //! Return the total mobile concentration scalar_t& get_total_mobile_concentration(index_t node, index_t component, Vector& concentrations) const { return concentrations.coeffRef(get_dof_component(node, component)); } //! Return the total mobile concentration scalar_t get_total_mobile_concentration(index_t node, index_t component, const Vector& concentrations) const { return concentrations.coeff(get_dof_component(node, component)); } //! Return the external flux Vector& external_flux() {return m_external_flux;} //! Return the external flux for 'component' at 'node' scalar_t& external_flux(index_t node, index_t component) {return m_external_flux(get_dof_component(node, component));} private: void element_residual_transport( index_t element, const Eigen::VectorXd& displacement, const Eigen::VectorXd& velocity, Eigen::VectorXd& residual); void element_residual_transport_component( index_t element, index_t component, const Vector& displacement, const Vector& velocity, Vector& element_residual); void element_jacobian_transport( index_t element, Vector& displacement, Vector& velocity, list_triplet_t& jacobian, scalar_t alphadt); index_t m_neq; Vector m_ideq; mesh::Mesh1DPtr m_mesh; database::RawDatabasePtr m_data; std::shared_ptr m_param; Vector 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/src/reactmicp/systems/saturated_diffusion/variables.cpp b/src/reactmicp/systems/saturated_diffusion/variables.cpp index 4a8c19f..c6b3c97 100644 --- a/src/reactmicp/systems/saturated_diffusion/variables.cpp +++ b/src/reactmicp/systems/saturated_diffusion/variables.cpp @@ -1,78 +1,95 @@ #include "variables.hpp" +#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp" +#include "reactmicp/meshes/mesh1d.hpp" namespace specmicp { namespace reactmicp { namespace systems { namespace siasaturated { +SIASaturatedVariables::SIASaturatedVariables( + RawDatabasePtr thedatabase, + std::shared_ptr themesh, + std::shared_ptr theparam + ): + SecondaryVariables(themesh->nb_nodes(), thedatabase), + m_data(thedatabase), + m_mesh(themesh), + m_param(theparam), + m_total_mobile_concentration( + Vector((thedatabase->nb_component-1)* themesh->nb_nodes())), + m_speciation_variables( + Matrix(thedatabase->nb_component+thedatabase->nb_mineral, themesh->nb_nodes())) +{} + void SIASaturatedVariables::update_composition(index_t node, const EquilibriumState& composition) { // speciation m_speciation_variables.col(node) = composition.main_variables(); m_secondary_concentration.col(node) = composition.molality_secondary(); m_secondary_variables.col(node) = composition.activity_coefficient(); m_ionic_strength(node) = composition.ionic_strength(); // transport for (index_t component: m_data->range_aqueous_component()) { total_mobile_concentration(node, component) = composition.total_aqueous_concentration_component(component); } m_speciation_variables(0, node) = mass_water(node); } scalar_t SIASaturatedVariables::immobile_total_amount(index_t node, index_t component) { return (m_data->nu_mineral.col(component).array()* m_speciation_variables.block(m_data->nb_component, node, m_data->nb_mineral, 1).array()).sum(); } void SIASaturatedVariables::nodal_update_total_amount(index_t node, Vector& total_amount) { total_amount.resize(m_data->nb_component); total_amount(0) = mass_water(node)/m_data->molar_mass_basis_si(0)+immobile_total_amount(node, 0); for (index_t component: m_data->range_aqueous_component()) { total_amount(component) = nodal_component_update_total_amount(node, component); } } scalar_t SIASaturatedVariables::nodal_component_update_total_amount(index_t node, index_t component) { return mass_water(node)*total_mobile_concentration(node,component)+immobile_total_amount(node, component); } scalar_t SIASaturatedVariables::mass_water(index_t node) { return m_param->density_water()*m_param->porosity(node)*m_mesh->get_volume_cell(node)*1e-3; } scalar_t SIASaturatedVariables::volume_minerals(index_t node) { scalar_t volume = 0; for (index_t m: m_data->range_mineral()) { volume += mineral_amount(node, m)*m_data->molar_volume_mineral(m); } return volume; } void SIASaturatedVariables::scale_volume() { for (index_t node: m_mesh->range_nodes()) { const scalar_t volume = volume_minerals(node); if (volume == 0) continue; const scalar_t vol_solid = (1-m_param->porosity(node))*m_mesh->get_volume_cell(node); const scalar_t factor = vol_solid/(volume*1e6); //std::cout << factor << std::endl; for (index_t mineral: m_data->range_mineral()) { mineral_amount(node, mineral) *= factor; } } } } // end namespace siasaturated } // end namespace systems } // end namespace reactmicp } // end namespace specmicp diff --git a/src/reactmicp/systems/saturated_diffusion/variables.hpp b/src/reactmicp/systems/saturated_diffusion/variables.hpp index cc88ca2..99eda00 100644 --- a/src/reactmicp/systems/saturated_diffusion/variables.hpp +++ b/src/reactmicp/systems/saturated_diffusion/variables.hpp @@ -1,147 +1,144 @@ #ifndef SPECMICP_REACTMICP_SATURATED_DIFFUSION_VARIABLES_HPP #define SPECMICP_REACTMICP_SATURATED_DIFFUSION_VARIABLES_HPP #include "reactmicp/systems/secondary_variables/secondary.hpp" #include "specmicp/equilibrium_data.hpp" -#include "reactmicp/meshes/mesh1d.hpp" -#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp" +#include "reactmicp/meshes/mesh1dfwd.hpp" + namespace specmicp { namespace reactmicp { + namespace systems { namespace siasaturated { +class SaturatedDiffusionTransportParameters; + + + //! Variables for the SIASaturatedReactiveTransportSolver class SIASaturatedVariables: public SecondaryVariables { public: SIASaturatedVariables( RawDatabasePtr thedatabase, std::shared_ptr themesh, std::shared_ptr theparam - ): - SecondaryVariables(themesh->nb_nodes(), thedatabase), - m_data(thedatabase), - m_mesh(themesh), - m_param(theparam), - m_total_mobile_concentration( - Vector((thedatabase->nb_component-1)* themesh->nb_nodes())), - m_speciation_variables( - Matrix(thedatabase->nb_component+thedatabase->nb_mineral, themesh->nb_nodes())) - {} + ); + index_t get_dof_total_concentration(index_t node, index_t component) const { specmicp_assert(component >= 1 and component < m_data->nb_component); specmicp_assert(node >= 0 and node < m_speciation_variables.cols()); return (component+offset_total_component())+(m_data->nb_component+offset_total_component())*node; } //! Return the total mobile concentration of 'component' at 'node' scalar_t total_mobile_concentration(index_t node, index_t component) const { return m_total_mobile_concentration.coeff(get_dof_total_concentration(node, component)); } //! Return a reference to the total mobile concentration of 'component' at 'node' scalar_t& total_mobile_concentration(index_t node, index_t component) { return m_total_mobile_concentration.coeffRef(get_dof_total_concentration(node, component)); } Vector& total_mobile_concentrations() { return m_total_mobile_concentration; } Eigen::Block total_mobile_concentrations(index_t node) { return m_total_mobile_concentration.segment(node*(m_data->nb_component+offset_total_component()), (m_data->nb_component+offset_total_component())); } //! Return the concentration of 'component' at 'node' scalar_t component_concentration(index_t node, index_t component) const { specmicp_assert(component >= 1 and component < m_data->nb_component); specmicp_assert(node >= 0 and node < m_speciation_variables.cols()); return pow10(m_speciation_variables.coeff(component+offset_component(), node)); } //! Return the concentration of 'component' at 'node' scalar_t log_component_concentration(index_t node, index_t component) const { specmicp_assert(component >= 1 and component < m_data->nb_component); specmicp_assert(node >= 0 and node < m_speciation_variables.cols()); return m_speciation_variables.coeff(component+offset_component(), node); } //! Return a reference to the concentration of 'component' at 'node' scalar_t& log_component_concentration(index_t node, index_t component) { specmicp_assert(component >= 1 and component < m_data->nb_component); specmicp_assert(node >= 0 and node < m_speciation_variables.cols()); return m_speciation_variables.coeffRef(component+offset_component(), node); } //! Return the number of moles of 'mineral' at 'node' scalar_t mineral_amount(index_t node, index_t mineral) const { specmicp_assert(mineral >= 0 and mineral < m_data->nb_mineral); specmicp_assert(node >= 0 and node < m_speciation_variables.cols()); return m_speciation_variables.coeff(offset_mineral()+mineral, node); } //! Return a reference to the number of moles of 'mineral' at 'node' scalar_t& mineral_amount(index_t node, index_t mineral) { specmicp_assert(mineral >= 0 and mineral < m_data->nb_mineral); specmicp_assert(node >= 0 and node < m_speciation_variables.cols()); return m_speciation_variables.coeffRef(offset_mineral()+mineral, node); } //! Reference to the speciation variables Matrix& speciation_variables() { return m_speciation_variables; } //! Eigen Expression to the speciation variables Matrix::ColXpr speciation_variables(index_t node) { return m_speciation_variables.col(node); } //! Return the chemical composition of a node - needed to initialize the speciation solver EquilibriumState equilibrium_composition(index_t node) { return EquilibriumState(m_speciation_variables.col(node), m_secondary_concentration.col(node), m_secondary_variables.col(node), m_ionic_strength(node), m_data ); } //! Update the composition at 'node' using the EquilibriumState instance 'composition' void update_composition(index_t node, const EquilibriumState &composition); //! Return the total moles number for 'component' at 'node' from the speciation scalar_t immobile_total_amount(index_t node, index_t component); //! Return the mass of water at 'node' scalar_t mass_water(index_t node); //! Compute the total moles number at 'node' using the transport mobile concentration void nodal_update_total_amount(index_t node, Vector& total_amount); // Used in the update FEM->speciation scalar_t nodal_component_update_total_amount(index_t node, index_t component); //! \brief Return the volume of minerals at node 'node' scalar_t volume_minerals(index_t node); //! \brief scale to volume the minerals void scale_volume(); private: // where the component informations are stored index_t offset_total_component() const {return -1;} // where the component informations are stored index_t offset_component() const {return 0;} // where the mineral informations are stored index_t offset_mineral() const {return m_data->nb_component+offset_component();} RawDatabasePtr m_data; mesh::Mesh1DPtr m_mesh; std::shared_ptr m_param; Vector m_total_mobile_concentration; Matrix m_speciation_variables; }; } // end namespace siasaturated } // end namespace systems } // end namespace reactmicp } // end namespace specmicp #endif // SPECMICP_REACTMICP_SATURATED_DIFFUSION_VARIABLES_HPP