Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F60142139
equilibrium_stagger.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
Sat, Apr 27, 19:51
Size
4 KB
Mime Type
text/x-c
Expires
Mon, Apr 29, 19:51 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17307768
Attached To
rSPECMICP SpecMiCP / ReactMiCP
equilibrium_stagger.cpp
View Options
#include "equilibrium_stagger.hpp"
#include "variables.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "utils/log.hpp"
#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
#ifdef SPECMICP_USE_OPENMP
#include <omp.h>
#endif
// SPECMICP_USE_OPENMP
namespace
specmicp
{
namespace
reactmicp
{
namespace
systems
{
namespace
satdiff
{
// EquilibriumStagger::
//! \brief Initialize the stagger at the beginning of an iteration
void
EquilibriumStagger
::
initialize_timestep
(
scalar_t
dt
,
VariablesBasePtr
var
)
{
m_dt
=
dt
;
SaturatedVariablesPtr
true_var
=
cast_var_from_base
(
var
);
// Initialize velocity using values from previous timestep
for
(
index_t
node
=
0
;
node
<
true_var
->
nb_nodes
();
++
node
)
{
if
(
true_var
->
is_fixed_composition
(
node
))
continue
;
scalar_t
alpha
=
1.0
;
for
(
index_t
component
;
component
<
true_var
->
nb_component
();
++
component
)
{
alpha
=
std
::
max
(
alpha
,
dt
*
true_var
->
solid_concentration
(
node
,
component
,
true_var
->
chemistry_rate
())
/
true_var
->
solid_concentration
(
node
,
component
,
true_var
->
displacement
())
);
}
auto
solid_velocity
=
true_var
->
velocity
().
segment
(
true_var
->
offset_node
(
node
)
+
true_var
->
offset_solid_concentration
(),
true_var
->
nb_component
());
auto
solid_chemistry_rate
=
true_var
->
chemistry_rate
().
segment
(
true_var
->
offset_node
(
node
)
+
true_var
->
offset_solid_concentration
(),
true_var
->
nb_component
());
solid_velocity
=
1
/
alpha
*
solid_chemistry_rate
;
}
}
//! \brief Solve the equation for the timestep
solver
::
StaggerReturnCode
EquilibriumStagger
::
restart_timestep
(
VariablesBasePtr
var
)
{
SaturatedVariablesPtr
true_var
=
cast_var_from_base
(
var
);
#ifdef SPECMICP_USE_OPENMP
#pragma omp parallel
{
#pragma omp for schedule(runtime)
#endif
// SPECMICP_USE_OPENMP
for
(
index_t
node
=
0
;
node
<
true_var
->
nb_component
();
++
node
)
{
if
(
true_var
->
is_fixed_composition
(
false
))
continue
;
solve_one_node
(
node
,
true_var
);
}
#ifdef SPECMICP_USE_OPENMP
}
#endif
// SPECMICP_USE_OPENMP
return
solver
::
StaggerReturnCode
::
NotConvergedYet
;
}
//!
void
EquilibriumStagger
::
solve_one_node
(
index_t
node
,
SaturatedVariablesPtr
true_var
)
{
AdimensionalSystemConstraints
constraints
(
m_constraints
);
constraints
.
total_concentrations
=
true_var
->
total_concentrations
(
node
);
AdimensionalSystemSolver
adim_solver
(
true_var
->
get_database
(),
constraints
,
true_var
->
equilibrium_solution
(
node
),
m_options
);
Vector
variables
(
true_var
->
equilibrium_solution
(
node
).
main_variables
);
micpsolver
::
MiCPPerformance
perf
=
adim_solver
.
solve
(
variables
);
micpsolver
::
MiCPSolverReturnCode
retcode
=
perf
.
return_code
;
if
(
retcode
<=
micpsolver
::
MiCPSolverReturnCode
::
NotConvergedYet
)
{
ERROR
<<
"Failed to solve chemistry problem at node "
<<
node
<<
", return code = "
<<
static_cast
<
int
>
(
retcode
);
}
true_var
->
equilibrium_solution
(
node
)
=
adim_solver
.
get_raw_solution
(
variables
);
AdimensionalSystemSolutionExtractor
extractor
(
true_var
->
equilibrium_solution
(
node
),
true_var
->
get_database
(),
m_options
.
units_set
);
for
(
index_t
component
=
0
;
component
<
true_var
->
nb_component
();
++
component
)
{
const
scalar_t
c_aq
=
extractor
.
density_water
()
*
extractor
.
total_aqueous_concentration
(
component
);
true_var
->
aqueous_concentration
(
node
,
component
,
true_var
->
displacement
())
=
c_aq
;
const
scalar_t
c_aq_0
=
true_var
->
aqueous_concentration
(
node
,
component
,
true_var
->
predictor
());
const
scalar_t
vel_aq
=
(
c_aq
-
c_aq_0
)
/
m_dt
;
true_var
->
aqueous_concentration
(
node
,
component
,
true_var
->
velocity
())
=
vel_aq
;
const
scalar_t
c_sol
=
extractor
.
total_solid_concentration
(
component
);
true_var
->
solid_concentration
(
node
,
component
,
true_var
->
displacement
())
=
c_sol
;
const
scalar_t
c_sol_0
=
true_var
->
solid_concentration
(
node
,
component
,
true_var
->
predictor
());
const
scalar_t
vel_sol
=
(
c_sol
-
c_sol_0
)
/
m_dt
;
true_var
->
solid_concentration
(
node
,
component
,
true_var
->
velocity
())
=
vel_sol
;
true_var
->
solid_concentration
(
node
,
component
,
true_var
->
chemistry_rate
())
=
vel_sol
;
}
}
}
// end namespace satdiff
}
// end namespace systems
}
// end namespace reactmicp
}
// end namespace specmicp
Event Timeline
Log In to Comment