Homec4science
Diffusion SpecMiCP / ReactMiCP (gasphase) Browse

specmicp/src/database/switch_basis.cpp
392b50f45f16gasphase

Last Change
Raw File

switch_basis.cpp
View Options

/*-------------------------------------------------------
- Module : database
- File : switch_basis.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "switch_basis.hpp"
#include <Eigen/Eigen>
#include "canonicalizer.hpp"
namespace specmicp {
namespace database {
void BasisSwitcher::change_basis(std::vector<int>& new_basis)
{
if (not is_database_canonical())
{
canonicalize_database(data);
}
Eigen::MatrixXd beta = Eigen::MatrixXd::Zero(data->nb_component, data->nb_component);
Eigen::VectorXd kswap(data->nb_component);
// swap coefficients, build beta matrix and kswap
for (int i=0; i<data->nb_component; ++i)
{
if (new_basis[i] >= data->nb_component)
{
int j = new_basis[i] - data->nb_component;
beta.row(i) = data->nu_aqueous.row(j);
data->nu_aqueous.row(j) = Eigen::VectorXd::Zero(data->nb_component);
data->nu_aqueous(j, i) = 1.0;
kswap(i) = data->logk_aqueous(j);
data->logk_aqueous(j) = 0.0;
}
else
{
beta(i, i) = 1.0;
kswap(i) = 0.0;
}
}
Eigen::MatrixXd betainv = beta.inverse();
data->molar_mass_basis = beta*data->molar_mass_basis;
data->nu_aqueous = data->nu_aqueous * betainv;
data->logk_aqueous -= data->nu_aqueous * kswap;
swap_aq_param(new_basis);
swap_labels(new_basis);
// now do the minerals
data->nu_mineral = data->nu_mineral*betainv;
data->logk_mineral -= data->nu_mineral * kswap;
data->nu_mineral_kinetic = data->nu_mineral_kinetic*betainv;
data->logk_mineral_kinetic -= data->nu_mineral_kinetic * kswap;
// and the gases
if (data->nb_gas != 0)
{
data->nu_gas = data->nu_gas*betainv;
data->logk_gas -= data->nu_gas * kswap;
}
}
void BasisSwitcher::swap_aq_param(std::vector<int> new_basis)
{
for (int i=0; i< data->nb_component; ++i)
{
if(new_basis[i] >= data->nb_component)
{
const int j = new_basis[i];
data->param_aq.row(i).swap(data->param_aq.row(j));
}
}
}
void BasisSwitcher::swap_labels(std::vector<int> new_basis)
{
for (int i=0; i<data->nb_component; ++i)
{
if (new_basis[i] >= data->nb_component)
{
std::swap(data->labels_basis[i], data->labels_aqueous[new_basis[i]-data->nb_component]);
}
}
}
} // end namespace database
} // end namespace specmicp
c4science · Help