Homec4science
Diffusion SpecMiCP / ReactMiCP (yaml/safe) Browse

specmicp/src/reactmicp/systems/saturated_diffusion/reactive_transport_solver.cpp
e5484eae1665yaml/safe

Last Change
Raw File

reactive_transport_solver.cpp
View Options

//#define EIGEN_DEFAULT_IO_FORMAT IOFormat(8)
#include "reactive_transport_solver.hpp"
#include "specmicp/reduced_system_solver.hpp"
#include <omp.h>
namespace specmicp {
namespace reactmicp {
namespace internals {
//! \brief Store residual information
struct SIAResiduals
{
int nb_iterations;
double transport_residual_0;
double transport_residual;
double transport_update;
Eigen::VectorXd speciation_residual;
Eigen::VectorXd speciation_update;
SIAResiduals(int nb_node):
nb_iterations(0),
speciation_residual(Eigen::VectorXd::Zero(nb_node)),
speciation_update(Eigen::VectorXd::Zero(nb_node))
{}
};
} // end namespace internals
void SIASaturatedReactiveTransportSolver::apply_boundary_conditions(SIABoundaryConditions bcs)
{
assert(bcs.bs_types.size() == (unsigned int) m_mesh->nnodes());
assert(bcs.initial_states.size() == (unsigned int) m_mesh->nnodes());
std::vector<int> list_fixed_nodes;
for (int node: m_mesh->range_nodes())
{
switch (bcs.bs_types[node]) {
case BoundaryConditionType::FixedComposition:
m_bcs(node) = 1;
list_fixed_nodes.push_back(node);
break;
default:
m_bcs(node) = 0;
break;
}
m_variables.update_composition(node, bcs.list_initial_states[bcs.initial_states[node]]);
m_transport_program.number_equations(list_fixed_nodes);
}
}
double SIASaturatedReactiveTransportSolver::residuals_transport()
{
Eigen::VectorXd residuals;
m_transport_solver.compute_residuals(m_variables.total_mobile_concentrations(), residuals);
//std::cout << "Residuals : " << std::endl << residuals << std::endl;
return residuals.norm();
}
SIASaturatedReactiveTransportSolverReturnCode convert_dfpm_retcode(dfpmsolver::ParabolicDriverReturnCode retcode)
{
switch (retcode) {
case dfpmsolver::ParabolicDriverReturnCode::MaxIterations:
return SIASaturatedReactiveTransportSolverReturnCode::MaxIterationsInTransport;
case dfpmsolver::ParabolicDriverReturnCode::StationaryPoint:
return SIASaturatedReactiveTransportSolverReturnCode::StationaryPointsInTransport;
default:
return SIASaturatedReactiveTransportSolverReturnCode::ErrorInTransport;
}
}
SIASaturatedReactiveTransportSolverReturnCode SIASaturatedReactiveTransportSolver::solve_timestep(double dt)
{
retcode_t return_code = retcode_t::NotConvergedYet;
internals::SIAResiduals residuals_storage(m_mesh->nnodes());
m_flag_compute_speciation = Eigen::VectorXi::Ones(m_mesh->nnodes());
// initialization transport program and transport solver
Eigen::VectorXd transport_variable(m_variables.total_mobile_concentrations());
m_transport_solver.get_velocity() = Eigen::VectorXd::Zero(transport_variable.rows());
m_transport_program.external_flux().setZero();
//
residuals_storage.transport_residual_0 = residuals_transport();
m_minerals_start_iteration = m_variables.speciation_variables().block(
m_database->nb_component, 0, m_database->nb_mineral, m_mesh->nnodes());
// solve transport
// ---------------
dfpmsolver::ParabolicDriverReturnCode retcode = m_transport_solver.solve_timestep(dt, transport_variable);
if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized)
{
ERROR << "Failed to solve 1st transport problem; return code " << (int) retcode;
return convert_dfpm_retcode(retcode);
}
residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual;
residuals_storage.transport_update = m_transport_solver.get_perfs().current_update;
// update
m_variables.total_mobile_concentrations() = transport_variable;
// Speciation
// ----------
restart_speciation_iterations(dt, residuals_storage);
// Convergence
// -----------
return_code = check_convergence(residuals_storage);
residuals_storage.nb_iterations +=1;
// If SNIA, we stop
if (get_options().max_iterations == 1) return retcode_t::Success;
// SIA algorithm
// -------------
while (return_code == retcode_t::NotConvergedYet)
{
//std::cout << m_flag_compute_speciation.sum() << std::endl;
transport_variable = m_variables.total_mobile_concentrations();
retcode = m_transport_solver.restart_timestep(transport_variable);
if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized)
{
if (retcode == dfpmsolver::ParabolicDriverReturnCode::MaxIterations and
m_transport_solver.get_perfs().current_residual/residuals_storage.transport_residual_0
< get_options().good_enough_transport_residual_tolerance)
{
restart_speciation_iterations(dt, residuals_storage);
return_code = retcode_t::Success;
break;
}
ERROR << "Failed to solve transport problem; return code " << (int) retcode;
return convert_dfpm_retcode(retcode);
}
residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual;
residuals_storage.transport_update = m_transport_solver.get_perfs().current_update;
m_variables.total_mobile_concentrations() = transport_variable;
residuals_storage.nb_iterations +=1;
restart_speciation_iterations(dt, residuals_storage);
return_code = check_convergence(residuals_storage);
}
std::cout << "End fixed point : return code " << (int) return_code
<< " - nb iter : " <<residuals_storage.nb_iterations
<< std::endl;
return return_code;
}
void SIASaturatedReactiveTransportSolver::restart_transport_iterations(
internals::SIAResiduals& residuals_storage
)
{
Eigen::VectorXd transport_variable(m_variables.total_mobile_concentrations());
dfpmsolver::ParabolicDriverReturnCode retcode = m_transport_solver.restart_timestep(transport_variable);
if (retcode != dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized)
{
ERROR << "Failed to solve transport problem; return code " << (int) retcode;
throw std::runtime_error("Failed to solve transport program");
}
residuals_storage.transport_residual = m_transport_solver.get_perfs().current_residual;
residuals_storage.transport_update = m_transport_solver.get_perfs().current_update;
}
void SIASaturatedReactiveTransportSolver::restart_speciation_iterations(
double dt,
internals::SIAResiduals& residuals_storage)
{
#ifdef _OPENMP
#pragma omp parallel
{
#pragma omp for schedule(static)
for (int node=0; node<m_mesh->nnodes(); ++node)
#else
for (int node: m_mesh->range_nodes())
#endif
{
if (m_bcs(node) == 0 and m_flag_compute_speciation(node))
{
solve_speciation_node(dt, node, residuals_storage);
}
}
#ifdef _OPENMP
}
#endif
}
void SIASaturatedReactiveTransportSolver::solve_speciation_node(double dt, int node, internals::SIAResiduals& residuals_storage)
{
Eigen::VectorXd totconc;
m_variables.nodal_update_total_amount(node, totconc);
//std::cout << "Node : " << node << " - Total concentration : " << std::endl << totconc << std::endl;
ReducedSystemSolver solver(m_database, totconc);
solver.get_options().solver_options.max_iter = 25;
Eigen::VectorXd variables = m_variables.speciation_variables(node);
micpsolver::MiCPPerformance perf = solver.solve(variables);
if (perf.return_code != micpsolver::MiCPSolverReturnCode::ResidualMinimized)
{
WARNING << "Failed to solve speciation solver for node " << node
<< "; return code = " << (int) perf.return_code;
}
residuals_storage.speciation_residual(node) = perf.current_residual;
residuals_storage.speciation_update(node) = perf.current_update;
if (perf.current_update < get_options().threshold_update_disable_speciation) m_flag_compute_speciation(node) = 0;
EquilibriumState solution = solver.get_solution(variables);
m_variables.update_composition(node, solution);
// mineral
const double denom = 1e-3*m_mesh->cell_volume(node)*dt;
for (int component: m_database->range_aqueous_component())
{
double sum = 0.0;
for (int mineral: m_database->range_mineral())
{
if (m_database->nu_mineral(mineral, component) == 0.0) continue;
sum -= m_database->nu_mineral(mineral, component)*(
m_variables.mineral_amount(node, mineral)
-m_minerals_start_iteration(mineral, node)
);
}
m_transport_program.external_flux(node, component) = sum/denom;
}
}
SIASaturatedReactiveTransportSolverReturnCode SIASaturatedReactiveTransportSolver::check_convergence(
internals::SIAResiduals& residuals_storage
)
{
// compute residuals of transport
const double norm_transport = residuals_transport();
retcode_t return_code = retcode_t::NotConvergedYet;
if (norm_transport/residuals_storage.transport_residual_0 < get_options().residual_tolerance)
{
return_code = retcode_t::Success;
}
else if (residuals_storage.transport_update < get_options().threshold_update_transport)
{
return_code = retcode_t::Success;
}
else if (residuals_storage.nb_iterations >= get_options().max_iterations)
{
WARNING << "Maximum number of iterations (" << get_options().max_iterations
<< ") reached during fixed-point iterations";
return_code = retcode_t::MaxIterations;
}
return return_code;
}
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help