Page MenuHomec4science
Files F63514466

equilibrium_stagger.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, May 20, 16:23

equilibrium_stagger.cpp
View Options

/*------------------------------------------------------------------------------
Copyright (c)2015 F. Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.
3. Neither the name of the copyright holder nor the names of its contributors
may be used to endorse or promote products derived from this software without
specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
------------------------------------------------------------------------------*/
#include "equilibrium_stagger.hpp"
#include "variables.hpp"
#include "equilibrium_constraints.hpp"
#include "../../solver/staggers_base/stagger_structs.hpp"
#include "../../../dfpm/meshes/mesh1d.hpp"
#include "../../../database/data_container.hpp"
#include "../../../utils/compat.hpp"
#include "../../../specmicp/adimensional/adimensional_system_solution.hpp"
#include "../../../specmicp/adimensional/adimensional_system_solution_extractor.hpp"
#include "../../../specmicp/adimensional/adimensional_system_solver.hpp"
#include "../../../utils/log.hpp"
namespace specmicp {
namespace reactmicp {
namespace systems {
namespace unsaturated {
// ============================================
//
// Declaration of implementation details
//
// ============================================
struct EquilibriumStagger::EquilibriumStaggerImpl
{
scalar_t m_dt {-1};
mesh::Mesh1DPtr m_mesh;
database::RawDatabasePtr m_database;
std::shared_ptr<EquilibriumConstraints> m_constraints;
EquilibriumStaggerImpl(
std::shared_ptr<UnsaturatedVariables> variables,
std::shared_ptr<EquilibriumConstraints> constraints
):
m_mesh(variables->get_mesh()),
m_database(variables->get_database()),
m_constraints(constraints)
{}
scalar_t dt() {return m_dt;}
int compute_one_node(index_t node, UnsaturatedVariables* vars);
void compute_total_concentrations(
index_t node,
Vector& total_concentrations,
UnsaturatedVariables* vars
);
void analyse_solution(
index_t node,
AdimensionalSystemSolution& solution,
UnsaturatedVariables* vars
);
void initialize_timestep_one_node(index_t node, UnsaturatedVariables* vars);
void initialize(UnsaturatedVariables* vars);
void initialize_timestep(scalar_t dt, UnsaturatedVariables* vars);
StaggerReturnCode restart_timestep(UnsaturatedVariables* vars);
units::UnitsSet& get_units() {
return m_constraints->get_options().units_set;
}
};
// ============================================
//
// Implementation
//
// ============================================
EquilibriumStagger::EquilibriumStagger(
std::shared_ptr<UnsaturatedVariables> variables,
std::shared_ptr<EquilibriumConstraints> constraints
):
m_impl(make_unique<EquilibriumStaggerImpl>(variables, constraints))
{
}
EquilibriumStagger::~EquilibriumStagger() = default;
//! \brief Initialize the stagger at the beginning of the computation
//!
//! \param var a shared_ptr to the variables
void EquilibriumStagger::initialize(VariablesBasePtr var)
{
UnsaturatedVariables* true_vars =
static_cast<UnsaturatedVariables*>(var.get());
m_impl->initialize(true_vars);
}
//! \brief Initialize the stagger at the beginning of an iteration
//!
//! This is where the predictor can be saved, the first trivial iteration done, ...
//!
//! \param dt the duration of the timestep
//! \param var a shared_ptr to the variables
void EquilibriumStagger::initialize_timestep(scalar_t dt, VariablesBasePtr var)
{
UnsaturatedVariables* true_vars =
static_cast<UnsaturatedVariables*>(var.get());
m_impl->initialize_timestep(dt, true_vars);
}
//! \brief Solve the equation for the timestep
//!
//! \param var a shared_ptr to the variables
EquilibriumStagger::StaggerReturnCode
EquilibriumStagger::restart_timestep(VariablesBasePtr var)
{
UnsaturatedVariables* true_vars =
static_cast<UnsaturatedVariables*>(var.get());
return m_impl->restart_timestep(true_vars);
}
int
EquilibriumStagger::EquilibriumStaggerImpl::compute_one_node(
index_t node,
UnsaturatedVariables* vars
)
{
AdimensionalSystemConstraints constraints = m_constraints->get_constraints(node);
compute_total_concentrations(node, constraints.total_concentrations, vars);
{
user_model_saturation_f callable = vars->get_vapor_pressure_model();
water_partial_pressure_f pv_model = [callable, node](scalar_t sat){
return callable(node, sat);
};
constraints.set_water_partial_pressure_model(pv_model);
}
// We use the current value of the saturation to avoid converging
// to the previous solution
// This is necessary because the update to the saturation is very small
AdimensionalSystemSolution& previous_solution = vars->get_adim_solution(node);
if (previous_solution.is_valid)
{
const auto porosity = vars->get_porosity()(node);
const auto saturation = vars->get_liquid_saturation()(node);
previous_solution.main_variables(0) = porosity*saturation;
}
//
AdimensionalSystemSolver the_solver(vars->get_database(),
constraints,
previous_solution,
m_constraints->get_options()
);
Vector x = previous_solution.main_variables;
micpsolver::MiCPPerformance perf = the_solver.solve(x);
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
ERROR << "Failed to solve equilibrium at node " << node
<< ". \n"
<< "Return code : " << (int) perf.return_code;
return -1;
}
AdimensionalSystemSolution solution = the_solver.get_raw_solution(x);
analyse_solution(node, solution, vars);
vars->set_adim_solution(node, solution);
return 0;
}
void
EquilibriumStagger::EquilibriumStaggerImpl::compute_total_concentrations(
index_t node,
Vector& total_concentrations,
UnsaturatedVariables* vars
)
{
total_concentrations.resize(m_database->nb_component());
const scalar_t saturation = vars->get_liquid_saturation()(node);
const scalar_t porosity = vars->get_porosity()(node);
const scalar_t rt = vars->get_rt();
// water
{
const scalar_t ctilde_w = vars->get_water_aqueous_concentration()(node);
const scalar_t cbar_w = vars->get_solid_concentration(0)(node);
scalar_t c_w = saturation*porosity*ctilde_w + cbar_w;
if (vars->component_has_gas(0))
{
const scalar_t pv_w = vars->get_pressure_main_variables(0)(node);
c_w += (1.0-saturation)*porosity*pv_w/rt;
}
total_concentrations(0) = c_w;
}
total_concentrations(1) = 0.0;
// aqueous components
for (index_t component: m_database->range_aqueous_component())
{
const scalar_t ctilde_i = vars->get_aqueous_concentration(component)(node);
const scalar_t cbar_i = vars->get_solid_concentration(component)(node);
scalar_t c_i = saturation*porosity*ctilde_i + cbar_i;
if (vars->component_has_gas(component))
{
const scalar_t pv_i = vars->get_pressure_main_variables(component)(node);
c_i += (1.0-saturation)*porosity*pv_i/rt;
}
total_concentrations(component) = c_i;
}
}
void EquilibriumStagger::EquilibriumStaggerImpl::analyse_solution(
index_t node,
AdimensionalSystemSolution& solution,
UnsaturatedVariables* vars
)
{
AdimensionalSystemSolutionExtractor extr(solution,
m_database,
get_units());
const scalar_t saturation = extr.saturation_water();
const scalar_t sat_g = 1-saturation;
const scalar_t rho_l = extr.density_water();
const scalar_t porosity = extr.porosity();
const scalar_t rt = vars->get_rt();
// porosity
SecondaryTransientVariable& porosity_vars = vars->get_porosity();
const scalar_t porosity_0 = porosity_vars.predictor(node);
const scalar_t porosity_vel = (porosity - porosity_0)/m_dt;
porosity_vars.variable(node) = porosity;
porosity_vars.velocity(node) = porosity_vel;
// water
// saturation
MainVariable& satvars = vars->get_liquid_saturation();
const scalar_t saturation_0 = satvars.predictor(node);
const scalar_t sat_vel = (saturation - saturation_0)/m_dt;
satvars.variable(node) = saturation;
satvars.velocity(node) = sat_vel;
{
// total aqueous concentration
SecondaryTransientVariable& cwtilde_vars = vars->get_water_aqueous_concentration();
const scalar_t cwtilde = rho_l*extr.total_aqueous_concentration(0);
const scalar_t cwtilde_0 = cwtilde_vars.predictor(node);
const scalar_t cwtilde_vel = (cwtilde - cwtilde_0)/m_dt;
cwtilde_vars.variable(node) = cwtilde;
cwtilde_vars.velocity(node) = cwtilde_vel;
// total solid concentration
MainVariable& solid_vars = vars->get_solid_concentration(0);
const scalar_t cwbar = extr.total_solid_concentration(0);
const scalar_t cwbar_0 = solid_vars.predictor(node);
const scalar_t cwbar_vel = (cwbar - cwbar_0)/m_dt;
solid_vars.variable(node) = cwbar;
solid_vars.velocity(node) = cwbar_vel;
solid_vars.chemistry_rate(node) = - cwbar_vel;
// vapor pressure
if (vars->component_has_gas(0))
{
MainVariable& pres_vars = vars->get_pressure_main_variables(0);
const scalar_t pv = vars->get_vapor_pressure_model()(node, saturation);
const scalar_t pv_0 = pres_vars.predictor(node);
const scalar_t pv_vel = (pv - pv_0)/m_dt;
pres_vars.variable(node) = pv;
pres_vars.velocity(node) = pv_vel;
const scalar_t transient = ( (pv*porosity*sat_g)
- (pv_0*porosity_0*(1-saturation_0))
)/
(rt*m_dt);
const scalar_t pv_chem_rate = - transient + pres_vars.transport_fluxes(node);
pres_vars.chemistry_rate(node) = pv_chem_rate;
}
}
// aqueous components
for (index_t component: m_database->range_aqueous_component())
{
// total aqueous concentration
MainVariable& aqconc = vars->get_aqueous_concentration(component);
const scalar_t ctildei = rho_l*extr.total_aqueous_concentration(component);
const scalar_t ctildei_0 = aqconc.predictor(node);
const scalar_t ctildei_vel = (ctildei - ctildei_0) / m_dt;
aqconc.variable(node) = ctildei;
aqconc.velocity(node) = ctildei_vel;
// total solid concentration
MainVariable& solconc = vars->get_solid_concentration(component);
const scalar_t cbari = extr.total_solid_concentration(component);
const scalar_t cbari_0 = solconc.predictor(node);
const scalar_t cbari_vel = (cbari - cbari_0) / m_dt;
solconc.variable(node) = cbari;
solconc.velocity(node) = cbari_vel;
solconc.chemistry_rate(node) = - cbari_vel;
// partial pressure
if (vars->component_has_gas(component))
{
MainVariable& pres_vars = vars->get_pressure_main_variables(component);
index_t id_g = vars->get_id_gas(component);
const scalar_t pi = extr.fugacity_gas(id_g)*vars->get_total_pressure();
const scalar_t pi_0 = pres_vars.predictor(node);
const scalar_t pi_vel = (pi - pi_0) / m_dt;
pres_vars.variable(node) = pi;
pres_vars.velocity(node) = pi_vel;
const scalar_t transient = (
(pi*porosity*sat_g)
- (pi_0*porosity_0*(1.0-saturation_0))
) / (rt * m_dt);
const scalar_t pi_chem_rate = - transient + pres_vars.transport_fluxes(node);
pres_vars.chemistry_rate(node) = pi_chem_rate;
}
}
}
void
EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep_one_node(
index_t node,
UnsaturatedVariables* vars
)
{
SecondaryTransientVariable& porosity = vars->get_porosity();
porosity.predictor = porosity.variable;
porosity.velocity.setZero();
}
EquilibriumStagger::StaggerReturnCode
EquilibriumStagger::EquilibriumStaggerImpl::restart_timestep(
UnsaturatedVariables* vars
)
{
int sum_retcode = 0;
for (index_t node: m_mesh->range_nodes())
{
sum_retcode += compute_one_node(node, vars);
}
if (sum_retcode != 0)
{
return StaggerReturnCode::UnknownError;
}
return StaggerReturnCode::ResidualMinimized;
}
void
EquilibriumStagger::EquilibriumStaggerImpl::initialize_timestep(
scalar_t dt,
UnsaturatedVariables* vars
)
{
m_dt = dt;
// initialize
//water
{
MainVariable& cbar_w = vars->get_solid_concentration(0);
cbar_w.predictor = cbar_w.variable;
cbar_w.velocity.setZero();
}
// aqueous component
for (index_t component: m_database->range_aqueous_component())
{
MainVariable& cbar_i = vars->get_solid_concentration(component);
cbar_i.predictor = cbar_i.variable;
cbar_i.velocity.setZero();
}
for (auto node: m_mesh->range_nodes())
{
initialize_timestep_one_node(node, vars);
}
}
void
EquilibriumStagger::EquilibriumStaggerImpl::initialize(
UnsaturatedVariables* vars
)
{
// do nothing by default
}
} //end namespace unsaturated
} //end namespace systems
} //end namespace reactmicp
} //end namespace specmicp

Event Timeline

c4science · Help