Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F120232646
reduced_system.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jul 2, 20:04
Size
11 KB
Mime Type
text/x-c
Expires
Fri, Jul 4, 20:04 (2 d)
Engine
blob
Format
Raw Data
Handle
27158952
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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
Event Timeline
Log In to Comment