Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61387576
adimensional_system.hpp
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
Mon, May 6, 08:41
Size
31 KB
Mime Type
text/x-c++
Expires
Wed, May 8, 08:41 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17504986
Attached To
rSPECMICP SpecMiCP / ReactMiCP
adimensional_system.hpp
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. *
============================================================================= */
#ifndef SPECMICP_ADIMENSIONALSYSTEM_HPP
#define SPECMICP_ADIMENSIONALSYSTEM_HPP
//! \file adimensional_system.hpp The MiCP program to solve speciation
#include "../../types.hpp"
#ifndef SPECMICP_DATABASE_HPP
#include "database.hpp"
#endif
#include "../../micpsolver/micpprog.hpp"
#include "../../physics/units.hpp"
#include "../../utils/options_handler.hpp"
#include "adimensional_system_numbering.hpp"
#include "adimensional_system_structs.hpp"
//! \namespace specmicp main namespace for the SpecMiCP solver
namespace
specmicp
{
struct
AdimensionalSystemSolution
;
//! \brief The equilibrium speciation problem
//!
//! Represent the equilibrium speciation problem as a Mixed Complementarity program
//!
//! The main variables are :
//! - \f$phi_w\f$ the volume fraction of water (total saturation)
//! - \f$b_i\f$ the molality of aqueous component
//! - \f$phi_m\f$ the volume fraction of mineral
//! - \f$s_f\f$ the concentration (molality) of free sorption sites
//!
//! The solid phase equilibrium is solved by using a complementarity formulation
//! \f[ S^t_m \geq 0, \; - \mathrm{SI}_m \geq 0 , \; \mathrm{and} \; - S^t_m \mathrm{SI}_m = 0 \f]
//!
//! The secondary variables (molalities of secondary aqueous species, gas fugacities, ionic strength, ...)
//! are solved using a fixed point iteration at the beginning of each iteration.
//!
//! This class should only be used through the AdimensionalSystemSolver.
class
SPECMICP_DLL_LOCAL
AdimensionalSystem
:
public
micpsolver
::
MiCPProg
<
AdimensionalSystem
>
,
public
AdimemsionalSystemNumbering
,
public
OptionsHandler
<
AdimensionalSystemOptions
>
,
private
units
::
UnitBaseClass
{
public
:
//! \brief Create a reduced system
//!
//! \param ptrdata Shared pointer to the thermodynamic database
//! \param constraints Constraints to apply to the system
//! \param options Numerical options, mainly related to the secondary variables
//! \param units_set The set of units
AdimensionalSystem
(
RawDatabasePtr
&
ptrdata
,
const
AdimensionalSystemConstraints
&
constraints
,
const
AdimensionalSystemOptions
&
options
=
AdimensionalSystemOptions
(),
const
units
::
UnitsSet
&
units_set
=
units
::
UnitsSet
()
);
//! \brief Create a reduced system, initialize the system using a previous solution
//!
//! \param ptrdata Shared pointer to the thermodynamic database
//! \param constraints Constraints to apply to the system
//! \param previous_solution The previous solution to use for initialization
//! \param options Numerical options, mainly related to the secondary variables
//! \param units_set The set of units
AdimensionalSystem
(
RawDatabasePtr
&
ptrdata
,
const
AdimensionalSystemConstraints
&
constraints
,
const
AdimensionalSystemSolution
&
previous_solution
,
const
AdimensionalSystemOptions
&
options
=
AdimensionalSystemOptions
(),
const
units
::
UnitsSet
&
units_set
=
units
::
UnitsSet
()
);
// Variables
// ==========
//! \name Variables
//! \brief How many ?
//!
//! @{
//! \brief Return the total number of variables
index_t
total_variables
()
const
{
return
m_equations
.
nb_tot_variables
;}
//! \brief Return the number of 'free' variables (not subject to complementarity conditions)
index_t
nb_free_variables
()
const
{
return
m_equations
.
nb_free_variables
;}
//! \brief Return the number of variables subjected to the complementarity conditions
index_t
nb_complementarity_variables
()
const
{
return
m_equations
.
nb_complementarity_variables
;}
//! @}
//! \brief Return a pointer to the database
RawDatabasePtr
get_database
()
{
return
m_data
;}
// The linear system
// ==================
//! \name Residuals and Jacobian
//! \brief Compute the residuals and the jacobian
//!
//! @{
//! \brief Return the residuals
//!
//! \param x Vector of the main variables
//! \param[out] residual Vector containing the residuals
void
get_residuals
(
const
Vector
&
x
,
Vector
&
residual
);
//! \brief Return the jacobian
//!
//! \param x Vector of the main variables
//! \param[out] jacobian Dense matrix representing the jacobian
void
get_jacobian
(
Vector
&
x
,
Matrix
&
jacobian
);
//! \brief Return the residual for the water conservation equation
//!
//! \param x Vector of the main variables
scalar_t
residual_water_conservation
(
const
Vector
&
x
)
const
;
//! \brief Return the residual for the water saturation equation
//!
//! \param x Vector of the main variables
scalar_t
residual_water_saturation
(
const
Vector
&
x
)
const
;
//! \brief Return the residual for the conservation equation of component (i)
//!
//! For water (i==0) use the method : 'residual_water'
//! \param x Vector of the main variables
//! \param i Index of the component (in the database)
scalar_t
residual_component
(
const
Vector
&
x
,
index_t
i
)
const
;
//! \brief Compute the residual for the surface sorption
//! \param x Vector of the main variables
scalar_t
residual_surface
(
const
Vector
&
x
)
const
;
//! \brief Equilibrium condition for the minerals
//! \param x Vector of the main variables
//! \param m Index of the mineral (in the database)
scalar_t
residual_mineral
(
const
Vector
&
x
,
index_t
m
)
const
;
//! \brief Return the residual of the charge conservation equation
//!
//! This equation may be used instead of a mass conservation equation to
//! explicitely enforce the electroneutrality.
//! The component corresponding to this equation is called the charge keeper
//! \param x Vector of the main variables
scalar_t
residual_charge_conservation
(
const
Vector
&
x
)
const
;
//! \brief Return the residual corresponding to a fixed fugacity equation
//! \param x Vector of the main variables
//! \param component Index of the corresponding component (in the database)
scalar_t
residual_fixed_fugacity
(
const
Vector
&
x
,
index_t
component
)
const
;
//! \brief Return the residual corresponding to a fixed activity equation
//!
//! If 'component ' is charged, then a charge keeper should be set.
//! \param x Vector of the main variables
//! \param component Index of the corresponding component (in the database)
scalar_t
residual_fixed_activity
(
const
Vector
&
x
,
index_t
component
)
const
;
//! \brief Return the residual corresponding to a fixed molality equation
//!
//! If 'component ' is charged, then a charge keeper should be set.
//! \param x Vector of the main variables
//! \param component Index of the corresponding component (in the database)
scalar_t
residual_fixed_molality
(
const
Vector
&
x
,
index_t
component
)
const
;
//! \brief Return the residual corresponding to the electron equation
scalar_t
residual_electron
(
const
Vector
&
x
)
const
;
//! \brief Compute the jacobian analytically
void
analytical_jacobian
(
Vector
&
x
,
Matrix
&
jacobian
);
//! \brief Compute the jacobian using finite difference
void
finite_difference_jacobian
(
Vector
&
x
,
Matrix
&
jacobian
);
//! \brief Return a scale factor to avoid negative mass during Newton's iteration
//!
//! \param x Vector of the main variables
//! \param update Update of the current iteration
double
max_lambda
(
const
Vector
&
x
,
const
Vector
&
update
);
//! @}
// Getters
// #######
// Equation id access
// ==================
//! \name Equation Id
//! \brief The equation id is the row of the equation in the residuals/jacobian
//!
//! The degree of freedom getter function are inherited from specmicp::AdimensionalSystemNumbering
//! @{
//! \brief Return the equation id
//! \param i Index of the main variable (in the set of the main variables)
index_t
ideq
(
index_t
i
)
const
{
return
m_equations
.
ideq
[
i
];}
//! \brief Return the equation number for conservation of water
index_t
ideq_w
()
const
{
return
m_equations
.
ideq
[
dof_water
()];}
//! \brief Return the equation number of component 'i'
//! \param i Index of the component (in the database)
index_t
ideq_paq
(
index_t
i
)
const
{
specmicp_assert
(
i
<
m_data
->
nb_component
());
return
m_equations
.
ideq
[
dof_component
(
i
)];
}
//! \brief Return the equation number of the electron equation
index_t
ideq_electron
()
const
{
return
m_equations
.
ideq
[
dof_electron
()];
}
//! \brief Return the equation number of the free surface concentration
index_t
ideq_surf
()
const
{
return
m_equations
.
ideq
[
dof_surface
()];
}
//! \brief Return the equation number of mineral 'm'
//! \param m Index of the mineral (in the database)
index_t
ideq_min
(
index_t
m
)
const
{
specmicp_assert
(
m
<
m_data
->
nb_mineral
());
return
m_equations
.
ideq
[
dof_mineral
(
m
)];
}
//! @}
//! \name Equation type
//! \brief Which equation is actually solved
//!
//! @{
//! \brief Return the type of equation used for the solvent
//!
//! Can be no equation, mass conservation, saturation of the system, ...
//!
//! \sa #WaterEquationType, #aqueous_component_equation_type, #AqueousComponentEquationType
WaterEquationType
water_equation_type
()
const
{
return
static_cast
<
WaterEquationType
>
(
m_equations
.
type_equation
(
dof_water
()));
}
//! \brief Return the type of equation for aqueous species 'component'
//!
//! For water used the method #water_equation_type
//! \param component Index of the aqueous component (in the database)
//!
//! \sa #AqueousComponentEquationType, water_equation_type, #WaterEquationType
AqueousComponentEquationType
aqueous_component_equation_type
(
index_t
component
)
const
{
specmicp_assert
(
component
>
0
and
component
<
m_data
->
nb_component
());
return
static_cast
<
AqueousComponentEquationType
>
(
m_equations
.
type_equation
(
dof_component
(
component
)));
}
//! \brief Return the type of equation for the electron
ElectronEquationType
electron_equation_type
()
const
{
return
static_cast
<
ElectronEquationType
>
(
m_equations
.
type_equation
(
dof_electron
()));
}
//! \brief Return true if an equation exist for this component
//!
//! \param component Index of the component (in the database)
bool
is_active_component
(
index_t
component
)
const
{
specmicp_assert
(
component
<
m_data
->
nb_component
()
and
component
>
no_species
);
return
m_equations
.
id_equation
(
dof_component
(
component
))
!=
no_equation
;
}
//! \brief Return the surface model
SurfaceEquationType
surface_model
()
const
{
if
(
ideq_surf
()
==
no_equation
)
return
SurfaceEquationType
::
NoEquation
;
return
SurfaceEquationType
::
Equilibrium
;
}
//! \brief Return true if 'component' is the charge keeper
//!
//! The charge keeper is the component added or removed to satisfy the charge balance
//! \param component Index of the aqueous component (in the database)
bool
is_charge_keeper
(
index_t
component
)
{
return
(
aqueous_component_equation_type
(
component
)
==
AqueousComponentEquationType
::
ChargeBalance
);
}
//! \brief Return the component that does not exist in the solution (inactive dof)
const
std
::
vector
<
index_t
>&
get_non_active_component
()
{
return
m_equations
.
nonactive_component
;}
//! @}
// Variables
// =========
//! \name2 Main variables
//! \brief Access to the main variables
//!
//! @{
//! \brief Return the volume fraction of water
//!
//! \param x Vector of the main variables
scalar_t
volume_fraction_water
(
const
Vector
&
x
)
const
;
//! \brief Return the density of water
scalar_t
density_water
()
const
;
//! \brief Return the volume fraction of a mineral
//!
//! \param x Vector of the main variables
//! \param mineral Index of the mineral (in the database)
scalar_t
volume_fraction_mineral
(
const
Vector
&
x
,
index_t
mineral
)
const
;
//! \brief Return the volume of minerals
//!
//! \param mineral Index of the mineral (in the database)
scalar_t
molar_volume_mineral
(
index_t
mineral
)
const
{
return
m_scaling_molar_volume
*
m_data
->
unsafe_molar_volume_mineral
(
mineral
);
}
//! \brief Return the sum of saturation of the minerals
//!
//! This corresponds to the volume fraction of the solid phases (minus the inert phase)
//! \param x Vector of the main variables
scalar_t
sum_volume_fraction_minerals
(
const
Vector
&
x
)
const
;
//! \brief Return the porosity
//!
//! Computed 'on-the-fly'
//!
//! \param x Vector of the main variables
scalar_t
porosity
(
const
Vector
&
x
)
const
{
return
1
-
sum_volume_fraction_minerals
(
x
)
-
m_inert_volume_fraction
;
}
//! \brief Return the log_10 of the component molality
//!
//! \param x Vector of the main variables
//! \param component Index of the aqueous component (in the database)
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 molality of 'component'
//!
//! \param x Vector of the main variables
//! \param component Index of the aqueous component (in the database)
scalar_t
component_molality
(
const
Vector
&
x
,
index_t
component
)
const
{
return
pow10
(
log_component_molality
(
x
,
component
));
}
//! \brief Return the concentration of free sorption site
scalar_t
free_sorption_site_concentration
(
const
Vector
&
x
)
const
{
return
pow10
(
log_free_sorption_site_concentration
(
x
));
}
//! \brief Return the log_10 of the free sorption site concentration
scalar_t
log_free_sorption_site_concentration
(
const
Vector
&
x
)
const
{
specmicp_assert
(
ideq_surf
()
!=
no_equation
);
return
x
(
ideq_surf
());
}
//! \brief log activity of the electron
scalar_t
log_activity_electron
(
const
Vector
&
x
)
const
{
specmicp_assert
(
ideq_electron
()
!=
no_equation
);
return
x
(
ideq_electron
());
}
//! \brief return the activity of the electro
scalar_t
activity_electron
(
const
Vector
&
x
)
const
{
return
pow10
(
log_activity_electron
(
x
));
}
//! @}
//! \name Secondary variables
//! \brief Access to the secondary variables
//!
//! @{
//! \brief Return the concentration of secondary species 'aqueous'
//!
//! \param aqueous Index of the secondary aqueous species (in the database)
scalar_t
secondary_molality
(
index_t
aqueous
)
const
{
if
(
not
is_aqueous_active
(
aqueous
))
return
0.0
;
return
m_second
.
secondary_molalities
(
dof_aqueous
(
aqueous
));
}
//! \brief Return true if 'aqueous' is active
//!
//! i.e. Return true if all its component are present in the system
//! \param aqueous Index of the secondary aqueous species (in the database)
bool
is_aqueous_active
(
index_t
aqueous
)
const
{
return
m_equations
.
active_aqueous
[
dof_aqueous
(
aqueous
)];
}
//! \brief Return log_10(γ_i) where i is a component
//!
//! γ is the activity coefficient
//! \param component Index of the aqueous component (in the database)
scalar_t
log_gamma_component
(
index_t
component
)
const
{
return
m_second
.
loggamma
(
dof_component_gamma
(
component
));
}
//! \brief Return log_10(γ_j) where j is a secondary aqueous species
//!
//! γ is the activity coefficient
//! \param aqueous Index of the secondary aqueous species (in the database)
scalar_t
log_gamma_secondary
(
index_t
aqueous
)
const
{
return
m_second
.
loggamma
(
dof_aqueous_gamma
(
aqueous
));
}
//! \brief Return the ionic strength
scalar_t
ionic_strength
()
const
noexcept
{
return
m_second
.
ionic_strength
;
}
//! \brief Return the total concentration of 'component'
//!
//! Only use this method if the mass is conserved for 'component'
//! \param component Index of the aqueous component (in the database)
scalar_t
total_concentration_bc
(
index_t
component
)
const
{
specmicp_assert
((
component
>
0
and
aqueous_component_equation_type
(
component
)
==
AqueousComponentEquationType
::
MassConservation
)
or
(
water_equation_type
()
==
WaterEquationType
::
MassConservation
));
return
m_fixed_values
(
component
);
}
//! \brief Return the fixed activity of 'component'
//!
//! Only use this method if 'component' has a fixed activity constraint
//! \param component Index of the aqueous component (in the database)
scalar_t
fixed_activity_bc
(
index_t
component
)
const
{
specmicp_assert
(
aqueous_component_equation_type
(
component
)
==
AqueousComponentEquationType
::
FixedActivity
);
return
m_fixed_values
(
component
);
}
//! \brief Return the fixed molality of 'component'
//!
//! Only use this method if 'component' has a fixed molality constraint
//! \param component Index of the aqueous component (in the database)
scalar_t
fixed_molality_bc
(
index_t
component
)
const
{
specmicp_assert
(
aqueous_component_equation_type
(
component
)
==
AqueousComponentEquationType
::
FixedMolality
);
return
m_fixed_values
(
component
);
}
//! \brief Return the fixed fugacity value for 'component'
scalar_t
fixed_fugacity_bc
(
index_t
component
)
const
{
specmicp_assert
(
aqueous_component_equation_type
(
component
)
==
AqueousComponentEquationType
::
FixedFugacity
);
return
m_fixed_values
(
component
);
}
//! \brief Return the total concentration for the electron
scalar_t
electron_total_concentration
()
const
{
return
0.0
;
}
// Gas
// ---
//! \brief Return the fugacity of 'gas'
//!
//! \param gas Index of the gas (in the database)
scalar_t
gas_fugacity
(
index_t
gas
)
const
{
return
m_second
.
gas_fugacity
(
gas
);
}
//! \brief Return the concentration of 'gas' in the system
//!
//! \param gas Index of the gas (in the database)
scalar_t
gas_concentration
(
index_t
gas
)
const
{
return
m_second
.
gas_concentration
(
gas
);
}
//! \brief Return true if gas is active
//!
//! \param gas Index of the gas (in the database)
bool
is_active_gas
(
index_t
gas
)
const
{
return
m_equations
.
active_gas
[
gas
];
}
//! \brief Return the volume fraction (total saturation) of the gas phase
scalar_t
volume_fraction_gas_phase
()
const
noexcept
{
return
m_second
.
volume_fraction_gas
;
}
// Sorbed species
// --------------
//! \brief Return the surface total concentration
const
scalar_t
&
surface_total_concentration
()
const
{
return
m_fixed_values
(
dof_surface
());
}
//! \brief Return true if 'sorbed' is an active species
//! \param sorbed Index of the sorbed species (in the database)
bool
is_active_sorbed
(
index_t
sorbed
)
const
{
return
m_equations
.
active_sorbed
[
sorbed
];
}
//! \brief Return the molality of the sorbed species 'sorbed'
//! \param sorbed Index of the sorbed species (in the database)
const
scalar_t
&
sorbed_species_concentration
(
index_t
sorbed
)
const
{
return
m_second
.
sorbed_concentrations
(
sorbed
);
}
// Pressure and temperature
// ------------------------
//! \brief Return the gas pressure
//!
//! This is a fixed value (1 atm)
scalar_t
gas_total_pressure
()
const
{
return
1.01325e5
;
// Note : all pressure are in pascal, unit conversion is done for the concentration
}
//! \brief Return the temperature
//!
//! This is a fixed value (25°C)
scalar_t
temperature
()
const
noexcept
{
return
273.16
+
25
;
}
//! @}
// Solution
// =================
//! \brief Return the equilibrium state of the system, the Solution of the speciation problem
//!
//! \param xtot the complete set of main variables
//! \param x the reduced set of main variables, only the variables with an equation
AdimensionalSystemSolution
get_solution
(
Vector
&
xtot
,
const
Vector
&
x
);
//! \name Secondary variables computation
//! \brief Algorithms to compute the secondary variables
//!
//! @{
//! \brief Compute the ionic strength
//!
//! \param x Vector of the main variables
void
set_ionic_strength
(
const
Vector
&
x
);
//! \brief Compute the activity coefficients
//!
//! \param x Vector of the main variables
void
compute_log_gamma
(
const
Vector
&
x
);
//! \brief Compute the secondary aqueous species molalities
//!
//! \param x Vector of the main variables
void
set_secondary_concentration
(
const
Vector
&
x
);
//! \brief Compute the secondary variables
//!
//! \param x Vector of the main variables
void
set_secondary_variables
(
const
Vector
&
x
);
//! \brief Compute the gas phase volume fraction
//!
//! \param x Vector of the main variables
void
set_volume_fraction_gas_phase
(
const
Vector
&
x
);
//! \brief Compute the fugacity for all the gas
//!
//! \param x Vector of the main variables
void
set_pressure_fugacity
(
const
Vector
&
x
);
//! \brief Compute the sorbed species concentrations
//!
//! \param x Vector of the main variables
void
set_sorbed_concentrations
(
const
Vector
&
x
);
//! \brief This function is called at the beginning of each iteration by the solver
//!
//! It does the fixed-point iteration to compute the secondary variables
//! \param x Vector of the main variables
//! \param norm_residual Norm of the current residuals
bool
hook_start_iteration
(
const
Vector
&
x
,
scalar_t
norm_residual
);
//! @}
//! \brief A reasonable (well... maybe...) starting guess
//!
//! \param xtot Vector of the main variables
//! \sa #reasonable_starting_guess
void
reasonable_starting_guess
(
Vector
&
xtot
);
//! \brief A reasonable (maybe...) restarting guess
//!
//! \param xtot Vector of the main variables
//! \sa #reasonable_restarting_guess
void
reasonable_restarting_guess
(
Vector
&
xtot
);
//! \brief Set the units
void
set_units
(
const
units
::
UnitsSet
&
unit_set
);
// Private Members and attributes
// ##############################
private
:
// scaling factors
void
compute_scaling_factors
();
//! \brief Sum of the aqueous molalities weighted by the stoichiometric coefficient of 'component'
scalar_t
weigthed_sum_aqueous
(
index_t
component
)
const
;
//! \brief Sum of the sorbed concentration weighted by the stoichiometric coefficient of 'component'
scalar_t
weigthed_sum_sorbed
(
index_t
component
)
const
;
//! \brief Sum of the mineral concentration weighted by the stoichiometric coefficient of 'component'
scalar_t
weigthed_sum_mineral
(
const
Vector
&
x
,
index_t
component
)
const
;
//! \brief Sum of the gas concentration weigthed by the stoichiometric coefficient of 'component'
scalar_t
weigthed_sum_gas
(
index_t
component
)
const
;
//! \brief Derivative of the aqueous weigthed sum with respect to 'diff_component'
//! \sa weigthed_sum_aqueous
scalar_t
diff_weigthed_sum_aqueous
(
index_t
diff_component
,
index_t
component
)
const
;
//! \brief Derivative of the sorbed weigthed sum with respect to 'diff_component'
//! \sa weigthed_sum_sorbed
scalar_t
diff_weigthed_sum_sorbed
(
index_t
diff_component
,
index_t
component
)
const
;
//! \brief Derivative of the sorbed weigthed sum with respect to the free surface concentration
//! \sa weigthed_sum_sorbed
scalar_t
diff_surface_weigthed_sum_sorbed
(
index_t
component
)
const
;
//! \brief Derivative of the gas weigthed sum with respect to 'diff_component'
//! \sa weigthed_sum_sorbed
scalar_t
diff_weigthed_sum_gas
(
index_t
diff_component
,
index_t
component
)
const
;
// Jacobian
// ########
//! \brief Compute the water equation contribution to the Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void
jacobian_water
(
Vector
&
x
,
Matrix
&
jacobian
);
//! \brief Compute the aqueous components equation contribution to the Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void
jacobian_aqueous_components
(
Vector
&
x
,
Matrix
&
jacobian
);
//! \brief Compute the mineral equations contribution to the Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void
jacobian_minerals
(
Vector
&
x
,
Matrix
&
jacobian
);
//! \brief Compute the contribution of the surface sorption equation to te Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void
jacobian_surface
(
Vector
&
x
,
Matrix
&
jacobian
);
//! \brief Compute the contribution of the electron equation to te Jacobian
//! \param x Reduced vector of the main variables
//! \param jacobian Dense matrix containing the jacobian
void
jacobian_electron
(
Vector
&
x
,
Matrix
&
jacobian
);
// Equation numbering
// ##################
//! \brief Number the equations
//! \param constraints the constraints to apply to the system
//! \sa number_eq_aqueous_component
void
number_eq
(
const
AdimensionalSystemConstraints
&
constraints
);
//! \brief Number the equations for the aqueous components
//! \param constraints the constraints to apply to the system
//! \param[in,out] neq the number of equations in the system
//! \sa number_eq
void
number_eq_aqueous_component
(
const
AdimensionalSystemConstraints
&
constraints
,
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_second
.
loggamma
(
dof_component_gamma
(
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_second
.
loggamma
(
dof_aqueous_gamma
(
aqueous
));
}
//! \brief Return a reference to the ionic strength
scalar_t
&
ionic_strength
()
{
return
m_second
.
ionic_strength
;
}
// Attributes
// ##########
//! \struct SecondaryVariables
//! \brief Contains information about the secondary variables
struct
SecondaryVariables
{
scalar_t
ionic_strength
{
0.0
};
//!< The ionic Strength
scalar_t
volume_fraction_gas
{
0.0
};
//!< The gas saturation
scalar_t
porosity
{
0.0
};
//!< The porosity
Vector
secondary_molalities
;
//!< The secondary molalities
Vector
loggamma
;
//!< The log of activity coefficients
Vector
gas_fugacity
;
//!< The gas fugacities
Vector
gas_concentration
;
//!< The gas concentrations
Vector
sorbed_concentrations
;
//!< The sorbed concentrations
//! \brief Initialization without a previous solution
SecondaryVariables
(
const
RawDatabasePtr
&
data
);
//! \brief Initialization with a previous solution
SecondaryVariables
(
const
AdimensionalSystemSolution
&
previous_solution
);
};
//! \struct IdEquations
//! \brief BookKeeper for the equations id and type
struct
IdEquations
{
bool
use_water_pressure_model
{
false
};
water_partial_pressure_f
water_pressure_model
{
nullptr
};
index_t
nb_tot_variables
{
0
};
index_t
nb_free_variables
{
0
};
index_t
nb_complementarity_variables
{
0
};
std
::
vector
<
index_t
>
ideq
;
std
::
vector
<
index_t
>
component_equation_type
;
std
::
vector
<
index_t
>
fixed_activity_species
;
std
::
vector
<
index_t
>
nonactive_component
{};
std
::
vector
<
bool
>
active_aqueous
;
std
::
vector
<
bool
>
active_gas
;
std
::
vector
<
bool
>
active_sorbed
;
//! \brief Initialize the data structure
IdEquations
(
index_t
nb_dofs
,
const
RawDatabasePtr
&
data
);
//! \brief Return the equation id
const
index_t
&
id_equation
(
index_t
dof
)
const
{
return
ideq
[
dof
];
}
//! \brief Return a reference to the equation id
index_t
&
id_equation
(
index_t
dof
)
{
return
ideq
[
dof
];
}
//! \brief Return a reference to the type of equation of 'dof'
index_t
&
type_equation
(
index_t
dof
)
{
return
component_equation_type
[
dof
];
}
//! \brief Return the equation type for 'dof'
index_t
type_equation
(
index_t
dof
)
const
{
return
component_equation_type
[
dof
];
}
//! \brief Return the related species for dof
//!
//! The exact meaning of this relation is dependant upon the equation type
index_t
&
related_species
(
index_t
dof
)
{
return
fixed_activity_species
[
dof
];
}
//! \brief Add a non active component to the system
void
add_non_active_component
(
index_t
id
)
{
nonactive_component
.
push_back
(
id
);
}
//! \brief Add an equation
void
add_equation
(
index_t
id
,
index_t
*
neq
)
{
ideq
[
id
]
=
*
neq
;
++
(
*
neq
);
}
//! \brief Set the active flag of aqueous species 'id' to 'is_active'
void
set_aqueous_active
(
index_t
id
,
bool
is_active
)
{
active_aqueous
[
id
]
=
is_active
;
}
//! \brief Set the active flag of gas 'id' to 'is_active'
void
set_gas_active
(
index_t
id
,
bool
is_active
)
{
active_gas
[
id
]
=
is_active
;
}
//! \brief Set the active flag of sorbed species 'id' to 'is_active'
void
set_sorbed_active
(
index_t
id
,
bool
is_active
)
{
active_sorbed
[
id
]
=
is_active
;
}
};
bool
not_in_linesearch
{
true
};
scalar_t
m_inert_volume_fraction
;
scalar_t
m_scaling_molar_volume
{
1.0
};
scalar_t
m_scaling_gas
{
1.0
};
Vector
m_fixed_values
{};
SecondaryVariables
m_second
;
IdEquations
m_equations
;
};
}
// end namespace specmicp
#endif
// SPECMICP_ADIMENSIONALSYSTEM_HPP
Event Timeline
Log In to Comment