Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82220782
reactive_transport_solver.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
Tue, Sep 10, 06:05
Size
9 KB
Mime Type
text/x-c
Expires
Thu, Sep 12, 06:05 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20664546
Attached To
rSPECMICP SpecMiCP / ReactMiCP
reactive_transport_solver.cpp
View Options
//#define EIGEN_DEFAULT_IO_FORMAT IOFormat(8)
#include "reactive_transport_solver.hpp"
#include "specmicp/reduced_system_solver.hpp"
#include <omp.h>
namespace
specmicp
{
namespace
reactmicp
{
namespace
internals
{
//! \brief Store residual information
struct
SIAResiduals
{
int
nb_iterations
;
double
transport_residual_0
;
double
transport_residual
;
double
transport_update
;
Eigen
::
VectorXd
speciation_residual
;
Eigen
::
VectorXd
speciation_update
;
SIAResiduals
(
int
nb_node
)
:
nb_iterations
(
0
),
speciation_residual
(
Eigen
::
VectorXd
::
Zero
(
nb_node
)),
speciation_update
(
Eigen
::
VectorXd
::
Zero
(
nb_node
))
{}
};
}
// end namespace internals
void
SIASaturatedReactiveTransportSolver
::
apply_boundary_conditions
(
SIABoundaryConditions
bcs
)
{
assert
(
bcs
.
bs_types
.
size
()
==
(
unsigned
int
)
m_mesh
->
nnodes
());
assert
(
bcs
.
initial_states
.
size
()
==
(
unsigned
int
)
m_mesh
->
nnodes
());
std
::
vector
<
int
>
list_fixed_nodes
;
for
(
int
node:
m_mesh
->
range_nodes
())
{
switch
(
bcs
.
bs_types
[
node
])
{
case
BoundaryConditionType
::
FixedComposition:
m_bcs
(
node
)
=
1
;
list_fixed_nodes
.
push_back
(
node
);
break
;
default
:
m_bcs
(
node
)
=
0
;
break
;
}
m_variables
.
update_composition
(
node
,
bcs
.
list_initial_states
[
bcs
.
initial_states
[
node
]]);
m_transport_program
.
number_equations
(
list_fixed_nodes
);
}
}
double
SIASaturatedReactiveTransportSolver
::
residuals_transport
()
{
Eigen
::
VectorXd
residuals
;
m_transport_solver
.
compute_residuals
(
m_variables
.
total_mobile_concentrations
(),
residuals
);
//std::cout << "Residuals : " << std::endl << residuals << std::endl;
return
residuals
.
norm
();
}
SIASaturatedReactiveTransportSolverReturnCode
convert_dfpm_retcode
(
dfpmsolver
::
ParabolicDriverReturnCode
retcode
)
{
switch
(
retcode
)
{
case
dfpmsolver
::
ParabolicDriverReturnCode
::
MaxIterations:
return
SIASaturatedReactiveTransportSolverReturnCode
::
MaxIterationsInTransport
;
case
dfpmsolver
::
ParabolicDriverReturnCode
::
StationaryPoint:
return
SIASaturatedReactiveTransportSolverReturnCode
::
StationaryPointsInTransport
;
default
:
return
SIASaturatedReactiveTransportSolverReturnCode
::
ErrorInTransport
;
}
}
SIASaturatedReactiveTransportSolverReturnCode
SIASaturatedReactiveTransportSolver
::
solve_timestep
(
double
dt
)
{
retcode_t
return_code
=
retcode_t
::
NotConvergedYet
;
internals
::
SIAResiduals
residuals_storage
(
m_mesh
->
nnodes
());
m_flag_compute_speciation
=
Eigen
::
VectorXi
::
Ones
(
m_mesh
->
nnodes
());
// initialization transport program and transport solver
Eigen
::
VectorXd
transport_variable
(
m_variables
.
total_mobile_concentrations
());
m_transport_solver
.
get_velocity
()
=
Eigen
::
VectorXd
::
Zero
(
transport_variable
.
rows
());
m_transport_program
.
external_flux
().
setZero
();
//
residuals_storage
.
transport_residual_0
=
residuals_transport
();
m_minerals_start_iteration
=
m_variables
.
speciation_variables
().
block
(
m_database
->
nb_component
,
0
,
m_database
->
nb_mineral
,
m_mesh
->
nnodes
());
// solve transport
// ---------------
dfpmsolver
::
ParabolicDriverReturnCode
retcode
=
m_transport_solver
.
solve_timestep
(
dt
,
transport_variable
);
if
(
retcode
!=
dfpmsolver
::
ParabolicDriverReturnCode
::
ResidualMinimized
)
{
ERROR
<<
"Failed to solve 1st transport problem; return code "
<<
(
int
)
retcode
;
return
convert_dfpm_retcode
(
retcode
);
}
residuals_storage
.
transport_residual
=
m_transport_solver
.
get_perfs
().
current_residual
;
residuals_storage
.
transport_update
=
m_transport_solver
.
get_perfs
().
current_update
;
// update
m_variables
.
total_mobile_concentrations
()
=
transport_variable
;
// Speciation
// ----------
restart_speciation_iterations
(
dt
,
residuals_storage
);
// Convergence
// -----------
return_code
=
check_convergence
(
residuals_storage
);
residuals_storage
.
nb_iterations
+=
1
;
// If SNIA, we stop
if
(
get_options
().
max_iterations
==
1
)
return
retcode_t
::
Success
;
// SIA algorithm
// -------------
while
(
return_code
==
retcode_t
::
NotConvergedYet
)
{
//std::cout << m_flag_compute_speciation.sum() << std::endl;
transport_variable
=
m_variables
.
total_mobile_concentrations
();
retcode
=
m_transport_solver
.
restart_timestep
(
transport_variable
);
if
(
retcode
!=
dfpmsolver
::
ParabolicDriverReturnCode
::
ResidualMinimized
)
{
if
(
retcode
==
dfpmsolver
::
ParabolicDriverReturnCode
::
MaxIterations
and
m_transport_solver
.
get_perfs
().
current_residual
/
residuals_storage
.
transport_residual_0
<
get_options
().
good_enough_transport_residual_tolerance
)
{
restart_speciation_iterations
(
dt
,
residuals_storage
);
return_code
=
retcode_t
::
Success
;
break
;
}
ERROR
<<
"Failed to solve transport problem; return code "
<<
(
int
)
retcode
;
return
convert_dfpm_retcode
(
retcode
);
}
residuals_storage
.
transport_residual
=
m_transport_solver
.
get_perfs
().
current_residual
;
residuals_storage
.
transport_update
=
m_transport_solver
.
get_perfs
().
current_update
;
m_variables
.
total_mobile_concentrations
()
=
transport_variable
;
residuals_storage
.
nb_iterations
+=
1
;
restart_speciation_iterations
(
dt
,
residuals_storage
);
return_code
=
check_convergence
(
residuals_storage
);
}
std
::
cout
<<
"End fixed point : return code "
<<
(
int
)
return_code
<<
" - nb iter : "
<<
residuals_storage
.
nb_iterations
<<
std
::
endl
;
return
return_code
;
}
void
SIASaturatedReactiveTransportSolver
::
restart_transport_iterations
(
internals
::
SIAResiduals
&
residuals_storage
)
{
Eigen
::
VectorXd
transport_variable
(
m_variables
.
total_mobile_concentrations
());
dfpmsolver
::
ParabolicDriverReturnCode
retcode
=
m_transport_solver
.
restart_timestep
(
transport_variable
);
if
(
retcode
!=
dfpmsolver
::
ParabolicDriverReturnCode
::
ResidualMinimized
)
{
ERROR
<<
"Failed to solve transport problem; return code "
<<
(
int
)
retcode
;
throw
std
::
runtime_error
(
"Failed to solve transport program"
);
}
residuals_storage
.
transport_residual
=
m_transport_solver
.
get_perfs
().
current_residual
;
residuals_storage
.
transport_update
=
m_transport_solver
.
get_perfs
().
current_update
;
}
void
SIASaturatedReactiveTransportSolver
::
restart_speciation_iterations
(
double
dt
,
internals
::
SIAResiduals
&
residuals_storage
)
{
#ifdef _OPENMP
#pragma omp parallel
{
#pragma omp for schedule(static)
for
(
int
node
=
0
;
node
<
m_mesh
->
nnodes
();
++
node
)
#else
for
(
int
node:
m_mesh
->
range_nodes
())
#endif
{
if
(
m_bcs
(
node
)
==
0
and
m_flag_compute_speciation
(
node
))
{
solve_speciation_node
(
dt
,
node
,
residuals_storage
);
}
}
#ifdef _OPENMP
}
#endif
}
void
SIASaturatedReactiveTransportSolver
::
solve_speciation_node
(
double
dt
,
int
node
,
internals
::
SIAResiduals
&
residuals_storage
)
{
Eigen
::
VectorXd
totconc
;
m_variables
.
nodal_update_total_amount
(
node
,
totconc
);
//std::cout << "Node : " << node << " - Total concentration : " << std::endl << totconc << std::endl;
ReducedSystemSolver
solver
(
m_database
,
totconc
);
solver
.
get_options
().
solver_options
.
max_iter
=
25
;
Eigen
::
VectorXd
variables
=
m_variables
.
speciation_variables
(
node
);
micpsolver
::
MiCPPerformance
perf
=
solver
.
solve
(
variables
);
if
(
perf
.
return_code
!=
micpsolver
::
MiCPSolverReturnCode
::
ResidualMinimized
)
{
WARNING
<<
"Failed to solve speciation solver for node "
<<
node
<<
"; return code = "
<<
(
int
)
perf
.
return_code
;
}
residuals_storage
.
speciation_residual
(
node
)
=
perf
.
current_residual
;
residuals_storage
.
speciation_update
(
node
)
=
perf
.
current_update
;
if
(
perf
.
current_update
<
get_options
().
threshold_update_disable_speciation
)
m_flag_compute_speciation
(
node
)
=
0
;
EquilibriumState
solution
=
solver
.
get_solution
(
variables
);
m_variables
.
update_composition
(
node
,
solution
);
// mineral
const
double
denom
=
1e-3
*
m_mesh
->
cell_volume
(
node
)
*
dt
;
for
(
int
component:
m_database
->
range_aqueous_component
())
{
double
sum
=
0.0
;
for
(
int
mineral:
m_database
->
range_mineral
())
{
if
(
m_database
->
nu_mineral
(
mineral
,
component
)
==
0.0
)
continue
;
sum
-=
m_database
->
nu_mineral
(
mineral
,
component
)
*
(
m_variables
.
mineral_amount
(
node
,
mineral
)
-
m_minerals_start_iteration
(
mineral
,
node
)
);
}
m_transport_program
.
external_flux
(
node
,
component
)
=
sum
/
denom
;
}
}
SIASaturatedReactiveTransportSolverReturnCode
SIASaturatedReactiveTransportSolver
::
check_convergence
(
internals
::
SIAResiduals
&
residuals_storage
)
{
// compute residuals of transport
const
double
norm_transport
=
residuals_transport
();
retcode_t
return_code
=
retcode_t
::
NotConvergedYet
;
if
(
norm_transport
/
residuals_storage
.
transport_residual_0
<
get_options
().
residual_tolerance
)
{
return_code
=
retcode_t
::
Success
;
}
else
if
(
residuals_storage
.
transport_update
<
get_options
().
threshold_update_transport
)
{
return_code
=
retcode_t
::
Success
;
}
else
if
(
residuals_storage
.
nb_iterations
>=
get_options
().
max_iterations
)
{
WARNING
<<
"Maximum number of iterations ("
<<
get_options
().
max_iterations
<<
") reached during fixed-point iterations"
;
return_code
=
retcode_t
::
MaxIterations
;
}
return
return_code
;
}
}
// end namespace reactmicp
}
// end namespace specmicp
Event Timeline
Log In to Comment