Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/reactmicp/systems/saturated_diffusion/transport_neutrality_program.cpp
998871ce7287reactivetransport

Last Change
Raw File

transport_neutrality_program.cpp
View Options

#include "transport_neutrality_program.hpp"
#include "transport_neutrality_parameters.hpp"
#include "transport_neutrality_variables.hpp"
#include "dfpm/meshes/mesh1d.hpp"
#include <iostream>
#define EPS_J 1e-8
namespace specmicp {
namespace reactmicp {
namespace systems {
namespace siasaturated {
index_t SaturatedNeutralityDiffusionProgram::get_tot_ndf() const
{
return m_mesh->nb_nodes()*get_ndf();
}
SaturatedNeutralityDiffusionProgram::SaturatedNeutralityDiffusionProgram(
std::shared_ptr<mesh::Mesh1D> the_mesh,
std::shared_ptr<database::DataContainer> the_data,
std::shared_ptr<SaturatedNeutralityDiffusionTransportParameters> the_param
):
m_variables(the_mesh->nb_nodes(), the_data),
m_mesh(the_mesh),
m_data(the_data),
m_param(the_param),
m_external_flux(Vector::Zero(the_mesh->nb_nodes()*(the_data->nb_component-1)))
{}
void SaturatedNeutralityDiffusionProgram::number_equations(std::vector<index_t> list_fixed_nodes
)
{
m_ideq = Vector::Zero(get_ndf()*m_mesh->nb_nodes());
for (auto it=list_fixed_nodes.begin(); it<list_fixed_nodes.end(); ++it)
{
for (index_t component: m_data->range_aqueous_component())
{
m_ideq(get_dof_component(*it, component)) = no_equation;
}
}
index_t neq = 0;
for (index_t dof=0; dof<m_ideq.rows(); ++dof)
{
if (m_ideq(dof) != no_equation)
{
m_ideq(dof) = neq;
++neq;
}
}
m_neq = neq;
}
void SaturatedNeutralityDiffusionProgram::copy_displacements_in_variable(
const Vector& displacement
)
{
//std::cout << "New displacements : " << std::endl;
//std::cout << displacement << std::endl;
m_variables.component_concentrations() = displacement;
//std::cout << m_variables.component_concentrations() << std::endl;
bool is_ok = m_variables.solve_secondary_variables();
if (not is_ok)
{
std::cout << " concentrations " << std::endl;
for (index_t component: m_data->range_aqueous_component())
{
std::cout << m_variables.log_component_concentration(1, component) << std::endl;
}
std::cout << "Secondary concentrations : " << std::endl;
for (index_t aqueous: m_data->range_aqueous())
{
std::cout << m_variables.secondary_concentration(1, aqueous) << std::endl;
}
std::cout << "Log gamma : " << std::endl;
for (index_t component: m_data->range_aqueous_component())
{
std::cout << m_variables.log_gamma_component(1, component) << std::endl;
}
for (index_t aqueous: m_data->range_aqueous())
{
std::cout << m_variables.log_gamma_secondary(1, aqueous) << std::endl;
}
throw std::runtime_error("secondary variables not ok");
}
}
void SaturatedNeutralityDiffusionProgram::compute_residuals(
const Vector& displacement,
const Vector& velocity,
Vector& residual
)
{
copy_displacements_in_variable(displacement);
residual.resize(get_neq());
residual.setZero();
//residual = Vector::Zero(get_neq());
for (index_t element: m_mesh->range_elements())
{
Vector eresiduals;
element_residual(element, velocity, eresiduals);
for (index_t enode: m_mesh->range_nen())
{
for (index_t component: m_data->range_aqueous_component())
{
const index_t node = m_mesh->get_node(element, enode);
const index_t id = m_ideq(get_dof_component(node, component));
if (id == no_equation) continue;
residual(id) += eresiduals(get_edof_component(enode, component));
}
}
}
}
void SaturatedNeutralityDiffusionProgram::element_residual(
index_t element,
const Vector& velocity,
Vector& eresiduals
)
{
eresiduals = Vector::Zero(get_ndf()*2);
element_residual_mass_fd(element, velocity, eresiduals);
//std::cout << "Mass : " << std::endl << eresiduals << std::endl;
element_residual_flux(element, eresiduals);
//std::cout << "Flux : " << std::endl << eresiduals << std::endl;
element_residual_source(element, eresiduals);
}
void SaturatedNeutralityDiffusionProgram::element_residual_source(
index_t element,
Vector& element_residual
)
{
for (index_t enode: m_mesh->range_nen())
{
index_t node = m_mesh->get_node(element, enode);
const scalar_t coeff_integ = ( m_mesh->get_volume_element(element) ) /2;
for (index_t component: m_data->range_aqueous_component())
{
element_residual(get_edof_component(enode, component)) +=
coeff_integ*m_external_flux(get_dof_component(node, component));
}
}
}
void SaturatedNeutralityDiffusionProgram::element_residual_flux(
index_t element,
Vector& element_residual
)
{
// Litchner (1996) pg 55-56
//
// Goal : compute the flux for each component (i)
//
// J_i = ψe_i ( Σ_l z_l ( Γ^d_l + Γ^γ_l ) )/( Σ_l z_l ψe_l ) - ( Γ^d_i + Γ^γ_i )
// \_________________________________________________/
// \/
// electroneutrality
//
Vector flux = Vector::Zero(m_data->nb_component); // Γ^d_i + Γ^γ_i
Vector potential_factor = Vector::Zero(m_data->nb_component); // ψe_i
for (index_t i: m_data->range_aqueous_component())
{
// flux : Γ^d_i
element_residual_transport_component_diffusion(element, i, flux);
// flux : Γ^d_i + Γ^γ_i
// potential_factor : ψe_i
element_residual_transport_component_nonideal(element, i, flux, potential_factor);
}
scalar_t sum_potential = 0.0; // Σ_l z_l ψe_l
scalar_t sum_flux_charge = 0.0; // Σ_l z_l ( Γ^d_i + Γ^γ_i )
for (index_t l: m_data->range_aqueous_component())
{
if (m_data->charge_component(l) == 0.0 ) continue;
sum_potential += m_data->charge_component(l)*potential_factor(l);
sum_flux_charge += m_data->charge_component(l)*flux(l);
}
// ( Σ_l z_l ( Γ^d_l + Γ^γ_l ) )/( Σ_l z_l ψe_l )
const scalar_t grad_pot = sum_flux_charge/sum_potential;
// integration
const index_t node_0 = m_mesh->get_node(element, 0);
const index_t node_1 = m_mesh->get_node(element, 1);
const scalar_t porosity = (m_param->porosity(node_0) +
m_param->porosity(node_1) )/2;
//const scalar_t porosity = 1.0/(0.5/m_param->porosity(node_0) +
// 0.5/m_param->porosity(node_1) );
const scalar_t coeff_integ = ( m_mesh->get_face_area(element)
* porosity
* m_param->density_water()
);
for (index_t component: m_data->range_aqueous_component())
{
const scalar_t total_flux = grad_pot*potential_factor(component) - flux(component); // J_i
//const scalar_t total_flux = - flux(component);
element_residual(get_edof_component(0, component)) += coeff_integ*total_flux; // ∫ J_i.n dS
element_residual(get_edof_component(1, component)) -= coeff_integ*total_flux; // ∫ J_i.n dS
}
}
void SaturatedNeutralityDiffusionProgram::element_residual_mass_fd(
index_t element,
const Vector& velocity,
Vector& element_residual
)
{
for (index_t enode: m_mesh->range_nen())
{
index_t node = m_mesh->get_node(element, enode);
const scalar_t integ_coeff = -( m_param->density_water()
* m_param->porosity(node)
* m_mesh->get_volume_element(element)/ 2.0
);
Vector totconc(m_data->nb_component);
m_variables.nodal_solve_secondary_variables(node);
for (index_t i: m_data->range_aqueous_component())
{
totconc(i) = m_variables.total_mobile_concentration(node, i);
}
for (index_t component: m_data->range_aqueous_component())
{
if (m_ideq(get_dof_component(node, component)) == no_equation) continue;
const scalar_t tmp_v = velocity(get_dof_component(node, component));
const scalar_t tmp_d = m_variables.log_component_concentration(node, component);
if (tmp_v == 0.0) continue;
scalar_t h = 1e-6*std::abs(tmp_d/tmp_v);
if (h == 0.0 or not std::isfinite(h)) h = 1e-2;
m_variables.log_component_concentration(node, component) = tmp_d + h*tmp_v;
h = (m_variables.log_component_concentration(node, component) - tmp_d)/tmp_v;
// std::cout << "before : " << tmp_d << " - after : "
// << m_variables.log_component_concentration(node, component)
// << " tmp_v " << tmp_v
// << " - h " << h << std::endl;
m_variables.nodal_solve_secondary_variables(node);
Vector totconc_perturb(m_data->nb_component);
for (index_t i: m_data->range_aqueous_component())
{
totconc_perturb(i) = m_variables.total_mobile_concentration(node, i);
}
m_variables.log_component_concentration(node, component) = tmp_d;
m_variables.nodal_solve_secondary_variables(node);
for (index_t k: m_data->range_aqueous_component())
{
if (m_ideq(get_dof_component(node, k)) == no_equation) continue;
//std::cout << "val : " << integ_coeff*(totconc_perturb(k) - totconc(k))/h << std::endl;
element_residual(get_edof_component(enode, k)) += integ_coeff*(totconc_perturb(k) - totconc(k))/h;
}
}
}
}
void SaturatedNeutralityDiffusionProgram::element_residual_transport_component_diffusion(
index_t element,
index_t component,
Vector& flux
)
{
// Litchner (1996)
//
// Diffusion for each component :
//
// Γ^d_i = D_i ∇ b_i + Σ_j ν_ji D_j ∇ b_j
//
const index_t node_0 = m_mesh->get_node(element, 0);
const index_t node_1 = m_mesh->get_node(element, 1);
// Component
{
const scalar_t grad_c = ( m_variables.component_concentration(node_0, component)
- m_variables.component_concentration(node_1, component)
) / m_mesh->get_dx(element);
const scalar_t diff_coeff = 1.0/(0.5/m_param->diffusion_component(node_0, component)
+0.5/m_param->diffusion_component(node_1, component));
const scalar_t flux_diffusion = diff_coeff*grad_c;
//std::cout << "flux c " << flux_diffusion << std::endl;
flux(component) += flux_diffusion;
}
// Secondary species
for (index_t sec: m_data->range_aqueous())
{
if (m_data->nu_aqueous(sec, component) == 0.0) continue;
const scalar_t grad_c = ( m_variables.secondary_concentration(node_0, sec)
- m_variables.secondary_concentration(node_1, sec)
) / m_mesh->get_dx(element);
const scalar_t diff_coeff = 1.0/(0.5/m_param->diffusion_secondary(node_0, sec)
+0.5/m_param->diffusion_secondary(node_1, sec));
const scalar_t flux_diffusion = m_data->nu_aqueous(sec, component)*diff_coeff*grad_c;
//std::cout << "flux s " << flux_diffusion << std::endl;
flux(component) += flux_diffusion;
}
}
void SaturatedNeutralityDiffusionProgram::element_residual_transport_component_nonideal(
index_t element,
index_t component,
Vector& flux,
Vector& potential_factor
)
{
// Lichtner (1996)
//
// Flux due to non-ideality :
//
// Γ^γ_i = D_i b_i ∇ log(γ_i) + Σ_j ν_ji D_j b_j ∇ log(γ_j)
//
// Contribution of the component to the potential :
//
// ψe_i = z_i D_i b_i + Σ_j ν_ji z_j D_j b_j
//
const index_t node_0 = m_mesh->get_node(element, 0);
const index_t node_1 = m_mesh->get_node(element, 1);
// Component
{
const scalar_t grad_log_gamma = ( m_variables.log_gamma_component(node_0, component)
- m_variables.log_gamma_component(node_1, component)
) / m_mesh->get_dx(element);
const scalar_t diff_coeff = 1.0/(0.5/m_param->diffusion_component(node_0, component)
+0.5/m_param->diffusion_component(node_1, component));
const scalar_t average_conc = ( m_variables.component_concentration(node_0, component)
+ m_variables.component_concentration(node_1, component)
) / 2.0;
const scalar_t flux_nonideal = diff_coeff*average_conc*grad_log_gamma;
flux(component) += flux_nonideal;
if (m_data->charge_component(component) != 0.0)
potential_factor(component) += m_data->charge_component(component)*diff_coeff*average_conc;
}
// Secondary species
for (index_t sec: m_data->range_aqueous())
{
if (m_data->nu_aqueous(sec, component) == 0.0) continue;
const scalar_t grad_log_gamma = ( m_variables.log_gamma_secondary(node_0, sec)
- m_variables.log_gamma_secondary(node_1, sec)
) / m_mesh->get_dx(element);
const scalar_t diff_coeff = 1.0/( 0.5/m_param->diffusion_secondary(node_0, sec)
+ 0.5/m_param->diffusion_secondary(node_1, sec));
const scalar_t average_conc = ( m_variables.secondary_concentration(node_0, sec)
+ m_variables.secondary_concentration(node_1, sec)
) / 2.0;
const scalar_t flux_nonideal = m_data->nu_aqueous(sec, component)
*diff_coeff*average_conc*grad_log_gamma;
flux(component) += flux_nonideal;
if (m_data->charge_aqueous(sec) == 0.0) continue;
potential_factor(component) += m_data->nu_aqueous(sec, component)
*m_data->charge_aqueous(sec)*diff_coeff*average_conc;
}
}
void SaturatedNeutralityDiffusionProgram::element_residual_mass(
index_t element,
const Vector& velocity,
Vector& element_residual
)
{
constexpr scalar_t logten = std::log(10.0);
const scalar_t integ_coeff = -(m_param->density_water() *
m_mesh->get_volume_element(element)/2.0);
for (index_t enode: m_mesh->range_nen())
{
index_t node = m_mesh->get_node(element, enode);
// dB_c/dt = db_c/dt + Σ_j ν_jc db_j/dt
// dB_c/dt = ln(10)*b_c*d(log(b_c))/dt + Σ_j ν_jc db_j/dt
for (index_t component: m_data->range_aqueous_component())
{
if (get_dof_component(node, component) == no_equation) continue;
scalar_t tot_vel_c = logten*m_variables.component_concentration(node, component)* velocity(get_dof_component(node, component));
//
// db_j/dt = 1/(K_j γ_j)* Σ_i [Π_k ( γ_k b_k )^{ν_jk} ] ν_ji ln(10) * d(log(b_i))/dt
//
for (index_t aqueous: m_data->range_aqueous())
{
if (m_data->nu_aqueous(aqueous, component) == 0.0 ) continue;
for (index_t i: m_data->range_aqueous_component())
{
if (m_data->nu_aqueous(aqueous, i) == 0.0 ) continue;
scalar_t product = m_data->nu_aqueous(aqueous, i)/ (
m_data->logk_aqueous(aqueous)
*pow10(m_variables.log_gamma_secondary(node, aqueous))
);
for (index_t k: m_data->range_aqueous_component())
{
if (m_data->nu_aqueous(aqueous, k) == 0.0) continue;
product *= pow10(m_data->nu_aqueous(aqueous, k)*
(m_variables.log_gamma_component(node, k)
+ m_variables.log_component_concentration(node, k)
));
}
product *= logten; ///m_variables.component_concentration(node, i);
tot_vel_c += m_data->nu_aqueous(aqueous, component)*product*velocity(get_dof_component(node, i));
}
}
element_residual(get_edof_component(enode, component)) += m_param->porosity(node)*integ_coeff*tot_vel_c;
}
}
}
void SaturatedNeutralityDiffusionProgram::compute_jacobian(
Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian,
scalar_t alphadt
)
{
copy_displacements_in_variable(displacement);
list_triplet_t jacob;
const index_t ncomp = m_data->nb_component-1;
const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen);
jacob.reserve(estimation);
for (index_t element: m_mesh->range_elements())
{
element_jacobian_transport(element, velocity, jacob, alphadt);
}
jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
jacobian.setFromTriplets(jacob.begin(), jacob.end());
}
// Transport
// =========
void SaturatedNeutralityDiffusionProgram::element_jacobian_transport(
index_t element,
Vector& velocity,
list_triplet_t& jacobian,
scalar_t alphadt)
{
for (index_t component: m_data->range_aqueous_component())
{
Vector element_residual_orig(Vector::Zero(2*get_ndf()));
element_residual(element, velocity, element_residual_orig);
for (index_t en: m_mesh->range_nen())
{
Vector element_residual_perturb(Vector::Zero(2*get_ndf()));
const index_t node = m_mesh->get_node(element, en);
const index_t dof = get_dof_component(node, component);
const index_t idc = m_ideq(dof);
if (idc == no_equation) continue;
const scalar_t tmp_v = velocity(dof);
const scalar_t tmp_d = m_variables.log_component_concentration(node, component);
scalar_t h = 0.0;
if (alphadt !=0) h = 1e-4/alphadt; //(1e-8)*(alphadt)*std::abs(tmp_v);
if (h == 0) h = EPS_J;
velocity(dof) = tmp_v + h;
h = velocity(dof) - tmp_v;
m_variables.log_component_concentration(node, component) = tmp_d + alphadt*h;
m_variables.nodal_solve_secondary_variables(node);
element_residual(element, velocity, element_residual_perturb);
velocity(dof) = tmp_v;
m_variables.log_component_concentration(node, component) = tmp_d;
m_variables.nodal_solve_secondary_variables(node);
for (index_t k: m_data->range_aqueous_component())
{
for (index_t enr: m_mesh->range_nen())
{
const index_t noder = m_mesh->get_node(element, enr);
const index_t idr = m_ideq(get_dof_component(noder, k));
if (idr == no_equation) continue;
const index_t edof = get_edof_component(enr, k);
jacobian.push_back(triplet_t(idr, idc,
(element_residual_perturb(edof) - element_residual_orig(edof))/h));
}
}
}
}
}
void SaturatedNeutralityDiffusionProgram::update_solution(const Vector& update,
scalar_t lambda,
scalar_t alpha_dt,
Vector& predictor,
Vector& displacement,
Vector& velocity
)
{
//std::cout << update << std::endl;
for (index_t node: m_mesh->range_nodes())
{
for (index_t component: m_data->range_aqueous_component())
{
const index_t dof = get_dof_component(node, component);
const index_t id = m_ideq(dof);
if (id == no_equation) continue;
velocity(dof) += lambda*update(id);
}
}
displacement = predictor + alpha_dt*velocity;
copy_displacements_in_variable(displacement);
}
void SaturatedNeutralityDiffusionProgram::get_total_mobile_concentrations(
const Vector& displacement,
Vector& total_concentrations)
{
copy_displacements_in_variable(displacement);
total_concentrations.resize(displacement.rows());
for (index_t node: m_mesh->range_nodes())
{
for (index_t component: m_data->range_aqueous_component())
{
total_concentrations(get_dof_component(node, component)) = m_variables.total_mobile_concentration(node, component);
}
}
}
void SaturatedNeutralityDiffusionProgram::extern_update_nodal_composition(
index_t node,
const EquilibriumState& composition
)
{
for (index_t component: m_data->range_aqueous_component())
{
m_variables.log_component_concentration(node, component) = composition.log_molality_component(component);
m_variables.log_gamma_component(node, component) = composition.loggamma_component(component);
}
for (index_t aqueous: m_data->range_aqueous())
{
m_variables.secondary_concentration(node, aqueous) = composition.loggamma_secondary(aqueous);
}
m_variables.nodal_solve_secondary_variables(node);
}
void SaturatedNeutralityDiffusionProgram::extern_update_nodal_composition(
index_t node,
const EquilibriumState& composition,
Vector& displacement
)
{
for (index_t component: m_data->range_aqueous_component())
{
displacement(get_dof_component(node, component)) = composition.log_molality_component(component);
m_variables.log_component_concentration(node, component) = composition.log_molality_component(component);
m_variables.log_gamma_component(node, component) = composition.loggamma_component(component);
}
for (index_t aqueous: m_data->range_aqueous())
{
m_variables.secondary_concentration(node, aqueous) = composition.loggamma_secondary(aqueous);
}
m_variables.nodal_solve_secondary_variables(node);
}
} // end namespace siasaturated
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help