Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79132690
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
Sat, Aug 24, 22:17
Size
5 KB
Mime Type
text/x-c
Expires
Mon, Aug 26, 22:17 (2 d)
Engine
blob
Format
Raw Data
Handle
20151015
Attached To
rSPECMICP SpecMiCP / ReactMiCP
reactive_transport_solver.cpp
View Options
#include "reactive_transport_solver.hpp"
#include "staggers_base/staggers_base.hpp"
#include "utils/log.hpp"
#include <iostream>
#include <chrono>
namespace
specmicp
{
namespace
reactmicp
{
namespace
solver
{
namespace
internal
{
// This contains internal information for the reactive transport solver
// It is used to check the convergence
struct
ReactiveTransportResiduals
{
index_t
nb_iterations
;
scalar_t
transport_residual_0
;
scalar_t
transport_residuals
;
scalar_t
update
;
ReactiveTransportResiduals
()
:
nb_iterations
(
0
),
transport_residual_0
(
-
1
),
transport_residuals
(
-
1
),
update
(
-
1
)
{}
};
}
// end namespace internal
// ReactiveTransportSolver::
// Solve a timestep
ReactiveTransportReturnCode
ReactiveTransportSolver
::
solve_timestep
(
scalar_t
timestep
,
VariablesBasePtr
variables
)
{
internal
::
ReactiveTransportResiduals
residuals
;
// initialization
// copy of variables // ?
m_transport_stagger
->
initialize_timestep
(
timestep
,
variables
);
m_chemistry_stagger
->
initialize_timestep
(
timestep
,
variables
);
m_upscaling_stagger
->
initialize_timestep
(
timestep
,
variables
);
//
residuals
.
transport_residual_0
=
m_transport_stagger
->
get_residual_0
(
variables
);
ReactiveTransportReturnCode
retcode
=
one_iteration
(
variables
,
residuals
);
// check for failure
if
(
retcode
<
ReactiveTransportReturnCode
::
StaggerFailure
)
{
ERROR
<<
"Failed to solve the iteration, return code : "
<<
(
int
)
retcode
;
return
retcode
;
}
retcode
=
check_convergence
(
variables
,
residuals
);
++
residuals
.
nb_iterations
;
// if sequential non-iterative algorithm
if
(
get_options
().
is_snia
())
{
goto
end
;
}
// else
while
(
retcode
==
ReactiveTransportReturnCode
::
NotConvergedYet
)
{
retcode
=
one_iteration
(
variables
,
residuals
);
// check for failure
if
(
retcode
<
ReactiveTransportReturnCode
::
StaggerFailure
)
{
ERROR
<<
"Failed to solve the iteration, return code : "
<<
(
int
)
retcode
;
return
retcode
;
}
retcode
=
check_convergence
(
variables
,
residuals
);
++
residuals
.
nb_iterations
;
}
// wrapup
end:
// upscaling, if needed
if
(
not
get_options
().
implicit_upscaling
)
{
std
::
chrono
::
time_point
<
std
::
chrono
::
system_clock
>
start
,
end
;
std
::
chrono
::
duration
<
scalar_t
>
elapsed_seconds
;
start
=
std
::
chrono
::
system_clock
::
now
();
m_upscaling_stagger
->
restart_timestep
(
variables
);
end
=
std
::
chrono
::
system_clock
::
now
();
elapsed_seconds
=
end
-
start
;
m_timer
.
upscaling_time
+=
elapsed_seconds
.
count
();
}
// record performance
get_perfs
().
nb_iterations
=
residuals
.
nb_iterations
;
get_perfs
().
return_code
=
retcode
;
get_perfs
().
residuals
=
residuals
.
transport_residuals
/
residuals
.
transport_residual_0
;
return
retcode
;
}
ReactiveTransportReturnCode
ReactiveTransportSolver
::
one_iteration
(
VariablesBasePtr
variables
,
internal
::
ReactiveTransportResiduals
&
residuals
)
{
// Transport
// ---------
std
::
chrono
::
time_point
<
std
::
chrono
::
system_clock
>
start
,
end
;
std
::
chrono
::
duration
<
scalar_t
>
elapsed_seconds
;
start
=
std
::
chrono
::
system_clock
::
now
();
StaggerReturnCode
transport_ret_code
=
m_transport_stagger
->
restart_timestep
(
variables
);
if
(
transport_ret_code
<=
StaggerReturnCode
::
NotConvergedYet
)
{
return
ReactiveTransportReturnCode
::
TransportFailure
;
}
end
=
std
::
chrono
::
system_clock
::
now
();
elapsed_seconds
=
end
-
start
;
m_timer
.
transport_time
+=
elapsed_seconds
.
count
();
// Chemistry
// ---------
start
=
std
::
chrono
::
system_clock
::
now
();
StaggerReturnCode
chemistry_ret_code
=
m_chemistry_stagger
->
restart_timestep
(
variables
);
if
(
chemistry_ret_code
<=
StaggerReturnCode
::
NotConvergedYet
)
{
return
ReactiveTransportReturnCode
::
ChemistryFailure
;
}
end
=
std
::
chrono
::
system_clock
::
now
();
elapsed_seconds
=
end
-
start
;
m_timer
.
chemistry_time
+=
elapsed_seconds
.
count
();
// Upscaling0
// ---------
if
(
get_options
().
implicit_upscaling
)
{
start
=
std
::
chrono
::
system_clock
::
now
();
StaggerReturnCode
upscaling_ret_code
=
m_upscaling_stagger
->
restart_timestep
(
variables
);
if
(
upscaling_ret_code
<=
StaggerReturnCode
::
NotConvergedYet
)
{
return
ReactiveTransportReturnCode
::
UpscalingFailure
;
}
end
=
std
::
chrono
::
system_clock
::
now
();
elapsed_seconds
=
end
-
start
;
m_timer
.
upscaling_time
+=
elapsed_seconds
.
count
();
}
// Final residuals
// ---------------
residuals
.
transport_residuals
=
m_transport_stagger
->
get_residual
(
variables
);
residuals
.
update
=
m_transport_stagger
->
get_update
(
variables
);
return
ReactiveTransportReturnCode
::
NotConvergedYet
;
}
// Check the convergence
ReactiveTransportReturnCode
ReactiveTransportSolver
::
check_convergence
(
VariablesBasePtr
_
,
const
internal
::
ReactiveTransportResiduals
&
residuals
)
{
// Residual
if
(
residuals
.
transport_residuals
/
residuals
.
transport_residual_0
<
get_options
().
residuals_tolerance
or
residuals
.
transport_residuals
<
get_options
().
absolute_residuals_tolerance
)
{
return
ReactiveTransportReturnCode
::
ResidualMinimized
;
}
// Step
if
(
residuals
.
update
<
get_options
().
step_tolerance
)
{
if
(
residuals
.
transport_residuals
/
residuals
.
transport_residual_0
<
get_options
().
good_enough_tolerance
)
{
return
ReactiveTransportReturnCode
::
ErrorMinimized
;
}
else
return
ReactiveTransportReturnCode
::
StationaryPoint
;
}
// Number of iterations
if
(
residuals
.
nb_iterations
>=
get_options
().
maximum_iterations
)
{
return
ReactiveTransportReturnCode
::
MaximumIterationsReached
;
}
return
ReactiveTransportReturnCode
::
NotConvergedYet
;
}
}
// end namespace solver
}
// end namespace reactmicp
}
// end namespace specmicp
Event Timeline
Log In to Comment