Homec4science
Diffusion SpecMiCP / ReactMiCP (database) Browse

specmicp/src/database/database.cpp
4d14a2cc5b7edatabase

Last Change
Raw File

database.cpp
View Options

/*-------------------------------------------------------
- Module : database
- File : database.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "database.hpp"
#include <Eigen/Eigen>
#include <iostream>
namespace specmicp {
namespace database {
void Database::make_aqueous_canonical()
{
for (int i=1; i<data->nb_aqueous; ++i)
{
for (int j=0; j<i; ++j)
{
const double stoech = data->nu_aqueous(i, data->nb_component+j);
if (stoech != 0)
{
data->nu_aqueous.row(i) += stoech*data->nu_aqueous.row(j);
data->logk_aqueous(i) += stoech*data->logk_aqueous(j);
data->nu_aqueous(i, data->nb_component+j) = 0;
}
}
}
data->nu_aqueous.conservativeResize(Eigen::NoChange_t(), data->nb_component);
}
void Database::make_mineral_canonical()
{
for (int i=1; i<data->nb_mineral; ++i)
{
for (int j=0; j<data->nb_aqueous; ++j)
{
const double stoech = data->nu_mineral(i, data->nb_component+j);
if (stoech != 0)
{
data->nu_mineral.block(i, 0, 1, data->nb_component) +=
stoech*data->nu_aqueous.row(j);
data->logk_mineral(i) += stoech*data->logk_aqueous(j);
data->nu_mineral(i, data->nb_component+j) = 0;
}
}
}
data->nu_mineral.conservativeResize(Eigen::NoChange_t(), data->nb_component);
}
void Database::change_basis(std::vector<int> new_basis)
{
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->nu_aqueous = data->nu_aqueous * betainv;
data->logk_aqueous -= data->nu_aqueous * kswap;
swap_aq_param(new_basis);
swap_labels(new_basis);
// no do the minerals
data->nu_mineral = data->nu_mineral*betainv;
data->logk_mineral -= data->nu_mineral * kswap;
}
void Database::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 Database::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