Homec4science
Diffusion SpecMiCP / ReactMiCP (specmicp_bin) Browse

specmicp/src/reactmicp/solver/reactive_transport_solver.cpp
9a4089e18a62specmicp_bin

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 <iostream>
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)
m_upscaling_stagger->restart_timestep(variables);
// 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
// ---------
StaggerReturnCode transport_ret_code = m_transport_stagger->restart_timestep(variables);
if (transport_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::TransportFailure;
}
// Chemistry
// ---------
StaggerReturnCode chemistry_ret_code = m_chemistry_stagger->restart_timestep(variables);
if (chemistry_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::ChemistryFailure;
}
// Upscaling
// ---------
if (get_options().implicit_upscaling)
{
StaggerReturnCode upscaling_ret_code = m_upscaling_stagger->restart_timestep(variables);
if (upscaling_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::UpscalingFailure;
}
}
// 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
c4science · Help