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] Implement Calculator-based workflow for Pw line element #12778

Merged
merged 52 commits into from
Nov 6, 2024

Conversation

rfaasse
Copy link
Contributor

@rfaasse rfaasse commented Oct 24, 2024

📝 Description
This PR implements a calculator-based workflow for the contributions added to the left hand side and right hand side in the Pw line element. In short, we now provide a list of contribution types to the element, which then creates the corresponding calculator and its input.

The advantage of this is that the actual calculation of the contributions is delegated to the calculators, meaning we can quite easily determine what contributions we would like in an element, without creating new elements in the class hierarchy.

If this workflow turns out to be satisfactory for the Pw line element, we can start taking steps to roll this out to our other (U)Pw elements as well as use it for the piping element.

@rfaasse rfaasse added the Refactor When code is moved or rewrote keeping the same behavior label Oct 25, 2024
@rfaasse rfaasse marked this pull request as ready for review October 25, 2024 07:53
Copy link
Contributor

@WPK4FEM WPK4FEM left a comment

Choose a reason for hiding this comment

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

I like the general applicability of the idea presented here. Some aspects of the implementation I have trouble understanding. Please explain a bit further to me.

namespace Kratos
{

enum class CalculationContribution { Permeability, Compressibility };
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 extended with mass, damping, stiffness, UPwCoupling and PwUCoupling when working on the UPw elements.
A similar thought process can perhaps be followed for the geometrical nonlinear contributions

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The GeometricStiffness matrix would be a very suitable one indeed, hadn't thought of that but when having a quick look, it seems we can probably remove the updated lagrangian elements from the hierarchy if we have a 'GeometricStiffnessCalculator' 🎉

For the mass and damping, I think it's a good idea to add the calculators, however since adding these terms is the responsibility of the scheme, this will require a bit more discussion.


virtual Matrix LHSContribution() = 0;
virtual Vector RHSContribution() = 0;
virtual std::pair<Matrix, Vector> CalculateLeftAndRightHandSide() = 0;
Copy link
Contributor

Choose a reason for hiding this comment

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

On first sight this looks a bit like CalculateLocalSystem,
Maybe include Contribution in the name ( CalculateLeftAndRightHandSideContirbution of CalculateLocalSystemContribution ).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I renamed it to LocalSystemContribution, to be in line with the LHS/RHS names, let me know whether this is clear enough!


Matrix CompressibilityCalculator::LHSContribution(const Matrix& rCompressibilityMatrix) const
{
return rCompressibilityMatrix * mInputProvider.GetDtPressureCoefficient();
Copy link
Contributor

Choose a reason for hiding this comment

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

The DtPressureCoefficient is determined by the time integration scheme. I wonder if we should do this multipication inside the calculator of leave this outside to avoid the mingling of functionality in the calculator.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

As discussed replaced it with the more general 'GetMatrixScalarFactor'

Comment on lines 375 to 378
static auto GetDtPressureCoefficientLambda(const ProcessInfo& rCurrentProcessInfo)
{
return [&rCurrentProcessInfo]() { return rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT]; };
}
Copy link
Contributor

Choose a reason for hiding this comment

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

This information belongs to the time integration scheme.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

As we discussed, I'll change this to GetMatrixScalarFactor, to be filled in by the scheme via the current process info object passed to the element.

Let me know if this is better

Comment on lines 72 to 89
if (Dimension == 1) {
BoundedMatrix<double, 1, 1> result_1d;
FillPermeabilityMatrix(result_1d, Prop);
return result_1d;
}
if (Dimension == 2) {
BoundedMatrix<double, 2, 2> result_2d;
FillPermeabilityMatrix(result_2d, Prop);
return result_2d;
}
if (Dimension == 3) {
BoundedMatrix<double, 3, 3> result_3d;
FillPermeabilityMatrix(result_3d, Prop);
return result_3d;
}

KRATOS_ERROR << "Dimension " << Dimension << " is not supported" << std::endl;
}
Copy link
Contributor

Choose a reason for hiding this comment

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

Looks like a case statement.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I had some trouble with the switch for some reason, maybe @avdg81 can help me out next week with this

Copy link
Contributor

@markelov208 markelov208 left a comment

Choose a reason for hiding this comment

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

Hi Richard, 'a small step to make major changes'.
I've got an impression that transientPwPipeElement is behind transientPwELement in terms of the functions' style. If this is true, shall it be another PR to make them consistent? Thank you.

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 is a very nice extension towards more generic handling of contributions to left hand side matrices and right hand side vectors. It will help us tremendously to reduce the explosion of element classes due to different combinations of features. Thank you so much for all the time and effort that you have put into this.
I have several suggestions, but overall this already looks quite good.

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.

I think I processed all comments, there might be a few improvements still possible, let's discuss next week!

Comment on lines 38 to 43
: mGetElementProperties(std::move(rElementProperties)),
mGetRetentionLaws(std::move(rRetentionLaws)),
mGetNContainer(std::move(rNContainer)),
mGetIntegrationCoefficients(std::move(rIntegrationCoefficients)),
mGetDtPressureCoefficient(std::move(DtPressureCoefficient)),
mGetNodalValues(std::move(rNodalValuesOfDtWaterPressure))
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Because the members are getter functions to be called when needed

Comment on lines 72 to 89
if (Dimension == 1) {
BoundedMatrix<double, 1, 1> result_1d;
FillPermeabilityMatrix(result_1d, Prop);
return result_1d;
}
if (Dimension == 2) {
BoundedMatrix<double, 2, 2> result_2d;
FillPermeabilityMatrix(result_2d, Prop);
return result_2d;
}
if (Dimension == 3) {
BoundedMatrix<double, 3, 3> result_3d;
FillPermeabilityMatrix(result_3d, Prop);
return result_3d;
}

KRATOS_ERROR << "Dimension " << Dimension << " is not supported" << std::endl;
}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

I had some trouble with the switch for some reason, maybe @avdg81 can help me out next week with this

WPK4FEM
WPK4FEM previously approved these changes Nov 4, 2024
Copy link
Contributor

@WPK4FEM WPK4FEM left a comment

Choose a reason for hiding this comment

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

I'm happy with this.

GeoElementUtilities::FillPermeabilityMatrix(r_properties, local_dimension);

const auto number_of_nodes = shape_function_gradients[0].size1();
auto result = Matrix{ZeroMatrix{number_of_nodes, number_of_nodes}};
Copy link
Contributor

Choose a reason for hiding this comment

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

This dimension is correct, as there is one D.O.F. per node. The more natural dimension is the number of Pw D.O.F. multiplied with itself.

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 is almost ready I think. I noticed there was one function renamed which you probably did not intend to... Please check whether my observation is correct. Apart from that I only have one more minor suggestion and a response to one of your questions.

Copy link
Contributor

@WPK4FEM WPK4FEM left a comment

Choose a reason for hiding this comment

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

This is a very scalable implementation that can be extended to our other elements.
Think it is ready to go.

@rfaasse rfaasse requested review from avdg81 and markelov208 November 5, 2024 14:29
Copy link
Contributor

@markelov208 markelov208 left a comment

Choose a reason for hiding this comment

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

Hi Richard, thank you very much for the introduction of calculators.

@rfaasse rfaasse merged commit 154939a into master Nov 6, 2024
11 checks passed
@rfaasse rfaasse deleted the geo/add-extendible-calculators-to-pw-line-element branch November 6, 2024 08:04
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
Development

Successfully merging this pull request may close these issues.

[GeoMechanicsApplication] Implement 'calculator-based' workflow for the Pw Line element
4 participants