Page MenuHomec4science
Files F74582517

reactive_transport_solver.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Jul 28, 14:02

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 {
// 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
)
{
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;
// 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, retcode);
++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;
}
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();
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;
get_perfs().update = residuals.update;
return retcode;
}
ReactiveTransportReturnCode ReactiveTransportSolver::one_iteration(
VariablesBasePtr variables,
internal::ReactiveTransportResiduals& residuals
)
{
return one_iteration(variables.get(), residuals);
}
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();
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 (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

Event Timeline

c4science · Help