Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F76756940
transport_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, Aug 10, 06:04
Size
23 KB
Mime Type
text/x-c++
Expires
Mon, Aug 12, 06:04 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
19761262
Attached To
rSPECMICP SpecMiCP / ReactMiCP
transport_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 "transport_stagger.hpp"
#include "../../solver/staggers_base/stagger_structs.hpp"
#include "variables.hpp"
#include "variables_box.hpp"
#include "saturation_equation.hpp"
#include "pressure_equation.hpp"
#include "aqueous_equation.hpp"
#include "../../../database/database_holder.hpp"
#include "../../../dfpmsolver/parabolic_driver.hpp"
#include "../../../utils/log.hpp"
namespace
specmicp
{
namespace
reactmicp
{
namespace
systems
{
namespace
unsaturated
{
//! \brief Type of the equation
//! \internal
enum
class
EquationType
{
Saturation
,
LiquidAqueous
,
Pressure
};
//! \brief Format type to a string (for message error)
//! \internal
static
std
::
string
to_string
(
EquationType
eq_type
);
//! \brief Format DFPM solver retcode to a string
static
std
::
string
to_string
(
dfpmsolver
::
ParabolicDriverReturnCode
retcode
);
// forward declaration of wrapper classes
// They are defined at the bootom of this file
class
TaskBase
;
class
SaturationTask
;
using
VariablesBasePtr
=
solver
::
VariablesBasePtr
;
// =================================== //
// //
// Declaration of the implementation //
// //
// =================================== //
//! \brief Implementation class for the Unsaturated transport stagger
//! \internal
//!
//! This class does all the work
//! It was designed to be OpenMP compatible
class
UnsaturatedTransportStagger
::
UnsaturatedTransportStaggerImpl
{
public
:
// Constructor
// -----------
UnsaturatedTransportStaggerImpl
(
UnsaturatedVariablesPtr
variables
,
const
TransportConstraints
&
constraints
);
// Main functions
// --------------
//! \brief Return the residual
scalar_t
get_residual
(
UnsaturatedVariables
*
vars
);
//! \brief Return the first residual (start of timestep)
scalar_t
get_residual_0
(
UnsaturatedVariables
*
vars
);
//! \brief Return the update (norm of velocity vector)
scalar_t
get_update
(
UnsaturatedVariables
*
vars
);
//! \brief Initialize timestep
void
initialize_timestep
(
scalar_t
dt
,
UnsaturatedVariables
*
vars
);
//! \brief Restart the timestep
solver
::
StaggerReturnCode
restart_timestep
(
UnsaturatedVariables
*
vars
);
// Timestep management
// -------------------
//! \brief Save the timestep
void
save_dt
(
scalar_t
dt
)
{
m_dt
=
dt
;}
//! \brief Return the timestep
scalar_t
get_dt
()
{
return
m_dt
;}
private
:
scalar_t
m_norm_0
{
1.0
};
// timestep management
scalar_t
m_dt
{
-
1.0
};
//! \brief The saturation equations
std
::
unique_ptr
<
SaturationTask
>
m_saturation_equation
;
// Equation and solver
std
::
vector
<
std
::
unique_ptr
<
TaskBase
>>
m_equation_list
;
};
// Main functions
// ===============
// call of the implementation
UnsaturatedTransportStagger
::
UnsaturatedTransportStagger
(
std
::
shared_ptr
<
UnsaturatedVariables
>
variables
,
const
TransportConstraints
&
constraints
)
:
m_impl
(
make_unique
<
UnsaturatedTransportStaggerImpl
>
(
variables
,
constraints
))
{
}
UnsaturatedTransportStagger
::~
UnsaturatedTransportStagger
()
=
default
;
static
UnsaturatedVariables
*
get_var
(
VariablesBasePtr
var
)
{
return
static_cast
<
UnsaturatedVariables
*>
(
var
.
get
());
}
void
UnsaturatedTransportStagger
::
initialize_timestep
(
scalar_t
dt
,
VariablesBasePtr
var
)
{
UnsaturatedVariables
*
true_var
=
get_var
(
var
);
m_impl
->
initialize_timestep
(
dt
,
true_var
);
}
solver
::
StaggerReturnCode
UnsaturatedTransportStagger
::
restart_timestep
(
VariablesBasePtr
var
)
{
UnsaturatedVariables
*
true_var
=
get_var
(
var
);
return
m_impl
->
restart_timestep
(
true_var
);
}
scalar_t
UnsaturatedTransportStagger
::
get_residual
(
VariablesBasePtr
var
)
{
UnsaturatedVariables
*
true_var
=
get_var
(
var
);
return
m_impl
->
get_residual
(
true_var
);
}
scalar_t
UnsaturatedTransportStagger
::
get_residual_0
(
VariablesBasePtr
var
)
{
UnsaturatedVariables
*
true_var
=
get_var
(
var
);
return
m_impl
->
get_residual_0
(
true_var
);
}
scalar_t
UnsaturatedTransportStagger
::
get_update
(
VariablesBasePtr
var
)
{
UnsaturatedVariables
*
true_var
=
get_var
(
var
);
return
m_impl
->
get_update
(
true_var
);
}
// =============================== //
// =============================== //
// //
// Implementation details //
// ---------------------- //
// //
// =============================== //
// =============================== //
// 2 main sections
// - Solver wrappers : wrapper around 1 couple equation/solver
// - UnsaturatedTransportStaggerImpl : call the wrappers
// =============================== //
// //
// Solver wrappers //
// //
// =============================== //
// This wrappers are used to abstract the residual computation
// and the methods to solve every equations
//! \brief Base class for a task
//!
//! A task is how we solve governing equations,
//! and obtain residuals and update
//!
//! \internal
class
SPECMICP_DLL_LOCAL
TaskBase
{
public
:
TaskBase
(
index_t
component
,
EquationType
eq_type
)
:
m_type
(
eq_type
),
m_component
(
component
)
{}
//! \brief Compute the squared residuals
virtual
scalar_t
compute_squared_residual
(
UnsaturatedVariables
*
vars
)
=
0
;
//! \brief Compute the squared residuals at the beginning of the timestep
virtual
scalar_t
compute_squared_residual_0
(
UnsaturatedVariables
*
vars
)
=
0
;
//! \brief Compute the squared update of the variables
virtual
scalar_t
compute_squared_update
(
UnsaturatedVariables
*
vars
)
=
0
;
//! \brief Initialize the timestep
virtual
void
initialize_timestep
(
scalar_t
dt
,
UnsaturatedVariables
*
vars
)
=
0
;
//! \brief Solve the governing equation
virtual
dfpmsolver
::
ParabolicDriverReturnCode
restart_timestep
(
UnsaturatedVariables
*
vars
)
=
0
;
//! \brief Return the component index (in the db)
index_t
component
()
{
return
m_component
;}
//! \brief The equation type
EquationType
equation_type
()
{
return
m_type
;}
private
:
EquationType
m_type
;
index_t
m_component
;
};
//! \brief Base class for a equation solver
//!
//! \internal
template
<
typename
Derived
,
typename
traits
>
class
SPECMICP_DLL_LOCAL
EquationTask:
public
TaskBase
{
public
:
using
EqT
=
typename
traits
::
EqT
;
using
SolverT
=
typename
traits
::
SolverT
;
EquationTask
(
index_t
component
,
mesh
::
Mesh1DPtr
the_mesh
,
typename
traits
::
VariableBoxT
&
variables
,
const
TransportConstraints
&
constraints
)
:
TaskBase
(
component
,
traits
::
equation_type
),
m_equation
(
the_mesh
,
variables
,
constraints
),
m_solver
(
m_equation
)
{}
EquationTask
(
index_t
component
,
mesh
::
Mesh1DPtr
the_mesh
,
typename
traits
::
VariableBoxT
&
variables
,
const
TransportConstraints
&
constraints
,
scalar_t
scaling
)
:
TaskBase
(
component
,
traits
::
equation_type
),
m_equation
(
the_mesh
,
variables
,
constraints
),
m_solver
(
m_equation
)
{
m_equation
.
set_scaling
(
scaling
);
}
Derived
*
derived
()
{
return
static_cast
<
Derived
*>
(
this
);}
// This function must be implemented by subclass
MainVariable
&
get_var
(
UnsaturatedVariables
*
vars
)
{
return
derived
()
->
get_var
(
vars
);
}
scalar_t
compute_squared_residual
(
UnsaturatedVariables
*
vars
)
override
{
Vector
residuals
;
const
MainVariable
&
main
=
get_var
(
vars
);
m_equation
.
compute_residuals
(
main
.
variable
,
main
.
velocity
,
residuals
,
true
);
return
residuals
.
squaredNorm
()
/
m_norm_square_0
;
}
scalar_t
compute_squared_residual_0
(
UnsaturatedVariables
*
vars
)
override
{
Vector
residuals
;
const
MainVariable
&
main
=
get_var
(
vars
);
m_equation
.
compute_residuals
(
main
.
variable
,
main
.
velocity
,
residuals
,
false
);
m_norm_square_0
=
residuals
.
squaredNorm
();
if
(
m_norm_square_0
<
1e-8
)
m_norm_square_0
=
1.0
;
return
1.0
;
}
scalar_t
compute_squared_update
(
UnsaturatedVariables
*
vars
)
override
{
const
MainVariable
&
main
=
get_var
(
vars
);
return
main
.
velocity
.
squaredNorm
();
}
void
initialize_timestep
(
scalar_t
dt
,
UnsaturatedVariables
*
vars
)
override
{
m_dt
=
dt
;
MainVariable
&
main
=
get_var
(
vars
);
main
.
predictor
=
main
.
variable
;
main
.
transport_fluxes
.
setZero
();
main
.
velocity
.
setZero
();
m_solver
.
initialize_timestep
(
dt
,
get_var
(
vars
).
variable
);
}
dfpmsolver
::
ParabolicDriverReturnCode
restart_timestep
(
UnsaturatedVariables
*
vars
)
override
{
MainVariable
&
main
=
get_var
(
vars
);
m_solver
.
velocity
()
=
main
.
velocity
;
dfpmsolver
::
ParabolicDriverReturnCode
retcode
=
m_solver
.
restart_timestep
(
main
.
variable
);
if
(
retcode
>
dfpmsolver
::
ParabolicDriverReturnCode
::
NotConvergedYet
)
{
main
.
velocity
=
m_solver
.
velocity
();
}
return
retcode
;
}
dfpmsolver
::
ParabolicDriverOptions
&
get_options
()
{
return
m_solver
.
get_options
();}
const
dfpmsolver
::
ParabolicDriverOptions
&
get_options
()
const
{
return
m_solver
.
get_options
();}
scalar_t
get_dt
()
{
return
m_dt
;}
protected
:
scalar_t
m_norm_square_0
{
-
1
};
scalar_t
m_dt
{
-
1
};
EqT
m_equation
;
SolverT
m_solver
;
};
//! \brief Traits struct for the SaturationTask class
//!
//! \internal
struct
SPECMICP_DLL_LOCAL
SaturationTaskTraits
{
using
EqT
=
SaturationEquation
;
using
SolverT
=
dfpmsolver
::
ParabolicDriver
<
EqT
>
;
using
VariableBoxT
=
SaturationVariableBox
;
static
constexpr
EquationType
equation_type
{
EquationType
::
Saturation
};
};
//! \brief Wrapper for the saturation equation solver
//!
//! \internal
class
SPECMICP_DLL_LOCAL
SaturationTask:
public
EquationTask
<
SaturationTask
,
SaturationTaskTraits
>
{
using
base
=
EquationTask
<
SaturationTask
,
SaturationTaskTraits
>
;
public
:
SaturationTask
(
index_t
component
,
mesh
::
Mesh1DPtr
the_mesh
,
SaturationTaskTraits
::
VariableBoxT
&&
variables
,
const
TransportConstraints
&
constraints
)
:
EquationTask
<
SaturationTask
,
SaturationTaskTraits
>
(
component
,
the_mesh
,
variables
,
constraints
)
{}
SaturationTask
(
index_t
component
,
mesh
::
Mesh1DPtr
the_mesh
,
SaturationTaskTraits
::
VariableBoxT
&&
variables
,
const
TransportConstraints
&
constraints
,
scalar_t
scaling
)
:
EquationTask
<
SaturationTask
,
SaturationTaskTraits
>
(
component
,
the_mesh
,
variables
,
constraints
,
scaling
)
{}
MainVariable
&
get_var
(
UnsaturatedVariables
*
vars
)
{
return
vars
->
get_liquid_saturation
();
}
// Also take into account the solid total concentration
scalar_t
compute_squared_update
(
UnsaturatedVariables
*
vars
)
override
{
scalar_t
solid_update
=
vars
->
get_solid_concentration
(
component
()).
velocity
.
squaredNorm
();
return
solid_update
+
base
::
compute_squared_update
(
vars
);
}
};
//! \brief Traits struct for the LiquidAqueousTask class
//!
//! \internal
struct
SPECMICP_DLL_LOCAL
LiquidAqueousTaskTraits
{
using
EqT
=
AqueousTransportEquation
;
using
SolverT
=
dfpmsolver
::
ParabolicDriver
<
EqT
>
;
using
VariableBoxT
=
LiquidAqueousComponentVariableBox
;
static
constexpr
EquationType
equation_type
{
EquationType
::
LiquidAqueous
};
};
//! \brief Wrapper for the liquid transport of aqueous component equation
//!
//! \internal
class
SPECMICP_DLL_LOCAL
LiquidAqueousTask:
public
EquationTask
<
LiquidAqueousTask
,
LiquidAqueousTaskTraits
>
{
using
base
=
EquationTask
<
LiquidAqueousTask
,
LiquidAqueousTaskTraits
>
;
public
:
LiquidAqueousTask
(
index_t
component
,
mesh
::
Mesh1DPtr
the_mesh
,
LiquidAqueousTaskTraits
::
VariableBoxT
&&
variables
,
const
TransportConstraints
&
constraints
)
:
EquationTask
<
LiquidAqueousTask
,
LiquidAqueousTaskTraits
>
(
component
,
the_mesh
,
variables
,
constraints
)
{}
LiquidAqueousTask
(
index_t
component
,
mesh
::
Mesh1DPtr
the_mesh
,
LiquidAqueousTaskTraits
::
VariableBoxT
&&
variables
,
const
TransportConstraints
&
constraints
,
scalar_t
scaling
)
:
EquationTask
<
LiquidAqueousTask
,
LiquidAqueousTaskTraits
>
(
component
,
the_mesh
,
variables
,
constraints
,
scaling
)
{}
MainVariable
&
get_var
(
UnsaturatedVariables
*
vars
)
{
return
vars
->
get_aqueous_concentration
(
component
());
}
// Also take into account the solid total concentration
scalar_t
compute_squared_update
(
UnsaturatedVariables
*
vars
)
override
{
scalar_t
solid_update
=
vars
->
get_solid_concentration
(
component
()).
velocity
.
squaredNorm
();
return
solid_update
+
base
::
compute_squared_update
(
vars
);
}
};
//! \brief Traits struct for the Pressure Task traits
//!
//! \internal
struct
SPECMICP_DLL_LOCAL
PressureTaskTraits
{
using
EqT
=
PressureEquation
;
using
SolverT
=
dfpmsolver
::
ParabolicDriver
<
EqT
>
;
using
VariableBoxT
=
PressureVariableBox
;
static
constexpr
EquationType
equation_type
{
EquationType
::
Pressure
};
};
//! \brief Wrapper for the pressure equation solver
//!
//! \internal
class
SPECMICP_DLL_LOCAL
PressureTask:
public
EquationTask
<
PressureTask
,
PressureTaskTraits
>
{
using
base
=
EquationTask
<
PressureTask
,
PressureTaskTraits
>
;
public
:
PressureTask
(
index_t
component
,
mesh
::
Mesh1DPtr
the_mesh
,
PressureTaskTraits
::
VariableBoxT
&&
variables
,
const
TransportConstraints
&
constraints
)
:
EquationTask
<
PressureTask
,
PressureTaskTraits
>
(
component
,
the_mesh
,
variables
,
constraints
)
{}
PressureTask
(
index_t
component
,
mesh
::
Mesh1DPtr
the_mesh
,
PressureTaskTraits
::
VariableBoxT
&&
variables
,
const
TransportConstraints
&
constraints
,
scalar_t
scaling
)
:
EquationTask
<
PressureTask
,
PressureTaskTraits
>
(
component
,
the_mesh
,
variables
,
constraints
,
scaling
)
{}
MainVariable
&
get_var
(
UnsaturatedVariables
*
vars
)
{
return
vars
->
get_pressure_main_variables
(
component
());
}
};
// ================================== //
// //
// UnsaturatedTransportStaggerImpl //
// //
// ================================== //
// constructor
// ===========
UnsaturatedTransportStagger
::
UnsaturatedTransportStaggerImpl
::
UnsaturatedTransportStaggerImpl
(
UnsaturatedVariablesPtr
variables
,
const
TransportConstraints
&
constraints
)
{
database
::
RawDatabasePtr
raw_db
=
variables
->
get_database
();
mesh
::
Mesh1DPtr
the_mesh
=
variables
->
get_mesh
();
// Saturation equations
m_saturation_equation
=
make_unique
<
SaturationTask
>
(
0
,
the_mesh
,
variables
->
get_saturation_variables
(),
constraints
,
variables
->
get_aqueous_scaling
(
0
)
);
dfpmsolver
::
ParabolicDriverOptions
&
opts
=
m_saturation_equation
->
get_options
();
opts
.
residuals_tolerance
=
1e-8
;
opts
.
absolute_tolerance
=
1e-12
;
opts
.
step_tolerance
=
1e-12
;
opts
.
threshold_stationary_point
=
1e-12
;
opts
.
maximum_iterations
=
500
;
opts
.
sparse_solver
=
sparse_solvers
::
SparseSolver
::
SparseLU
;
opts
.
quasi_newton
=
3
;
opts
.
maximum_step_length
=
0.1
;
opts
.
max_iterations_at_max_length
=
5000
;
const
index_t
size
=
raw_db
->
nb_aqueous_components
()
+
variables
->
nb_gas
();
m_equation_list
.
reserve
(
size
);
// Liquid aqueous diffusion-advection
for
(
index_t
id:
raw_db
->
range_aqueous_component
())
{
m_equation_list
.
emplace_back
(
make_unique
<
LiquidAqueousTask
>
(
id
,
the_mesh
,
variables
->
get_liquid_aqueous_component_variables
(
id
),
constraints
,
variables
->
get_aqueous_scaling
(
id
))
);
dfpmsolver
::
ParabolicDriverOptions
&
opts
=
static_cast
<
LiquidAqueousTask
*>
(
m_equation_list
[
m_equation_list
.
size
()
-
1
].
get
())
->
get_options
();
opts
.
maximum_iterations
=
500
;
opts
.
residuals_tolerance
=
1e-8
;
opts
.
absolute_tolerance
=
1e-12
;
opts
.
step_tolerance
=
1e-12
;
opts
.
threshold_stationary_point
=
1e-12
;
opts
.
sparse_solver
=
sparse_solvers
::
SparseSolver
::
SparseLU
;
opts
.
maximum_step_length
=
1e6
;
}
// Gaseous diffusion
for
(
index_t
id:
raw_db
->
range_component
())
{
if
(
variables
->
component_has_gas
(
id
))
{
m_equation_list
.
emplace_back
(
make_unique
<
PressureTask
>
(
id
,
the_mesh
,
variables
->
get_pressure_variables
(
id
),
constraints
,
variables
->
get_gaseous_scaling
(
id
)
));
dfpmsolver
::
ParabolicDriverOptions
&
opts
=
static_cast
<
PressureTask
*>
(
m_equation_list
[
m_equation_list
.
size
()
-
1
].
get
())
->
get_options
();
opts
.
maximum_iterations
=
500
;
opts
.
residuals_tolerance
=
1e-8
;
opts
.
absolute_tolerance
=
1e-12
;
opts
.
step_tolerance
=
1e-12
;
opts
.
threshold_stationary_point
=
1e-12
;
opts
.
sparse_solver
=
sparse_solvers
::
SparseSolver
::
SparseLU
;
opts
.
maximum_step_length
=
1e6
;
}
}
}
// Residuals
// =========
scalar_t
UnsaturatedTransportStagger
::
UnsaturatedTransportStaggerImpl
::
get_residual
(
UnsaturatedVariables
*
vars
)
{
scalar_t
sum
=
m_saturation_equation
->
compute_squared_residual
(
vars
);
for
(
auto
&
task:
m_equation_list
)
{
sum
+=
task
->
compute_squared_residual
(
vars
);
}
auto
norm
=
std
::
sqrt
(
sum
);
return
norm
;
}
scalar_t
UnsaturatedTransportStagger
::
UnsaturatedTransportStaggerImpl
::
get_residual_0
(
UnsaturatedVariables
*
vars
)
{
return
m_norm_0
;
}
scalar_t
UnsaturatedTransportStagger
::
UnsaturatedTransportStaggerImpl
::
get_update
(
UnsaturatedVariables
*
vars
)
{
scalar_t
sum
=
m_saturation_equation
->
compute_squared_update
(
vars
);
for
(
auto
&
task:
m_equation_list
)
{
sum
+=
task
->
compute_squared_update
(
vars
);
}
return
std
::
sqrt
(
sum
);
}
// Solving the equations
// =====================
void
UnsaturatedTransportStagger
::
UnsaturatedTransportStaggerImpl
::
initialize_timestep
(
scalar_t
dt
,
UnsaturatedVariables
*
vars
)
{
save_dt
(
dt
);
vars
->
get_advection_flux
().
set_zero
();
m_saturation_equation
->
initialize_timestep
(
dt
,
vars
);
scalar_t
sum
=
m_saturation_equation
->
compute_squared_residual_0
(
vars
);
for
(
auto
&
task:
m_equation_list
)
{
task
->
initialize_timestep
(
dt
,
vars
);
sum
+=
task
->
compute_squared_residual_0
(
vars
);
}
m_norm_0
=
1.0
;
//std::sqrt(sum);
}
solver
::
StaggerReturnCode
UnsaturatedTransportStagger
::
UnsaturatedTransportStaggerImpl
::
restart_timestep
(
UnsaturatedVariables
*
vars
)
{
dfpmsolver
::
ParabolicDriverReturnCode
retcode
;
bool
flag_fail
=
false
;
bool
flag_error_minimized
=
false
;
// Saturation equation
retcode
=
m_saturation_equation
->
restart_timestep
(
vars
);
if
(
dfpmsolver
::
has_failed
(
retcode
))
{
WARNING
<<
"Failed to solve saturation equation, return code :"
<<
to_string
(
retcode
);
flag_fail
=
true
;
}
if
(
retcode
==
dfpmsolver
::
ParabolicDriverReturnCode
::
ErrorMinimized
)
flag_error_minimized
=
true
;
// Other equation
if
(
not
(
flag_fail
or
flag_error_minimized
))
{
vars
->
set_relative_variables
();
for
(
auto
&
task:
m_equation_list
)
{
retcode
=
task
->
restart_timestep
(
vars
);
if
(
dfpmsolver
::
has_failed
(
retcode
))
{
WARNING
<<
"Equation of type '"
<<
to_string
(
task
->
equation_type
())
<<
"' for component "
<<
task
->
component
()
<<
" has failed with return code : "
<<
to_string
(
retcode
);
flag_fail
=
true
;
}
flag_error_minimized
=
flag_error_minimized
and
(
retcode
==
dfpmsolver
::
ParabolicDriverReturnCode
::
ErrorMinimized
);
}
}
// Return code
solver
::
StaggerReturnCode
return_code
=
solver
::
StaggerReturnCode
::
ResidualMinimized
;
if
(
flag_error_minimized
)
{
return_code
=
solver
::
StaggerReturnCode
::
ErrorMinimized
;
}
if
(
flag_fail
)
return_code
=
solver
::
StaggerReturnCode
::
UnknownError
;
return
return_code
;
}
// ================================ //
// //
// Helper functions //
// //
// ================================ //
static
std
::
string
to_string
(
EquationType
eq_type
)
{
std
::
string
str
;
switch
(
eq_type
)
{
case
EquationType
::
LiquidAqueous:
str
=
"Liquid advection-diffusion"
;
break
;
case
EquationType
::
Pressure:
str
=
"Gaseous diffusion"
;
break
;
case
EquationType
::
Saturation:
str
=
"Saturation equation"
;
break
;
default
:
str
=
"Unknown equation"
;
break
;
}
return
str
;
}
std
::
string
to_string
(
dfpmsolver
::
ParabolicDriverReturnCode
retcode
)
{
using
RetCode
=
dfpmsolver
::
ParabolicDriverReturnCode
;
std
::
string
str
;
switch
(
retcode
)
{
case
RetCode
::
ResidualMinimized:
str
=
"ResidualMinimized"
;
break
;
case
RetCode
::
ErrorMinimized:
str
=
"ErrorMinimized"
;
break
;
case
RetCode
::
NotConvergedYet:
str
=
"NotConvergedYet"
;
break
;
case
RetCode
::
StationaryPoint:
str
=
"StationaryPoint"
;
break
;
case
RetCode
::
ErrorLinearSystem:
str
=
"ErrorLinearSystem"
;
break
;
case
RetCode
::
MaxStepTakenTooManyTimes:
str
=
"MaxStepTakenTooManyTimes"
;
break
;
case
RetCode
::
MaxIterations:
str
=
"MaxIterations"
;
break
;
case
RetCode
::
LinesearchFailed:
str
=
"LinesearchFailed"
;
default
:
str
=
"Unknown return code"
;
break
;
}
return
str
;
}
}
//end namespace unsaturated
}
//end namespace systems
}
//end namespace reactmicp
}
//end namespace specmicp
Event Timeline
Log In to Comment