Homec4science
Diffusion SpecMiCP / ReactMiCP (bin/unsaturated) Browse

specmicp/tests/reactmicp/binary/unsaturated/acc_carbo_ups.cpp
85239b62bc3cbin/unsaturated

Last Change
Raw File

acc_carbo_ups.cpp
View Options

#include "reactmicp_unsaturated.hpp"
#include "physics/unsaturated_laws.hpp"
#include "utils/plugins/plugin_interface.h"
#include "utils/plugins/plugin_base.hpp"
#include "utils/compat.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() {}
~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;
//! \brief Solve the equation for the timestep
//!
//! \param var a shared_ptr to the variables
virtual StaggerReturnCode restart_timestep(
VariablesBase * const var) override;
// 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);
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};
};
class SPECMICP_DLL_PUBLIC AccCarboUpscalingFactory: public UpscalingStaggerFactory
{
public:
static AccCarboUpscalingFactory* get() {
return new AccCarboUpscalingFactory();
}
virtual std::shared_ptr<UpscalingStaggerBase> get_upscaling_stagger(
SimulationData& data,
const io::RORichYAMLNode& configuration
) override
{
return std::make_shared<UpscalingStagger>();
}
};
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(
"upscaling_factory", "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);
// 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=1; 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);
}
}
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));
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
)
{
}
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=1; 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;
}
c4science · Help