Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F84426367
reduced_system.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
Sun, Sep 22, 19:46
Size
19 KB
Mime Type
text/x-c
Expires
Tue, Sep 24, 19:46 (2 d)
Engine
blob
Format
Raw Data
Handle
20895555
Attached To
rSPECMICP SpecMiCP / ReactMiCP
reduced_system.cpp
View Options
/*-------------------------------------------------------
- Module : specmicp
- File : specmicp.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 "reduced_system.hpp"
#include "utils/log.hpp"
#include "equilibrium_data.hpp"
#include "physics/constants.hpp"
#include "physics/laws.hpp"
namespace
specmicp
{
constexpr
scalar_t
log10
=
std
::
log
(
10.0
);
inline
double
pow10
(
double
x
)
{
return
std
::
pow
(
10.0
,
x
);
}
ReducedSystem
::
ReducedSystem
(
RawDatabasePtr
ptrdata
,
const
Vector
&
totconc
,
const
ReducedSystemOptions
&
options
)
:
OptionsHandler
<
ReducedSystemOptions
>
(
options
),
m_data
(
ptrdata
),
m_tot_conc
(
totconc
),
m_secondary_conc
(
ptrdata
->
nb_aqueous
),
m_loggamma
(
ptrdata
->
nb_component
+
ptrdata
->
nb_aqueous
)
{
if
(
not
m_simulbox
.
is_fixed_temperature
()
or
m_simulbox
.
get_temperature
()
!=
units
::
celsius
(
25
))
{
throw
std
::
runtime_error
(
"Only a fixed temperature of 25°C is supported"
);
}
number_eq
();
m_secondary_conc
.
setZero
();
m_loggamma
.
setZero
();
}
ReducedSystem
::
ReducedSystem
(
RawDatabasePtr
ptrdata
,
const
Vector
&
totconc
,
const
SimulationBox
&
simul_box
,
const
ReducedSystemOptions
&
options
)
:
OptionsHandler
<
ReducedSystemOptions
>
(
options
),
m_simulbox
(
simul_box
),
m_data
(
ptrdata
),
m_tot_conc
(
totconc
),
m_secondary_conc
(
ptrdata
->
nb_aqueous
),
m_loggamma
(
ptrdata
->
nb_component
+
ptrdata
->
nb_aqueous
)
{
if
(
not
m_simulbox
.
is_fixed_temperature
()
or
m_simulbox
.
get_temperature
()
!=
units
::
celsius
(
25
))
{
throw
std
::
runtime_error
(
"Only a fixed temperature of 25°C is supported"
);
}
number_eq
();
m_secondary_conc
.
setZero
();
m_loggamma
.
setZero
();
}
ReducedSystem
::
ReducedSystem
(
RawDatabasePtr
ptrdata
,
const
Vector
&
totconc
,
const
SimulationBox
&
simul_box
,
const
EquilibriumState
&
previous_solution
,
const
ReducedSystemOptions
&
options
)
:
OptionsHandler
<
ReducedSystemOptions
>
(
options
),
m_simulbox
(
simul_box
),
m_data
(
ptrdata
),
m_tot_conc
(
totconc
),
m_secondary_conc
(
previous_solution
.
molality_secondary
()),
m_loggamma
(
previous_solution
.
activity_coefficient
())
{
if
(
not
m_simulbox
.
is_fixed_temperature
()
or
m_simulbox
.
get_temperature
()
!=
units
::
celsius
(
25
))
{
throw
std
::
runtime_error
(
"Only a fixed temperature of 25°C is supported"
);
}
number_eq
();
}
void
ReducedSystem
::
number_eq
()
{
index_t
neq
=
0
;
m_ideq
.
reserve
(
m_data
->
nb_component
+
m_data
->
nb_mineral
);
if
(
get_options
().
conservation_water
)
{
m_ideq
.
push_back
(
neq
);
++
neq
;
}
else
{
m_ideq
.
push_back
(
no_equation
);
}
for
(
index_t
i:
m_data
->
range_aqueous_component
())
{
// Remove components that does not exist
if
(
m_tot_conc
(
i
)
==
0
and
i
!=
1
)
{
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
);
continue
;
}
else
{
m_ideq
.
push_back
(
neq
);
++
neq
;
}
}
m_nb_free_vars
=
neq
;
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
);
++
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
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
);
}
}
scalar_t
ReducedSystem
::
residual_water
(
const
Vector
&
x
)
const
{
scalar_t
res
=
m_tot_conc
(
0
)
-
mass_water
(
x
)
/
m_data
->
molar_mass_basis_si
(
0
);
for
(
index_t
j:
m_data
->
range_aqueous
())
{
if
(
m_data
->
nu_aqueous
(
j
,
0
)
==
0
)
continue
;
res
-=
mass_water
(
x
)
*
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
)
continue
;
res
-=
m_data
->
nu_mineral
(
m
,
0
)
*
x
(
ideq_min
(
m
));
}
return
res
;
}
double
ReducedSystem
::
residual_component
(
const
Vector
&
x
,
index_t
i
)
const
{
const
index_t
idp
=
ideq_paq
(
i
);
scalar_t
res
=
m_tot_conc
(
i
)
-
mass_water
(
x
)
*
pow10
(
x
(
idp
));
for
(
index_t
j:
m_data
->
range_aqueous
())
{
if
(
not
m_active_aqueous
[
j
])
continue
;
res
-=
mass_water
(
x
)
*
m_data
->
nu_aqueous
(
j
,
i
)
*
m_secondary_conc
(
j
);
}
for
(
index_t
m:
m_data
->
range_mineral
())
{
if
(
ideq_min
(
m
)
==
no_equation
)
continue
;
res
-=
m_data
->
nu_mineral
(
m
,
i
)
*
x
(
ideq_min
(
m
));
}
return
res
;
}
double
ReducedSystem
::
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
)
*
(
x
(
ideq_paq
(
i
))
+
m_loggamma
(
i
));
}
return
res
;
}
double
ReducedSystem
::
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_concentration
(
x
,
i
);
}
for
(
index_t
j:
m_data
->
range_aqueous
())
{
if
(
m_data
->
charge_aqueous
(
j
)
==
0
and
not
m_active_aqueous
[
j
])
continue
;
res
+=
m_data
->
charge_aqueous
(
j
)
*
m_secondary_conc
(
j
);
}
return
res
;
}
void
ReducedSystem
::
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
())
{
if
(
ideq_paq
(
i
)
==
no_equation
)
continue
;
if
(
is_charge_keeper
(
i
))
residual
(
ideq_paq
(
i
))
=
residual_charge_conservation
(
x
);
else
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
);
}
}
//non-optimized Finite difference, for test only
//void ReducedSystem::get_jacobian(Eigen::VectorXd& x,
// Eigen::MatrixXd& 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;
//}
void
ReducedSystem
::
jacobian_water
(
Vector
&
x
,
Matrix
&
jacobian
)
{
if
(
ideq_w
()
!=
no_equation
)
{
scalar_t
tmp
=
-
1
/
m_data
->
molar_mass_basis_si
(
0
);
for
(
index_t
j:
m_data
->
range_aqueous
())
{
if
(
m_data
->
nu_aqueous
(
j
,
0
)
==
0
)
continue
;
tmp
-=
m_data
->
nu_aqueous
(
j
,
0
)
*
m_secondary_conc
(
j
);
}
jacobian
(
0
,
0
)
=
tmp
;
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
)
*
m_secondary_conc
(
j
);
}
jacobian
(
0
,
ideq_paq
(
k
))
=
x
(
ideq_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
);
}
}
}
void
ReducedSystem
::
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
;
if
(
is_charge_keeper
(
i
))
{
// 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
m_active_aqueous
[
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
*=
m_secondary_conc
(
j
)
/
component_concentration
(
x
,
k
);
tmp_drdb
+=
tmp_value
;
}
jacobian
(
idp
,
idc
)
+=
component_concentration
(
x
,
k
)
*
log10
*
tmp_drdb
;
}
}
else
{
// 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
-=
pow10
(
x
(
ideq_paq
(
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
)
*
m_secondary_conc
(
j
);
}
jacobian
(
idp
,
ideq_paq
(
k
))
=
mass_water
(
x
)
*
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
);
}
if
(
ideq_w
()
!=
no_equation
)
// Water
{
scalar_t
tmp_iw
=
-
pow10
(
x
(
idp
));
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
)
*
m_secondary_conc
(
j
);
}
jacobian
(
idp
,
ideq_w
())
=
tmp_iw
;
}
}
}
}
void
ReducedSystem
::
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
ReducedSystem
::
get_jacobian
(
Vector
&
x
,
Matrix
&
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
);
}
void
ReducedSystem
::
set_secondary_concentration
(
const
Vector
&
x
)
{
for
(
index_t
j:
m_data
->
range_aqueous
())
{
if
(
m_active_aqueous
[
j
]
==
false
)
{
m_secondary_conc
(
j
)
=
0
;
continue
;
}
scalar_t
conc
=
-
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
;
conc
+=
m_data
->
nu_aqueous
(
j
,
k
)
*
(
x
(
ideq_paq
(
k
))
+
m_loggamma
(
k
));
}
m_secondary_conc
(
j
)
=
pow10
(
conc
);
}
}
void
ReducedSystem
::
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
+=
pow10
(
x
(
ideq_paq
(
i
)))
*
std
::
pow
(
m_data
->
charge_component
(
i
),
2
);
}
for
(
index_t
j:
m_data
->
range_aqueous
())
{
if
(
not
m_active_aqueous
[
j
]
or
m_data
->
charge_aqueous
(
j
)
==
0
)
continue
;
ionic
+=
m_secondary_conc
(
j
)
*
std
::
pow
(
m_data
->
charge_aqueous
(
j
),
2
);
}
m_ionic_strength
=
ionic
/
2
;
}
void
ReducedSystem
::
compute_log_gamma
(
const
Vector
&
x
)
{
set_ionic_strength
(
x
);
if
(
get_options
().
non_ideality
)
{
const
scalar_t
sqrti
=
std
::
sqrt
(
m_ionic_strength
);
for
(
index_t
i:
m_data
->
range_aqueous_component
())
{
if
(
ideq_paq
(
i
)
==
no_equation
)
{
m_loggamma
(
i
)
=
0
;
continue
;
}
m_loggamma
(
i
)
=
laws
::
extended_debye_huckel
(
m_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
m_active_aqueous
[
j
])
{
m_loggamma
(
m_data
->
nb_component
+
j
)
=
0
;
continue
;
}
m_loggamma
(
m_data
->
nb_component
+
j
)
=
laws
::
extended_debye_huckel
(
m_ionic_strength
,
sqrti
,
m_data
->
charge_aqueous
(
j
),
m_data
->
a_debye_aqueous
(
j
),
m_data
->
b_debye_aqueous
(
j
)
);
}
}
}
bool
ReducedSystem
::
hook_start_iteration
(
const
Vector
&
x
,
scalar_t
norm_residual
)
{
// if we consider ideal solution, we break
if
(
not
get_options
().
non_ideality
)
return
true
;
// else we solve for non ideality by fixed-point iterations
//
// initialisation
not_in_linesearch
=
true
;
scalar_t
previous_norm
=
m_loggamma
.
norm
();
bool
may_have_converged
=
false
;
// too much iteration is not a failure
// solve
if
(
norm_residual
<
nb_free_variables
())
{
// Use fix point iterations for non-ideality
for
(
int
i
=
0
;
i
<
get_options
().
non_ideality_max_iter
;
++
i
)
{
set_secondary_concentration
(
x
);
compute_log_gamma
(
x
);
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
ReducedSystem
::
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
ReducedSystem
::
reasonable_starting_guess
(
Vector
&
x
,
bool
for_min
)
{
x
.
resize
(
m_data
->
nb_component
+
m_data
->
nb_mineral
);
x
(
0
)
=
2
*
m_tot_conc
(
0
)
*
m_data
->
molar_mass_basis_si
(
0
);
for
(
index_t
i:
m_data
->
range_aqueous_component
())
{
x
(
i
)
=
-
4.0
;
}
if
(
for_min
)
x
.
block
(
m_data
->
nb_component
,
0
,
m_data
->
nb_mineral
,
1
).
setZero
();
else
x
.
block
(
m_data
->
nb_component
,
0
,
m_data
->
nb_mineral
,
1
).
setZero
();
m_loggamma
.
setZero
();
m_secondary_conc
.
setZero
();
}
void
ReducedSystem
::
reasonable_restarting_guess
(
Vector
&
x
)
{
x
(
0
)
=
m_tot_conc
(
0
)
*
m_data
->
molar_mass_basis_si
(
0
);
for
(
index_t
i
=
m_data
->
nb_component
;
i
<
m_data
->
nb_component
+
m_data
->
nb_mineral
;
++
i
)
{
x
(
i
)
=
0.0
;
}
m_loggamma
.
setZero
();
m_secondary_conc
.
setZero
();
}
EquilibriumState
ReducedSystem
::
get_solution
(
const
Vector
&
xtot
,
const
Vector
&
x
)
{
double
previous_norm
=
m_loggamma
.
norm
();
set_secondary_concentration
(
x
);
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
EquilibriumState
(
xtot
,
m_secondary_conc
,
m_loggamma
,
m_ionic_strength
,
m_data
);
}
void
ReducedSystem
::
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
ReducedSystem
::
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
Event Timeline
Log In to Comment