Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F74723391
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
Mon, Jul 29, 08:24
Size
9 KB
Mime Type
text/x-c
Expires
Wed, Jul 31, 08:24 (2 d)
Engine
blob
Format
Raw Data
Handle
19429726
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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
{
index_t nb_iterations;
scalar_t transport_residual_0;
scalar_t transport_residual;
scalar_t transport_update;
Vector speciation_residual;
Vector speciation_update;
SIAResiduals(index_t nb_node):
nb_iterations(0),
speciation_residual(Vector::Zero(nb_node)),
speciation_update(Vector::Zero(nb_node))
{}
};
} // end namespace internals
void SIASaturatedReactiveTransportSolver::apply_boundary_conditions(SIABoundaryConditions bcs)
{
specmicp_assert(bcs.bs_types.size() == (unsigned int) m_mesh->nnodes());
specmicp_assert(bcs.initial_states.size() == (unsigned int) m_mesh->nnodes());
std::vector<index_t> list_fixed_nodes;
for (index_t 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);
}
}
scalar_t SIASaturatedReactiveTransportSolver::residuals_transport()
{
Eigen::VectorXd residuals;
m_transport_solver.compute_residuals(m_variables.total_mobile_concentrations(), residuals);
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(scalar_t 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(
scalar_t dt,
internals::SIAResiduals& residuals_storage
)
{
#ifdef _OPENMP
#pragma omp parallel
{
#pragma omp for schedule(static)
for (index_t node=0; node<m_mesh->nnodes(); ++node)
#else
for (index_t 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(
scalar_t dt,
index_t node,
internals::SIAResiduals& residuals_storage
)
{
Eigen::VectorXd totconc;
m_variables.nodal_update_total_amount(node, totconc);
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 scalar_t denom = 1e-3*m_mesh->cell_volume(node)*dt;
for (index_t component: m_database->range_aqueous_component())
{
scalar_t sum = 0.0;
for (index_t 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 scalar_t 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
Event Timeline
Log In to Comment