Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F70915345
adimensional_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
Mon, Jul 8, 06:07
Size
19 KB
Mime Type
text/x-c
Expires
Wed, Jul 10, 06:07 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18885895
Attached To
rSPECMICP SpecMiCP / ReactMiCP
adimensional_system.cpp
View Options
/*-------------------------------------------------------
- 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"
namespace
specmicp
{
constexpr
scalar_t
log10
=
std
::
log
(
10.0
);
AdimensionalSystem
::
AdimensionalSystem
(
RawDatabasePtr
ptrdata
,
const
Vector
&
totconc
,
bool
conservation_water
,
index_t
charge_keeper
)
:
m_data
(
ptrdata
),
m_tot_conc
(
totconc
),
m_secondary_conc
(
ptrdata
->
nb_aqueous
),
m_loggamma
(
ptrdata
->
nb_component
+
ptrdata
->
nb_aqueous
),
m_charge_keeper
(
charge_keeper
)
{
number_eq
(
conservation_water
);
m_secondary_conc
.
setZero
();
m_loggamma
.
setZero
();
}
AdimensionalSystem
::
AdimensionalSystem
(
RawDatabasePtr
ptrdata
,
const
Vector
&
totconc
,
const
AdimensionalSystemSolution
&
previous_solution
,
bool
conservation_water
,
index_t
charge_keeper
)
:
m_data
(
ptrdata
),
m_tot_conc
(
totconc
),
m_secondary_conc
(
previous_solution
.
secondary_molalities
),
m_loggamma
(
previous_solution
.
log_gamma
),
m_ionic_strength
(
previous_solution
.
ionic_strength
),
m_charge_keeper
(
charge_keeper
)
{
number_eq
(
conservation_water
);
}
void
AdimensionalSystem
::
number_eq
(
bool
water_equation
)
{
index_t
neq
=
0
;
m_ideq
.
reserve
(
m_data
->
nb_component
+
m_data
->
nb_mineral
);
if
(
water_equation
)
{
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
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
(
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
);
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
);
}
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
);
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
);
}
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
())
{
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 AdimensionalSystem::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
AdimensionalSystem
::
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
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
// ----------------------
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
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
;
}
}
// Mass balance equation
// ---------------------
else
{
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_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
)
{
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
<
10
;
++
i
)
{
set_secondary_concentration
(
x
);
compute_log_gamma
(
x
);
if
(
std
::
abs
(
previous_norm
-
m_loggamma
.
norm
())
/
previous_norm
<
1e-10
)
{
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
)
/
(
0.9
*
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_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
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
Event Timeline
Log In to Comment