Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/reactmicp/systems/diffusion/diffusion.hpp
7009d71a1296reactivetransport

Last Change
Raw File

diffusion.hpp
View Options

/*-------------------------------------------------------
- Module : reactmicp/systems/diffusion
- File : diffusion.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#ifndef SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
#define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
#include <vector>
#include <memory>
#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<DiffusionProgram>
{
public:
using Variables = ParabolicVariables;
const int no_equation = -1;
DiffusionProgram(std::shared_ptr<mesh::Mesh1D> themesh,
std::shared_ptr<database::DataContainer> thedatabase,
std::shared_ptr<DiffusionParameter> 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<mesh::Mesh1D> themesh,
std::shared_ptr<database::DataContainer> thedatabase,
std::shared_ptr<DiffusionParameter> 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<bool>& 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<mesh::Mesh1D> m_mesh;
std::shared_ptr<database::DataContainer> m_data;
std::shared_ptr<DiffusionParameter> m_param;
diffusion::DiffusionNumbering m_ideq;
diffusion::DiffusionSecondaryVariables m_second;
std::vector<bool> m_is_micp;
//SparseVector m_boundary_condition;
};
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
#endif // SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
c4science · Help