Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79242896
reactive_transport_solver.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Aug 25, 11:09
Size
5 KB
Mime Type
text/x-c
Expires
Tue, Aug 27, 11:09 (2 d)
Engine
blob
Format
Raw Data
Handle
20151015
Attached To
rSPECMICP SpecMiCP / ReactMiCP
reactive_transport_solver.cpp
View Options
#include "reactive_transport_solver.hpp"
#include "staggers_base/staggers_base.hpp"
#include "utils/log.hpp"
#include <iostream>
#include <chrono>
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)
{
std::chrono::time_point<std::chrono::system_clock> start, end;
std::chrono::duration<scalar_t> elapsed_seconds;
start = std::chrono::system_clock::now();
m_upscaling_stagger->restart_timestep(variables);
end = std::chrono::system_clock::now();
elapsed_seconds = end-start;
m_timer.upscaling_time += elapsed_seconds.count();
}
// 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
// ---------
std::chrono::time_point<std::chrono::system_clock> start, end;
std::chrono::duration<scalar_t> elapsed_seconds;
start = std::chrono::system_clock::now();
StaggerReturnCode transport_ret_code = m_transport_stagger->restart_timestep(variables);
if (transport_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::TransportFailure;
}
end = std::chrono::system_clock::now();
elapsed_seconds = end-start;
m_timer.transport_time += elapsed_seconds.count();
// Chemistry
// ---------
start = std::chrono::system_clock::now();
StaggerReturnCode chemistry_ret_code = m_chemistry_stagger->restart_timestep(variables);
if (chemistry_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::ChemistryFailure;
}
end = std::chrono::system_clock::now();
elapsed_seconds = end-start;
m_timer.chemistry_time += elapsed_seconds.count();
// Upscaling0
// ---------
if (get_options().implicit_upscaling)
{
start = std::chrono::system_clock::now();
StaggerReturnCode upscaling_ret_code = m_upscaling_stagger->restart_timestep(variables);
if (upscaling_ret_code <= StaggerReturnCode::NotConvergedYet)
{
return ReactiveTransportReturnCode::UpscalingFailure;
}
end = std::chrono::system_clock::now();
elapsed_seconds = end-start;
m_timer.upscaling_time += elapsed_seconds.count();
}
// 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
Event Timeline
Log In to Comment