Homec4science
Diffusion SpecMiCP / ReactMiCP (master) Browse

specmicp/src/reactmicp/solver/reactive_transport_solver.cpp
921a58bdaf61master

Last Change
Raw File

reactive_transport_solver.cpp
View Options

/* =============================================================================
Copyright (c) 2014 - 2016
F. Georget <fabieng@princeton.edu> Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
3. Neither the name of the copyright holder nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
============================================================================= */
#include "reactive_transport_solver.hpp"
#include "staggers_base/staggers_base.hpp"
#include "specmicp_common/log.hpp"
#include "specmicp_common/timer.hpp"
namespace specmicp {
namespace reactmicp {
namespace solver {
namespace internal {
//! \brief Residuals info to check convergence in coupling algorithm
// This contains internal information for the reactive transport solver
// It is used to check the convergence
struct ReactiveTransportResiduals
{
index_t nb_iterations; //!< The number of iteratons
scalar_t transport_residual_0; //!< The residuals at the beginning of the timestep
scalar_t transport_residuals; //!< The current residuals
scalar_t update; //!< Variables 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
)
{
return solve_timestep(timestep, variables.get());
}
ReactiveTransportReturnCode ReactiveTransportSolver::solve_timestep(
scalar_t timestep,
VariablesBase* variables
)
{
// start the timer
Timer tot_timer;
tot_timer.start();
// initialization
internal::ReactiveTransportResiduals residuals;
reset_perfs();
get_perfs().timestep = timestep;
m_transport_stagger->initialize_timestep(timestep, variables);
m_chemistry_stagger->initialize_timestep(timestep, variables);
m_upscaling_stagger->initialize_timestep(timestep, variables);
// get convergence criterion from transport stagger
residuals.transport_residual_0 = m_transport_stagger->get_residual_0(variables);
// run the first iteration
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, retcode);
++residuals.nb_iterations;
// if sequential non-iterative algorithm
if (get_options().is_snia())
{
goto end;
}
// else run the other iterations
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;
}
retcode = check_convergence(variables, residuals, retcode);
++residuals.nb_iterations;
}
// wrapup
end:
// upscaling, if needed
if (not get_options().implicit_upscaling)
{
Timer timer;
timer.start();
StaggerReturnCode upscaling_ret_code = m_upscaling_stagger->restart_timestep(variables);
if (upscaling_ret_code <= StaggerReturnCode::NotConvergedYet)
{
retcode = ReactiveTransportReturnCode::UpscalingFailure;
goto set_return;
}
else if (upscaling_ret_code == StaggerReturnCode::UserTermination)
retcode = ReactiveTransportReturnCode::UserTermination;
timer.stop();
const scalar_t utime = timer.elapsed_time();
m_timer.upscaling_time += utime;
get_perfs().upscaling_time += utime;
}
// record performance and return code
set_return:
if (retcode <= ReactiveTransportReturnCode::NotConvergedYet)
{
ERROR << "Current Residual : "
<< residuals.transport_residuals/residuals.transport_residual_0;
m_transport_stagger->print_debug_information(variables);
}
else if (retcode == ReactiveTransportReturnCode::GoodEnough)
{
WARNING << "Good enough convergence was triggered";
m_transport_stagger->print_debug_information(variables);
}
// perfs informations
get_perfs().nb_iterations = residuals.nb_iterations;
get_perfs().return_code = retcode;
get_perfs().residuals = residuals.transport_residuals/residuals.transport_residual_0;
get_perfs().update = residuals.update;
// at the end, stop timer, and record time
tot_timer.stop();
get_perfs().total_time = tot_timer.elapsed_time();
// timer for staggers are already registered
return retcode;
}
ReactiveTransportReturnCode ReactiveTransportSolver::one_iteration(
VariablesBasePtr variables,
internal::ReactiveTransportResiduals& residuals
)
{
return one_iteration(variables.get(), residuals);
}
// Set the clock
void ReactiveTransportSolver::set_clock(scalar_t clock_time)
{
m_upscaling_stagger->set_clock(clock_time);
}
ReactiveTransportReturnCode ReactiveTransportSolver::one_iteration(
VariablesBase* variables,
internal::ReactiveTransportResiduals& residuals
)
{
Timer timer;
bool bypass = false;
bool user_termination = false;
// Transport
// ---------
timer.start();
StaggerReturnCode transport_ret_code = m_transport_stagger->restart_timestep(variables);
if (transport_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::TransportFailure;
}
else 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;
// Upscaling
// ---------
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;
}
else if (upscaling_ret_code == StaggerReturnCode::UserTermination)
{
user_termination = true;
}
timer.stop();
const scalar_t utime = timer.elapsed_time();
get_perfs().upscaling_time += utime;
m_timer.upscaling_time += utime;
}
// Final residuals
// ---------------
residuals.transport_residuals = m_transport_stagger->get_residual(variables);
residuals.update = m_transport_stagger->get_update(variables);
if (user_termination)
return ReactiveTransportReturnCode::UserTermination;
else if (bypass)
return ReactiveTransportReturnCode::TransportBypass;
else
return ReactiveTransportReturnCode::NotConvergedYet;
}
// Check the convergence
ReactiveTransportReturnCode ReactiveTransportSolver::check_convergence(
VariablesBase* _,
const internal::ReactiveTransportResiduals& residuals,
ReactiveTransportReturnCode iteration_return_code
)
{
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 (iteration_return_code == ReactiveTransportReturnCode::TransportBypass)
{
return ReactiveTransportReturnCode::TransportBypass;
}
else if (iteration_return_code == ReactiveTransportReturnCode::UserTermination)
{
return ReactiveTransportReturnCode::UserTermination;
}
// 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