Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77578385
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
Thu, Aug 15, 07:08
Size
15 KB
Mime Type
text/x-c
Expires
Sat, Aug 17, 07:08 (2 d)
Engine
blob
Format
Raw Data
Handle
19891064
Attached To
rSPECMICP SpecMiCP / ReactMiCP
equilibrium_stagger.cpp
View Options
/*------------------------------------------------------------------------------
Copyright (c)2015 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.
------------------------------------------------------------------------------*/
#include "equilibrium_stagger.hpp"
#include "variables.hpp"
#include "equilibrium_constraints.hpp"
#include "../../solver/staggers_base/stagger_structs.hpp"
#include "../../../dfpm/meshes/mesh1d.hpp"
#include "../../../database/data_container.hpp"
#include "../../../utils/compat.hpp"
#include "../../../specmicp/adimensional/adimensional_system_solution.hpp"
#include "../../../specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "../../../specmicp/adimensional/adimensional_system_solver.hpp"
#include "../../../utils/log.hpp"
namespace
specmicp
{
namespace
reactmicp
{
namespace
systems
{
namespace
unsaturated
{
// ============================================
//
// Declaration of implementation details
//
// ============================================
struct
EquilibriumStagger
::
EquilibriumStaggerImpl
{
scalar_t
m_dt
{
-
1
};
mesh
::
Mesh1DPtr
m_mesh
;
database
::
RawDatabasePtr
m_database
;
std
::
shared_ptr
<
EquilibriumConstraints
>
m_constraints
;
EquilibriumStaggerImpl
(
std
::
shared_ptr
<
UnsaturatedVariables
>
variables
,
std
::
shared_ptr
<
EquilibriumConstraints
>
constraints
)
:
m_mesh
(
variables
->
get_mesh
()),
m_database
(
variables
->
get_database
()),
m_constraints
(
constraints
)
{}
scalar_t
dt
()
{
return
m_dt
;}
int
compute_one_node
(
index_t
node
,
UnsaturatedVariables
*
vars
);
void
compute_total_concentrations
(
index_t
node
,
Vector
&
total_concentrations
,
UnsaturatedVariables
*
vars
);
void
analyse_solution
(
index_t
node
,
AdimensionalSystemSolution
&
solution
,
UnsaturatedVariables
*
vars
);
void
initialize_timestep_one_node
(
index_t
node
,
UnsaturatedVariables
*
vars
);
void
initialize
(
UnsaturatedVariables
*
vars
);
void
initialize_timestep
(
scalar_t
dt
,
UnsaturatedVariables
*
vars
);
StaggerReturnCode
restart_timestep
(
UnsaturatedVariables
*
vars
);
units
::
UnitsSet
&
get_units
()
{
return
m_constraints
->
get_options
().
units_set
;
}
};
// ============================================
//
// Implementation
//
// ============================================
EquilibriumStagger
::
EquilibriumStagger
(
std
::
shared_ptr
<
UnsaturatedVariables
>
variables
,
std
::
shared_ptr
<
EquilibriumConstraints
>
constraints
)
:
m_impl
(
make_unique
<
EquilibriumStaggerImpl
>
(
variables
,
constraints
))
{
}
EquilibriumStagger
::~
EquilibriumStagger
()
=
default
;
//! \brief Initialize the stagger at the beginning of the computation
//!
//! \param var a shared_ptr to the variables
void
EquilibriumStagger
::
initialize
(
VariablesBasePtr
var
)
{
UnsaturatedVariables
*
true_vars
=
static_cast
<
UnsaturatedVariables
*>
(
var
.
get
());
m_impl
->
initialize
(
true_vars
);
}
//! \brief Initialize the stagger at the beginning of an iteration
//!
//! This is where the predictor can be saved, the first trivial iteration done, ...
//!
//! \param dt the duration of the timestep
//! \param var a shared_ptr to the variables
void
EquilibriumStagger
::
initialize_timestep
(
scalar_t
dt
,
VariablesBasePtr
var
)
{
UnsaturatedVariables
*
true_vars
=
static_cast
<
UnsaturatedVariables
*>
(
var
.
get
());
m_impl
->
initialize_timestep
(
dt
,
true_vars
);
}
//! \brief Solve the equation for the timestep
//!
//! \param var a shared_ptr to the variables
EquilibriumStagger
::
StaggerReturnCode
EquilibriumStagger
::
restart_timestep
(
VariablesBasePtr
var
)
{
UnsaturatedVariables
*
true_vars
=
static_cast
<
UnsaturatedVariables
*>
(
var
.
get
());
return
m_impl
->
restart_timestep
(
true_vars
);
}
int
EquilibriumStagger
::
EquilibriumStaggerImpl
::
compute_one_node
(
index_t
node
,
UnsaturatedVariables
*
vars
)
{
AdimensionalSystemConstraints
constraints
=
m_constraints
->
get_constraints
(
node
);
compute_total_concentrations
(
node
,
constraints
.
total_concentrations
,
vars
);
{
user_model_saturation_f
callable
=
vars
->
get_vapor_pressure_model
();
water_partial_pressure_f
pv_model
=
[
callable
,
node
](
scalar_t
sat
){
return
callable
(
node
,
sat
);
};
constraints
.
set_water_partial_pressure_model
(
pv_model
);
}
// We use the current value of the saturation to avoid converging
// to the previous solution
// This is necessary because the update to the saturation is very small
AdimensionalSystemSolution
&
previous_solution
=
vars
->
get_adim_solution
(
node
);
if
(
previous_solution
.
is_valid
)
{
const
auto
porosity
=
vars
->
get_porosity
()(
node
);
const
auto
saturation
=
vars
->
get_liquid_saturation
()(
node
);
previous_solution
.
main_variables
(
0
)
=
porosity
*
saturation
;
}
//
AdimensionalSystemSolver
the_solver
(
vars
->
get_database
(),
constraints
,
previous_solution
,
m_constraints
->
get_options
()
);
Vector
x
=
previous_solution
.
main_variables
;
micpsolver
::
MiCPPerformance
perf
=
the_solver
.
solve
(
x
);
if
(
perf
.
return_code
<
micpsolver
::
MiCPSolverReturnCode
::
Success
)
{
ERROR
<<
"Failed to solve equilibrium at node "
<<
node
<<
".
\n
"
<<
"Return code : "
<<
(
int
)
perf
.
return_code
;
return
-
1
;
}
AdimensionalSystemSolution
solution
=
the_solver
.
get_raw_solution
(
x
);
analyse_solution
(
node
,
solution
,
vars
);
vars
->
set_adim_solution
(
node
,
solution
);
return
0
;
}
void
EquilibriumStagger
::
EquilibriumStaggerImpl
::
compute_total_concentrations
(
index_t
node
,
Vector
&
total_concentrations
,
UnsaturatedVariables
*
vars
)
{
total_concentrations
.
resize
(
m_database
->
nb_component
());
const
scalar_t
saturation
=
vars
->
get_liquid_saturation
()(
node
);
const
scalar_t
porosity
=
vars
->
get_porosity
()(
node
);
const
scalar_t
rt
=
vars
->
get_rt
();
// water
{
const
scalar_t
ctilde_w
=
vars
->
get_water_aqueous_concentration
()(
node
);
const
scalar_t
cbar_w
=
vars
->
get_solid_concentration
(
0
)(
node
);
scalar_t
c_w
=
saturation
*
porosity
*
ctilde_w
+
cbar_w
;
if
(
vars
->
component_has_gas
(
0
))
{
const
scalar_t
pv_w
=
vars
->
get_pressure_main_variables
(
0
)(
node
);
c_w
+=
(
1.0
-
saturation
)
*
porosity
*
pv_w
/
rt
;
}
total_concentrations
(
0
)
=
c_w
;
}
total_concentrations
(
1
)
=
0.0
;
// aqueous components
for
(
index_t
component:
m_database
->
range_aqueous_component
())
{
const
scalar_t
ctilde_i
=
vars
->
get_aqueous_concentration
(
component
)(
node
);
const
scalar_t
cbar_i
=
vars
->
get_solid_concentration
(
component
)(
node
);
scalar_t
c_i
=
saturation
*
porosity
*
ctilde_i
+
cbar_i
;
if
(
vars
->
component_has_gas
(
component
))
{
const
scalar_t
pv_i
=
vars
->
get_pressure_main_variables
(
component
)(
node
);
c_i
+=
(
1.0
-
saturation
)
*
porosity
*
pv_i
/
rt
;
}
total_concentrations
(
component
)
=
c_i
;
}
}
void
EquilibriumStagger
::
EquilibriumStaggerImpl
::
analyse_solution
(
index_t
node
,
AdimensionalSystemSolution
&
solution
,
UnsaturatedVariables
*
vars
)
{
AdimensionalSystemSolutionExtractor
extr
(
solution
,
m_database
,
get_units
());
const
scalar_t
saturation
=
extr
.
saturation_water
();
const
scalar_t
sat_g
=
1
-
saturation
;
const
scalar_t
rho_l
=
extr
.
density_water
();
const
scalar_t
porosity
=
extr
.
porosity
();
const
scalar_t
rt
=
vars
->
get_rt
();
// porosity
SecondaryTransientVariable
&
porosity_vars
=
vars
->
get_porosity
();
const
scalar_t
porosity_0
=
porosity_vars
.
predictor
(
node
);
const
scalar_t
porosity_vel
=
(
porosity
-
porosity_0
)
/
m_dt
;
porosity_vars
.
variable
(
node
)
=
porosity
;
porosity_vars
.
velocity
(
node
)
=
porosity_vel
;
// water
// saturation
MainVariable
&
satvars
=
vars
->
get_liquid_saturation
();
const
scalar_t
saturation_0
=
satvars
.
predictor
(
node
);
const
scalar_t
sat_vel
=
(
saturation
-
saturation_0
)
/
m_dt
;
satvars
.
variable
(
node
)
=
saturation
;
satvars
.
velocity
(
node
)
=
sat_vel
;
{
// total aqueous concentration
SecondaryTransientVariable
&
cwtilde_vars
=
vars
->
get_water_aqueous_concentration
();
const
scalar_t
cwtilde
=
rho_l
*
extr
.
total_aqueous_concentration
(
0
);
const
scalar_t
cwtilde_0
=
cwtilde_vars
.
predictor
(
node
);
const
scalar_t
cwtilde_vel
=
(
cwtilde
-
cwtilde_0
)
/
m_dt
;
cwtilde_vars
.
variable
(
node
)
=
cwtilde
;
cwtilde_vars
.
velocity
(
node
)
=
cwtilde_vel
;
// total solid concentration
MainVariable
&
solid_vars
=
vars
->
get_solid_concentration
(
0
);
const
scalar_t
cwbar
=
extr
.
total_solid_concentration
(
0
);
const
scalar_t
cwbar_0
=
solid_vars
.
predictor
(
node
);
const
scalar_t
cwbar_vel
=
(
cwbar
-
cwbar_0
)
/
m_dt
;
solid_vars
.
variable
(
node
)
=
cwbar
;
solid_vars
.
velocity
(
node
)
=
cwbar_vel
;
solid_vars
.
chemistry_rate
(
node
)
=
-
cwbar_vel
;
// vapor pressure
if
(
vars
->
component_has_gas
(
0
))
{
MainVariable
&
pres_vars
=
vars
->
get_pressure_main_variables
(
0
);
const
scalar_t
pv
=
vars
->
get_vapor_pressure_model
()(
node
,
saturation
);
const
scalar_t
pv_0
=
pres_vars
.
predictor
(
node
);
const
scalar_t
pv_vel
=
(
pv
-
pv_0
)
/
m_dt
;
pres_vars
.
variable
(
node
)
=
pv
;
pres_vars
.
velocity
(
node
)
=
pv_vel
;
const
scalar_t
transient
=
(
(
pv
*
porosity
*
sat_g
)
-
(
pv_0
*
porosity_0
*
(
1
-
saturation_0
))
)
/
(
rt
*
m_dt
);
const
scalar_t
pv_chem_rate
=
-
transient
+
pres_vars
.
transport_fluxes
(
node
);
pres_vars
.
chemistry_rate
(
node
)
=
pv_chem_rate
;
}
}
// aqueous components
for
(
index_t
component:
m_database
->
range_aqueous_component
())
{
// total aqueous concentration
MainVariable
&
aqconc
=
vars
->
get_aqueous_concentration
(
component
);
const
scalar_t
ctildei
=
rho_l
*
extr
.
total_aqueous_concentration
(
component
);
const
scalar_t
ctildei_0
=
aqconc
.
predictor
(
node
);
const
scalar_t
ctildei_vel
=
(
ctildei
-
ctildei_0
)
/
m_dt
;
aqconc
.
variable
(
node
)
=
ctildei
;
aqconc
.
velocity
(
node
)
=
ctildei_vel
;
// total solid concentration
MainVariable
&
solconc
=
vars
->
get_solid_concentration
(
component
);
const
scalar_t
cbari
=
extr
.
total_solid_concentration
(
component
);
const
scalar_t
cbari_0
=
solconc
.
predictor
(
node
);
const
scalar_t
cbari_vel
=
(
cbari
-
cbari_0
)
/
m_dt
;
solconc
.
variable
(
node
)
=
cbari
;
solconc
.
velocity
(
node
)
=
cbari_vel
;
solconc
.
chemistry_rate
(
node
)
=
-
cbari_vel
;
// partial pressure
if
(
vars
->
component_has_gas
(
component
))
{
MainVariable
&
pres_vars
=
vars
->
get_pressure_main_variables
(
component
);
index_t
id_g
=
vars
->
get_id_gas
(
component
);
const
scalar_t
pi
=
extr
.
fugacity_gas
(
id_g
)
*
vars
->
get_total_pressure
();
const
scalar_t
pi_0
=
pres_vars
.
predictor
(
node
);
const
scalar_t
pi_vel
=
(
pi
-
pi_0
)
/
m_dt
;
pres_vars
.
variable
(
node
)
=
pi
;
pres_vars
.
velocity
(
node
)
=
pi_vel
;
const
scalar_t
transient
=
(
(
pi
*
porosity
*
sat_g
)
-
(
pi_0
*
porosity_0
*
(
1.0
-
saturation_0
))
)
/
(
rt
*
m_dt
);
const
scalar_t
pi_chem_rate
=
-
transient
+
pres_vars
.
transport_fluxes
(
node
);
pres_vars
.
chemistry_rate
(
node
)
=
pi_chem_rate
;
}
}
}
void
EquilibriumStagger
::
EquilibriumStaggerImpl
::
initialize_timestep_one_node
(
index_t
node
,
UnsaturatedVariables
*
vars
)
{
SecondaryTransientVariable
&
porosity
=
vars
->
get_porosity
();
porosity
.
predictor
=
porosity
.
variable
;
porosity
.
velocity
.
setZero
();
}
EquilibriumStagger
::
StaggerReturnCode
EquilibriumStagger
::
EquilibriumStaggerImpl
::
restart_timestep
(
UnsaturatedVariables
*
vars
)
{
int
sum_retcode
=
0
;
#if SPECMICP_USE_OPENMP
#pragma omp parallel for \
reduction(+: sum_retcode)
#endif
for
(
index_t
node
=
0
;
node
<
m_mesh
->
nb_nodes
();
++
node
)
{
if
(
m_constraints
->
is_fixed_node
(
node
))
continue
;
sum_retcode
+=
compute_one_node
(
node
,
vars
);
}
if
(
sum_retcode
!=
0
)
{
return
StaggerReturnCode
::
UnknownError
;
}
return
StaggerReturnCode
::
ResidualMinimized
;
}
void
EquilibriumStagger
::
EquilibriumStaggerImpl
::
initialize_timestep
(
scalar_t
dt
,
UnsaturatedVariables
*
vars
)
{
m_dt
=
dt
;
// initialize
//water
{
MainVariable
&
cbar_w
=
vars
->
get_solid_concentration
(
0
);
cbar_w
.
predictor
=
cbar_w
.
variable
;
cbar_w
.
velocity
.
setZero
();
}
// aqueous component
for
(
index_t
component:
m_database
->
range_aqueous_component
())
{
MainVariable
&
cbar_i
=
vars
->
get_solid_concentration
(
component
);
cbar_i
.
predictor
=
cbar_i
.
variable
;
cbar_i
.
velocity
.
setZero
();
}
for
(
auto
node:
m_mesh
->
range_nodes
())
{
initialize_timestep_one_node
(
node
,
vars
);
}
}
void
EquilibriumStagger
::
EquilibriumStaggerImpl
::
initialize
(
UnsaturatedVariables
*
vars
)
{
// do nothing by default
}
}
//end namespace unsaturated
}
//end namespace systems
}
//end namespace reactmicp
}
//end namespace specmicp
Event Timeline
Log In to Comment