diff --git a/src/specmicp/adimensional/adimensional_system.cpp b/src/specmicp/adimensional/adimensional_system.cpp
index 3cd2da3..ae889d9 100644
--- a/src/specmicp/adimensional/adimensional_system.cpp
+++ b/src/specmicp/adimensional/adimensional_system.cpp
@@ -1,687 +1,707 @@
/*-------------------------------------------------------
- Module : specmicp
- File : adim_system.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
//#include <iostream>
#include <cmath>
#include "adimensional_system.hpp"
#include "utils/log.hpp"
#include "../equilibrium_data.hpp"
#include "physics/constants.hpp"
#include "physics/laws.hpp"
#include "adimensional_system_solution.hpp"
//#define SPECMICP_DEBUG_EQUATION_FD_JACOBIAN
namespace specmicp {
constexpr scalar_t log10 = std::log(10.0);
-AdimensionalSystem::AdimensionalSystem(
- RawDatabasePtr ptrdata,
- const Vector &totconc,
+AdimensionalSystem::AdimensionalSystem(RawDatabasePtr ptrdata,
+ const AdimensionalSystemBC& conditions,
const AdimensionalSystemOptions& options)
:
OptionsHandler<AdimensionalSystemOptions>(options),
m_data(ptrdata),
- m_tot_conc(totconc),
+ m_fixed_values(conditions.total_concentrations),
m_secondary_conc(ptrdata->nb_aqueous),
m_loggamma(ptrdata->nb_component+ptrdata->nb_aqueous),
m_saturation_gas(0),
m_gas_fugacity(ptrdata->nb_gas)
{
- number_eq();
+ number_eq(conditions);
m_secondary_conc.setZero();
m_gas_fugacity.setZero();
m_loggamma.setZero();
}
AdimensionalSystem::AdimensionalSystem(
RawDatabasePtr ptrdata,
- const Vector& totconc,
+ const AdimensionalSystemBC& conditions,
const AdimensionalSystemSolution& previous_solution,
const AdimensionalSystemOptions& options
):
OptionsHandler<AdimensionalSystemOptions>(options),
m_data(ptrdata),
- m_tot_conc(totconc),
+ m_fixed_values(conditions.total_concentrations),
m_secondary_conc(previous_solution.secondary_molalities),
m_loggamma(previous_solution.log_gamma),
m_ionic_strength(previous_solution.ionic_strength),
m_saturation_gas(0),
m_gas_fugacity(ptrdata->nb_gas)
{
- number_eq();
+ number_eq(conditions);
m_gas_fugacity.setZero();
}
-
-void AdimensionalSystem::number_eq()
+void AdimensionalSystem::number_eq_aqueous_component(
+ const AdimensionalSystemBC& conditions,
+ index_t& neq
+ )
{
- index_t neq = 0;
- m_ideq.reserve(m_data->nb_component+m_data->nb_mineral);
- m_aqueous_component_equation_type.reserve(m_data->nb_component);
- if (get_options().conservation_water)
+ using EqT = AqueousComponentEquationType;
+ // First set the charge keeper
+ if (conditions.charge_keeper != no_species)
{
- m_ideq.push_back(neq);
- m_aqueous_component_equation_type.push_back(AqueousComponentEquationType::MassConservation);
- ++neq;
- } else {
- m_ideq.push_back(no_equation);
- m_aqueous_component_equation_type.push_back(AqueousComponentEquationType::NoEquation);
+ if (conditions.charge_keeper == 0 or conditions.charge_keeper > m_data->nb_component)
+ {
+ throw std::invalid_argument("The charge keeper must be an aqueous component : "
+ + m_data->labels_basis[conditions.charge_keeper]);
+ }
+ m_aqueous_component_equation_type[conditions.charge_keeper] = EqT::ChargeBalance;
}
- for (index_t i: m_data->range_aqueous_component())
+ // Then go over fix fugacity gas
+ for (auto it=conditions.fixed_fugacity_bc.begin(); it!=conditions.fixed_fugacity_bc.end(); ++it)
{
- // Remove components that does not exist
- if (m_tot_conc(i) == 0 and i!= 1)
+ if (m_aqueous_component_equation_type[it->id_component] != EqT::NoEquation)
{
- INFO << "Component " << m_data->labels_basis[i]
- << "has total concentration equal to zero, we desactivate it" << std::endl;
- m_nonactive_component.push_back(i);
- m_ideq.push_back(no_equation);
- m_aqueous_component_equation_type.push_back(AqueousComponentEquationType::NoEquation);
- continue;
+ throw std::invalid_argument("Component '" + m_data->labels_basis[it->id_component]
+ + "' is already restrained");
}
- else
+ m_aqueous_component_equation_type[it->id_component] = EqT::FixedFugacity;
+ m_fixed_values(it->id_component) = it->value;
+ m_fixed_activity_species[it->id_component] = it->id_component;
+ }
+ // Then over the fix activity species
+ for (auto it=conditions.fixed_activity_bc.begin(); it!=conditions.fixed_activity_bc.end(); ++it)
+ {
+ if (m_aqueous_component_equation_type[it->id_component] != EqT::NoEquation)
{
- m_ideq.push_back(neq);
- if (get_options().charge_keeper == i) {
- m_aqueous_component_equation_type.push_back(AqueousComponentEquationType::ChargeBalance);
- }
- else {
- m_aqueous_component_equation_type.push_back(AqueousComponentEquationType::MassConservation);
- }
+ throw std::invalid_argument("Component '" + m_data->labels_basis[it->id_component]
+ + "' is already restrained");
+ }
+ m_aqueous_component_equation_type[it->id_component] = EqT::FixedActivity;
+ m_fixed_values(it->id_component) = it->value;
+ }
+ // Finally number the equations
+ for (index_t component: m_data->range_aqueous_component())
+ {
+ // Mass is conserved for this component
+ if (m_aqueous_component_equation_type[component] == EqT::NoEquation
+ and m_fixed_values(component) != 0.0)
+ {
+ m_aqueous_component_equation_type[component] = EqT::MassConservation;
+ m_ideq[component] = neq;
+ ++neq;
+ }
+ // Or this is another type of equation
+ else if (m_aqueous_component_equation_type[component] != EqT::NoEquation)
+ {
+ m_ideq[component] = neq;
++neq;
}
+ // else add component to the nonactive component list
+ else
+ {
+ m_nonactive_component.push_back(component);
+ }
+ }
+}
+
+void AdimensionalSystem::number_eq(
+ const AdimensionalSystemBC& conditions
+ )
+{
+ index_t neq = 0;
+ m_ideq =std::vector<index_t>(m_data->nb_component+m_data->nb_mineral, no_equation);
+ m_aqueous_component_equation_type = std::vector<AqueousComponentEquationType>(
+ m_data->nb_component, AqueousComponentEquationType::NoEquation);
+ m_fixed_activity_species = std::vector<index_t>(m_data->nb_component, no_species);
+ // Water
+ // =====
+ if (conditions.conservation_water)
+ {
+ m_ideq[0] = neq;
+ m_aqueous_component_equation_type[0] = AqueousComponentEquationType::MassConservation;
+ ++neq;
}
+ // Aqueous components
+ // ==================
+ number_eq_aqueous_component(conditions, neq);
m_nb_free_vars = neq;
+ // Minerals
+ // ========
for (index_t m: m_data->range_mineral())
{
bool can_precipitate = true;
// Remove minerals that cannot precipitate
for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it)
{
if (m_data->nu_mineral(m, *it) != 0)
{
can_precipitate = false;
break; // this is not a mineral that can precipitate
}
}
if (can_precipitate)
{
- m_ideq.push_back(neq);
+ m_ideq[m_data->nb_component+m] = 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
+ // Secondary species
+ // =================
+ // Secondary aqueous species
+ // -------------------------
m_active_aqueous.reserve(m_data->nb_aqueous);
for (index_t j: m_data->range_aqueous())
{
bool can_exist = true;
for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it)
{
if (m_data->nu_aqueous(j,*it) != 0)
{
can_exist = false;
}
}
m_active_aqueous.push_back(can_exist);
}
m_active_gas.reserve(m_data->nb_gas);
+ // Gas
+ // ---
for (index_t k: m_data->range_gas())
{
bool can_exist = true;
for (auto it=m_nonactive_component.begin(); it!=m_nonactive_component.end(); ++it)
{
if (m_data->nu_gas(k,*it) != 0)
{
can_exist = false;
}
}
m_active_gas.push_back(can_exist);
}
}
scalar_t AdimensionalSystem::saturation_water(const Vector& x) const
{
if (ideq_w() != no_equation)
return x(ideq_w());
else
return 1.0 - sum_saturation_minerals(x);
}
scalar_t AdimensionalSystem::saturation_mineral(const Vector& x, index_t mineral) const
{
specmicp_assert(mineral >= 0 and mineral < m_data->nb_mineral);
if (ideq_min(mineral) == no_equation) return 0.0;
else return x(ideq_min(mineral));
}
scalar_t AdimensionalSystem::sum_saturation_minerals(const Vector& x) const
{
scalar_t sum_saturations = 0.0;
for (index_t mineral: m_data->range_mineral())
{
sum_saturations += saturation_mineral(x, mineral);
}
return sum_saturations;
}
scalar_t AdimensionalSystem::density_water() const
{
return laws::density_water(units::celsius(25.0), length_unit(), mass_unit());
}
scalar_t AdimensionalSystem::molar_volume_mineral(index_t mineral) const
{
return m_data->molar_volume_mineral(mineral, length_unit());
}
scalar_t AdimensionalSystem::residual_water(const Vector& x) const
{
const scalar_t conc_w = density_water()*saturation_water(x);
- scalar_t res = m_tot_conc(0) - conc_w/m_data->molar_mass_basis_si(0);
+ scalar_t res = total_concentration_bc(0) - conc_w/m_data->molar_mass_basis_si(0);
for (index_t j: m_data->range_aqueous())
{
if (not is_aqueous_active(j)) continue;
res -= conc_w*m_data->nu_aqueous(j, 0)*m_secondary_conc(j);
}
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation or m_data->nu_mineral(m, 0) == 0.0) continue;
res -= m_data->nu_mineral(m, 0)*saturation_mineral(x, m)/molar_volume_mineral(m);
}
for (index_t k: m_data->range_gas())
{
if (m_data->nu_gas(k, 0) == 0.0) continue;
res -= m_data->nu_gas(k, 0)*saturation_gas_phase()*gas_fugacity(k)/(
constants::gas_constant*gas_total_pressure()*temperature());
}
return res;
}
double AdimensionalSystem::residual_component(const Vector &x, index_t component) const
{
const scalar_t conc_w = density_water()*saturation_water(x);
- scalar_t res = m_tot_conc(component) - conc_w*component_molality(x, component);
+ scalar_t res = total_concentration_bc(component) - conc_w*component_molality(x, component);
for (index_t j: m_data->range_aqueous())
{
if (not is_aqueous_active(j)) continue;
res -= conc_w*m_data->nu_aqueous(j, component)*secondary_molality(j);
}
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation) continue;
res -= m_data->nu_mineral(m, component)*saturation_mineral(x, m)/molar_volume_mineral(m);
}
for (index_t k: m_data->range_gas())
{
if (m_data->nu_gas(k, component) == 0.0) continue;
res -= m_data->nu_gas(k, component)*saturation_gas_phase()*gas_fugacity(k)/(
constants::gas_constant*gas_total_pressure()*temperature());
}
return res;
}
double AdimensionalSystem::residual_mineral(const Vector& x, index_t m) const
{
scalar_t res = m_data->logk_mineral(m);
for (index_t i: m_data->range_aqueous_component())
{
if (m_data->nu_mineral(m, i) != 0)
res -= m_data->nu_mineral(m, i)*(log_component_molality(x, i) + log_gamma_component(i));
}
return res;
}
double AdimensionalSystem::residual_charge_conservation(const Vector& x) const
{
scalar_t res = 0.0;
for (index_t i: m_data->range_aqueous_component())
{
if (m_data->charge_component(i) != 0 and ideq_paq(i) != no_equation)
res += m_data->charge_component(i)*component_molality(x, i);
}
for (index_t j: m_data->range_aqueous())
{
if (m_data->charge_aqueous(j) == 0 and not is_aqueous_active(j)) continue;
res += m_data->charge_aqueous(j)*secondary_molality(j);
}
return res;
}
void AdimensionalSystem::get_residuals(const Vector& x, Vector& residual)
{
set_secondary_concentration(x);
if (ideq_w() != no_equation) residual(ideq_w()) = residual_water(x);
for (index_t i: m_data->range_aqueous_component())
{
switch (aqueous_component_equation_type(i))
{
case AqueousComponentEquationType::NoEquation:
continue;
case AqueousComponentEquationType::ChargeBalance:
residual(ideq_paq(i)) = residual_charge_conservation(x);
break;
default:
+ specmicp_assert(m_aqueous_component_equation_type[i]
+ == AqueousComponentEquationType::MassConservation);
residual(ideq_paq(i)) = residual_component(x, i);
}
}
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) != no_equation) residual(ideq_min(m)) = residual_mineral(x, m);
}
}
void AdimensionalSystem::get_jacobian(Vector& x, Matrix& jacobian)
// //non-optimized Finite difference, for test only
#ifdef SPECMICP_DEBUG_EQUATION_FD_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;
}
#else // analytical jacobian
{
jacobian.resize(total_variables(), total_variables());
jacobian.setZero();
// water
jacobian_water(x, jacobian);
// aqueous component
jacobian_aqueous_components(x, jacobian);
// mineral equilibrium
jacobian_minerals(x, jacobian);
}
#endif
void AdimensionalSystem::jacobian_water(Vector& x, Matrix& jacobian)
{
if (ideq_w() != no_equation)
{
scalar_t tmp = -density_water()/m_data->molar_mass_basis_si(0);
for (index_t j: m_data->range_aqueous())
{
if ( m_data->nu_aqueous(j, 0) == 0 and not is_aqueous_active(j)) continue;
tmp -= m_data->nu_aqueous(j, 0)*secondary_molality(j);
}
jacobian(0, 0) = tmp;
const scalar_t conc_w = density_water()*saturation_water(x);
for (index_t k: m_data->range_aqueous_component())
{
if (ideq_paq(k) == no_equation) continue;
scalar_t tmp = 0;
for (index_t j: m_data->range_aqueous())
{
tmp -= m_data->nu_aqueous(j,0)*m_data->nu_aqueous(j,k)*secondary_molality(j);
}
jacobian(0, ideq_paq(k)) = conc_w*log10 * tmp;
}
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation) continue;
jacobian(0, ideq_min(m)) = -m_data->nu_mineral(m, 0)/molar_volume_mineral(m);
}
}
}
void AdimensionalSystem::jacobian_aqueous_components(Vector& x, Matrix& jacobian)
{
for (index_t i: m_data->range_aqueous_component())
{
const index_t idp = ideq_paq(i);
if (idp == no_equation) continue;
// Charge balance equation
// ----------------------
switch (aqueous_component_equation_type(i))
{
case AqueousComponentEquationType::NoEquation:
continue;
case AqueousComponentEquationType::ChargeBalance:
// Aqueous components
for (index_t k: m_data->range_aqueous_component())
{
const index_t idc = ideq_paq(k);
if (idc == no_equation) continue;
scalar_t tmp_drdb = 0.0;
if (m_data->charge_component(k) != 0.0)
tmp_drdb = m_data->charge_component(k);
// Secondary species
for (index_t j: m_data->range_aqueous())
{
if ( not is_aqueous_active(j)
or m_data->nu_aqueous(j, k) == 0.0
or m_data->charge_aqueous(j) == 0.0
)
continue;
scalar_t tmp_value = m_data->nu_aqueous(j, k)*m_data->charge_aqueous(j);
tmp_value *= secondary_molality(j)/component_molality(x, k);
tmp_drdb += tmp_value;
}
jacobian(idp, idc) += component_molality(x,k)*log10*tmp_drdb;
}
break;
// Mass balance equation
// ---------------------
default:
const scalar_t conc_w = density_water()*saturation_water(x);
// Aqueous components
for (index_t k: m_data->range_aqueous_component())
{
if (ideq_paq(k) == no_equation) continue;
scalar_t tmp_iip = 0;
if (k == i) tmp_iip -= component_molality(x, i)*log10;
for (index_t j: m_data->range_aqueous())
{
tmp_iip -= log10*m_data->nu_aqueous(j, i)*m_data->nu_aqueous(j, k)*secondary_molality(j);
}
jacobian(idp, ideq_paq(k)) = conc_w*tmp_iip;
}
// Minerals
for (index_t m: m_data->range_mineral())
{
if (ideq_min(m) == no_equation) continue;
jacobian(idp, ideq_min(m)) = - m_data->nu_mineral(m, i)/molar_volume_mineral(m);
}
// Water
if (ideq_w() != no_equation)
{
scalar_t tmp_iw = -component_molality(x, i);
for (index_t j: m_data->range_aqueous())
{
if ( m_data->nu_aqueous(j, i) == 0 ) continue;
tmp_iw -= m_data->nu_aqueous(j, i)*secondary_molality(j);
}
jacobian(idp, ideq_w()) = density_water()*tmp_iw;
}
}
}
}
void AdimensionalSystem::jacobian_minerals(Vector& x, Matrix& jacobian)
{
for (index_t m: m_data->range_mineral())
{
const index_t idm = ideq_min(m);
if (idm == no_equation) continue;
for (index_t i: m_data->range_aqueous_component())
{
if (ideq_paq(i) == no_equation) continue;
jacobian(idm, ideq_paq(i)) = -m_data->nu_mineral(m, i);
}
}
}
void AdimensionalSystem::set_secondary_variables(const Vector& x)
{
set_saturation_gas_phase(x);
set_pressure_fugacity(x);
set_secondary_concentration(x);
compute_log_gamma(x);
}
void AdimensionalSystem::set_saturation_gas_phase(const Vector& x)
{
m_saturation_gas = saturation_water(x) - sum_saturation_minerals(x);
}
void AdimensionalSystem::set_pressure_fugacity(const Vector& x)
{
for (index_t k: m_data->range_gas())
{
if (not is_active_gas(k)) continue;
scalar_t logp = m_data->logk_gas(k);
for (index_t i: m_data->range_aqueous_component())
{
if (m_data->nu_gas(k, i) == 0.0) continue;
logp += m_data->nu_gas(k, i) * (log_component_molality(x, i) + log_gamma_component(i));
}
m_gas_fugacity(k) = pow10(logp);
}
}
void AdimensionalSystem::set_secondary_concentration(const Vector& x)
{
for (index_t j: m_data->range_aqueous())
{
if (not is_aqueous_active(j))
{
m_secondary_conc(j) = 0;
continue;
}
scalar_t logconc = -m_data->logk_aqueous(j) - m_loggamma(j+m_data->nb_component);
for (index_t k: m_data->range_aqueous_component())
{
if (m_data->nu_aqueous(j, k) == 0) continue;
logconc += m_data->nu_aqueous(j, k)*(log_component_molality(x, k) + m_loggamma(k));
}
m_secondary_conc(j) = pow10(logconc);
}
}
void AdimensionalSystem::set_ionic_strength(const Vector& x)
{
scalar_t ionic = 0;
for (index_t i: m_data->range_aqueous_component())
{
if (ideq_paq(i) == no_equation or m_data->charge_component(i) == 0) continue;
ionic += component_molality(x, i)*std::pow(m_data->charge_component(i),2);
}
for (index_t j: m_data->range_aqueous())
{
if (not is_aqueous_active(j) or m_data->charge_aqueous(j) == 0) continue;
ionic += secondary_molality(j)*std::pow(m_data->charge_aqueous(j),2);
}
ionic_strength() = ionic/2;
}
void AdimensionalSystem::compute_log_gamma(const Vector& x)
{
set_ionic_strength(x);
const scalar_t sqrti = std::sqrt(ionic_strength());
for (index_t i: m_data->range_aqueous_component())
{
if (ideq_paq(i) == no_equation)
{
log_gamma_component(i) = 0.0;
continue;
}
log_gamma_component(i) = laws::extended_debye_huckel(
ionic_strength(), sqrti,
m_data->charge_component(i),
m_data->a_debye_component(i),
m_data->b_debye_component(i)
);
}
for (index_t j: m_data->range_aqueous())
{
if (not is_aqueous_active(j))
{
log_gamma_secondary(j) = 0.0;
continue;
}
log_gamma_secondary(j) = laws::extended_debye_huckel(
ionic_strength(), sqrti,
m_data->charge_aqueous(j),
m_data->a_debye_aqueous(j),
m_data->b_debye_aqueous(j)
);
}
}
bool AdimensionalSystem::hook_start_iteration(const Vector& x, scalar_t norm_residual)
{
if (not get_options().non_ideality) return true;
not_in_linesearch = true;
scalar_t previous_norm = m_loggamma.norm();
bool may_have_converged = false;
if (norm_residual < nb_free_variables())
{
// Use fixed point iterations for non-ideality
for (int i=0; i<get_options().non_ideality_max_iter; ++i)
{
set_secondary_variables(x);
compute_log_gamma(x);
// convergence check
if (std::abs(previous_norm - m_loggamma.norm())/previous_norm <
get_options().non_ideality_tolerance) {
may_have_converged = true;
break;
}
previous_norm = m_loggamma.norm();
}
}
return may_have_converged;
}
double AdimensionalSystem::max_lambda(const Vector& x, const Vector& update)
{
if (ideq_w() != no_equation)
{
return 1.0/std::max(1.0, -update(0)/(get_options().under_relaxation_factor*x(0)));
}
else
{
return 1.0;
}
}
void AdimensionalSystem::reasonable_starting_guess(Vector &x)
{
x.resize(m_data->nb_component+ m_data->nb_mineral);
x(0) = 1.0;
for (index_t i: m_data->range_aqueous_component())
{
x(i) = -4.0;
}
x.block(m_data->nb_component, 0, m_data->nb_mineral, 1).setZero();
m_loggamma.setZero();
m_secondary_conc.setZero();
}
void AdimensionalSystem::reasonable_restarting_guess(Vector& x)
{
x(0) = 0.5;
for (index_t i: m_data->range_aqueous_component())
{
if (x(i) > 0 or x(i) < -9)
x(i) = -3;
}
x.segment(m_data->nb_component, m_data->nb_mineral).setZero();
m_loggamma.setZero();
m_secondary_conc.setZero();
}
AdimensionalSystemSolution AdimensionalSystem::get_solution(const Vector& xtot, const Vector& x)
{
double previous_norm = m_loggamma.norm();
set_saturation_gas_phase(x);
set_pressure_fugacity(x);
set_secondary_concentration(x);
if (get_options().non_ideality)
{
compute_log_gamma(x);
if (std::abs(previous_norm - m_loggamma.norm()) > 1e-6)
{
WARNING << "Activity coefficient have not converged !" << std::endl
<< "output can not be trusted\n Difference : "
+std::to_string(std::abs(previous_norm - m_loggamma.norm()));
}
}
return AdimensionalSystemSolution(xtot, m_secondary_conc, m_loggamma, m_ionic_strength);
}
-void AdimensionalSystem::reduce_mineral_problem(Vector& true_var)
-{
- for (index_t mineral: m_data->range_mineral())
- {
- if (ideq_min(mineral) == no_equation) continue;
- if (m_data->stability_mineral(mineral) == database::MineralStabilityClass::cannot_dissolve)
- {
- if (true_var(ideq_min(mineral)) == 0.0) continue;
- for (index_t component: m_data->range_component())
- {
- if (m_data->nu_mineral(mineral, component) == 0.0) continue;
- m_tot_conc(component) -= m_data->nu_mineral(mineral, component)*true_var(ideq_min(mineral));
- }
- true_var(ideq_min(mineral)) = 0.0;
- }
- }
-}
-
-void AdimensionalSystem::reset_mineral_system(Vector& true_var, Vector& output_var)
-{
- for (index_t mineral: m_data->range_mineral())
- {
- if (ideq_min(mineral) == no_equation) continue;
- if (m_data->stability_mineral(mineral) == database::MineralStabilityClass::cannot_dissolve)
- {
- output_var(mineral) += true_var(ideq_min(mineral));
- }
- }
-}
} // end namespace specmicp
diff --git a/src/specmicp/adimensional/adimensional_system.hpp b/src/specmicp/adimensional/adimensional_system.hpp
index 4e50eec..a7f16e2 100644
--- a/src/specmicp/adimensional/adimensional_system.hpp
+++ b/src/specmicp/adimensional/adimensional_system.hpp
@@ -1,356 +1,366 @@
/*-------------------------------------------------------
- Module : specmicp
- File : adimensional_system.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#ifndef SPECMICP_ADIMENSIONALSYSTEM_HPP
#define SPECMICP_ADIMENSIONALSYSTEM_HPP
//! \file reduced_system.hpp The MiCP program to solve speciation
#include "common.hpp"
#include "database.hpp"
#include "micpsolver/micpprog.hpp"
#include "../simulation_box.hpp"
#include "physics/units.hpp"
#include "utils/options_handler.hpp"
#include "adimensional_system_structs.hpp"
//! \namespace specmicp main namespace for the SpecMiCP solver
namespace specmicp {
class AdimensionalSystemSolution;
//! \brief The equilibrium speciation problem
//!
//! Represent the equilibrium speciation problem as a Mixed Complementarity program
//! It is called reduced becaused it does not solve explicitely for the secondary species.
//!
class AdimensionalSystem :
public micpsolver::MiCPProg<AdimensionalSystem>,
public OptionsHandler<AdimensionalSystemOptions>,
public units::UnitBaseClass
{
public:
//! \brief Create a reduced system
//!
//! \param ptrthermo Shared pointer to the thermodynamic data
//! \param totconc Vector containing the total concentrations (in mol/volume) for each components
AdimensionalSystem(
RawDatabasePtr ptrdata,
- const Vector& total_concentration,
+ const AdimensionalSystemBC& conditions,
const AdimensionalSystemOptions& options=AdimensionalSystemOptions()
);
AdimensionalSystem(
RawDatabasePtr ptrdata,
- const Vector& totconc,
+ const AdimensionalSystemBC& conditions,
const AdimensionalSystemSolution& previous_solution,
- const AdimensionalSystemOptions& options=AdimensionalSystemOptions()
+ const AdimensionalSystemOptions& options
);
// Variables
// ==========
//! \brief Return the total number of variables
index_t total_variables() const {return m_nb_tot_vars;}
//! \brief Return the number of 'free' variables (not subject to complementarity conditions)
index_t nb_free_variables() const {return m_nb_free_vars;}
//! \brief Return the number of variables subjected to the complementarity conditions
index_t nb_complementarity_variables() const {return m_nb_compl_vars;}
// The linear system
// ==================
//! \brief Return the residual for the water conservation equation
scalar_t residual_water(const Vector& x) const;
//! \brief Return the residual for the conservation equation of component (i)
//!
//! For (i==0), water use the method : 'residual_water'
double residual_component(const Vector& x, index_t i) const;
//! \brief Equilibrium condition for the minerals \f$ 1 - \mathrm{SI}_l \f$
double residual_mineral(const Vector& x, index_t m) const;
//! \brief Charge conservation
double residual_charge_conservation(const Vector& x) const;
//! \brief Return the residuals
void get_residuals(const Vector& x, Vector& residual);
//! \brief Return the jacobian
void get_jacobian(Vector& x,
Matrix& jacobian);
// Getters
// #######
// Equation id access
// ==================
//! \brief Return the equation id
index_t ideq(index_t i) const {return m_ideq[i];}
//! \brief Return the equation number for conservation of water
index_t ideq_w() const {return m_ideq[0];}
//! \brief Return the equation number of component 'i'
index_t ideq_paq(index_t i) const {
specmicp_assert(i < m_data->nb_component);
return m_ideq[i];
}
//! \brief Return the equation number of mineral 'm'
index_t ideq_min(index_t m) const {
specmicp_assert(m < m_data->nb_mineral);
return m_ideq[m+m_data->nb_component];
}
//! \brief Return the type of equation for aqueous species 'component'
AqueousComponentEquationType aqueous_component_equation_type(index_t component) const {
specmicp_assert(component >= 0 and component < m_data->nb_component);
return m_aqueous_component_equation_type[component];
}
+ //! Return true if an equation exist for this component
+ bool is_active_component(index_t component) const {
+ specmicp_assert(component < m_data->nb_component);
+ return m_ideq[component] != no_equation;
+ }
+
// Variables
// =========
// Primary
// -------
//! \brief Return the saturation of water
scalar_t saturation_water(const Vector& x) const;
//! \brief Return the density of water
scalar_t density_water() const;
//! \brief Return the saturation (mol/volume) of a mineral
scalar_t saturation_mineral(const Vector& x, index_t mineral) const;
//! \brief Return the volume of minerals
scalar_t molar_volume_mineral(index_t mineral) const;
//! \brief Return the sum of saturation of the minerals
scalar_t sum_saturation_minerals(const Vector& x) const;
//! \brief Return the log_10 of component concentration
scalar_t log_component_molality(const Vector& x, index_t component) const
{
specmicp_assert(ideq_paq(component) != no_equation and component < m_data->nb_component);
return x(ideq_paq(component));
}
//! \brief Return the concentration of 'component'
scalar_t component_molality(const Vector& x, index_t component) const
{
return pow10(log_component_molality(x, component));
}
// Secondary
// ---------
//! \brief Return the concentration of secondary specis 'aqueous'
scalar_t secondary_molality(index_t aqueous) const
{
if (not m_active_aqueous[aqueous]) return 0.0;
return m_secondary_conc(aqueous);
}
//! \brief Return true if 'aqueous' is active
//!
//! i.e. Return true if all its component are present in the system
bool is_aqueous_active(index_t aqueous) const
{
return m_active_aqueous[aqueous];
}
//! \brief Return log_10(γ_i) where i is a component
scalar_t log_gamma_component(index_t component) const {
return m_loggamma(component);
}
//! \brief Return log_10(γ_j) where j is a secondary aqueous species
scalar_t log_gamma_secondary(index_t aqueous) const {
return m_loggamma(m_data->nb_component + aqueous);
}
//! \brief Return the ionic strength
scalar_t ionic_strength() const {
return m_ionic_strength;
}
//! \brief Return the total concentration of 'component'
- scalar_t total_concentration(index_t component) const {return m_tot_conc(component);}
+ scalar_t total_concentration_bc(index_t component) const {
+ specmicp_assert(m_aqueous_component_equation_type[component]
+ == AqueousComponentEquationType::MassConservation);
+ return m_fixed_values(component);
+ }
+ //! \brief Return the total concentration of 'component'
+ scalar_t fixed_activity_bc(index_t component) const {
+ specmicp_assert(m_aqueous_component_equation_type[component]
+ == AqueousComponentEquationType::FixedActivity);
+ return m_fixed_values(component);
+ }
+ //! \brief Return the total concentration of 'component'
+ scalar_t fixed_fugacity_bc(index_t component) const {
+ specmicp_assert(m_aqueous_component_equation_type[component]
+ == AqueousComponentEquationType::FixedFugacity);
+ return m_fixed_values(component);
+ }
//! \brief Return the fugacity of 'gas'
scalar_t gas_fugacity(index_t gas) const {
return m_gas_fugacity(gas);
}
//! \brief Return true if gas is active
bool is_active_gas(index_t gas) const {
return m_active_gas[gas];
}
//! \brief Return the saturation of the gas phase
scalar_t saturation_gas_phase() const {
return m_saturation_gas;
}
//! \brief Return the gas pressure
scalar_t gas_total_pressure() const {
return 1.01325e5;
}
//! \brief Return the temperature
scalar_t temperature() const {
return 273.16+25;
}
// Equilibrium State
// =================
//! \brief Return the equilibrium state of the system, the 'Solution' of the speciation problem
AdimensionalSystemSolution get_solution(const Vector& xtot, const Vector& x);
// Algorithm
// #########
//! \brief Compute the ionic strength
void set_ionic_strength(const Vector& x);
//! \brief Compute the activity coefficients
void compute_log_gamma(const Vector& x);
// For the normal algorithm
void set_secondary_concentration(const Vector& x);
//! \brief Compute the secondary variables
void set_secondary_variables(const Vector& x);
//! \brief Set the gas phase saturation
void set_saturation_gas_phase(const Vector& x);
//! \brief Set the fugacity pressure
void set_pressure_fugacity(const Vector& x);
//! \brief Compute the activity model, called by the solver at the beginning of an iteration
bool hook_start_iteration(const Vector &x, scalar_t norm_residual);
//! \brief Return a scale factor to avoid negative mass during Newton's iteration
double max_lambda(const Vector &x, const Vector &update);
//! \brief Return the component that does not exist in the solution (inactive dof)
const std::vector<index_t>& get_non_active_component() {return m_nonactive_component;}
//! \brief A reasonable (well... maybe...) starting guess
void reasonable_starting_guess(Vector& x);
//! \brief A reasonable (maybe...) restarting guess
void reasonable_restarting_guess(Vector& x);
- //! \brief Reduce the mineral system if one cannot dissolve
- void reduce_mineral_problem(Vector& true_var);
- //! \brief Reset the mineral system to the true problem
- void reset_mineral_system(Vector& true_var, Vector& output_var);
-
- //! \brief set the charge keeper to 'component'
+ //! \brief Return true if 'component' is the charge keeper
//!
- //! The equation for the charge keeper is replaced by the charge balance
- //!
- //! Set component to 'no_species' to solve charge conservation implicitely.
- //! In this case, the total concentration must conserve the charge.
- void set_charge_keeper(index_t component)
- {
- specmicp_assert(component >= 1 and component < m_data->nb_component);
- get_options().charge_keeper = component;
- }
-
- //! \brief Return true if 'component' is thecharge keeper
+ //! The charge keeper is the component added or removed to satisfy the charge balance
bool is_charge_keeper(index_t component)
{
- return ( get_options().charge_keeper != no_species
- and component == get_options().charge_keeper
- );
+ return (m_aqueous_component_equation_type[component] == AqueousComponentEquationType::ChargeBalance);
}
// Private Members and attributes
// ##############################
private:
// Jacobian
// ########
void jacobian_water(Vector& x, Matrix& jacobian);
void jacobian_aqueous_components(Vector& x, Matrix& jacobian);
void jacobian_minerals(Vector& x, Matrix& jacobian);
// Equation numbering
// ##################
- void number_eq();
+ void number_eq(const AdimensionalSystemBC& condition);
+ void number_eq_aqueous_component(
+ const AdimensionalSystemBC& conditions,
+ index_t& neq
+ );
// Setter
// ######
// main variables
// --------------
// they are set by the solver
// Secondary variables
// -------------------
//! \brief Return a reference to log_10(γ_i) where i is a component
scalar_t& log_gamma_component(index_t component) {
return m_loggamma(component);
}
//! \brief Return a reference to log_10(γ_j) where j is a secondary aqueous species
scalar_t& log_gamma_secondary(index_t aqueous) {
return m_loggamma(m_data->nb_component + aqueous);
}
//! \brief Return a reference to the ionic strength
scalar_t& ionic_strength() {
return m_ionic_strength;
}
// Attributes
// ##########
RawDatabasePtr m_data;
- Vector m_tot_conc;
+ Vector m_fixed_values;
index_t m_nb_tot_vars;
index_t m_nb_free_vars;
index_t m_nb_compl_vars;
std::vector<index_t> m_ideq;
std::vector<AqueousComponentEquationType> m_aqueous_component_equation_type;
+ std::vector<index_t> m_fixed_activity_species;
Vector m_secondary_conc;
Vector m_loggamma;
scalar_t m_ionic_strength;
scalar_t m_saturation_gas;
Vector m_gas_fugacity;
bool not_in_linesearch;
std::vector<index_t> m_nonactive_component;
std::vector<bool> m_active_aqueous;
std::vector<bool> m_active_gas;
};
} // end namespace specmicp
#endif // SPECMICP_ADIMENSIONALSYSTEM_HPP
diff --git a/src/specmicp/adimensional/adimensional_system_solver.cpp b/src/specmicp/adimensional/adimensional_system_solver.cpp
index 553b2d4..9c5ccae 100644
--- a/src/specmicp/adimensional/adimensional_system_solver.cpp
+++ b/src/specmicp/adimensional/adimensional_system_solver.cpp
@@ -1,213 +1,210 @@
/*-------------------------------------------------------
- Module : specmicp
- File : reduced_system_solver
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#include "adimensional_system_solver.hpp"
#include "micpsolver/micpsolver.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
namespace specmicp {
-AdimensionalSystemSolver::AdimensionalSystemSolver(
- RawDatabasePtr data,
- const Vector& totconc,
+AdimensionalSystemSolver::AdimensionalSystemSolver(RawDatabasePtr data,
+ const AdimensionalSystemBC& conditons,
AdimensionalSystemSolverOptions options
):
OptionsHandler<AdimensionalSystemSolverOptions>(options),
m_data(data),
m_system(std::make_shared<AdimensionalSystem>(
data,
- totconc,
+ conditons,
options.system_options
)),
m_var(Vector::Zero(data->nb_component+data->nb_mineral))
{}
AdimensionalSystemSolver::AdimensionalSystemSolver(RawDatabasePtr data,
- const Vector& totconc,
+ const AdimensionalSystemBC& conditions,
const AdimensionalSystemSolution& previous_solution,
AdimensionalSystemSolverOptions options
):
OptionsHandler<AdimensionalSystemSolverOptions>(options),
m_data(data),
m_system(std::make_shared<AdimensionalSystem>(
data,
- totconc,
+ conditions,
previous_solution,
options.system_options
)),
m_var(Vector::Zero(data->nb_component+data->nb_mineral))
{}
AdimensionalSystemSolution AdimensionalSystemSolver::get_raw_solution(const Eigen::VectorXd& x)
{
set_true_variable_vector(x);
return m_system->get_solution(x, m_var);
}
micpsolver::MiCPPerformance AdimensionalSystemSolver::solve(Vector& x, bool init)
{
m_system->set_units(get_options().units_set);
if (init)
{
m_system->reasonable_starting_guess(x);
}
set_true_variable_vector(x);
micpsolver::MiCPPerformance perf = solve();
if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized and get_options().allow_restart)
{
WARNING << "Failed to solve the system ! - We shake it up and start again";
const scalar_t save_penalization_factor = get_options().solver_options.penalization_factor;
const scalar_t save_factor_descent = get_options().solver_options.factor_descent_condition;
get_options().solver_options.factor_descent_condition *= 1e-2;
if (save_penalization_factor == 1)
get_options().solver_options.penalization_factor = 0.8;
set_return_vector(x);
m_system->reasonable_restarting_guess(x);
set_true_variable_vector(x);
micpsolver::MiCPPerformance perf2 = solve();
get_options().solver_options.penalization_factor = save_penalization_factor;
get_options().solver_options.factor_descent_condition = save_factor_descent;
perf += perf2;
}
if (perf.return_code > micpsolver::MiCPSolverReturnCode::NotConvergedYet) set_return_vector(x);
return perf;
}
void AdimensionalSystemSolver::set_true_variable_vector(const Vector& x)
{
const std::vector<index_t>& non_active_component = m_system->get_non_active_component();
if (non_active_component.size() == 0)
{
// we still copy the data, if we failed to solve the problem, we can restart
- if (get_options().system_options.conservation_water)
+ if (m_system->is_active_component(0))
m_var = x; // direct copy
else
m_var = x.block(1, 0, x.rows()-1, 1);
for (int i=0; i<m_var.rows(); ++i)
{
if (m_var(i) == -HUGE_VAL) // check for previously undefined value
{
m_var(i) = -4;
}
}
}
else // remove the dof that are not part of the problem
{
uindex_t new_i = 0;
- if (get_options().system_options.conservation_water)
+ if (m_system->is_active_component(0))
{
m_var(0) = x(0);
++new_i;
}
for (index_t i: m_data->range_aqueous_component())
{
auto it = std::find(non_active_component.begin(), non_active_component.end(),i);
if (it != non_active_component.end()) continue;
m_var(new_i) = x(i);
++new_i;
}
for (index_t m: m_data->range_mineral())
{
for (auto it=non_active_component.begin(); it!=non_active_component.end(); ++it)
{
if (m_data->nu_mineral(m, *it) != 0) continue;
m_var(new_i) = x(m_data->nb_component+m);
++new_i;
}
}
m_var.conservativeResize(new_i);
}
- m_system->reduce_mineral_problem(m_var);
}
micpsolver::MiCPPerformance AdimensionalSystemSolver::solve()
{
micpsolver::MiCPSolver<AdimensionalSystem> solver(m_system);
solver.set_options(get_options().solver_options);
solver.solve(m_var);
return solver.get_performance();
}
void AdimensionalSystemSolver::set_return_vector(Vector& x)
{
const std::vector<index_t>& non_active_component = m_system->get_non_active_component();
- if (non_active_component.size() == 0)
+ if (non_active_component.size() == 0) // shortcut
{
- if (get_options().system_options.conservation_water)
+ if (m_system->is_active_component(0))
x = m_var; //direct copy
else
x.block(1, 0, x.rows()-1, 1) = m_var;
// at that point we should have the correct solution
}
else
{
uindex_t new_i = 0;
- if (get_options().system_options.conservation_water)
+ if (m_system->is_active_component(0))
{
x(0) = m_var(new_i);
++new_i;
}
for (index_t i: m_data->range_aqueous_component())
{
auto it = std::find(non_active_component.begin(), non_active_component.end(),i);
if (it != non_active_component.end())
{
x(i) = -HUGE_VAL;
continue;
}
x(i) = m_var(new_i) ;
++new_i;
}
for (index_t m: m_data->range_mineral())
{
for (auto it=non_active_component.begin(); it!=non_active_component.end(); ++it)
{
if (m_data->nu_mineral(m, *it) != 0)
{
x(m_data->nb_component+m) = 0;
continue;
}
x(m_data->nb_component+m) =m_var(new_i);
++new_i;
}
}
}
- m_system->reset_mineral_system(m_var, x);
}
} // end namespace specmicp
diff --git a/src/specmicp/adimensional/adimensional_system_solver.hpp b/src/specmicp/adimensional/adimensional_system_solver.hpp
index 9612aae..87fc255 100644
--- a/src/specmicp/adimensional/adimensional_system_solver.hpp
+++ b/src/specmicp/adimensional/adimensional_system_solver.hpp
@@ -1,111 +1,111 @@
/*-------------------------------------------------------
- Module : specmicp
- File : reduced_system_solver.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP
#define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP
#include "database.hpp"
#include "adimensional_system.hpp"
#include "adimensional_system_solver_structs.hpp"
#include "utils/options_handler.hpp"
//! \file reduced_system_solver.hpp Solve a reduced system
namespace specmicp {
// forward declaration
class AdimensionalSystemSolution;
//! \brief Solver a reduced system
//!
//! Take care of possibly non-existing component in the system
class AdimensionalSystemSolver: public OptionsHandler<AdimensionalSystemSolverOptions>
{
public:
using SystemPtr = std::shared_ptr<AdimensionalSystem>;
//! Default constructor - you need to use another method to initialize it correctly
AdimensionalSystemSolver() {}
AdimensionalSystemSolver(RawDatabasePtr data,
- const Vector& totconc,
+ const AdimensionalSystemBC& conditons,
AdimensionalSystemSolverOptions options=AdimensionalSystemSolverOptions()
);
AdimensionalSystemSolver(RawDatabasePtr data,
- const Vector& totconc,
+ const AdimensionalSystemBC& conditions,
const AdimensionalSystemSolution& previous_solution,
AdimensionalSystemSolverOptions options=AdimensionalSystemSolverOptions()
);
//! \brief solve the problem using initial guess x
//!
//! \param[in,out] x in -> initial guess, out -> solution
//! \param init if true, the algorithm guess a starting point
micpsolver::MiCPPerformance solve(Vector& x, bool init=false);
//! \brief Return the system used for the computation
SystemPtr get_system() {return m_system;}
//! \brief Return the solution in a manageable form
AdimensionalSystemSolution get_raw_solution(const Eigen::VectorXd& x);
private:
//! \brief set up the true variable vector
void set_true_variable_vector(const Vector& x);
//! \brief set up the true solution vector
//!
//! add zero components
void set_return_vector(Vector& x);
//! \brief solve the problem
micpsolver::MiCPPerformance solve();
RawDatabasePtr m_data; //! The raw database
SystemPtr m_system; //! The system to solve
Vector m_var; //! Copy of the solution vector (necessary in case of failing)
};
//! \brief Solve a reduced system, function provided for convenience
inline micpsolver::MiCPPerformance solve_reduced_system(
std::shared_ptr<database::DataContainer> data,
Vector& totconc,
Vector& x
)
{
return AdimensionalSystemSolver(data, totconc).solve(x);
}
} // end namespace specmicp
#endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP
diff --git a/src/specmicp/adimensional/adimensional_system_structs.hpp b/src/specmicp/adimensional/adimensional_system_structs.hpp
index 8364fd4..eee9477 100644
--- a/src/specmicp/adimensional/adimensional_system_structs.hpp
+++ b/src/specmicp/adimensional/adimensional_system_structs.hpp
@@ -1,42 +1,94 @@
#ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP
#define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP
#include "common.hpp"
namespace specmicp {
//! \brief Options for the Adimensional Systems
struct AdimensionalSystemOptions
{
- bool conservation_water; //!< Solve the conservation of water
- index_t charge_keeper; //!< Index of species used for charge balance
bool non_ideality; //!< Solve for non ideality
scalar_t non_ideality_tolerance; //!< Tolerance for non ideality
index_t non_ideality_max_iter; //!< Max iterations fornon ideality
scalar_t under_relaxation_factor; //!< Under relaxation factor for the conservation of water
AdimensionalSystemOptions():
- conservation_water(true),
- charge_keeper(no_species),
non_ideality(true),
non_ideality_tolerance(1e-8),
non_ideality_max_iter(10),
under_relaxation_factor(0.9)
{}
};
//! \brief Type of an aqueous component equation
enum class AqueousComponentEquationType
{
NoEquation, //!< Not an equation, component is not present in the system
MassConservation, //!< Mass balance
ChargeBalance, //!< M.B. replaced by charge balance
FixedFugacity, //!< M.B. replaced by a fixed fugacity equation
FixedActivity //!< M.B. replaced by a fixed activity equation
};
+//! \brief Struct to contain information needed to solve a fix fugacity problem
+struct FixedFugacityBC
+{
+ index_t id_gas;
+ index_t id_component;
+ scalar_t value;
+};
+
+//! \brief Struct to contain information needed to solve a fix activity problem.
+struct FixedActivityBC
+{
+ index_t id_component;
+ scalar_t value;
+};
+
+//! \brief Struct to contains the "Boundary conditions" for
+struct AdimensionalSystemBC
+{
+ Vector total_concentrations; //!< Total concentrations
+ bool conservation_water; //!< True if conservation of water is solved
+ index_t charge_keeper; //!< The equation for this component is replace by the charge balance
+ std::vector<FixedFugacityBC> fixed_fugacity_bc; //!< Contains information about fixed fugacity gas
+ std::vector<FixedActivityBC> fixed_activity_bc; //!< Contains information about fixed activity component
+
+
+ AdimensionalSystemBC():
+ conservation_water(true),
+ charge_keeper(no_species)
+ {}
+
+ AdimensionalSystemBC(const Vector& total_concs):
+ total_concentrations(total_concs),
+ conservation_water(true),
+ charge_keeper(no_species)
+ {}
+
+ //! \brief Enable the conservation of water
+ void enable_conservation_water() {conservation_water = true;}
+ //! \brief Disable the conservation of water
+ void disable_conservation_water() {conservation_water = false;}
+
+ //! \brief Set the charge keeper to 'component'
+ void set_charge_keeper(index_t component) {
+ charge_keeper = component;
+ }
+
+ //! Add a fixed fugacity gas
+ void add_fixed_fugacity_gas(const FixedFugacityBC& bc) {
+ fixed_fugacity_bc.push_back(bc);
+ }
+
+ //! Add a fixed activity component
+ void add_fixed_activity_component(const FixedActivityBC& bc) {
+ fixed_activity_bc.push_back(bc);
+ }
+};
} // end namespace specmicp
#endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSTRUCTS_HPP