Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/src/specmicp/io/configuration.cpp
235eaf0b5b35master

Last Change
Raw File

configuration.cpp
View Options

/* =============================================================================
Copyright (c) 2014 - 2016
F. 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:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. 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.
3. Neither the name of the copyright holder 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 HOLDER 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 "configuration.hpp"
#include "../../utils/io/yaml.hpp"
#include "../adimensional/adimensional_system_solver_structs.hpp"
#include "../adimensional/config_default_options_solver.h"
#include "../problem_solver/formulation.hpp"
#include "../../database.hpp"
#include "../../physics/units.hpp"
#include "../../physics/laws.hpp"
#define S_UNITS "units"
#define S_UNITS_A_LENGTH "length"
#define S_SPECMICP "specmicp_options"
#define S_SPECMICP_A_MAX_ITER "maximum_iteration"
#define S_SPECMICP_A_RES_TOL "residual_tolerance"
#define S_SPECMICP_A_STEP_TOL "step_tolerance"
#define S_SPECMICP_A_MAX_STEP_LENGTH "maximum_step_length"
#define S_SPECMICP_A_MAX_STEP_MAX_ITER "maximum_step_max_iteration"
#define S_SPECMICP_A_SCALING "enable_scaling"
#define S_SPECMICP_A_NONMONOTONE "enable_nonmonotone_linesearch"
#define S_SPECMICP_A_DESCENT_DIRECTION "factor_descent_direction"
#define S_SPECMICP_A_COND_CHECK "threshold_condition_check"
#define S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH "threshold_cycling_linesearch"
#define S_SPECMICP_A_NONIDEAL_TOL "non_ideality_tolerance"
#define S_SPECMICP_A_CUTOFF_TOT_CONC "cutoff_total_concentration"
#define S_FORMULATION "formulation"
#define S_FORMULATION_A_SOLUTION "solution"
#define S_FORMULATION_A_AQUEOUS "aqueous"
#define S_FORMULATION_A_MINERALS "solid_phases"
#define S_FORMULATION_A_AMOUNT "amount"
#define S_FORMULATION_A_UNIT "unit"
#define S_FORMULATION_A_LABEL "label"
#define S_CONSTRAINTS "constraints"
#define S_CONSTRAINTS_A_CHARGEKEEPER "charge_keeper"
#define S_CONSTRAINTS_A_FIXEDACTIVITY "fixed_activity"
#define S_CONSTRAINTS_A_FIXEDMOLALITY "fixed_molality"
#define S_CONSTRAINTS_A_FIXEDFUGACITY "fixed_fugacity"
#define S_CONSTRAINTS_A_SATURATED_SYSTEM "saturated_system"
#define S_CONSTRAINTS_A_CONSERVATION_WATER "conservation_water"
#define S_CONSTRAINTS_A_SURFACE_MODEL "surface_model"
#define S_CONSTRAINTS_A_SURFACE_CONCENTRATION "surface_site_concentration"
#define S_CONSTRAINTS_A_LABEL "label"
#define S_CONSTRAINTS_A_AMOUNT "amount"
#define S_CONSTRAINTS_A_AMOUNTLOG "amount_log"
#define S_CONSTRAINTS_A_LABEL_COMPONENT "label_component"
#define S_CONSTRAINTS_A_LABEL_GAS "label_gas"
namespace specmicp {
namespace io {
// Additional declarations
// =======================
void configure_specmicp_constraints_fixed_activity(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints_fixed_activity,
const RawDatabasePtr& raw_db
);
void configure_specmicp_constraints_fixed_fugacity(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints_fixed_fugacity,
const RawDatabasePtr& raw_db
);
void configure_specmicp_constraints_fixed_molality(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints_fixed_molality,
const RawDatabasePtr& raw_db
);
scalar_t molar_mass_aqueous_unknown_class(const std::string& label,
const RawDatabasePtr& raw_db);
scalar_t molar_mass_solid_unknown_class(const std::string& label,
const RawDatabasePtr& raw_db
);
scalar_t molar_volume_solid_unknown_class(const std::string& label,
const RawDatabasePtr& raw_db
);
// Implementation
// ==============
void configure_specmicp_options(
AdimensionalSystemSolverOptions& options,
const YAML::Node& configuration
)
{
check_mandatory_yaml_node(configuration, S_SPECMICP, "__main__");
const YAML::Node& conf = configuration[S_SPECMICP];
// fixme units
//options.units_set = the_units_set;
options.solver_options.fvectol = get_yaml_optional<double>(conf,
S_SPECMICP_A_RES_TOL, S_SPECMICP, SPECMICP_DEFAULT_RES_TOL);
options.solver_options.steptol = get_yaml_optional<double>(conf,
S_SPECMICP_A_STEP_TOL, S_SPECMICP, SPECMICP_DEFAULT_STEP_TOL);
options.solver_options.max_iter = get_yaml_optional<int>(conf,
S_SPECMICP_A_MAX_ITER, S_SPECMICP, SPECMICP_DEFAULT_MAX_ITER);
options.solver_options.maxstep = get_yaml_optional<double>(conf,
S_SPECMICP_A_MAX_STEP_LENGTH, S_SPECMICP, SPECMICP_DEFAULT_MAX_STEP);
options.solver_options.maxiter_maxstep = get_yaml_optional<int>(conf,
S_SPECMICP_A_MAX_STEP_MAX_ITER, S_SPECMICP,
std::min(SPECMICP_DEFAULT_MAX_ITER_MAXSTEP, options.solver_options.max_iter)
);
options.solver_options.use_scaling = get_yaml_optional<bool>(conf,
S_SPECMICP_A_SCALING, S_SPECMICP, SPECMICP_DEFAULT_USE_SCALING);
options.solver_options.non_monotone_linesearch = get_yaml_optional<bool>(conf,
S_SPECMICP_A_NONMONOTONE, S_SPECMICP, SPECMICP_DEFAULT_USE_NONMONOTONE_LSEARCH);
options.solver_options.factor_descent_condition = get_yaml_optional<double>(conf,
S_SPECMICP_A_DESCENT_DIRECTION, S_SPECMICP, SPECMICP_DEFAULT_FACTOR_DESC_COND
);
options.solver_options.condition_limit = get_yaml_optional<double>(conf,
S_SPECMICP_A_COND_CHECK, S_SPECMICP, SPECMICP_DEFAULT_COND_THRESHOLD);
options.solver_options.threshold_cycling_linesearch = get_yaml_optional<double>(conf,
S_SPECMICP_A_TRSHOLD_CYCLING_LSEARCH, S_SPECMICP, SPECMICP_DEFAULT_TRHSOLD_CYCLING_LSEARCH
);
options.system_options.non_ideality_tolerance = get_yaml_optional<double>(conf,
S_SPECMICP_A_NONIDEAL_TOL, S_SPECMICP, SPECMICP_DEFAULT_NONIDEAL_TOL);
options.system_options.cutoff_total_concentration = get_yaml_optional<double>(conf,
S_SPECMICP_A_CUTOFF_TOT_CONC, S_SPECMICP, SPECMICP_DEFAULT_CUTOFF_TOT_CONC);
if (configuration[S_UNITS])
{
const YAML::Node& conf = configuration[S_UNITS];
if (conf[S_UNITS_A_LENGTH])
{
const std::string length_u = conf[S_UNITS_A_LENGTH].as<std::string>();
if (length_u == "centimeter")
options.units_set.length = units::LengthUnit::centimeter;
else if (length_u == "decimeter")
options.units_set.length = units::LengthUnit::decimeter;
else if (length_u == "meter")
options.units_set.length = units::LengthUnit::meter;
else
throw std::invalid_argument("Unknown length unit : '"+length_u+"'.\n" +
"Valid options : 'centimeter', 'decimeter', 'meter'.");
}
}
}
void configure_specmicp_formulation(
Formulation& formulation,
const YAML::Node& conf_formulation,
const RawDatabasePtr& raw_db,
const units::UnitsSet& units_sets
)
{
configure_specmicp_formulation_solution(formulation, conf_formulation, units_sets);
configure_specmicp_formulation_aqueous(formulation, conf_formulation, raw_db, units_sets);
configure_specmicp_formulation_solid(formulation, conf_formulation, raw_db, units_sets);
}
void configure_specmicp_formulation_solution(
Formulation& formulation,
const YAML::Node& conf_formulation,
const units::UnitsSet& _)
{
check_mandatory_yaml_node(conf_formulation, S_FORMULATION_A_SOLUTION, S_FORMULATION);
const YAML::Node& sol_conf = conf_formulation[S_FORMULATION_A_SOLUTION];
check_mandatory_yaml_node(sol_conf, S_FORMULATION_A_AMOUNT, S_FORMULATION_A_SOLUTION);
scalar_t factor=1.0;
if (sol_conf[S_FORMULATION_A_UNIT])
{
units::AmountUnit unit;
units::parse_amount_unit(sol_conf[S_FORMULATION_A_UNIT].as<std::string>(), unit);
if (unit.type != units::AmountUnitType::Mass )
{
throw std::invalid_argument("'" + sol_conf[S_FORMULATION_A_UNIT].as<std::string>() +
"'' is not a mass unit (in the solution configuration) !");
}
factor = unit.factor_si;
}
formulation.set_mass_solution(sol_conf[S_FORMULATION_A_AMOUNT].as<double>()*factor);
}
void configure_specmicp_formulation_aqueous(
Formulation& formulation,
const YAML::Node& conf_formulation,
const RawDatabasePtr& raw_db,
const units::UnitsSet& _)
{
check_mandatory_yaml_node(conf_formulation, S_FORMULATION_A_AQUEOUS, S_FORMULATION);
const YAML::Node& aq_conf = conf_formulation[S_FORMULATION_A_AQUEOUS];
for (auto aqueous_conf: aq_conf)
{
check_mandatory_yaml_node(aqueous_conf, S_FORMULATION_A_LABEL, S_FORMULATION_A_AQUEOUS);
check_mandatory_yaml_node(aqueous_conf, S_FORMULATION_A_AMOUNT, S_FORMULATION_A_AQUEOUS);
check_mandatory_yaml_node(aqueous_conf, S_FORMULATION_A_UNIT, S_FORMULATION_A_AQUEOUS);
std::string label = aqueous_conf[S_FORMULATION_A_LABEL].as<std::string>();
scalar_t factor = 1.0;
units::AmountUnit unit;
units::parse_amount_unit(aqueous_conf[S_FORMULATION_A_UNIT].as<std::string>(), unit);
if (unit.type == units::AmountUnitType::Molality)
{
factor = unit.factor_si;
}
else if (unit.type == units::AmountUnitType::MoleConcentration)
{
factor = unit.factor_si/laws::density_water(units::celsius(25.0), units::LengthUnit::meter, units::MassUnit::kilogram);
}
else if (unit.type == units::AmountUnitType::MassConcentration)
{
const scalar_t molar_mass = molar_mass_aqueous_unknown_class(label, raw_db);
factor = unit.factor_si / molar_mass
/ laws::density_water(units::celsius(25.0), units::LengthUnit::meter, units::MassUnit::kilogram);
}
else if (unit.type == units::AmountUnitType::Mass)
{
const scalar_t molar_mass = molar_mass_aqueous_unknown_class(label, raw_db);
factor = unit.factor_si / molar_mass / formulation.mass_solution;
}
else if (unit.type == units::AmountUnitType::NumberOfMoles)
{
factor = unit.factor_si / formulation.mass_solution;
}
else
{
throw std::invalid_argument("'"+aqueous_conf[S_FORMULATION_A_UNIT].as<std::string>() +
"' is not a valid unit for an aqueous species");
}
const scalar_t molality = factor*aqueous_conf[S_FORMULATION_A_AMOUNT].as<scalar_t>();
formulation.add_aqueous_species(label, molality);
}
}
void configure_specmicp_formulation_solid(
Formulation& formulation,
const YAML::Node& conf_formulation,
const RawDatabasePtr& raw_db,
const units::UnitsSet& unit_set)
{
check_mandatory_yaml_node(conf_formulation, S_FORMULATION_A_MINERALS, S_FORMULATION);
const YAML::Node& minerals_conf = conf_formulation[S_FORMULATION_A_MINERALS];
for (auto solid_conf : minerals_conf)
{
check_mandatory_yaml_node(solid_conf, S_FORMULATION_A_LABEL, S_FORMULATION_A_AQUEOUS);
check_mandatory_yaml_node(solid_conf, S_FORMULATION_A_AMOUNT, S_FORMULATION_A_AQUEOUS);
check_mandatory_yaml_node(solid_conf, S_FORMULATION_A_UNIT, S_FORMULATION_A_AQUEOUS);
const std::string label = solid_conf[S_FORMULATION_A_LABEL].as<std::string>();
scalar_t factor = 1.0;
units::AmountUnit unit;
units::parse_amount_unit(solid_conf[S_FORMULATION_A_UNIT].as<std::string>(), unit);
// convert to mole concentration is SI
if (unit.type == units::AmountUnitType::MoleConcentration or
unit.type == units::AmountUnitType::NumberOfMoles)
{
factor = unit.factor_si;
}
else if (unit.type == units::AmountUnitType::MassConcentration or
unit.type == units::AmountUnitType::Mass
)
{
const scalar_t molar_mass = molar_mass_solid_unknown_class(label, raw_db);
factor = unit.factor_si / molar_mass;
}
else if (unit.type == units::AmountUnitType::Volume)
{
const scalar_t molar_volume = molar_volume_solid_unknown_class(label, raw_db); // may raise an error
factor = unit.factor_si / molar_volume;
}
else
{
throw std::invalid_argument("'"+solid_conf[S_FORMULATION_A_UNIT].as<std::string>() +
"' is not a valid unit for a solid phase");
}
// set to the correct units
switch (unit_set.length) {
case units::LengthUnit::decimeter:
factor *= 1e-3;
break;
case units::LengthUnit::centimeter:
factor *= 1e-6;
break;
case units::LengthUnit::meter:
break; // just to remove warning
}
const scalar_t mole_conc = factor*solid_conf[S_FORMULATION_A_AMOUNT].as<scalar_t>();
formulation.add_mineral(label, mole_conc);
}
}
void configure_specmicp_constraints(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints,
const RawDatabasePtr& raw_db
)
{
// Component equations
// ===================
// Charge keeper
// -------------
if (conf_constraints[S_CONSTRAINTS_A_CHARGEKEEPER])
{
const std::string label = conf_constraints[S_CONSTRAINTS_A_CHARGEKEEPER].as<std::string>();
const index_t id = raw_db->get_id_component(label);
if (id == no_species)
{
throw std::invalid_argument("'"+label+"' is not a valid component (charge_keeper).");
}
constraints.set_charge_keeper(id);
}
// Fixed activity
// --------------
if (conf_constraints[S_CONSTRAINTS_A_FIXEDACTIVITY])
{
configure_specmicp_constraints_fixed_activity(constraints, conf_constraints[S_CONSTRAINTS_A_FIXEDACTIVITY], raw_db);
}
// Fixed fugacity
// --------------
if (conf_constraints[S_CONSTRAINTS_A_FIXEDFUGACITY])
{
configure_specmicp_constraints_fixed_fugacity(constraints, conf_constraints[S_CONSTRAINTS_A_FIXEDFUGACITY], raw_db);
}
// Fixed molality
// --------------
if (conf_constraints[S_CONSTRAINTS_A_FIXEDMOLALITY])
{
configure_specmicp_constraints_fixed_molality(constraints, conf_constraints[S_CONSTRAINTS_A_FIXEDMOLALITY], raw_db);
}
// Water
// =====
if (conf_constraints[S_CONSTRAINTS_A_CONSERVATION_WATER])
{
if (conf_constraints[S_CONSTRAINTS_A_SATURATED_SYSTEM] and
conf_constraints[S_CONSTRAINTS_A_SATURATED_SYSTEM].as<bool>())
{
throw std::invalid_argument("Options '" S_CONSTRAINTS_A_CONSERVATION_WATER "' and '"
S_CONSTRAINTS_A_SATURATED_SYSTEM "' cannot be used at the same time");
}
const bool value = conf_constraints[S_CONSTRAINTS_A_CONSERVATION_WATER].as<bool>();
if (value)
constraints.enable_conservation_water();
else
constraints.disable_conservation_water();
}
if (conf_constraints[S_CONSTRAINTS_A_SATURATED_SYSTEM])
{
const bool value = conf_constraints[S_CONSTRAINTS_A_SATURATED_SYSTEM].as<bool>();
if (value)
constraints.set_saturated_system();
}
// Surface model
// =============
if (conf_constraints[S_CONSTRAINTS_A_SURFACE_MODEL])
{
const bool is_enabled = conf_constraints[S_CONSTRAINTS_A_SURFACE_MODEL].as<bool>();
if (is_enabled)
{
check_mandatory_yaml_node(conf_constraints, S_CONSTRAINTS_A_SURFACE_CONCENTRATION, S_CONSTRAINTS);
const scalar_t value = conf_constraints[S_CONSTRAINTS_A_SURFACE_CONCENTRATION].as<scalar_t>();
constraints.enable_surface_model(value);
}
else
{
constraints.disable_surface_model();
}
}
}
void configure_specmicp_constraints_fixed_activity(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints_fixed_activity,
const RawDatabasePtr& raw_db
)
{
for (auto fixact_conf: conf_constraints_fixed_activity)
{
check_mandatory_yaml_node(fixact_conf, S_CONSTRAINTS_A_LABEL, S_CONSTRAINTS_A_FIXEDACTIVITY);
const std::string label = fixact_conf[S_CONSTRAINTS_A_LABEL].as<std::string>();
const index_t id = raw_db->get_id_component(label);
if (id == no_species)
{
throw std::invalid_argument("'"+label+"' is not a valid component (fixed_activity).");
}
scalar_t value;
if (fixact_conf[S_CONSTRAINTS_A_AMOUNT])
{
value = std::log10(fixact_conf[S_CONSTRAINTS_A_AMOUNT].as<scalar_t>());
}
else if (fixact_conf[S_CONSTRAINTS_A_AMOUNTLOG])
{
value = fixact_conf[S_CONSTRAINTS_A_AMOUNTLOG].as<scalar_t>();
}
else
{
throw std::invalid_argument("Either '" S_CONSTRAINTS_A_AMOUNT "' or '" S_CONSTRAINTS_A_AMOUNTLOG
"' are required in section '" S_CONSTRAINTS_A_FIXEDACTIVITY "'.");
}
constraints.add_fixed_activity_component(id, value);
}
}
void configure_specmicp_constraints_fixed_fugacity(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints_fixed_fugacity,
const RawDatabasePtr& raw_db
)
{
for(auto fixfug_conf: conf_constraints_fixed_fugacity)
{
check_mandatory_yaml_node(fixfug_conf, S_CONSTRAINTS_A_LABEL_COMPONENT, S_CONSTRAINTS_A_FIXEDFUGACITY);
check_mandatory_yaml_node(fixfug_conf, S_CONSTRAINTS_A_LABEL_GAS, S_CONSTRAINTS_A_FIXEDFUGACITY);
const std::string& label_c = fixfug_conf[S_CONSTRAINTS_A_LABEL_COMPONENT].as<std::string>();
const index_t id_c = raw_db->get_id_component(label_c);
if (id_c == no_species)
{
throw std::invalid_argument("'"+label_c+"' is not a valid component (fixed_fugacity)");
}
const std::string& label_g = fixfug_conf[S_CONSTRAINTS_A_LABEL_GAS].as<std::string>();
const index_t id_g = raw_db->get_id_gas(label_g);
if (id_g == no_species)
{
throw std::invalid_argument("'"+label_g+"' is not a valid gas (fixed_fugacity)");
}
scalar_t value;
if (fixfug_conf[S_CONSTRAINTS_A_AMOUNT])
{
value = std::log10(fixfug_conf[S_CONSTRAINTS_A_AMOUNT].as<scalar_t>());
}
else if (fixfug_conf[S_CONSTRAINTS_A_AMOUNTLOG])
{
value = fixfug_conf[S_CONSTRAINTS_A_AMOUNTLOG].as<scalar_t>();
}
else
{
throw std::invalid_argument("Either '" S_CONSTRAINTS_A_AMOUNT "' or '" S_CONSTRAINTS_A_AMOUNTLOG
"' are required in section '" S_CONSTRAINTS_A_FIXEDFUGACITY "'.");
}
constraints.add_fixed_fugacity_gas(id_g, id_c, value);
}
}
void configure_specmicp_constraints_fixed_molality(
AdimensionalSystemConstraints& constraints,
const YAML::Node& conf_constraints_fixed_molality,
const RawDatabasePtr& raw_db
)
{
for (auto fixmol_conf: conf_constraints_fixed_molality)
{
check_mandatory_yaml_node(fixmol_conf, S_CONSTRAINTS_A_LABEL, S_CONSTRAINTS_A_FIXEDMOLALITY);
const std::string& label = fixmol_conf[S_CONSTRAINTS_A_LABEL].as<std::string>();
const index_t id = raw_db->get_id_component(label);
if (id == no_species)
{
throw std::invalid_argument("'"+label+"' is not a valid component (fixed_molality)");
}
scalar_t value;
if (fixmol_conf[S_CONSTRAINTS_A_AMOUNT])
{
value = std::log10(fixmol_conf[S_CONSTRAINTS_A_AMOUNT].as<scalar_t>());
}
else if (fixmol_conf[S_CONSTRAINTS_A_AMOUNTLOG])
{
value = fixmol_conf[S_CONSTRAINTS_A_AMOUNTLOG].as<scalar_t>();
}
else
{
throw std::invalid_argument("Either '" S_CONSTRAINTS_A_AMOUNT "' or '" S_CONSTRAINTS_A_AMOUNTLOG
"' are required in section '" S_CONSTRAINTS_A_FIXEDMOLALITY "'.");
}
constraints.add_fixed_molality_component(id, value);
}
}
scalar_t molar_mass_aqueous_unknown_class(const std::string& label,
const RawDatabasePtr& raw_db)
{
scalar_t molar_mass;
if (raw_db->get_id_component(label) != no_species)
{
molar_mass = raw_db->molar_mass_basis(raw_db->get_id_component(label), units::MassUnit::kilogram);
}
else if (raw_db->get_id_aqueous(label) != no_species)
{
molar_mass = raw_db->molar_mass_aqueous(raw_db->get_id_aqueous(label), units::MassUnit::kilogram);
}
else if (raw_db->get_id_compound(label) != no_species)
{
molar_mass = raw_db->molar_mass_compound(raw_db->get_id_compound(label), units::MassUnit::kilogram);
}
else
{
throw std::invalid_argument("'"+label+"is not a valid species");
}
return molar_mass;
}
scalar_t molar_mass_solid_unknown_class(const std::string& label,
const RawDatabasePtr& raw_db
)
{
scalar_t molar_mass;
if (raw_db->get_id_mineral(label) != no_species)
{
molar_mass = raw_db->molar_mass_mineral(raw_db->get_id_mineral(label), units::MassUnit::kilogram);
}
else if (raw_db->get_id_mineral_kinetic(label) != no_species)
{
molar_mass = raw_db->molar_mass_mineral_kinetic(raw_db->get_id_mineral_kinetic(label), units::MassUnit::kilogram);
}
else
{
throw std::invalid_argument("'"+label+"' is not a valid solid phase.");
}
return molar_mass;
}
scalar_t molar_volume_solid_unknown_class(const std::string& label,
const RawDatabasePtr& raw_db
)
{
scalar_t molar_volume = -1;
if (raw_db->get_id_mineral(label) != no_species)
{
molar_volume = raw_db->molar_volume_mineral(raw_db->get_id_mineral(label), units::LengthUnit::meter);
}
else if (raw_db->get_id_mineral_kinetic(label) != no_species)
{
molar_volume = raw_db->molar_volume_mineral_kinetic(raw_db->get_id_mineral_kinetic(label), units::LengthUnit::meter);
}
else
{
throw std::invalid_argument("'"+label+"' is not a valid solid phase.");
}
if (molar_volume < 0)
{
throw std::runtime_error("The molar volume for solid phase '"+label+"' is not defined");
}
return molar_volume;
}
} //end namespace io
} //end namespace specmicp
c4science · Help