Homec4science
Diffusion SpecMiCP / ReactMiCP (reaction_path) Browse

specmicp/src/specmicp/reduced_system.cpp
a17e5379d078reaction_path

Last Change
Raw File

reduced_system.cpp
View Options

/*-------------------------------------------------------
- Module : specmicp
- File : specmicp.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include <iostream>
#include "reduced_system.hpp"
#include "utils/log.hpp"
namespace specmicp {
inline double pow10(double x)
{
return std::pow(10.0, x);
}
const double Adebye = 0.5092;
const double Bdebye = 0.3283;
void ReducedSystem::number_eq(bool water_equation)
{
int neq = 0;
m_ideq.reserve(m_thermo->nb_component()+m_thermo->nb_mineral());
if (water_equation)
{
m_ideq.push_back(neq);
++neq;
} else {
m_ideq.push_back(no_equation);
}
for (int i=1;i<m_thermo->nb_component();++i)
{
// Remve components that does not exist
if (m_tot_conc(i) == 0 and i!= 1)
{
INFO << "Component " << m_thermo->label_aqueous(i)
<< "has total concentration equal to zero, we desactivate it" << std::endl;
m_nonactive_component.push_back(i);
m_ideq.push_back(no_equation);
continue;
}
else
{
m_ideq.push_back(neq);
++neq;
}
}
m_nb_free_vars = neq;
for (int m=0; m<m_thermo->nb_mineral();++m)
{
bool can_precipitate = true;
// Remove minerals that cannot precipitate
for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it)
{
if (m_thermo->nu_min(m, *it) != 0)
{
can_precipitate = false;
break; // this is not a mineral that can precipitate
}
}
if (can_precipitate)
{
m_ideq.push_back(neq);
++neq;
}
else
{
m_ideq.push_back(no_equation);
}
}
m_nb_tot_vars = neq;
m_nb_compl_vars = m_nb_tot_vars - m_nb_free_vars;
// aqueous species
m_active_aqueous.reserve(m_thermo->nb_aqueous());
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
bool can_exist = true;
for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it)
{
if (m_thermo->nu_aq(j,*it) != 0)
{
can_exist = false;
}
}
m_active_aqueous.push_back(can_exist);
}
}
double ReducedSystem::residual_water(const Eigen::VectorXd& x) const
{
double res = m_tot_conc(0) - mass_water(x)/molar_mass_water;
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
res -= mass_water(x)*m_thermo->nu_aq(j, 0)*m_secondary_conc(j);
}
for (int m=0; m<m_thermo->nb_mineral(); ++m)
{
if (ideq_min(m) == no_equation) continue;
res -= m_thermo->nu_min(m, 0)*x(ideq_min(m));
}
return res;
}
double ReducedSystem::residual_component(const Eigen::VectorXd& x, int i) const
{
const int idp = ideq_paq(i);
double res = m_tot_conc(i) - mass_water(x)*pow10(x(idp));
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
if (not m_active_aqueous[j]) continue;
res -= mass_water(x)*m_thermo->nu_aq(j, i)*m_secondary_conc(j);
}
for (int m=0; m<m_thermo->nb_mineral(); ++m)
{
if (ideq_min(m) == no_equation) continue;
res -= m_thermo->nu_min(m, i)*x(ideq_min(m));
}
return res;
}
double ReducedSystem::residual_mineral(const Eigen::VectorXd& x, int m) const
{
double res = m_thermo->logkr_min(m);
for (int i=1; i<m_thermo->nb_component(); ++i)
{
if (m_thermo->nu_min(m, i) != 0)
res -= m_thermo->nu_min(m, i)*(x(ideq_paq(i)) + m_loggamma(i));
}
return res;
}
void ReducedSystem::get_residuals(const Eigen::VectorXd& x,
Eigen::VectorXd& residual)
{
set_secondary_concentration(x);
if (ideq_w() != no_equation) residual(ideq_w()) = residual_water(x);
for (int i=1; i<m_thermo->nb_component(); ++i)
{
if (ideq_paq(i) != no_equation) residual(ideq_paq(i)) = residual_component(x, i);
}
for (int m=0; m<m_thermo->nb_mineral(); ++m)
{
if (ideq_min(m) != no_equation) residual(ideq_min(m)) = residual_mineral(x, m);
}
}
//non-optimized Finite difference, for test only
//void ReducedSystem::get_jacobian(Eigen::VectorXd& x,
// Eigen::MatrixXd& jacobian)
//{
// const int neq = total_variables();
// Eigen::VectorXd res(total_variables());
// Eigen::VectorXd perturbed_res(total_variables());
// get_residuals(x, res);
// for (int j=0; j<neq; ++j)
// {
// double h = 1e-8*std::abs(x(j));
// //h = std::copysign(h, x(j));
// if (h==0) h = 1e-8;
// double tmp = x(j);
// x(j) += h;
// h = x(j) - tmp;
// get_residuals(x, perturbed_res);
// for (int i=0; i<neq; ++i)
// {
// jacobian(i, j) = (perturbed_res(i) - res(i))/h;
// }
// x(j) = tmp;
// }
// std::cout << jacobian << std::endl;
// return;
//}
void ReducedSystem::get_jacobian(Eigen::VectorXd& x,
Eigen::MatrixXd& jacobian)
{
const double log10 = std::log(10.0);
jacobian.setZero();
// water
if (ideq_w() != no_equation)
{
double tmp = -1/molar_mass_water;
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
if ( m_thermo->nu_aq(j, 0) == 0 ) continue;
tmp -= m_thermo->nu_aq(j, 0)*m_secondary_conc(j);
}
jacobian(0, 0) = tmp;
for (int k=1; k<m_thermo->nb_component(); ++k)
{
if (ideq_paq(k) == no_equation) continue;
double tmp = 0;
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
tmp -= m_thermo->nu_aq(j,0)*m_thermo->nu_aq(j,k)*m_secondary_conc(j);
}
jacobian(0, ideq_paq(k)) = x(ideq_w())*log10 * tmp;
}
for (int m=0; m<m_thermo->nb_mineral(); ++m)
{
if (ideq_min(m) == no_equation) continue;
jacobian(0, ideq_min(m)) = -m_thermo->nu_min(m, 0);
}
}
// aqueous component
for (int i=1; i<m_thermo->nb_component(); ++i)
{
if (ideq_paq(i) == no_equation) continue;
const int idp = ideq_paq(i);
// aqueous components
for (int k=1; k<m_thermo->nb_component(); ++k)
{
if (ideq_paq(k) == no_equation) continue;
double tmp_iip = 0;
if (k == i) tmp_iip -= pow10(x(ideq_paq(i)))*log10;
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
tmp_iip -= log10*m_thermo->nu_aq(j, i)*m_thermo->nu_aq(j, k)*m_secondary_conc(j);
}
jacobian(idp, ideq_paq(k)) = mass_water(x)*tmp_iip;
}
// minerals
for (int m=0; m<m_thermo->nb_mineral(); ++m)
{
if (ideq_min(m) == no_equation) continue;
jacobian(idp, ideq_min(m)) = - m_thermo->nu_min(m, i);
}
if (ideq_w() != no_equation) // Water
{
double tmp_iw = -std::pow(10.0, x(idp));
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
if ( m_thermo->nu_aq(j, i) == 0 ) continue;
tmp_iw -= m_thermo->nu_aq(j, i)*m_secondary_conc(j);
}
jacobian(idp, ideq_w()) = tmp_iw;
}
}
// mineral equilibrium
for (int m=0; m<m_thermo->nb_mineral(); ++m)
{
const int idm = ideq_min(m);
if (idm == no_equation) continue;
for (int i=1; i<m_thermo->nb_component(); ++i)
{
if (ideq_paq(i) == no_equation) continue;
jacobian(idm, ideq_paq(i)) = -m_thermo->nu_min(m, i);
}
}
//std::cout << jacobian << std::endl;
}
void ReducedSystem::set_secondary_concentration(const Eigen::VectorXd& x)
{
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
if (m_active_aqueous[j] == false)
{
m_secondary_conc(j) = 0;
continue;
}
double conc = -m_thermo->logkr_aq(j) - m_loggamma(j+m_thermo->nb_component());
for (int k=1; k<m_thermo->nb_component(); ++k)
{
if (m_thermo->nu_aq(j, k) == 0) continue;
conc += m_thermo->nu_aq(j, k)*(x(ideq_paq(k)) + m_loggamma(k));
}
m_secondary_conc(j) = pow10(conc);
}
}
void ReducedSystem::set_ionic_strength(const Eigen::VectorXd& x)
{
double ionic = 0;
for (int i=1; i<m_thermo->nb_component(); ++i)
{
if (ideq_paq(i) == no_equation or m_thermo->charge_paq(i) == 0) continue;
ionic += pow10(x(ideq_paq(i)))*std::pow(m_thermo->charge_paq(i),2);
}
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
if (not m_active_aqueous[j] or m_thermo->charge_saq(j) == 0) continue;
ionic += m_secondary_conc(j)*std::pow(m_thermo->charge_saq(j),2);
}
m_ionic_strength = ionic;
}
void ReducedSystem::compute_log_gamma(const Eigen::VectorXd& x)
{
set_ionic_strength(x);
const double sqrti = std::sqrt(m_ionic_strength);
for (int i=0; i<m_thermo->nb_component(); ++i)
{
if (ideq_paq(i) == no_equation)
{
m_loggamma(i) = 0;
continue;
}
double tmp = 0;
if (m_thermo->charge(i) != 0)
{
tmp = - Adebye*std::pow(2, m_thermo->charge(i))*sqrti;
tmp /= (1 + m_thermo->ao_debye(i)*Bdebye*sqrti);
}
tmp += m_thermo->bdot_debye(i)*m_ionic_strength;
m_loggamma(i) = tmp;
}
for (int j=0; j<m_thermo->nb_aqueous(); ++j)
{
if (not m_active_aqueous[j])
{
m_loggamma(m_thermo->nb_component() + j) = 0;
continue;
}
double tmp = 0;
if (m_thermo->charge_saq(j) != 0)
{
tmp = - Adebye*std::pow(2, m_thermo->charge_saq(j))*sqrti;
tmp /= (1 + m_thermo->ao_debye_saq(j)*Bdebye*sqrti);
}
tmp += m_thermo->bdot_debye_saq(j)*m_ionic_strength;
m_loggamma(m_thermo->nb_component() + j) = tmp;
}
}
void ReducedSystem::aqueous_tot_conc(const Eigen::VectorXd& x,
Eigen::VectorXd& aq_totconc)
{
assert(aq_totconc.rows() == m_thermo->nb_component());
set_secondary_concentration(x); // just to be sure
for (int i=1; i<m_thermo->nb_component();++i)
{
if (ideq_paq(i) != no_equation)
{
double tmp = pow10(x(ideq_paq(i)));
for (int j=0; j<m_thermo->nb_aqueous();++j)
{
if (m_thermo->nu_aq(j, i) == 0) continue;
tmp += m_thermo->nu_aq(j, i)*m_secondary_conc(j);
}
aq_totconc(i) = tmp;
}
else
{
aq_totconc(i) = 0;
}
}
}
double ReducedSystem::max_lambda(const Eigen::VectorXd& x, const Eigen::VectorXd& update)
{
//return 1.0/std::max(1.0,(-update.block(0,0,5,1).cwiseQuotient(0.95*x.block(0,0,5,1))).maxCoeff());
if (ideq_w() != no_equation)
{
// std::cout << "max lambda : " << 1.0/std::max(1.0, -update(0)/(0.95*x(0))) << std::endl;
return 1.0/std::max(1.0, -update(0)/(0.5*x(0)));
}
else
{
return 1.0;
}
}
} // end namespace specmicp
c4science · Help