Homec4science
Diffusion SpecMiCP / ReactMiCP (bin/unsaturated) Browse

specmicp/src/reactmicp/solver/reactive_transport_solver.cpp
a1bb082f067bbin/unsaturated

Last Change
Raw File

reactive_transport_solver.cpp
View Options

#include "reactive_transport_solver.hpp"
#include "staggers_base/staggers_base.hpp"
#include "utils/log.hpp"
#include "utils/timer.hpp"
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
//
// Attention : This function uses goto to handle return code and performance
ReactiveTransportReturnCode ReactiveTransportSolver::solve_timestep(
scalar_t timestep,
VariablesBasePtr variables
)
{
// start the timer
Timer tot_timer;
tot_timer.start();
// initialization
internal::ReactiveTransportResiduals residuals;
reset_perfs();
get_perfs().timestep = timestep;
// 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;
goto set_return;
}
retcode = check_convergence(variables, residuals, false);
++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;
goto set_return;
}
if (retcode != ReactiveTransportReturnCode::TransportBypass)
retcode = check_convergence(variables, residuals, true);
retcode = check_convergence(variables, residuals, false);
++residuals.nb_iterations;
}
// wrapup
end:
// upscaling, if needed
if (not get_options().implicit_upscaling)
{
Timer timer;
timer.start();
m_upscaling_stagger->restart_timestep(variables);
timer.stop();
m_timer.upscaling_time += timer.elapsed_time();
}
// record performance and return code
set_return:
tot_timer.stop();
get_perfs().total_time = tot_timer.elapsed_time();
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
)
{
Timer timer;
bool bypass = false;
// Transport
// ---------
timer.start();
StaggerReturnCode transport_ret_code = m_transport_stagger->restart_timestep(variables);
if (transport_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::TransportFailure;
}
if (transport_ret_code == StaggerReturnCode::ErrorMinimized)
{
bypass = true;
}
timer.stop();
const scalar_t ttime = timer.elapsed_time();
m_timer.transport_time += ttime;
get_perfs().transport_time += ttime;
// Chemistry
// ---------
timer.start();
StaggerReturnCode chemistry_ret_code = m_chemistry_stagger->restart_timestep(variables);
if (chemistry_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::ChemistryFailure;
}
timer.stop();
const scalar_t ctime = timer.elapsed_time();
m_timer.chemistry_time += ctime;
get_perfs().chemistry_time += ctime;
// Upscaling0
// ---------
if (get_options().implicit_upscaling)
{
timer.start();
StaggerReturnCode upscaling_ret_code = m_upscaling_stagger->restart_timestep(variables);
if (upscaling_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::UpscalingFailure;
}
timer.stop();
m_timer.upscaling_time += timer.elapsed_time();
}
// Final residuals
// ---------------
residuals.transport_residuals = m_transport_stagger->get_residual(variables);
residuals.update = m_transport_stagger->get_update(variables);
if (bypass)
return ReactiveTransportReturnCode::TransportBypass;
return ReactiveTransportReturnCode::NotConvergedYet;
}
// Check the convergence
ReactiveTransportReturnCode ReactiveTransportSolver::check_convergence(
VariablesBasePtr _,
const internal::ReactiveTransportResiduals& residuals,
bool bypass
)
{
const scalar_t relative_residual = residuals.transport_residuals/residuals.transport_residual_0;
// Residual
if (relative_residual < get_options().residuals_tolerance
or residuals.transport_residuals < get_options().absolute_residuals_tolerance)
{
return ReactiveTransportReturnCode::ResidualMinimized;
}
// Step
else if (residuals.update < get_options().step_tolerance)
{
if (relative_residual < get_options().good_enough_tolerance)
{
return ReactiveTransportReturnCode::ErrorMinimized;
}
else
return ReactiveTransportReturnCode::StationaryPoint;
}
else if (bypass)
{
return ReactiveTransportReturnCode::TransportBypass;
}
// Number of iterations
else if (residuals.nb_iterations >= get_options().maximum_iterations)
{
if (relative_residual < get_options().good_enough_tolerance)
{
return ReactiveTransportReturnCode::GoodEnough;
}
return ReactiveTransportReturnCode::MaximumIterationsReached;
}
return ReactiveTransportReturnCode::NotConvergedYet;
}
} // end namespace solver
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help