Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/reactmicp/systems/diffusion/diffusion.cpp
c46a22fd727breactivetransport

Last Change
Raw File

diffusion.cpp
View Options

/*-------------------------------------------------------
- Module : reactmicp/systems/diffusion
- File : diffusion.cpp
- 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.
---------------------------------------------------------*/
#include "diffusion.hpp"
#include "physics/laws.hpp"
#include "utils/log.hpp"
#include <iostream>
#define EPS_J 1e-8
namespace specmicp {
namespace reactmicp {
namespace systems {
DiffusionProgram::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):
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<bool>(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->crosssection()*m_mesh->get_dx(element)/2);
mass << 1, 0, 0, 1;
mass *= mass_coeff;
double flux_coeff = -( m_mesh->crosssection() / m_mesh->get_dx(element)
* m_param->porosity(0)
* m_param->diffusion_coefficient(element)
* m_param->density_water()
);
jacob << 1, -1, -1, 1;
jacob *= flux_coeff;
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->crosssection()*m_mesh->get_dx(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->nnodes()*(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->crosssection()*m_mesh->get_dx(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; edof<get_ndf(); ++edof)
// {
// if (m_ideq.id_equation(node, edof) == no_equation) continue;
// const ind_t dof = m_ideq.get_dof(node, edof);
// const ind_t id_eq = m_ideq.id_equation(node, edof);
// x.velocity(dof) += factor*update(id_eq);
// x.displacement(dof) = predictor(dof) + alpha_dt * x.velocity(dof);
// }
// }
for (ind_t node: m_mesh->range_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->nnodes()*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->cell_volume(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
c4science · Help