Page MenuHomec4science
Files F65495407

acc_carbo_ups.cpp
No OneTemporary
Actions

File Metadata

Created
Tue, Jun 4, 04:48

acc_carbo_ups.cpp
View Options

#include "reactmicp/reactmicp_unsaturated.hpp"
#include "specmicp_common/physics/unsaturated_laws.hpp"
#include "specmicp_common/plugins/plugin_interface.h"
#include "specmicp_common/plugins/plugin_base.hpp"
#include "specmicp_common/compat.hpp"
#include "specmicp_common/io/safe_config.hpp"
#include <iostream>
using namespace specmicp;
using namespace specmicp::reactmicp::solver;
using namespace specmicp::reactmicp::systems::unsaturated;
class SPECMICP_DLL_PUBLIC UpscalingStagger:
public specmicp::reactmicp::solver::UpscalingStaggerBase
{
public:
UpscalingStagger(
SimulationData& data,
io::YAMLConfigHandle&& conf
);
~UpscalingStagger() {}
//! \brief Initialize the stagger at the beginning of the computation
//!
//! \param var a shared_ptr to the variables
virtual void initialize(VariablesBase * const var) override;
//! \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
virtual void initialize_timestep(
scalar_t dt,
VariablesBase * const var
) override;
//! \param var a shared_ptr to the variables
virtual StaggerReturnCode restart_timestep(
VariablesBase * const var) override;
virtual void register_chemistry_stagger(
std::shared_ptr<ChemistryStaggerBase> chem_stagger
) override {
EquilibriumStagger* stagger = static_cast<EquilibriumStagger*>(chem_stagger.get());
m_chem_opts = stagger->get_options();
}
// User models
scalar_t capillary_pressure(index_t node, scalar_t saturation);
scalar_t vapor_pressure(index_t node, scalar_t saturation);
scalar_t relative_liquid_diffusivity(index_t node, scalar_t saturation);
scalar_t relative_liquid_permeability(index_t node, scalar_t saturation);
//! \brief Solve the equation for the timestep
//!
scalar_t relative_gas_diffusivity(index_t node, scalar_t saturation);
scalar_t intrinsic_permeability(scalar_t porosity);
scalar_t intrinsic_liquid_diffusivity(scalar_t porosity);
scalar_t resistance_gas_diffusivity(scalar_t porosity);
private:
scalar_t m_a {37.5479e6};
scalar_t m_b {2.1684};
scalar_t m_alpha {7.27};
scalar_t m_beta {0.440};
scalar_t m_pv_sat {3170};
scalar_t m_exp_r_l_d {3.00};
scalar_t m_exp_i_g_d {2.74};
scalar_t m_exp_r_g_d {4.20};
scalar_t m_exp_r_l_k {0.4191};
scalar_t m_k_0 {1e-21*1e4};
scalar_t m_d_0 {2.3e-13*1e4};
scalar_t m_phi_0 {0.2};
std::shared_ptr<EquilibriumOptionsVector> m_chem_opts;
std::shared_ptr<BoundaryConditions> m_bcs;
};
UpscalingStagger::UpscalingStagger(
SimulationData& data,
io::YAMLConfigHandle&& conf
)
{
m_bcs = data.bcs; // register boundary conditions
auto param = conf.get_section("parameters");
m_a = param.get_required_attribute<scalar_t>("capillary_a");
m_b = param.get_required_attribute<scalar_t>("capillary_b");
m_pv_sat = param.get_required_attribute<scalar_t>("pv_sat");
m_exp_r_l_d = param.get_required_attribute<scalar_t>("exp_r_l_d");
m_exp_i_g_d = param.get_required_attribute<scalar_t>("exp_i_g_d");
m_exp_r_g_d = param.get_required_attribute<scalar_t>("exp_r_g_d");
m_exp_r_l_k = param.get_required_attribute<scalar_t>("exp_r_l_k");
m_k_0 = param.get_required_attribute<scalar_t>("k_0");
m_d_0 = param.get_required_attribute<scalar_t>("d_0");
m_phi_0 = param.get_required_attribute<scalar_t>("phi_0");
}
class SPECMICP_DLL_PUBLIC AccCarboUpscalingFactory: public UpscalingStaggerFactory
{
public:
static AccCarboUpscalingFactory* get() {
return new AccCarboUpscalingFactory();
}
virtual std::shared_ptr<UpscalingStaggerBase> get_upscaling_stagger(
SimulationData& data,
io::YAMLConfigHandle&& configuration
) override
{
//auto bcs = data.bcs;
//bcs->set_fixed_flux_liquid_dof(1, 0, 2.2e-8*1e4);
return std::make_shared<UpscalingStagger>(data, std::move(configuration));
}
};
class SPECMICP_DLL_PUBLIC AccCarboUpscalingPlugin:
public specmicp::plugins::PluginBase
{
public:
AccCarboUpscalingPlugin():
PluginBase()
{
set_api_version(0, 1, 0);
}
bool initialize(const specmicp::plugins::PluginManagerServices& services) override
{
auto retcode = services.register_object(
"reactmicp/unsaturated/upscaling", "acc_carbo_ups",
&AccCarboUpscalingFactory::get
);
if (not retcode) {
return false;
}
return true;
}
};
SPECMICP_PLUGIN(AccCarboUpscalingPlugin)
// implementation
void UpscalingStagger::initialize(VariablesBase * const var)
{
UnsaturatedVariables * const true_var = static_cast<UnsaturatedVariables * const>(var);
m_bcs->add_fixed_flux_gas_dof(0, 0);
// set relative models
true_var->set_capillary_pressure_model(
[this](index_t x, scalar_t sat){
return this->capillary_pressure(x, sat);
});
true_var->set_relative_liquid_permeability_model(
[this](index_t x, scalar_t sat){
return this->relative_liquid_permeability(x, sat);
});
true_var->set_relative_liquid_diffusivity_model(
[this](index_t x, scalar_t sat){
return this->relative_liquid_diffusivity(x, sat);
});
true_var->set_relative_gas_diffusivity_model(
[this](index_t x, scalar_t sat){
return this->relative_gas_diffusivity(x, sat);
});
true_var->set_vapor_pressure_model(
[this](index_t x, scalar_t sat){
return this->vapor_pressure(x, sat);
});
// initialization
SecondaryTransientVariable& porosity = true_var->get_porosity();
SecondaryVariable& permeability = true_var->get_liquid_permeability();
SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity();
SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity();
//true_var->set_relative_variables();
//gas_diffusivity(0) = resistance_gas_diffusivity(1.0);
//true_var->get_relative_gas_diffusivity()(0) = relative_gas_diffusivity(0, 1.0);
for (index_t node=0; node<true_var->get_mesh()->nb_nodes(); ++node)
{
scalar_t phi = porosity(node);
permeability(node) = intrinsic_permeability(phi);
liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi);
gas_diffusivity(node) = resistance_gas_diffusivity(phi);
true_var->set_relative_variables(node);
}
// true_var->set_relative_variables(0);
}
scalar_t UpscalingStagger::capillary_pressure(index_t node, scalar_t saturation)
{
if (saturation == 0.0) return 0.0;
return laws::capillary_pressure_BaroghelBouny(saturation, m_a, m_b);
}
scalar_t UpscalingStagger::vapor_pressure(index_t node, scalar_t saturation)
{
scalar_t pv;
if (saturation == 0.0) pv = 0.0;
else if (saturation >= 1.0) pv = m_pv_sat;
else pv = m_pv_sat*laws::invert_kelvin_equation(
capillary_pressure(node, saturation));
//if (node == 1)
//std::cout << "saturation " << saturation << " - pv " << pv << "\n";
return pv;
}
scalar_t UpscalingStagger::relative_liquid_diffusivity(
index_t node,
scalar_t saturation
)
{
if (saturation == 0.0) return 0.0;
return std::pow(saturation, m_exp_r_l_d);
}
scalar_t UpscalingStagger::relative_liquid_permeability(
index_t node,
scalar_t saturation
)
{
if (saturation == 0.0) return 0.0;
return laws::relative_liquid_permeability_Mualem(saturation, 1.0/m_b);
}
scalar_t UpscalingStagger::relative_gas_diffusivity(index_t node, scalar_t saturation)
{
if (saturation == 0.0) return 1.0;
return std::pow(1.0 - saturation, m_exp_r_g_d);
}
scalar_t UpscalingStagger::intrinsic_permeability(scalar_t porosity)
{
scalar_t tmp_1 = porosity / m_phi_0;
scalar_t tmp_2 = (1.0 - m_phi_0) / (1.0 - porosity);
tmp_1 = std::pow(tmp_1, 3);
tmp_2 = std::pow(tmp_2, 2);
return m_k_0 * tmp_1 * tmp_2;
}
scalar_t UpscalingStagger::intrinsic_liquid_diffusivity(scalar_t porosity)
{
return m_d_0 * std::exp(-9.95*porosity);
}
scalar_t UpscalingStagger::resistance_gas_diffusivity(scalar_t porosity)
{
return std::pow(porosity, m_exp_i_g_d);
}
void UpscalingStagger::initialize_timestep(
scalar_t dt,
VariablesBase * const var
)
{
UnsaturatedVariables* true_var = static_cast<UnsaturatedVariables*>(var);
scalar_t flux = - 1.22e-7 * true_var->get_rt() * (
true_var->get_pressure_main_variables(0).variable(0) - 2000);
m_bcs->set_fixed_flux_gas_dof(0, 0, flux);
// const auto& the_mesh = true_var->get_mesh();
// const auto& raw_data = true_var->get_database();
// auto& liq_sat = true_var->get_liquid_saturation();
// auto& ca_conc = true_var->get_solid_concentration(raw_data->get_id_component("Ca[2+]"));
// auto& co2_conc = true_var->get_solid_concentration(raw_data->get_id_component("CO2"));
// for (index_t node: the_mesh->range_nodes()) {
// if (not liq_sat(node) > 0) continue;
// if (co2_conc(node) >= ca_conc(node)) {
// WARNING << "Node " << node << "is carbonated";
// m_chem_opts->set(node, "carbonated");
// }
// }
}
StaggerReturnCode UpscalingStagger::restart_timestep(VariablesBase * const var)
{
UnsaturatedVariables* true_var = static_cast<UnsaturatedVariables*>(var);
SecondaryTransientVariable& porosity = true_var->get_porosity();
SecondaryVariable& permeability = true_var->get_liquid_permeability();
SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity();
SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity();
for (index_t node=0; node<true_var->get_mesh()->nb_nodes(); ++node)
{
scalar_t phi = porosity(node);
permeability(node) = intrinsic_permeability(phi);
liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi);
gas_diffusivity(node) = resistance_gas_diffusivity(phi);
true_var->set_relative_variables(node);
}
return StaggerReturnCode::ResidualMinimized;
}

Event Timeline

c4science · Help