Commit 56383707 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

[MPL] Replace liquid_saturation_rate.

Now using the previous timestep's values.
parent 206e22fd
......@@ -55,13 +55,14 @@ void SaturationDependentSwelling::checkScale() const
PropertyDataType SaturationDependentSwelling::value(
VariableArray const& variable_array,
VariableArray const& variable_array_prev,
ParameterLib::SpatialPosition const& pos, double const /*t*/,
double const dt) const
{
auto const S_L = std::get<double>(
variable_array[static_cast<int>(Variable::liquid_saturation)]);
auto const S_L_dot = std::get<double>(
variable_array[static_cast<int>(Variable::liquid_saturation_rate)]);
auto const S_L_prev = std::get<double>(
variable_array_prev[static_cast<int>(Variable::liquid_saturation)]);
Eigen::Matrix<double, 3, 3> const e =
local_coordinate_system_ == nullptr
......@@ -76,8 +77,6 @@ PropertyDataType SaturationDependentSwelling::value(
return delta_sigma_sw; // still being zero.
}
double const S_L_prev = S_L - S_L_dot * dt;
double const S_eff = std::clamp((S_L - S_min_) / (S_max_ - S_min_), 0., 1.);
double const S_eff_prev =
std::clamp((S_L_prev - S_min_) / (S_max_ - S_min_), 0., 1.);
......@@ -100,9 +99,10 @@ PropertyDataType SaturationDependentSwelling::value(
}
PropertyDataType SaturationDependentSwelling::dValue(
VariableArray const& variable_array, Variable const primary_variable,
VariableArray const& variable_array,
VariableArray const& variable_array_prev, Variable const primary_variable,
ParameterLib::SpatialPosition const& pos, double const /*t*/,
double const dt) const
double const /*dt*/) const
{
(void)primary_variable;
assert((primary_variable == Variable::liquid_saturation) &&
......@@ -111,8 +111,8 @@ PropertyDataType SaturationDependentSwelling::dValue(
auto const S_L = std::get<double>(
variable_array[static_cast<int>(Variable::liquid_saturation)]);
auto const S_L_dot = std::get<double>(
variable_array[static_cast<int>(Variable::liquid_saturation_rate)]);
auto const S_L_prev = std::get<double>(
variable_array_prev[static_cast<int>(Variable::liquid_saturation)]);
Eigen::Matrix<double, 3, 3> const e =
local_coordinate_system_ == nullptr
......@@ -127,8 +127,6 @@ PropertyDataType SaturationDependentSwelling::dValue(
return delta_sigma_sw; // still being zero.
}
double const S_L_prev = S_L - S_L_dot * dt;
double const S_eff = std::clamp((S_L - S_min_) / (S_max_ - S_min_), 0., 1.);
double const S_eff_prev =
std::clamp((S_L_prev - S_min_) / (S_max_ - S_min_), 0., 1.);
......
......@@ -56,9 +56,11 @@ public:
void checkScale() const override;
PropertyDataType value(VariableArray const& variable_array,
VariableArray const& variable_array_prev,
ParameterLib::SpatialPosition const& pos,
double const t, double const dt) const override;
PropertyDataType dValue(VariableArray const& variable_array,
VariableArray const& variable_array_prev,
Variable const variable,
ParameterLib::SpatialPosition const& pos,
double const t, double const dt) const override;
......
......@@ -49,7 +49,6 @@ enum class Variable : int
effective_pore_pressure_rate,
grain_compressibility,
liquid_saturation,
liquid_saturation_rate,
phase_pressure,
porosity,
solid_grain_pressure,
......
......@@ -362,8 +362,8 @@ void RichardsMechanicsLocalAssembler<
S_L = medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
variables[static_cast<int>(MPL::Variable::liquid_saturation_rate)] =
(S_L - S_L_prev) / dt;
variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
S_L_prev;
double const dS_L_dp_cap =
medium->property(MPL::PropertyType::saturation)
......@@ -417,8 +417,8 @@ void RichardsMechanicsLocalAssembler<
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
.template value<DimMatrix>(variables, x_position, t,
dt));
.template value<DimMatrix>(
variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
......@@ -692,8 +692,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
S_L = medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
variables[static_cast<int>(MPL::Variable::liquid_saturation_rate)] =
(S_L - S_L_prev) / dt;
variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
S_L_prev;
double const dS_L_dp_cap =
medium->property(MPL::PropertyType::saturation)
......@@ -748,8 +748,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
.template value<DimMatrix>(variables, x_position, t,
dt));
.template value<DimMatrix>(
variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
variables[static_cast<int>(
......@@ -859,8 +859,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
.template dValue<DimMatrix>(
variables, MPL::Variable::liquid_saturation,
x_position, t, dt));
variables, variables_prev,
MPL::Variable::liquid_saturation, x_position, t,
dt));
local_Jac
.template block<displacement_size, pressure_size>(
displacement_index, pressure_index)
......@@ -1412,8 +1413,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
S_L = medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
variables[static_cast<int>(MPL::Variable::liquid_saturation_rate)] =
(S_L - S_L_prev) / dt;
variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
S_L_prev;
auto const chi = [medium, x_position, t, dt](double const S_L) {
MPL::VariableArray variables;
......@@ -1476,8 +1477,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
.template value<DimMatrix>(variables, x_position, t,
dt));
.template value<DimMatrix>(
variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
variables[static_cast<int>(
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment