diff --git a/src/reactmicp/systems/diffusion/diffusion.cpp b/src/reactmicp/systems/diffusion/diffusion.cpp deleted file mode 100644 index b0ec6de..0000000 --- a/src/reactmicp/systems/diffusion/diffusion.cpp +++ /dev/null @@ -1,623 +0,0 @@ -/*------------------------------------------------------- - - - Module : reactmicp/systems/diffusion - - File : diffusion.cpp - - Author : Fabien Georget - -Copyright (c) 2014, Fabien Georget , 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. - ----------------------------------------------------------*/ - -#include "diffusion.hpp" -#include "physics/laws.hpp" -#include "utils/log.hpp" -#include - -#define EPS_J 1e-8 - -namespace specmicp { -namespace reactmicp { -namespace systems { - - -DiffusionProgram::DiffusionProgram(std::shared_ptr themesh, - std::shared_ptr thedatabase, - std::shared_ptr parameters, - diffusion::ListBoundaryCondition& list_bcs, - const EquilibriumState& initialstate, - Variables& var): - DiffusionProgram(themesh, thedatabase, parameters, list_bcs, var) -{ - initialize_no_bc_with_init(initialstate, var); -} - -// ================================== // -// // -// Residuals // -// // -// ================================== // - -void DiffusionProgram::get_residuals(const Variables& variable, - Vector& residual) -{ - m_is_micp = std::vector(get_neq(), false); - residual = Eigen::VectorXd::Zero(get_neq()); - residual_transport(variable, residual); - residual_speciation(variable, residual); -} - -// Transport -// ========= - -void DiffusionProgram::residual_transport(const Variables& variable, - Vector& residual) -{ - for (ind_t element: m_mesh->range_elements()) - { - element_residual_transport(element, variable, residual); - } -} - -void DiffusionProgram::element_residual_transport_component( - ind_t element, - int component, - const Variables& variable, - Vector& element_residual) -{ - - Eigen::Matrix2d mass, jacob; - Eigen::Vector2d velocity, displacement; - double mass_coeff = -(m_param->density_water() * - m_param->porosity(0) * - m_mesh->get_volume_element(element)/2); - mass << 1, 0, 0, 1; - mass *= mass_coeff; - double flux_coeff = -( m_mesh->get_face_area(element) / 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; - velocity << variable.velocity(m_ideq.get_dof_diffusion(m_mesh->get_node(element,0), component)), - variable.velocity(m_ideq.get_dof_diffusion(m_mesh->get_node(element,1), component)); - displacement << m_ideq.mobile_total_concentration(m_mesh->get_node(element,0), component, variable), - m_ideq.mobile_total_concentration(m_mesh->get_node(element,1), component, variable); - - element_residual = mass*velocity+jacob*displacement; - - for (int en: m_mesh->range_nen()) - { - element_residual(en) += m_mesh->get_volume_element(element) * - nodal_mineral_transient_term_transport( - m_mesh->get_node(element, en), component, variable)/2.0; - } -} - -void DiffusionProgram::element_residual_transport(ind_t element, - const Variables& variable, - Vector& residual) -{ - Eigen::VectorXd element_residual(2); - for (ind_t component: m_data->range_aqueous_component()) - { - element_residual_transport_component(element, component, variable, 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.id_equation_diffusion(node, component); - if (id != no_equation) {residual.coeffRef(id) += element_residual(en);} - } - } -} - -double DiffusionProgram::nodal_mineral_transient_term_transport(ind_t node, - ind_t component, - const Variables& variable) -{ - double mineral_term = 0; - for (ind_t mineral: m_data->range_mineral()) - { - if (m_data->nu_mineral(mineral, component) == 0 ) continue; - mineral_term -= ( - m_data->nu_mineral(mineral, component) * - variable.velocity(m_ideq.get_dof_mineral(node, mineral)) - ); - } - return mineral_term; -} - -// Speciation -// ========== - -void DiffusionProgram::residual_speciation(const Variables& variable, - Vector& residual) -{ - for (ind_t node: m_mesh->range_nodes()) - { - nodal_residual_massbalance(node, variable, residual); - nodal_residual_mineral(node, variable, residual); - } -} - -void DiffusionProgram::nodal_residual_massbalance(ind_t node, - const Variables& variable, - Vector& residual) -{ - for (ind_t component: m_data->range_aqueous_component()) - { - ind_t id = m_ideq.id_equation_massbalance(node, component); - if (id != no_equation) - { - residual(id) += nodal_component_residual_massbalance(node, component, variable); - } - } -} - -double DiffusionProgram::nodal_component_residual_massbalance(ind_t node, - int component, - const Variables& variable) -{ - double sum = pow10(m_ideq.component_concentration(node, component, variable)); - for (ind_t aqueous: m_data->range_aqueous()) - { - if (m_data->nu_aqueous(aqueous, component) != 0) - { - sum += m_data->nu_aqueous(aqueous, component)*m_second.secondary_concentration(node, aqueous); - } - } - return sum - m_ideq.mobile_total_concentration(node, component, variable); -} - -void DiffusionProgram::nodal_residual_mineral(ind_t node, - const Variables& variable, - Vector& residual) -{ - for (ind_t mineral: m_data->range_mineral()) - { - ind_t id = m_ideq.id_equation_mineral(node, mineral); - if (id != no_equation) - { - - residual(id) = nodal_mineral_residual_mineral(node, mineral, variable); - m_is_micp[id] = true; - } - } -} - -double DiffusionProgram::nodal_mineral_residual_mineral(ind_t node, - int mineral, - const Variables& variable) -{ - double res = m_data->logk_mineral(mineral); - for (ind_t component: m_data->range_aqueous_component()) - { - res -= m_data->nu_mineral(mineral, component) * ( - m_ideq.component_concentration(node, component, variable) - + m_second.loggamma_component(node, component)); - } - return res; -} - -// ================================== // -// // -// Jacobian // -// // -// ================================== // - -void DiffusionProgram::get_jacobian(Variables &variable, - Vector &residual, - SparseMatrix &jacobian, - double alphadt) -{ - list_triplet_t jacob; - const ind_t ncomp = m_data->nb_component-1; - const ind_t nmin = m_data->nb_mineral; - const ind_t estimation = m_mesh->nb_nodes()*(ncomp*(3+ncomp) + ncomp*(ncomp+nmin)+nmin+ncomp*2); - jacob.reserve(estimation); - jacobian_transport(variable, residual, jacob, alphadt); - jacobian_speciation(variable, jacob, alphadt); - jacobian = SparseMatrix(get_neq(), get_neq()); - jacobian.setFromTriplets(jacob.begin(), jacob.end()); -} - -// Transport -// ========= - -void DiffusionProgram::jacobian_transport(Variables& variable, - Vector& residual, - list_triplet_t& jacobian, - double alphadt) -{ - for (ind_t element: m_mesh->range_elements()) - { - element_jacobian_transport(element, variable, residual, jacobian, alphadt); - } -} - - -void DiffusionProgram::element_jacobian_transport(ind_t element, - Variables& variable, - Vector& residual, - 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, variable, 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 ind_t idc = m_ideq.id_equation_diffusion(node, component); - const ind_t dof = m_ideq.get_dof_diffusion(node, component); - - if (idc == no_equation) continue; - - const double tmp_v = variable.velocity(dof); - const double tmp_d = variable.displacement(dof); - - double h = EPS_J*std::abs(tmp_v); - if (h == 0) h = EPS_J; - variable.velocity(dof) = tmp_v + h; - h = variable.velocity(dof) - tmp_v; - - variable.displacement(dof) = tmp_d + alphadt*h; - - element_residual_transport_component(element, component, variable, element_residual); - variable.velocity(dof) = tmp_v; - variable.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.id_equation_diffusion(noder, component); - - if (idr == no_equation) continue; - jacobian.push_back(triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h)); - } - // mineral -> not using finite difference - for (ind_t mineral: m_data->range_mineral()) - { - const ind_t idm = m_ideq.id_equation_mineral(node, mineral); - if (idm != no_equation) - jacobian.push_back(triplet_t(idc, idm, - -alphadt*m_mesh->get_volume_element(element)* - m_data->nu_mineral(mineral, component)/(2.0))); - } - } - - } -} - -// Speciation -// ========== - -void DiffusionProgram::jacobian_speciation(Variables& variable, - list_triplet_t& jacobian, - double alphadt) -{ - for (ind_t node: m_mesh->range_nodes()) - { - nodal_jacobian_speciation_fd(node, variable, jacobian, alphadt); - //nodal_jacobian_massbalance(node, variable, jacobian, alphadt); - //nodal_jacobian_mineral(node, variable, jacobian, alphadt); - - } -} - - -void DiffusionProgram::nodal_jacobian_speciation_fd( - ind_t node, - Variables& variable, - list_triplet_t& jacobian, - double alphadt) -{ - m_second.nodal_solve_secondary_variables(node, variable); - // Residuals - Eigen::VectorXd residual_massbalance_orig(m_data->nb_component); - for (int component: m_data->range_aqueous_component()) - { - //std::cout << nodal_component_residual_massbalance(node, component, variable)<< std::endl; - residual_massbalance_orig(component) = nodal_component_residual_massbalance(node, component, variable); - } - Eigen::VectorXd residual_mineral_orig(m_data->nb_mineral); - for (int mineral: m_data->range_mineral()) - { - residual_mineral_orig(mineral) = nodal_mineral_residual_mineral(node, mineral, variable); - } - - for (int component: m_data->range_aqueous_component()) - { - const ind_t id_c = m_ideq.id_equation_massbalance(node, component); - if (id_c == no_equation) continue; - const ind_t dof_c = m_ideq.get_dof_massbalance(node, component); - double tmp_v = variable.velocity(dof_c); - double tmp_d = variable.displacement(dof_c); - - double h = std::abs(tmp_v)*EPS_J; - if (h < EPS_J*1e-2) h = EPS_J; - variable.velocity(dof_c) += h; - h = variable.velocity(dof_c) - tmp_v; - variable.displacement(dof_c) = update_massbalance(variable.velocity(dof_c), tmp_d, alphadt); - - m_second.nodal_secondary_concentrations(node, variable); - //m_second.nodal_solve_secondary_variables(node, variable); - - // Aqueous mass balance equations - for (int k: m_data->range_aqueous_component()) - { - const ind_t id_k = m_ideq.id_equation_massbalance(node, k); - if (id_k == no_equation) continue; - const double residual = nodal_component_residual_massbalance(node, k, variable); - jacobian.push_back(triplet_t(id_k, id_c, - (residual-residual_massbalance_orig(k))/h)); - // contribution of total mobile concentration - if (k == component) - { - const ind_t id_tc = m_ideq.id_equation_diffusion(node, component); - if (id_tc != no_equation) - { - jacobian.push_back(triplet_t(id_c, id_tc, -alphadt)); - } - } - } - // Mineral stability - for (int mineral: m_data->range_mineral()) - { - const ind_t id_m = m_ideq.id_equation_mineral(node, mineral); - if (id_m == no_equation) continue; - double residual = nodal_mineral_residual_mineral(node, mineral, variable); - jacobian.push_back(triplet_t(id_m, id_c, - (residual-residual_mineral_orig(mineral))/h)); - } - - variable.velocity(dof_c) = tmp_v; - variable.displacement(dof_c) = tmp_d; - m_second.nodal_secondary_concentrations(node, variable); - //m_second.nodal_solve_secondary_variables(node, variable); - } - - // add dummy diagonal element for mineral => needed for reformulation - for (int mineral: m_data->range_mineral()) - { - const ind_t id_m = m_ideq.id_equation_mineral(node, mineral); - if (id_m == no_equation) continue; - jacobian.push_back(triplet_t(id_m, id_m, 0)); - } -} - - -void DiffusionProgram::nodal_jacobian_massbalance(ind_t node, - const Variables& variable, - list_triplet_t& jacobian, - double alphadt) -{ - const double logten = std::log(10.0); - for (ind_t component: m_data->range_aqueous_component()) - { - const ind_t idp = m_ideq.id_equation_massbalance(node, component); - if (idp == no_equation) continue; - // total concentration - const ind_t idc = m_ideq.id_equation_diffusion(node, component); - if (idc != no_equation) - { - jacobian.push_back(triplet_t(idp, idc, -alphadt)); - } - for (ind_t k: m_data->range_aqueous_component()) - { - // concentration - const ind_t ids = m_ideq.id_equation_massbalance(node, k); - if (ids == no_equation) continue; - double tmp_iip = 0; - if (k == component) tmp_iip += logten*pow10(m_ideq.component_concentration(node, component, variable)); - for (ind_t aqueous: m_data->range_aqueous()) - { - if (m_data->nu_aqueous(aqueous, component) == 0.0) continue; - tmp_iip += (logten*m_data->nu_aqueous(aqueous, component)* - m_data->nu_aqueous(aqueous, k)*m_second.secondary_concentration(node, aqueous)); - // - } - jacobian.push_back(triplet_t(idp, ids, alphadt*tmp_iip)); - } - } -} - -void DiffusionProgram::nodal_jacobian_mineral(ind_t node, - const Variables& variable, - list_triplet_t& jacobian, - double alphadt) -{ - for (ind_t mineral: m_data->range_mineral()) - { - const ind_t idm = m_ideq.id_equation_mineral(node, mineral); - if (idm == no_equation) continue; - for (ind_t component: m_data->range_aqueous_component()) - { - const ind_t idc = m_ideq.id_equation_massbalance(node, component); - if (idc == no_equation) continue; - jacobian.push_back(triplet_t(idm, idc, -alphadt*m_data->nu_mineral(mineral, component))); - } - // make place for an element - jacobian.push_back(triplet_t(idm, idm, 0)); - } -} - -// Variables -// ========= - -void DiffusionProgram::update_variables(ParabolicVariables& x, - const Eigen::VectorXd& predictor, - const Eigen::VectorXd& update, - double factor, - double alpha_dt - ) -{ -// for (ind_t node: m_mesh->range_nodes()) -// { -// for (ind_t edof=0; edofrange_nodes()) - { - for (int component: m_data->range_aqueous_component()) - { - const ind_t id_eq_d = m_ideq.id_equation_diffusion(node, component); - if (id_eq_d != no_equation) - { - const ind_t dof_d = m_ideq.get_dof_diffusion(node, component); - x.velocity(dof_d) += factor*update(id_eq_d); - x.displacement(dof_d) = update_diffusion(x.velocity(dof_d), predictor(dof_d), alpha_dt); - } - const ind_t id_eq_b = m_ideq.id_equation_massbalance(node, component); - if (id_eq_b != no_equation) - { - const ind_t dof_b = m_ideq.get_dof_massbalance(node, component); - x.velocity(dof_b) += factor*update(id_eq_b); - x.displacement(dof_b) = update_massbalance(x.velocity(dof_b), predictor(dof_b), alpha_dt); - } - } - for (int mineral: m_data->range_mineral()) - { - const ind_t id_eq_m = m_ideq.id_equation_mineral(node, mineral); - if (id_eq_m != no_equation) - { - const ind_t dof_m = m_ideq.get_dof_mineral(node, mineral); - x.velocity(dof_m) += factor*update(id_eq_m); - x.displacement(dof_m) = update_mineral(x.velocity(dof_m), predictor(dof_m), alpha_dt); - } - } - } - m_second.solve_secondary_variables(x); - -} - -void DiffusionProgram::set_predictor(ParabolicVariables& x, - Eigen::VectorXd& predictor, - double alpha, - double dt) -{ - predictor.resize(m_mesh->nb_nodes()*get_ndf()); - predictor = x.displacement + (1-alpha)*dt*x.velocity; - x.velocity.setZero(); - WARNING << "Predictor size : " << predictor.rows(); -} - -// compute secondary conc -bool DiffusionProgram::hook_start_iteration(const Variables& x, double norm_residual) { - return m_second.solve_secondary_variables(x); -} - - -double DiffusionProgram::maximum_lambda(const ParabolicVariables& x, const Eigen::VectorXd& update, double alphadt) -{ - double inv_maximum_lambda = 1.0; -// for (ind_t node: m_mesh->range_nodes()) -// { -// for (int mineral: m_data->range_mineral()) -// { -// if (m_ideq.id_equation_mineral(node, mineral) == no_equation -// or x.displacement(m_ideq.get_dof_mineral(node, mineral)) == 0) continue; -// inv_maximum_lambda = std::max(inv_maximum_lambda, -// -alphadt*update(m_ideq.id_equation_mineral(node, mineral))/ -// (0.9*x.displacement(m_ideq.get_dof_mineral(node, mineral)))); -// } -// } - return 1.0/inv_maximum_lambda; -} - -void DiffusionProgram::initialize_no_bc_with_init(const EquilibriumState& initialstate, - ParabolicVariables& var) -{ - for (ind_t node: m_mesh->range_nodes()) - { - if (not m_ideq.node_has_bc(node)) - { - for (ind_t component: m_data->range_aqueous_component()) - { - m_ideq.component_concentration(node, component, var) = - std::log10(initialstate.molality_component(component)); - m_ideq.mobile_total_concentration(node, component, var) = - initialstate.total_aqueous_concentration_component(component); - m_second.loggamma_component(node, component) = initialstate.loggamma_component(component); - } - for (ind_t aqueous: m_data->range_aqueous()) - { - m_second.secondary_concentration(node, aqueous) = initialstate.molality_secondary()[aqueous]; - m_second.loggamma_secondary(node, aqueous) = initialstate.loggamma_secondary(aqueous); - } - for (ind_t mineral: m_data->range_mineral()) - { - m_ideq.mole_mineral(node, mineral, var) = initialstate.moles_mineral(mineral)/m_mesh->get_volume_cell(1); - } - } - - } -} - -void DiffusionProgram::copy_cp_variables(const ParabolicVariables& x, - Eigen::VectorXd &moles_mineral) -{ - moles_mineral.resize(get_neq()); - moles_mineral.setZero(); - for (ind_t node: m_mesh->range_nodes()) - { - for (int mineral: m_data->range_mineral()) - { - const int id_eq = m_ideq.id_equation_mineral(node, mineral); - if (id_eq == no_equation) continue; - moles_mineral(id_eq) = m_ideq.mole_mineral(node, mineral, x); - } - } -} - -void DiffusionProgram::projection(ParabolicVariables& x) -{ - for (int node: m_mesh->range_nodes()) - { - for (int mineral: m_data->range_mineral()) - { - if (m_ideq.mole_mineral(node, mineral, x) < 1e-10) - { - m_ideq.mole_mineral(node, mineral, x) = 0.0; - } - } - } -} - -} // end namespace systems -} // end namespace reactmicp -} // end namespace specmicp diff --git a/src/reactmicp/systems/diffusion/diffusion.hpp b/src/reactmicp/systems/diffusion/diffusion.hpp deleted file mode 100644 index 5a610d9..0000000 --- a/src/reactmicp/systems/diffusion/diffusion.hpp +++ /dev/null @@ -1,306 +0,0 @@ -/*------------------------------------------------------- - - - Module : reactmicp/systems/diffusion - - File : diffusion.hpp - - Author : Fabien Georget - -Copyright (c) 2014, Fabien Georget , 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_SYSTEMS_DIFFUSION_HPP -#define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP - -#include -#include - -#include "database/data_container.hpp" - -#include "reactmicp/common.hpp" -#include "reactmicp/micpsolvers/parabolicprogram.hpp" -#include "reactmicp/meshes/mesh1d.hpp" - -#include "diffusion_numbering.hpp" -#include "diffusion_secondary.hpp" - - -namespace specmicp { -namespace reactmicp { -namespace systems { - -//! \brief Parameters for the diffusion system -class DiffusionParameter -{ -public: - - DiffusionParameter(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(ind_t element) {return m_diff;} - - //! \brief Porosity (at a node (function of composition)) - double porosity(ind_t node) {return m_poro;} - -private: - double m_diff; - double m_poro; -}; - -//! \brief Saturated diffusion in reactive porous media -class DiffusionProgram: public solvers::ParabolicProgram -{ - -public: - using Variables = ParabolicVariables; - const int no_equation = -1; - - DiffusionProgram(std::shared_ptr themesh, - std::shared_ptr thedatabase, - std::shared_ptr parameters, - diffusion::ListBoundaryCondition& list_bcs, - Variables& var): - m_mesh(themesh), m_data(thedatabase), m_param(parameters), - m_ideq(themesh->nb_nodes(), thedatabase, list_bcs, var), - m_second(themesh->nb_nodes(), thedatabase, m_ideq) - { - - } - - DiffusionProgram(std::shared_ptr themesh, - std::shared_ptr thedatabase, - std::shared_ptr parameters, - diffusion::ListBoundaryCondition& list_bcs, - const EquilibriumState &initialstate, - Variables& var); - - - //! \brief Return the number of degree of freedoms - ind_t get_neq() const {return m_ideq.get_neq();} - ind_t get_ndf() const {return m_ideq.get_ndf();} - - //! \brief Return the residuals - void get_residuals(const ParabolicVariables& variable, - Vector& residual); - - //! brief Return the jacobian - void get_jacobian(ParabolicVariables& variable, - Vector& residual, - SparseMatrix& jacobian, - double alphadt); - - //! \brief Return the index of the MiCP equations - const std::vector& micp_index() const { return m_is_micp;} - - - //! \brief Update the variables - void update_variables(ParabolicVariables& x, - const Eigen::VectorXd& predictor, - const Eigen::VectorXd& update, - double factor, - double alpha_dt - ); - //! \brief set the predictor - void set_predictor(ParabolicVariables& x, - Eigen::VectorXd& predictor, - double alpha, - double dt); - - //! \brief compute secondary concentrations - bool hook_start_iteration(const ParabolicVariables& x, double norm_residual); - - //! \brief Return the maximum lambda allowed for the linesearch - double maximum_lambda(const ParabolicVariables& x, const Eigen::VectorXd& update, double alphadt); - - //! \brief copy mineral amount from variables - void copy_cp_variables(const ParabolicVariables& x, - Eigen::VectorXd& moles_mineral); - - //! \brief projection step - void projection(ParabolicVariables& x); - - // todo or not........ -// double update_equation(int id_eq, -// double update, -// double predictor, -// double alphadt -// ) -// { -// // find _id eq .... -// } - -private: - - // Residuals - // ========= - - // transport - // --------- - - //! \brief Compute the residuals for the diffusion equations - void residual_transport(const Variables& variable, - Vector& residual); - - //! \brief Contribution of the diffusion operator to the residuals - void element_residual_transport(ind_t element, - const Variables& variable, - Vector& residual); - - //! \brief Contribution of the minerals to the residuals of the diffusion equations - double nodal_mineral_transient_term_transport(ind_t node, - ind_t component, - const Variables& variable); - - void element_residual_transport_component( - ind_t element, - int component, - const Variables& variable, - Vector& element_residual); - - // speciation - // ---------- - - //! \brief Compute the residual for the speciation problem - void residual_speciation(const Variables& variable, - Vector& residual); - - //! \brief Compute the residuals of the aqueous mass balances - void nodal_residual_massbalance(ind_t node, - const Variables& variable, - Vector& residual); - - //! \brief compute the residual for the aqueous mass balance of 'component' at 'node' - double nodal_component_residual_massbalance(ind_t node, - int component, - const Variables& variable); - - //! \brief Compute the residuals of the mineral speciation - void nodal_residual_mineral(ind_t node, - const Variables& variable, - Vector& residual); - - //! \brief compute the residual for the stability equation of 'mineral' at 'node' - double nodal_mineral_residual_mineral(ind_t node, - int mineral, - const Variables& variable); - - // Jacobian - // ======== - - // transport - // --------- - - //! \brief Contribution of the diffusion equation to the jacobian - void jacobian_transport(Variables& variable, - Vector& residual, - list_triplet_t& jacobian, - double alphadt); - - //! \brief Element by element assembly procedure - void element_jacobian_transport(ind_t element, - Variables& variable, - Vector& residual, - list_triplet_t& jacobian, - double alphadt); - // speciation - // ---------- - - //! \brief Contribution of the speciation problem to the mass balance - void jacobian_speciation(Variables &variable, - list_triplet_t& jacobian, - double alphadt); - - //! \brief Contribution of the mass balance at 'node' to the jacobian - void nodal_jacobian_massbalance(ind_t node, - const Variables& variable, - list_triplet_t& jacobian, - double alphadt); - - //! \brief Contribution of the mineral at equilibrium problem to the jacobian - void nodal_jacobian_mineral(ind_t node, - const Variables& variable, - list_triplet_t& jacobian, - double alphadt); - - //! \brief Finite difference jacobian for speciation - void nodal_jacobian_speciation_fd(ind_t node, - Variables &variable, - list_triplet_t& jacobian, - double alphadt); - - // Initialization - // ============== - - void initialize_no_bc_with_init(const EquilibriumState& initialstate, - ParabolicVariables& var); - - // Update - // ====== - - double update_diffusion(double update, - double predictor, - double alphadt) - { - return predictor + alphadt*update; - } - double update_massbalance(double update, - double predictor, - double alphadt) - { - return predictor + alphadt*update; - } - double update_mineral(double update, - double predictor, - double alphadt) - { - return predictor + alphadt*update; - } - - // Member - // ======= - - std::shared_ptr m_mesh; - std::shared_ptr m_data; - std::shared_ptr m_param; - - diffusion::DiffusionNumbering m_ideq; - diffusion::DiffusionSecondaryVariables m_second; - - std::vector m_is_micp; - - //SparseVector m_boundary_condition; - -}; - -} // end namespace systems -} // end namespace reactmicp -} // end namespace specmicp - -#endif // SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP diff --git a/src/reactmicp/systems/diffusion/diffusion_numbering.cpp b/src/reactmicp/systems/diffusion/diffusion_numbering.cpp deleted file mode 100644 index 144e0d5..0000000 --- a/src/reactmicp/systems/diffusion/diffusion_numbering.cpp +++ /dev/null @@ -1,167 +0,0 @@ -/*------------------------------------------------------- - - - Module : reactmicp/systems/diffusion - - File : diffusion_numbering.cpp - - Author : Fabien Georget - -Copyright (c) 2014, Fabien Georget , 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. - ----------------------------------------------------------*/ - -#include "diffusion_numbering.hpp" -#include "reactmicp/meshes/mesh1d.hpp" - -namespace specmicp { -namespace reactmicp { -namespace systems { -namespace diffusion { - - -// Boundary Condition - -void BoundaryCondition::fix_all() -{ - if (not composition.is_valid()) - throw std::runtime_error("The composition has not been correctly initialized"); - - std::shared_ptr data = composition.get_database(); - component.reserve(data->nb_component - 1); - tot_mobile_concentration.reserve(data->nb_component -1); - mineral.reserve(data->nb_mineral); - for (int comp: data->range_aqueous_component()) - { - component.push_back(comp); - tot_mobile_concentration.push_back(comp); - } - for (int min: data->range_mineral()) - { - mineral.push_back(min); - } -} - -// Diffusion numbering - -DiffusionNumbering::DiffusionNumbering(ind_t number_nodes, - ind_t number_aqueous_component, - ind_t number_mineral, const ListBoundaryCondition &list_bc, Variables &variable): - nb_aq_component(number_aqueous_component), nb_mineral(number_mineral), - nb_ndf(2*number_aqueous_component+number_mineral), - m_ideq(nb_ndf, number_nodes) -{ - number_equations(list_bc, variable); -} - - -void DiffusionNumbering::number_equations(const ListBoundaryCondition& list_bc, - Variables& variable) -{ - initialize_variable(variable); - // First we parse the bcs, to find the correct bc to apply - parse_bcs(list_bc, variable); - // Then we number the equations - do_numbering(); -} - - -void DiffusionNumbering::initialize_variable(Variables& variable) -{ - variable.displacement = Vector::Zero(m_ideq.rows()*m_ideq.cols()); - variable.velocity = Vector::Zero(m_ideq.rows()*m_ideq.cols()); -} - -void DiffusionNumbering::parse_bcs(const ListBoundaryCondition& list_bc, - Variables& variable) -{ - variable.displacement.setZero(); - variable.velocity.setZero(); - m_ideq.setZero(); - m_nodes_with_bc.reserve(list_bc.size()); - for (BoundaryCondition bc: list_bc) - { - const ind_t node = bc.node; - m_nodes_with_bc.push_back(node); - for (ind_t component: bc.tot_mobile_concentration) - { - m_ideq(get_ndof_diffusion(component), node) = no_equation; - mobile_total_concentration(node, component, variable) = bc.composition.total_aqueous_concentration_component(component); - } - for (ind_t component: bc.component) - { - m_ideq(get_ndof_massbalance(component), node) = no_equation; - component_concentration(node, component, variable) = std::log10(bc.composition.molality_component(component)); - } - for (ind_t mineral: bc.mineral) - { - m_ideq(get_ndof_mineral(mineral), node) = no_equation; - mole_mineral(node, mineral, variable) = bc.composition.moles_mineral(mineral); - } - } - - // - for (int node=0; node, 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_SYSTEMS_DIFFUSION_NUMBERING_HPP -#define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_NUMBERING_HPP - - -#include -#include "reactmicp/common.hpp" -#include "specmicp/equilibrium_data.hpp" - -namespace specmicp { -namespace reactmicp { - -namespace mesh { -// forward declaration -class Mesh1D; -} // end namespace mesh - -namespace systems { - -namespace diffusion { - -struct BoundaryCondition -{ - int node; - EquilibriumState composition; - std::vector tot_mobile_concentration; - std::vector component; - std::vector mineral; - - void fix_all(); -}; - -using ListBoundaryCondition = std::vector; - -//! \brief Equation numbering scheme -class DiffusionNumbering -{ -public: - const int no_equation = -1; - using Variables = ParabolicVariables; - - DiffusionNumbering(ind_t number_nodes, - std::shared_ptr data, - const ListBoundaryCondition& list_bc, - Variables& variable): - DiffusionNumbering(number_nodes, data->nb_component-1, data->nb_mineral, list_bc, variable) - {} - - DiffusionNumbering(ind_t number_nodes, - ind_t number_aqueous_component, - ind_t number_mineral, - const ListBoundaryCondition& list_bc, - Variables& variable); - - //! \brief Return the number of equations - ind_t get_neq() const {return nb_neq;} - //! \brief Return the number of degree of freedom (per node) - ind_t get_ndf() const {return nb_ndf;} - - // Degree of freedom numbering - // =========================== - - //! \brief Return the number of the dof - ind_t get_dof(ind_t node, ind_t elemental_dof) const { - return elemental_dof + nb_ndf*node;} - - //! \brief Return the ndof of a diffusion equation - ind_t get_ndof_diffusion(ind_t ind_component) const { - return ind_component-1;} - - //! \brief Return the ndof of a mass balance equation - ind_t get_ndof_massbalance(ind_t ind_component) const { - return nb_aq_component + ind_component - 1;} - - //! \brief Return the ndof of a mineral equation - ind_t get_ndof_mineral(ind_t ind_mineral) const { - return 2*nb_aq_component + ind_mineral;} - - //! \brief Return the dof of a diffusion equation - ind_t get_dof_diffusion(ind_t node, ind_t ind_component) const { - return get_dof(node, get_ndof_diffusion(ind_component));} - - //! \brief Return the dof of a mass balance equation - ind_t get_dof_massbalance(ind_t node, ind_t ind_component) const { - return get_dof(node, get_ndof_massbalance(ind_component));} - - //! \brief Return the dof of a mineral equation - ind_t get_dof_mineral(ind_t node, ind_t ind_mineral) const { - return get_dof(node, get_ndof_mineral(ind_mineral));} - - //! \brief Return the index of the equation - //! - //! Return no_equation if there is no equation for this degree of freedom - ind_t id_equation(ind_t node, ind_t elemental_dof) const {return m_ideq(elemental_dof, node);} - - //! \brief Return the equation's index of a diffusion equation - ind_t id_equation_diffusion(ind_t node, ind_t ind_component) const { - return id_equation(node, get_ndof_diffusion(ind_component));} - - //! \brief Return the equation's index of a mass balance equation - ind_t id_equation_massbalance(ind_t node, ind_t ind_component) const { - return id_equation(node, get_ndof_massbalance(ind_component));} - - //! \brief Return the equation's index of a mineral equation - ind_t id_equation_mineral(ind_t node, ind_t ind_mineral) const { - return id_equation(node, get_ndof_mineral(ind_mineral));} - - // Access primary variable - // ======================= - - //! \brief Return the total mobile concentration for 'component' at 'node' - double mobile_total_concentration(ind_t node, ind_t component, const Variables& var) const { - return var.displacement(get_dof_diffusion(node, component)); - } - //! \brief Return a reference to the total mobile concentration for 'component' at 'node' - double& mobile_total_concentration(ind_t node, ind_t component, Variables& var) { - return var.displacement.coeffRef(get_dof_diffusion(node, component)); - } - - - //! \brief Return the component concentration for 'component' at 'node' - double component_concentration(ind_t node, ind_t component, const Variables& var) const { - return var.displacement(get_dof_massbalance(node, component)); - } - //! \brief Return a reference to the component concentration for 'component' at 'node' - double& component_concentration(ind_t node, ind_t component, Variables& var) { - return var.displacement.coeffRef(get_dof_massbalance(node, component)); - } - - //! \brief Return the mineral mole number for 'mineral' at 'node' - double mole_mineral(ind_t node, ind_t mineral, const Variables& var) const { - return var.displacement(get_dof_mineral(node, mineral)); - } - //! \brief Return a reference to the mineral mole number for 'mineral' at 'node' - double& mole_mineral(ind_t node, ind_t mineral, Variables& var) { - return var.displacement.coeffRef(get_dof_mineral(node, mineral)); - } - - //! \brief Return true if the node has a boundary condition - double node_has_bc(ind_t node); - -private: - //! \brief Initialize the variable - void initialize_variable(Variables& variable); - - //! \brief Parse the BCs and number the equations - void number_equations(const ListBoundaryCondition& list_bc, - Variables& variable); - //! \brief Parse the BCS, set the corresponding values in variable - void parse_bcs(const ListBoundaryCondition& list_bc, - Variables& variable); - //! \brief number the equations - void do_numbering(); - - ind_t nb_aq_component; - ind_t nb_mineral; - - ind_t nb_ndf; - ind_t nb_neq; - - Eigen::MatrixXi m_ideq; - std::vector m_nodes_with_bc; -}; - -} // end namespace diffusion -} // end namespace systems -} // end namespace reactmicp -} // end namespace specmicp - -#endif // SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_NUMBERING_HPP diff --git a/src/reactmicp/systems/diffusion/diffusion_secondary.cpp b/src/reactmicp/systems/diffusion/diffusion_secondary.cpp deleted file mode 100644 index 107fabc..0000000 --- a/src/reactmicp/systems/diffusion/diffusion_secondary.cpp +++ /dev/null @@ -1,216 +0,0 @@ -/*------------------------------------------------------- - - - Module : reactmicp/systems/diffusion - - File : diffusion_secondary.cpp - - Author : Fabien Georget - -Copyright (c) 2014, Fabien Georget , 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. - ----------------------------------------------------------*/ - -#include "diffusion_secondary.hpp" -#include "physics/laws.hpp" - -namespace specmicp { -namespace reactmicp { -namespace systems { -namespace diffusion { - - -// ================================== // -// // -// Secondary variables wibbly-timey // -// // -// ================================== // - -//! \brief Set secondary concentrations -void DiffusionSecondaryVariables::compute_secondary_concentrations(const Variables& variable) -{ - for (ind_t node=0; noderange_aqueous()) - { - double sum = -m_data->logk_aqueous(species) - loggamma_secondary(node, species); - for (int component: m_data->range_aqueous_component()) - { - sum += m_data->nu_aqueous(species, component)*( - loggamma_component(node, component) - + m_ideq.component_concentration(node, component, variable) - ); - } - secondary_concentration(node, species) = pow10(sum); - } -} - -double DiffusionSecondaryVariables::norm_secondary_variables() -{ - return m_secondary_variables.block( - 1, 0, - m_data->nb_aqueous+m_data->nb_component, - m_secondary_concentration.cols() - ).norm(); -} - -double DiffusionSecondaryVariables::nodal_norm_secondary_variables(ind_t node) -{ - return m_secondary_variables.block( - 1, node, - m_data->nb_aqueous+m_data->nb_component, - 1 - ).norm(); -} - - -//! \brief Solve for secondary variables -int DiffusionSecondaryVariables::solve_secondary_variables(const Variables& variable) -{ - bool may_have_converged = false; - double previous_norm = norm_secondary_variables(); - for (int i=0; i<6; ++i) - { - compute_secondary_concentrations(variable); - compute_secondary_variables(variable); - double new_norm = norm_secondary_variables(); - if (std::abs(previous_norm - new_norm)/previous_norm < 1e-6) {may_have_converged=true; break;} - previous_norm = new_norm; - } - return may_have_converged; -} - - -//! \brief Solve for secondary variables -int DiffusionSecondaryVariables::nodal_solve_secondary_variables(ind_t node, const Variables& variable) -{ - bool may_have_converged = false; - double previous_norm = nodal_norm_secondary_variables(node); - for (int i=0; i<6; ++i) - { - nodal_secondary_concentrations(node, variable); - nodal_secondary_variables(node, variable); - double new_norm = nodal_norm_secondary_variables(node); - if (std::abs(previous_norm - new_norm)/previous_norm < 1e-6) {may_have_converged=true; break;} - previous_norm = new_norm; - } - return may_have_converged; -} - -//! \brief Compute secondary variables -void DiffusionSecondaryVariables::compute_secondary_variables(const Variables& variable) -{ - for (ind_t node=0; noderange_aqueous_component()) - { - if (m_data->param_aq(component, 0) == 0) continue; - ionics += std::pow(m_data->param_aq(component, 0), 2)*pow10( - m_ideq.component_concentration(node, component, variable)); - } - for (int aqueous: m_data->range_aqueous()) - { - ind_t idaq = m_data->nb_component + aqueous; - if (m_data->param_aq(idaq, 0) == 0) continue; - ionics += std::pow(m_data->param_aq(idaq, 0), 2)*secondary_concentration(node, aqueous); - } - ionic_strength(node) = ionics/2; -} - -void DiffusionSecondaryVariables::compute_loggamma() -{ - for (ind_t node=0; noderange_aqueous_component()) - { - loggamma_component(node, component) = laws::extended_debye_huckel(ionic_strength(node), sqrti, - m_data->param_aq(component, 0), - m_data->param_aq(component, 1), - m_data->param_aq(component, 2)); - } -} - -//! \brief compute the logarithm of the activity coefficients for node 'node' -void DiffusionSecondaryVariables::nodal_loggamma_aqueous(ind_t node) -{ - const double sqrti = std::sqrt(ionic_strength(node)); - for (int aqueous: m_data->range_aqueous()) - { - ind_t idaq = m_data->nb_component+aqueous; - loggamma_secondary(node, aqueous) = laws::extended_debye_huckel(ionic_strength(node), sqrti, - m_data->param_aq(idaq, 0), - m_data->param_aq(idaq, 1), - m_data->param_aq(idaq, 2)); - } -} - -} // end namespace diffusion -} // end namespace systems -} // end namespace reactmicp -} // end namespace specmicp diff --git a/src/reactmicp/systems/diffusion/diffusion_secondary.hpp b/src/reactmicp/systems/diffusion/diffusion_secondary.hpp deleted file mode 100644 index b5bd5f2..0000000 --- a/src/reactmicp/systems/diffusion/diffusion_secondary.hpp +++ /dev/null @@ -1,180 +0,0 @@ -/*------------------------------------------------------- - - - Module : reactmicp/systems/diffusion - - File : diffusion_secondary.hpp - - Author : Fabien Georget - -Copyright (c) 2014, Fabien Georget , 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_SYSTEMS_DIFFUSION_SECONDARY_HPP -#define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_SECONDARY_HPP - - -#include - -#include "database/data_container.hpp" - -#include "reactmicp/common.hpp" -#include "diffusion_numbering.hpp" - -namespace specmicp { -namespace reactmicp { -namespace systems { -namespace diffusion { - -//! \brief handle secondary variables for the diffusion system -class DiffusionSecondaryVariables -{ -public: - using Variables = ParabolicVariables; - - DiffusionSecondaryVariables(int nb_nodes, - std::shared_ptr thedatabase, - DiffusionNumbering& numberer): - m_data(thedatabase), m_ideq(numberer), - m_secondary_variables(Eigen::MatrixXd(thedatabase->nb_component+thedatabase->nb_aqueous+1, nb_nodes)), - m_secondary_concentration(Eigen::MatrixXd(thedatabase->nb_aqueous, nb_nodes)) - { - m_secondary_variables.setZero(); - m_secondary_concentration.setZero(); - } - - // Secondary variables - // =================== - - // Algorithms - // ---------- - - //! \brief Set secondary concentrations - void compute_secondary_concentrations(const Variables& variable); - - //! \brief Set secondary concentrations for node 'node' - void nodal_secondary_concentrations(ind_t node, const Variables& variable); - - //! \brief Solve for secondary variables - int solve_secondary_variables(const Variables& variable); - - //! \brief Solve for secondary variables in one node - int nodal_solve_secondary_variables(ind_t node, const Variables& variable); - - //! \brief Compute secondary variables - void compute_secondary_variables(const Variables& variable); - - //! \brief Compute the secondary variables at node 'node' - void nodal_secondary_variables(ind_t node, const Variables& variable); - - //! \brief compute the ionic strengths - void compute_ionic_strength(const Variables& variable); - - //! \brief compute the ionic strength for node 'node' - void nodal_ionic_strength(ind_t node, const Variables& variable); - - //! \brief compute the activity coefficients - void compute_loggamma(); - - //! \brief compute the logarithm of the activity coefficients for node 'node' - void nodal_loggamma(ind_t node); - - //! \brief compute the activity coefficient for the components of node 'node' - void nodal_loggamma_component(ind_t node); - - //! \brief compute the activity coefficient for the aqueous species of node 'node' - void nodal_loggamma_aqueous(ind_t node); - - // Getters and setters - // ------------------- - - //! \brief Secondary species concentration - double& secondary_concentration(ind_t node, ind_t aqueous) { - assert(aqueous >= 0 and aqueous < m_data->nb_aqueous); - return m_secondary_concentration.coeffRef(aqueous, node); - } - //! \brief Secondary species concentration - double secondary_concentration(ind_t node, ind_t aqueous) const { - assert(aqueous >= 0 and aqueous < m_data->nb_aqueous); - return m_secondary_concentration.coeff(aqueous, node); - } - - //! \brief Return the loggamma of a component - double& loggamma_component(ind_t node, ind_t component) { - assert(component >= 1 and component < m_data->nb_component); - return m_secondary_variables.coeffRef(component, node); - } - //! \brief Return the loggamma of a component - double loggamma_component(ind_t node, ind_t component) const { - assert(component >= 1 and component < m_data->nb_component); - return m_secondary_variables.coeff(component, node); - } - - //! \brief Return the loggamma of a secondary species - double& loggamma_secondary(ind_t node, ind_t aqueous) { - assert(aqueous >= 0 and aqueous < m_data->nb_aqueous); - return m_secondary_variables.coeffRef(m_data->nb_component + aqueous, node); - } - //! \brief Return the loggamma of a secondary species - double loggamma_secondary(ind_t node, ind_t aqueous) const { - assert(aqueous >= 0 and aqueous < m_data->nb_aqueous); - return m_secondary_variables.coeff(m_data->nb_component + aqueous, node); - } - - //! \brief Return ionic strength - double& ionic_strength(ind_t node) { - return m_secondary_variables.coeffRef(m_data->nb_component+m_data->nb_aqueous, node); - } - //! \brief Return ionic strength - double ionic_strength(ind_t node) const { - return m_secondary_variables.coeff(m_data->nb_component+m_data->nb_aqueous, node); - } - - - //! \brief Return norm of the loggamma - double loggamma_norm() const { - return m_secondary_variables.norm(); - } - -private: - // Return the norm used to check the convergence of the fixed-point iterations - double norm_secondary_variables(); - double nodal_norm_secondary_variables(ind_t node); - - std::shared_ptr m_data; - DiffusionNumbering& m_ideq; - - Eigen::MatrixXd m_secondary_variables; - Eigen::MatrixXd m_secondary_concentration; - -}; - - -} // end namespace diffusion -} // end namespace systems -} // end namespace reactmicp -} // end namespace specmicp - -#endif // SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_SECONDARY_HPP -