Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/reactmicp/systems/diffusion/diffusion.cpp
68b7dab0dd3dreactmicp/unsaturated

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>
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) += nodal_mineral_transient_term_transport(
m_mesh->get_node(element, en), component, variable)/2;
}
}
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))
);
}
//mineral_term *= 1.0/2.0; //(2.0/m_mesh->cell_volume(node));
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)
{
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);
}
}
residual(id) += sum - m_ideq.mobile_total_concentration(node, component, variable);///m_param->density_water();
}
}
}
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)
{
double sum = m_data->logk_mineral(mineral);
for (ind_t component: m_data->range_aqueous_component())
{
sum -= m_data->nu_mineral(mineral, component) * (
m_ideq.component_concentration(node, component, variable)
+ m_second.loggamma_component(node, component));
}
residual(id) = sum;
m_is_micp[id] = true;
}
}
}
// ================================== //
// //
// 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_fd(element, variable, residual, jacobian, alphadt);
}
}
void DiffusionProgram::element_jacobian_transport_fd(ind_t element,
Variables& variable,
Vector& residual,
list_triplet_t& jacobian,
double alphadt)
{
const double eps = 1e-6;
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);
//ind_t idl = m_ideq.id_equation_diffusion(m_mesh->get_node(element, 0), component);
//ind_t idr = m_ideq.id_equation_diffusion(m_mesh->get_node(element, 1), component);
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*std::abs(tmp_v);
if (h == 0) h = eps;
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, -m_data->nu_mineral(mineral, component)/2.0));
}
}
}
}
void DiffusionProgram::element_jacobian_transport(ind_t element,
const Variables& variable,
Vector& residual,
list_triplet_t& jacobian,
double alphadt)
{
double deriv = -alphadt*( m_mesh->crosssection() /m_mesh->get_dx(element)
* m_param->diffusion_coefficient(element)*m_param->porosity(0)
* m_param->density_water());
for (ind_t component: m_data->range_aqueous_component())
{
ind_t idl = m_ideq.id_equation_diffusion(m_mesh->get_node(element, 0), component);
ind_t idr = m_ideq.id_equation_diffusion(m_mesh->get_node(element, 1), component);
if (idl != no_equation)
{
jacobian.push_back(triplet_t(idl, idl, deriv));
if (idr != no_equation) jacobian.push_back(triplet_t(idl, idr, -deriv));
}
if (idr != no_equation)
{
jacobian.push_back(triplet_t(idr, idr, deriv));
if (idl != no_equation) jacobian.push_back(triplet_t(idr, idl, -deriv));
}
// Mass Matrix
double mass_coeff = -(m_param->porosity(m_mesh->get_node(element, 0)));
mass_coeff *= m_param->density_water()*m_mesh->crosssection()*m_mesh->get_dx(element)/2;
if (idl != no_equation) jacobian.push_back(triplet_t(idl, idl, mass_coeff));
if (idr != no_equation) jacobian.push_back(triplet_t(idr, idr, mass_coeff));
// mineral
// for (ind_t mineral: m_data->range_mineral())
// {
// ind_t idm_0 = m_ideq.id_equation_mineral(m_mesh->get_node(element, 0), mineral);
// ind_t idm_1 = m_ideq.id_equation_mineral(m_mesh->get_node(element, 1), mineral);
// if (idm_0 != no_equation and idl != no_equation)
// jacobian.push_back(triplet_t(idl, idm_0,
// -m_data->nu_mineral(mineral, component)/2.0));
// if (idm_1 != no_equation and idr != no_equation)
// jacobian.push_back(triplet_t(idr, idm_1,
// -m_data->nu_mineral(mineral, component)/2.0));
// }
}
}
// Speciation
// ==========
void DiffusionProgram::jacobian_speciation(const Variables& variable,
list_triplet_t& jacobian,
double alphadt)
{
for (ind_t node: m_mesh->range_nodes())
{
nodal_jacobian_massbalance(node, variable, jacobian, alphadt);
nodal_jacobian_mineral(node, variable, jacobian, alphadt);
}
}
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)); // /m_param->density_water()));
}
for (ind_t k: m_data->range_aqueous_component())
{
// concentration
const ind_t ids = m_ideq.id_equation_massbalance(node, component);
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 j: m_data->range_aqueous())
{
tmp_iip += (logten*m_data->nu_aqueous(j, component)*
m_data->nu_aqueous(j, k)*m_second.secondary_concentration(node, j));
}
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);
}
}
}
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;
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);
}
}
}
}
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);
}
}
}
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help