Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[GeoMechanicsApplication] Simplified fluid flux for pw element #12338

Merged
merged 8 commits into from
May 6, 2024
Merged

[GeoMechanicsApplication] Simplified fluid flux for pw element #12338

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
3564929
Copied fluid flux vector calculation to pw element
rfaasse May 2, 2024
4c02764
Extract `CalculateFluidFluxes` which takes permeability update factor…
rfaasse May 2, 2024
4a5056e
Use `CalculateFluidFluxes` in U_Pw_small_strain_element.cpp
rfaasse May 2, 2024
14f8280
Cleaning up and simplifying CalculateFluidFluxes
rfaasse May 2, 2024
e3d770c
Minor formatting commit
rfaasse May 2, 2024
f4e155a
Revert "Cleaning up and simplifying CalculateFluidFluxes"
rfaasse May 3, 2024
0c56f55
Merge remote-tracking branch 'origin/master' into geo/12319_simplifie…
rfaasse May 3, 2024
c25296e
Processing review comments
rfaasse May 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 49 additions & 28 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -620,42 +620,21 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(
if (rVariable == FLUID_FLUX_VECTOR) {
ElementVariables Variables;
this->InitializeElementVariables(Variables, rCurrentProcessInfo);

const PropertiesType& rProp = this->GetProperties();

array_1d<double, TDim> GradPressureTerm;
array_1d<double, TDim> FluidFlux;

// Create general parameters of retention law
RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);


std::vector<double> permeability_update_factors;
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
// Compute Np, GradNpT, B and StrainVector
this->CalculateKinematics(Variables, GPoint);

// Compute strain, needed for update permeability
this->CalculateStrain(Variables, GPoint);
this->CalculatePermeabilityUpdateFactor(Variables);
permeability_update_factors.push_back(Variables.PermeabilityUpdateFactor);
}

GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint);
Variables.FluidPressure = CalculateFluidPressure(Variables);

RetentionParameters.SetFluidPressure(Variables.FluidPressure);

Variables.RelativePermeability =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters);

noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector);

noalias(GradPressureTerm) +=
PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration;

noalias(FluidFlux) = PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse *
Variables.RelativePermeability * Variables.PermeabilityUpdateFactor *
prod(Variables.PermeabilityMatrix, GradPressureTerm);

GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], FluidFlux);
const auto fluid_fluxes = CalculateFluidFluxes(permeability_update_factors, rCurrentProcessInfo);
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], fluid_fluxes[GPoint]);
}
} else {
if (rOutput.size() != mConstitutiveLawVector.size())
Expand Down Expand Up @@ -1112,6 +1091,48 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeBiotCoefficients(ElementV
KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
std::vector<array_1d<double, TDim>> UPwSmallStrainElement<TDim, TNumNodes>::CalculateFluidFluxes(
const std::vector<double>& permeability_update_factors, const ProcessInfo& rCurrentProcessInfo)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To be in line with the Kratos Style Guide, we should rename the first parameter:

Suggested change
const std::vector<double>& permeability_update_factors, const ProcessInfo& rCurrentProcessInfo)
const std::vector<double>& rPermeabilityUpdateFactors, const ProcessInfo& rCurrentProcessInfo)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

{
const GeometryType& rGeom = this->GetGeometry();
const IndexType NumGPoints = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod());

std::vector<array_1d<double, TDim>> FluidFluxes;
ElementVariables Variables;
this->InitializeElementVariables(Variables, rCurrentProcessInfo);

const PropertiesType& rProp = this->GetProperties();

array_1d<double, TDim> GradPressureTerm;

// Create general parameters of retention law
RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);

for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
this->CalculateKinematics(Variables, GPoint);
Variables.PermeabilityUpdateFactor = permeability_update_factors[GPoint];

GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint);
Variables.FluidPressure = CalculateFluidPressure(Variables);
RetentionParameters.SetFluidPressure(Variables.FluidPressure);

Variables.RelativePermeability =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters);

noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector);

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

lines 1124 and 1126 are constructing the contents of the same variable. This interspacing contradicts that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed blank line

noalias(GradPressureTerm) += PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration;

FluidFluxes.push_back(PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse *
Variables.RelativePermeability * Variables.PermeabilityUpdateFactor *
prod(Variables.PermeabilityMatrix, GradPressureTerm));
}

return FluidFluxes;
}

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainElement<TDim, TNumNodes>::CalculatePermeabilityUpdateFactor(ElementVariables& rVariables)
{
Expand Down
5 changes: 4 additions & 1 deletion applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -240,7 +240,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa

virtual void CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables);

virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, const ElementVariables& rVariables);
virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix,
const ElementVariables& rVariables);

virtual void CalculateAndAddRHS(VectorType& rRightHandSideVector, ElementVariables& rVariables, unsigned int GPoint);

Expand Down Expand Up @@ -283,6 +284,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa

void InitializeProperties(ElementVariables& rVariables);
double CalculateFluidPressure(const ElementVariables& rVariables);
std::vector<array_1d<double, TDim>> CalculateFluidFluxes(const std::vector<double>& permeability_update_factors,
const ProcessInfo& rCurrentProcessInfo);

void CalculateRetentionResponse(ElementVariables& rVariables,
RetentionLaw::Parameters& rRetentionParameters,
Expand Down
11 changes: 10 additions & 1 deletion applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -383,8 +383,17 @@ void TransientPwElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const Var
{
KRATOS_TRY

const GeometryType& rGeom = this->GetGeometry();
const IndexType NumGPoints = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod());
if (rOutput.size() != NumGPoints) rOutput.resize(NumGPoints);

if (rVariable == FLUID_FLUX_VECTOR) {
BaseType::CalculateOnIntegrationPoints(rVariable, rOutput, rCurrentProcessInfo);
std::vector<double> permeability_update_factors(NumGPoints, 1.0);
const auto fluid_fluxes = this->CalculateFluidFluxes(permeability_update_factors, rCurrentProcessInfo);

for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], fluid_fluxes[GPoint]);
}
} else {
if (rOutput.size() != mRetentionLawVector.size())
rOutput.resize(mRetentionLawVector.size());
Expand Down
Loading