Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F63790889
configuration.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, May 22, 12:21
Size
23 KB
Mime Type
text/x-c
Expires
Fri, May 24, 12:21 (2 d)
Engine
blob
Format
Raw Data
Handle
17819712
Attached To
rSPECMICP SpecMiCP / ReactMiCP
configuration.cpp
View Options
/*-------------------------------------------------------------------------------
Copyright (c) 2015 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 "../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_fugactiy"
#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_COND_THRESHOLD
);
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
;
}
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_component
(
label_c
);
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
Event Timeline
Log In to Comment