Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78286893
equilibrium_stagger.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, Aug 19, 15:09
Size
5 KB
Mime Type
text/x-c
Expires
Wed, Aug 21, 15:09 (2 d)
Engine
blob
Format
Raw Data
Handle
20018881
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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;
ERROR << "Total concentration : \n" << constraints.total_concentrations;
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
Event Timeline
Log In to Comment