Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/reactmicp/systems/saturated_react/equilibrium_stagger.cpp
f7c3cd7d6b9breactmicp/unsaturated

Last Change
Raw File

equilibrium_stagger.cpp
View Options

#include "equilibrium_stagger.hpp"
#include "variables.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "utils/log.hpp"
#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
#ifdef SPECMICP_USE_OPENMP
#include <omp.h>
#endif // SPECMICP_USE_OPENMP
namespace specmicp {
namespace reactmicp {
namespace systems {
namespace satdiff {
// EquilibriumStagger::
//! \brief Initialize the stagger at the beginning of an iteration
void EquilibriumStagger::initialize_timestep(scalar_t dt, VariablesBasePtr var)
{
m_dt = dt;
SaturatedVariablesPtr true_var = cast_var_from_base(var);
// Initialize velocity using values from previous timestep
for (index_t node=0; node<true_var->nb_nodes(); ++node)
{
if (true_var->is_fixed_composition(node)) continue;
scalar_t alpha = 1.0;
for (index_t component=1; component<true_var->nb_component(); ++component)
{
alpha = std::max(alpha,
dt*true_var->solid_concentration(node, component, true_var->chemistry_rate())
/true_var->solid_concentration(node, component, true_var->displacement())
);
}
auto solid_velocity = true_var->velocity().segment(
true_var->offset_node(node)+true_var->offset_solid_concentration(),
true_var->nb_component());
auto solid_chemistry_rate = true_var->chemistry_rate().segment(
true_var->offset_node(node)+true_var->offset_solid_concentration(),
true_var->nb_component());
solid_velocity = 1/alpha * solid_chemistry_rate;
}
}
//! \brief Solve the equation for the timestep
solver::StaggerReturnCode EquilibriumStagger::restart_timestep(VariablesBasePtr var)
{
SaturatedVariablesPtr true_var = cast_var_from_base(var);
int failed_chemistry = 0;
#ifdef SPECMICP_USE_OPENMP
#pragma omp parallel
{
#pragma omp for schedule(runtime)
for (index_t node=0; node<true_var->nb_nodes(); ++node)
{
if (true_var->is_fixed_composition(node)) continue;
int retcode = solve_one_node(node, true_var);
if (retcode > 0)
{
#pragma omp atomic
++failed_chemistry;
}
}
}
#else
{
for (index_t node=0; node<true_var->nb_nodes(); ++node)
{
if (true_var->is_fixed_composition(node)) continue;
int retcode = solve_one_node(node, true_var);
if (retcode > 0)
{
++failed_chemistry;
break;
}
}
}
#endif // SPECMICP_USE_OPENMP
if (failed_chemistry > 0)
return solver::StaggerReturnCode::UnknownError;
return solver::StaggerReturnCode::ResidualMinimized;
}
//!
int EquilibriumStagger::solve_one_node(index_t node, SaturatedVariablesPtr true_var)
{
AdimensionalSystemConstraints constraints(get_constraints(node));
constraints.total_concentrations = true_var->total_concentrations(node);
AdimensionalSystemSolver adim_solver(true_var->get_database(),
constraints,
true_var->equilibrium_solution(node),
m_options);
Vector variables(true_var->equilibrium_solution(node).main_variables);
micpsolver::MiCPPerformance perf = adim_solver.solve(variables);
micpsolver::MiCPSolverReturnCode retcode = perf.return_code;
if (retcode <= micpsolver::MiCPSolverReturnCode::NotConvergedYet)
{
ERROR << "Failed to solve chemistry problem at node " << node
<< ", return code = " << static_cast<int>(retcode)
<< ", residual = " << perf.current_residual;
return 1;
}
true_var->equilibrium_solution(node) = adim_solver.get_raw_solution(variables);
AdimensionalSystemSolutionExtractor extractor(true_var->equilibrium_solution(node),
true_var->get_database(),
m_options.units_set);
for (index_t component=0; component<true_var->nb_component(); ++component)
{
const scalar_t c_aq = extractor.density_water()*extractor.total_aqueous_concentration(component);
true_var->aqueous_concentration(node, component, true_var->displacement()) = c_aq;
const scalar_t c_aq_0 = true_var->aqueous_concentration(node, component, true_var->predictor());
const scalar_t vel_aq = (c_aq - c_aq_0)/m_dt;
true_var->aqueous_concentration(node, component, true_var->velocity()) = vel_aq;
const scalar_t c_sol = extractor.total_immobile_concentration(component);
true_var->solid_concentration(node, component, true_var->displacement()) = c_sol;
const scalar_t c_sol_0 = true_var->solid_concentration(node, component, true_var->predictor());
const scalar_t vel_sol = (c_sol - c_sol_0)/m_dt;
true_var->solid_concentration(node, component, true_var->velocity()) = vel_sol;
true_var->solid_concentration(node, component, true_var->chemistry_rate()) = vel_sol;
}
return 0;
}
} // end namespace satdiff
} // end namespace systems
} // end namespace reactmicp
} // end namespace specmicp
c4science · Help