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.

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] Extract biot modulus inverse list #12402

Merged
merged 20 commits into from
May 29, 2024

Conversation

rfaasse
Copy link
Contributor

@rfaasse rfaasse commented May 24, 2024

📝 Description
Similar to the related PRs of extracting pre-calculated lists with clear in- and outputs, this extracts lists for the biot modulus inverse, the degree of saturation and the derivative of the saturation.

In this PR, also some entries from the ElementVariables bucket could be removed, since dependencies became more clear and turned out to be very local.

Comment on lines 1024 to 1052
template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDerivativesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateDerivativeOfSaturation(RetentionParameters);
});

return result;
}

template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDegreesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateSaturation(RetentionParameters);
});

return result;
}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

It's a shame this is a 100% duplication with diff order. We could extract a utility, but we'd need to supply the vector of retention laws. It's still open for discussion though, I'm doubting between the two possibilities

Copy link
Contributor

Choose a reason for hiding this comment

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

I've had the exact same struggle for the calculation of the relative permeability values. I have also opted for keeping the calculation in the element itself, since you need some information that is maintained by the element. Since our intention is to ultimately remove the different branches in the element hierarchy, I'm willing to accept the technical debt that arises from the code duplication.

Copy link
Contributor

Choose a reason for hiding this comment

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

As you mentioned, this part is also duplication, and it has to be shifted to a utility.

@rfaasse rfaasse changed the title Geo/12319 extract biot modulus inverse list [GeoMechanicsApplication] Extract biot modulus inverse list May 27, 2024
@rfaasse rfaasse marked this pull request as ready for review May 27, 2024 08:48
@rfaasse rfaasse requested review from avdg81 and mnabideltares May 27, 2024 08:49
Copy link
Contributor

@mnabideltares mnabideltares left a comment

Choose a reason for hiding this comment

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

@rfaasse Thank you for this nice work and making the code more readible.
I have few comments, mainly about the duplications. In my opinion the rest of implimentation is done nicely.

const auto derivatives_of_saturation =
this->CalculateDerivativesOfSaturation(fluid_pressures, RetentionParameters);
const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli(
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, Prop);

Copy link
Contributor

Choose a reason for hiding this comment

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

This part is duplicated several times. It has to be shifted to a utilility. I think it can be done in a seperate PR.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Agreed, but as discussed it's better to do it in a separate PR (same for the other comments)

const auto derivatives_of_saturation =
CalculateDerivativesOfSaturation(fluid_pressures, RetentionParameters);
const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli(
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, rProp);

Copy link
Contributor

Choose a reason for hiding this comment

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

Same as above comment. Can be shifted to utility.

Comment on lines 1024 to 1052
template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDerivativesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateDerivativeOfSaturation(RetentionParameters);
});

return result;
}

template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDegreesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateSaturation(RetentionParameters);
});

return result;
}
Copy link
Contributor

Choose a reason for hiding this comment

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

As you mentioned, this part is also duplication, and it has to be shifted to a utility.

const auto derivatives_of_saturation =
this->CalculateDerivativesOfSaturation(fluid_pressures, RetentionParameters);
const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli(
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, Prop);

Copy link
Contributor

Choose a reason for hiding this comment

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

Duplication, as mentioned above


return result;
}

void SmallStrainUPwDiffOrderElement::CalculateMaterialStiffnessMatrix(MatrixType& rStiffnessMatrix,
Copy link
Contributor

Choose a reason for hiding this comment

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

Duplication, as mentioned

@rfaasse rfaasse self-assigned this May 28, 2024
@rfaasse rfaasse added Refactor When code is moved or rewrote keeping the same behavior GeoMechanics Issues related to the GeoMechanicsApplication labels May 28, 2024
Copy link
Contributor

@avdg81 avdg81 left a comment

Choose a reason for hiding this comment

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

It's great to see that we are almost there with precalculating lists of input parameters. This PR looks good to me; I only have several minor remarks. You may also need to resolve some merge conflicts after PR #12408 has been merged. Just so you know ;-)

@@ -1015,6 +1021,36 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLe
KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDerivativesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
Copy link
Contributor

Choose a reason for hiding this comment

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

I would suggest to make RetentionParameters a local variable of this function, to keep its scope as small as possible. To create such an object, with my recent changes, all you need is a Properties instance.

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


template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDegreesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
Copy link
Contributor

Choose a reason for hiding this comment

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

Also here, I would suggest to change RetentionParameters to a local variable rather than it being a function parameter (to reduce its scope as much as possible).

Also, please be aware that GeoTransportEquationUtilities::CalculateDegreesSaturation implements the same functionality, although it necessarily has more function parameters than your solution. I would suggest to replace that one by your new member function, since it seems more natural to have it in the element itself.

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, as per your suggestion I removed the GeoTransportEquationUtilities::CalculateDegreesSaturation and used the new function instead. Although there are then 2 versions of the function, we hope this duplication will be removed soon by merging the diff-order/non-diff order branch.

FYI @markelov208 since you separated out the utility. Please let me know if you have any reservations with this approach

Comment on lines 1024 to 1052
template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDerivativesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateDerivativeOfSaturation(RetentionParameters);
});

return result;
}

template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDegreesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateSaturation(RetentionParameters);
});

return result;
}
Copy link
Contributor

Choose a reason for hiding this comment

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

I've had the exact same struggle for the calculation of the relative permeability values. I have also opted for keeping the calculation in the element itself, since you need some information that is maintained by the element. Since our intention is to ultimately remove the different branches in the element hierarchy, I'm willing to accept the technical debt that arises from the code duplication.

@@ -1312,6 +1313,34 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi
KRATOS_CATCH("")
}

std::vector<double> SmallStrainUPwDiffOrderElement::CalculateDerivativesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
Copy link
Contributor

Choose a reason for hiding this comment

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

To comply with the Kratos Style Guide, we need to rename the first parameter:

Suggested change
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
const std::vector<double>& rFluidPressures, RetentionLaw::Parameters& RetentionParameters)

And as mentioned before, I would make RetentionParameters a local variable, to limit its scope as much as possible.

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 std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
Copy link
Contributor

Choose a reason for hiding this comment

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

You might consider to raise an error when the number of fluid pressure values does not match the number of constitutive law pointers.

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

}

std::vector<double> SmallStrainUPwDiffOrderElement::CalculateDegreesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
Copy link
Contributor

Choose a reason for hiding this comment

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

The above comments also apply to this member function.

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

Copy link
Contributor

Choose a reason for hiding this comment

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

I believe you may now also be able to remove data member DerivativeOfSaturation of nested class ElementVariables (just as you did for class UPwSmallStrainElement::ElementVariables).

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

Comment on lines 316 to 319
std::vector<double> CalculateDegreesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);
std::vector<double> CalculateDerivativesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);
Copy link
Contributor

Choose a reason for hiding this comment

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

Is there any reason why these member functions can't be const? And perhaps you may want to add attribute [[nodiscard]] to them?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, the function calls in retention law that we call in these are not const. Since that's a member, I assumed we couldn't make this one const. Turns out these members are ::Pointers, meaning we can do it.

However, I think it's quite confusing to have a const function which calls a non-const function on an object that this element has ownership over. What do you think in this case?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

But indeed we can add [[nodiscard]] regardless

@rfaasse rfaasse requested review from avdg81 and mnabideltares May 29, 2024 10:59
Copy link
Contributor Author

@rfaasse rfaasse left a comment

Choose a reason for hiding this comment

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

Processed the review comments, thank you for your thorough reviews!

std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Clang-tidy complains about this because it thinks the retention law is copied (it doesn't recognize this is a pointer). I could make it const auto&, but not sure if that makes it more clear.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Used auto to remove the clang-tidy issue

Comment on lines 316 to 319
std::vector<double> CalculateDegreesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);
std::vector<double> CalculateDerivativesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, the function calls in retention law that we call in these are not const. Since that's a member, I assumed we couldn't make this one const. Turns out these members are ::Pointers, meaning we can do it.

However, I think it's quite confusing to have a const function which calls a non-const function on an object that this element has ownership over. What do you think in this case?

Comment on lines 316 to 319
std::vector<double> CalculateDegreesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);
std::vector<double> CalculateDerivativesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);
Copy link
Contributor Author

Choose a reason for hiding this comment

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

But indeed we can add [[nodiscard]] regardless

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

Comment on lines +35 to +38
class StressStatePolicy;

template <class T, std::size_t N>
class array_1d;
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Added these as forward declare. I would like to discuss this, especially the one from core as e.g. the google style guide discourages forward declares (however, mainly focused on external libraries) https://google.github.io/styleguide/cppguide.html#Forward_Declarations. I didn't find anything in the core guidelines

@@ -1015,6 +1021,36 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLe
KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDerivativesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
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


template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDegreesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
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, as per your suggestion I removed the GeoTransportEquationUtilities::CalculateDegreesSaturation and used the new function instead. Although there are then 2 versions of the function, we hope this duplication will be removed soon by merging the diff-order/non-diff order branch.

FYI @markelov208 since you separated out the utility. Please let me know if you have any reservations with this approach

@@ -1312,6 +1313,34 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi
KRATOS_CATCH("")
}

std::vector<double> SmallStrainUPwDiffOrderElement::CalculateDerivativesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
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

}

std::vector<double> SmallStrainUPwDiffOrderElement::CalculateDegreesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
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 std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
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

Copy link
Contributor

@avdg81 avdg81 left a comment

Choose a reason for hiding this comment

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

This one looks good to go for me. I have only two minor comments, but those can be addressed later.

@@ -463,6 +463,12 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, Prop);
const auto relative_permeability_values = this->CalculateRelativePermeabilityValues(
GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector));
const auto fluid_pressures = GeoTransportEquationUtilities::CalculateFluidPressures(
Copy link
Contributor

Choose a reason for hiding this comment

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

A very minor thing, which we can also do later: if you would move the calculation of the fluid pressures one line up, you could also pass it to CalculateRelativePermeabilityValues.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Will do in the next PR

@@ -976,6 +975,12 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLe
strain_vectors, mStressVector, constitutive_matrices);
const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients(
constitutive_matrices, this->GetProperties());
const auto fluid_pressures = GeoTransportEquationUtilities::CalculateFluidPressures(
Copy link
Contributor

Choose a reason for hiding this comment

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

Also here the calculated fluid_pressures could be passed to CalculateRelativePermeabilityValues.

@rfaasse rfaasse enabled auto-merge (squash) May 29, 2024 13:12
Copy link
Contributor

@mnabideltares mnabideltares left a comment

Choose a reason for hiding this comment

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

Thank you @rfaasse for this work. It looks good to me.

@rfaasse rfaasse merged commit 1ddd4f1 into master May 29, 2024
11 checks passed
@rfaasse rfaasse deleted the geo/12319-extract-biot-modulus-inverse-list branch May 29, 2024 13:35
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
GeoMechanics Issues related to the GeoMechanicsApplication Refactor When code is moved or rewrote keeping the same behavior
Projects
None yet
3 participants