Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65334300
reaction_path.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, Jun 3, 01:17
Size
2 KB
Mime Type
text/x-c
Expires
Wed, Jun 5, 01:17 (2 d)
Engine
blob
Format
Raw Data
Handle
18040123
Attached To
rSPECMICP SpecMiCP / ReactMiCP
reaction_path.cpp
View Options
/*-------------------------------------------------------
- Module : specmicp
- File : reaction_path.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "reaction_path.hpp"
#include "thermodata.hpp"
#include "micpsolver/micpsolver.hpp"
#include "utils/log.hpp"
#include <iostream>
namespace
specmicp
{
void
ReactionPathDriver
::
read_database
()
{
m_database
.
parse_database
(
m_model
->
database_path
);
m_database
.
make_canonical
();
m_data
=
m_database
.
get_database
();
}
int
ReactionPathDriver
::
label_to_id
(
const
std
::
string
&
label
,
const
std
::
vector
<
std
::
string
>&
list_labels
)
{
auto
it
=
std
::
find
(
list_labels
.
begin
(),
list_labels
.
end
(),
label
);
if
(
it
==
list_labels
.
end
())
{
throw
std
::
invalid_argument
(
"Unknown species : "
+
label
);
}
return
it
-
list_labels
.
begin
();
}
void
ReactionPathDriver
::
dissolve_to_components
()
{
m_tot_conc
=
Eigen
::
VectorXd
::
Zero
(
m_data
->
nb_component
);
m_tot_conc_increase
=
Eigen
::
VectorXd
::
Zero
(
m_data
->
nb_component
);
// everything is dissolved ...
for
(
auto
it
=
m_model
->
amount_components
.
begin
();
it
!=
m_model
->
amount_components
.
end
();
++
it
)
{
const
int
id
=
label_to_id
(
it
->
first
,
m_data
->
labels_basis
);
m_tot_conc
(
id
)
=
it
->
second
.
first
;
m_tot_conc_increase
(
id
)
=
it
->
second
.
second
;
}
// Fixme : add contribution for minerals and aqueous species
std
::
vector
<
int
>
to_remove
;
to_remove
.
reserve
(
10
);
int
new_i
=
0
;
for
(
int
i
=
0
;
i
<
m_data
->
nb_component
;
++
i
)
{
if
(
m_tot_conc
(
i
)
==
0
and
m_tot_conc_increase
(
i
)
==
0
)
{
to_remove
.
push_back
(
i
);
continue
;
}
m_tot_conc
(
new_i
)
=
m_tot_conc
(
i
);
m_tot_conc_increase
(
new_i
)
=
m_tot_conc_increase
(
i
);
++
new_i
;
}
m_tot_conc
.
conservativeResize
(
new_i
);
m_tot_conc_increase
.
conservativeResize
(
new_i
);
if
(
to_remove
.
size
()
>
0
)
m_database
.
remove_components
(
to_remove
);
}
micpsolver
::
MiCPPerformance
ReactionPathDriver
::
one_step
(
Eigen
::
VectorXd
&
x
,
bool
init
)
{
m_current_solver
=
ReducedSystemSolver
(
m_data
,
m_tot_conc
);
m_current_solver
.
get_options
()
=
m_solver_options
;
micpsolver
::
MiCPPerformance
perf
=
m_current_solver
.
solve
(
x
,
init
);
if
(
perf
.
return_code
!=
micpsolver
::
MiCPSolverReturnCode
::
ResidualMinimized
)
{
ERROR
<<
"Failed to solve the system !"
;
}
m_tot_conc
+=
m_tot_conc_increase
;
return
perf
;
}
void
ReactionPathDriver
::
run_all
()
{
Eigen
::
VectorXd
x
;
// first step
one_step
(
x
,
true
);
for
(
int
i
=
1
;
i
<
m_model
->
nb_step
;
++
i
)
{
one_step
(
x
);
save_data
(
x
);
}
}
void
ReactionPathDriver
::
save_data
(
Eigen
::
VectorXd
&
x
)
{
}
}
// end namespace specmicp
Event Timeline
Log In to Comment