Homec4science
Diffusion SpecMiCP / ReactMiCP (yaml/safe) Browse

specmicp/src/reactmicp/solver/reactive_transport_solver.cpp
6f21451dd340yaml/safe

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 {
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())
{
return retcode;
}
// else
while (retcode == ReactiveTransportReturnCode::NotConvergedYet)
{
std::cout << "iter : " << residuals.nb_iterations << std::endl;
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;
}
return retcode;
}
ReactiveTransportReturnCode ReactiveTransportSolver::one_iteration(
VariablesBasePtr variables,
internal::ReactiveTransportResiduals& residuals
)
{
StaggerReturnCode transport_ret_code = m_transport_stagger->restart_timestep(variables);
if (transport_ret_code <= StaggerReturnCode::NotconvergedYet)
{
return ReactiveTransportReturnCode::TransportFailure;
}
StaggerReturnCode chemistry_ret_code = m_chemistry_stagger->restart_timestep(variables);
if (chemistry_ret_code <= StaggerReturnCode::NotconvergedYet)
{
return ReactiveTransportReturnCode::ChemistryFailure;
}
if (get_options().implicit_upscaling)
{
StaggerReturnCode upscaling_ret_code = m_upscaling_stagger->restart_timestep(variables);
if (upscaling_ret_code <= StaggerReturnCode::NotconvergedYet)
{
return ReactiveTransportReturnCode::UpscalingFailure;
}
}
residuals.transport_residuals = m_transport_stagger->get_residual(variables);
residuals.update = m_transport_stagger->get_update(variables);
}
// Check the convergence
ReactiveTransportReturnCode ReactiveTransportSolver::check_convergence(
VariablesBasePtr _,
const internal::ReactiveTransportResiduals& residuals
)
{
if (residuals.transport_residuals/residuals.transport_residual_0 < get_options().residuals_tolerance
or residuals.transport_residuals < get_options().absolute_residuals_tolerance)
{
return ReactiveTransportReturnCode::ResidualMinimized;
}
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;
}
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