[GeoMechanicsApplication] Implement Calculator-based workflow for Pw line element #12778
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,19 @@ | ||
| // KRATOS___ | ||
| // // ) ) | ||
| // // ___ ___ | ||
| // // ____ //___) ) // ) ) | ||
| // // / / // // / / | ||
| // ((____/ / ((____ ((___/ / MECHANICS | ||
| // | ||
| // License: geo_mechanics_application/license.txt | ||
| // | ||
| // Main authors: Richard Faasse | ||
| // | ||
| #pragma once | ||
|
|
||
| namespace Kratos | ||
| { | ||
|
|
||
| enum class CalculationContribution { Permeability, Compressibility }; | ||
|
|
||
| } | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,113 @@ | ||
| // KRATOS___ | ||
| // // ) ) | ||
| // // ___ ___ | ||
| // // ____ //___) ) // ) ) | ||
| // // / / // // / / | ||
| // ((____/ / ((____ ((___/ / MECHANICS | ||
| // | ||
| // License: geo_mechanics_application/license.txt | ||
| // | ||
| // Main authors: Richard Faasse | ||
| // | ||
| #include "compressibility_calculator.h" | ||
| #include "custom_utilities/transport_equation_utilities.hpp" | ||
| #include "geo_mechanics_application_variables.h" | ||
|
|
||
| namespace Kratos | ||
| { | ||
|
|
||
| CompressibilityCalculator::CompressibilityCalculator(InputProvider rInputProvider) | ||
| : mInputProvider(std::move(rInputProvider)) | ||
| { | ||
| } | ||
|
|
||
| Matrix CompressibilityCalculator::LHSContribution() | ||
| { | ||
| return LHSContribution(CalculateCompressibilityMatrix()); | ||
| } | ||
|
|
||
| Vector CompressibilityCalculator::RHSContribution() | ||
| { | ||
| return RHSContribution(CalculateCompressibilityMatrix()); | ||
| } | ||
|
|
||
| Vector CompressibilityCalculator::RHSContribution(const Matrix& rCompressibilityMatrix) const | ||
| { | ||
| return -prod(rCompressibilityMatrix, mInputProvider.GetNodalValues(DT_WATER_PRESSURE)); | ||
| } | ||
|
|
||
| Matrix CompressibilityCalculator::LHSContribution(const Matrix& rCompressibilityMatrix) const | ||
| { | ||
| return mInputProvider.GetMatrixScalarFactor()* rCompressibilityMatrix; | ||
| } | ||
|
|
||
| std::pair<Matrix, Vector> CompressibilityCalculator::LocalSystemContribution() | ||
| { | ||
| const auto compressibility_matrix = CalculateCompressibilityMatrix(); | ||
| return {LHSContribution(compressibility_matrix), RHSContribution(compressibility_matrix)}; | ||
| } | ||
|
|
||
| Matrix CompressibilityCalculator::CalculateCompressibilityMatrix() const | ||
| { | ||
| const auto& r_N_container = mInputProvider.GetNContainer(); | ||
| const auto& integration_coefficients = mInputProvider.GetIntegrationCoefficients(); | ||
| auto result = Matrix{ZeroMatrix{r_N_container.size2(), r_N_container.size2()}}; | ||
| for (unsigned int integration_point_index = 0; | ||
rfaasse marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| integration_point_index < integration_coefficients.size(); ++integration_point_index) { | ||
| const auto N = Vector{row(r_N_container, integration_point_index)}; | ||
| result += GeoTransportEquationUtilities::CalculateCompressibilityMatrix( | ||
| N, CalculateBiotModulusInverse(mInputProvider.GetRetentionLaws()[integration_point_index]), | ||
| integration_coefficients[integration_point_index]); | ||
| } | ||
| return result; | ||
| } | ||
|
|
||
| double CompressibilityCalculator::CalculateBiotModulusInverse(const RetentionLaw::Pointer& rRetentionLaw) const | ||
| { | ||
| const auto& r_properties = mInputProvider.GetElementProperties(); | ||
| const double biot_coefficient = r_properties[BIOT_COEFFICIENT]; | ||
|
|
||
| double bulk_fluid = TINY; | ||
| if (!r_properties[IGNORE_UNDRAINED]) { | ||
markelov208 marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| bulk_fluid = r_properties[BULK_MODULUS_FLUID]; | ||
| } | ||
| double result = (biot_coefficient - r_properties[POROSITY]) / r_properties[BULK_MODULUS_SOLID] + | ||
| r_properties[POROSITY] / bulk_fluid; | ||
WPK4FEM marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
|
||
| RetentionLaw::Parameters retention_parameters(r_properties); | ||
| result *= rRetentionLaw->CalculateSaturation(retention_parameters); | ||
| result -= rRetentionLaw->CalculateDerivativeOfSaturation(retention_parameters) * r_properties[POROSITY]; | ||
| return result; | ||
| } | ||
|
|
||
| const Properties& CompressibilityCalculator::InputProvider::GetElementProperties() const | ||
| { | ||
| return mGetElementProperties(); | ||
| } | ||
|
|
||
| const std::vector<RetentionLaw::Pointer>& CompressibilityCalculator::InputProvider::GetRetentionLaws() const | ||
| { | ||
| return mGetRetentionLaws(); | ||
| } | ||
|
|
||
| const Matrix& CompressibilityCalculator::InputProvider::GetNContainer() const | ||
| { | ||
| return mGetNContainer(); | ||
| } | ||
|
|
||
| Vector CompressibilityCalculator::InputProvider::GetIntegrationCoefficients() const | ||
| { | ||
| return mGetIntegrationCoefficients(); | ||
| } | ||
|
|
||
| double CompressibilityCalculator::InputProvider::GetMatrixScalarFactor() const | ||
| { | ||
| return mGetMatrixScalarFactor(); | ||
| } | ||
|
|
||
| Vector CompressibilityCalculator::InputProvider::GetNodalValues(const Variable<double>& rVariable) const | ||
| { | ||
| return mGetNodalValues(rVariable); | ||
| } | ||
|
|
||
| } // namespace Kratos | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,78 @@ | ||
| // KRATOS___ | ||
| // // ) ) | ||
| // // ___ ___ | ||
| // // ____ //___) ) // ) ) | ||
| // // / / // // / / | ||
| // ((____/ / ((____ ((___/ / MECHANICS | ||
| // | ||
| // License: geo_mechanics_application/license.txt | ||
| // | ||
| // Main authors: Richard Faasse | ||
| // | ||
|
|
||
| #pragma once | ||
|
|
||
| #include "contribution_calculator.h" | ||
| #include "custom_retention/retention_law.h" | ||
| #include "includes/properties.h" | ||
| #include "includes/ublas_interface.h" | ||
|
|
||
| #include <utility> | ||
| #include <vector> | ||
|
|
||
| namespace Kratos | ||
| { | ||
|
|
||
| class CompressibilityCalculator : public ContributionCalculator | ||
| { | ||
| public: | ||
| class InputProvider | ||
| { | ||
| public: | ||
| InputProvider(std::function<const Properties&()> GetElementProperties, | ||
| std::function<const std::vector<RetentionLaw::Pointer>&()> GetRetentionLaws, | ||
| std::function<const Matrix&()> GetNContainer, | ||
| std::function<Vector()> GetIntegrationCoefficients, | ||
| std::function<double()> GetMatrixScalarFactor, | ||
| std::function<Vector(const Variable<double>&)> GetNodalValuesOf) | ||
| : mGetElementProperties(std::move(GetElementProperties)), | ||
| mGetRetentionLaws(std::move(GetRetentionLaws)), | ||
| mGetNContainer(std::move(GetNContainer)), | ||
| mGetIntegrationCoefficients(std::move(GetIntegrationCoefficients)), | ||
| mGetMatrixScalarFactor(std::move(GetMatrixScalarFactor)), | ||
| mGetNodalValues(std::move(GetNodalValuesOf)) | ||
| { | ||
| } | ||
|
|
||
| [[nodiscard]] const Properties& GetElementProperties() const; | ||
| [[nodiscard]] const std::vector<RetentionLaw::Pointer>& GetRetentionLaws() const; | ||
| [[nodiscard]] const Matrix& GetNContainer() const; | ||
| [[nodiscard]] Vector GetIntegrationCoefficients() const; | ||
| [[nodiscard]] double GetMatrixScalarFactor() const; | ||
| [[nodiscard]] Vector GetNodalValues(const Variable<double>& rVariable) const; | ||
|
|
||
| private: | ||
| std::function<const Properties&()> mGetElementProperties; | ||
| std::function<const std::vector<RetentionLaw::Pointer>&()> mGetRetentionLaws; | ||
| std::function<const Matrix&()> mGetNContainer; | ||
| std::function<Vector()> mGetIntegrationCoefficients; | ||
| std::function<double()> mGetMatrixScalarFactor; | ||
| std::function<Vector(const Variable<double>&)> mGetNodalValues; | ||
| }; | ||
|
|
||
| explicit CompressibilityCalculator(InputProvider rInputProvider); | ||
|
|
||
| Matrix LHSContribution() override; | ||
| Vector RHSContribution() override; | ||
| std::pair<Matrix, Vector> LocalSystemContribution() override; | ||
|
|
||
| private: | ||
| [[nodiscard]] Matrix CalculateCompressibilityMatrix() const; | ||
| [[nodiscard]] double CalculateBiotModulusInverse(const RetentionLaw::Pointer& rRetentionLaw) const; | ||
| [[nodiscard]] Vector RHSContribution(const Matrix& rCompressibilityMatrix) const; | ||
| [[nodiscard]] Matrix LHSContribution(const Matrix& rCompressibilityMatrix) const; | ||
|
|
||
| InputProvider mInputProvider; | ||
| }; | ||
|
|
||
| } // namespace Kratos |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,29 @@ | ||
| // KRATOS___ | ||
| // // ) ) | ||
| // // ___ ___ | ||
| // // ____ //___) ) // ) ) | ||
| // // / / // // / / | ||
| // ((____/ / ((____ ((___/ / MECHANICS | ||
| // | ||
| // License: geo_mechanics_application/license.txt | ||
| // | ||
| // Main authors: Richard Faasse | ||
| // | ||
| #pragma once | ||
|
|
||
| #include "includes/ublas_interface.h" | ||
|
|
||
| namespace Kratos | ||
| { | ||
|
|
||
| class ContributionCalculator | ||
| { | ||
| public: | ||
| virtual ~ContributionCalculator() = default; | ||
|
|
||
| virtual Matrix LHSContribution() = 0; | ||
| virtual Vector RHSContribution() = 0; | ||
| virtual std::pair<Matrix, Vector> LocalSystemContribution() = 0; | ||
| }; | ||
|
|
||
| } // namespace Kratos |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,93 @@ | ||
| // KRATOS___ | ||
| // // ) ) | ||
| // // ___ ___ | ||
| // // ____ //___) ) // ) ) | ||
| // // / / // // / / | ||
| // ((____/ / ((____ ((___/ / MECHANICS | ||
| // | ||
| // License: geo_mechanics_application/license.txt | ||
| // | ||
| // Main authors: Richard Faasse | ||
| // | ||
| #include "permeability_calculator.h" | ||
| #include "custom_retention/retention_law.h" | ||
| #include "custom_utilities/element_utilities.hpp" | ||
| #include "custom_utilities/transport_equation_utilities.hpp" | ||
| #include "includes/cfd_variables.h" | ||
|
|
||
| namespace Kratos | ||
| { | ||
|
|
||
| PermeabilityCalculator::PermeabilityCalculator(InputProvider InputProvider) | ||
| : mInputProvider{std::move(InputProvider)} | ||
| { | ||
| } | ||
|
|
||
| Matrix PermeabilityCalculator::LHSContribution() { return CalculatePermeabilityMatrix(); } | ||
|
|
||
| Vector PermeabilityCalculator::RHSContribution() | ||
| { | ||
| return RHSContribution(CalculatePermeabilityMatrix()); | ||
| } | ||
|
|
||
| std::pair<Matrix, Vector> PermeabilityCalculator::LocalSystemContribution() | ||
| { | ||
| const auto permeability_matrix = CalculatePermeabilityMatrix(); | ||
| return {permeability_matrix, RHSContribution(permeability_matrix)}; | ||
| } | ||
|
|
||
| Vector PermeabilityCalculator::RHSContribution(const Matrix& rPermeabilityMatrix) const | ||
| { | ||
| return -prod(rPermeabilityMatrix, mInputProvider.GetNodalValues(WATER_PRESSURE)); | ||
| } | ||
|
|
||
| Matrix PermeabilityCalculator::CalculatePermeabilityMatrix() const | ||
| { | ||
| RetentionLaw::Parameters retention_parameters(mInputProvider.GetElementProperties()); | ||
| const auto& r_properties = mInputProvider.GetElementProperties(); | ||
| auto integration_coefficients = mInputProvider.GetIntegrationCoefficients(); | ||
| const auto shape_function_gradients = mInputProvider.GetShapeFunctionGradients(); | ||
| const auto local_dimension = shape_function_gradients[0].size2(); | ||
| const Matrix constitutive_matrix = | ||
| 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}}; | ||
|
Contributor
There was a problem hiding this comment. 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. |
||
| const double dynamic_viscosity_inverse = 1.0 / r_properties[DYNAMIC_VISCOSITY]; | ||
| for (unsigned int integration_point_index = 0; | ||
| integration_point_index < integration_coefficients.size(); ++integration_point_index) { | ||
| const double relative_permeability = | ||
| mInputProvider.GetRetentionLaws()[integration_point_index]->CalculateRelativePermeability(retention_parameters); | ||
| result += GeoTransportEquationUtilities::CalculatePermeabilityMatrix( | ||
| shape_function_gradients[integration_point_index], dynamic_viscosity_inverse, constitutive_matrix, | ||
| relative_permeability, integration_coefficients[integration_point_index]); | ||
| } | ||
| return result; | ||
| } | ||
|
|
||
| const Properties& PermeabilityCalculator::InputProvider::GetElementProperties() const | ||
| { | ||
| return mGetElementProperties(); | ||
| } | ||
|
|
||
| const std::vector<RetentionLaw::Pointer>& PermeabilityCalculator::InputProvider::GetRetentionLaws() const | ||
| { | ||
| return mGetRetentionLaws(); | ||
| } | ||
|
|
||
| Vector PermeabilityCalculator::InputProvider::GetIntegrationCoefficients() const | ||
| { | ||
| return mGetIntegrationCoefficients(); | ||
| } | ||
|
|
||
| Vector PermeabilityCalculator::InputProvider::GetNodalValues(const Variable<double>& rVariable) const | ||
| { | ||
| return mGetNodalValues(rVariable); | ||
| } | ||
|
|
||
| Geometry<Node>::ShapeFunctionsGradientsType PermeabilityCalculator::InputProvider::GetShapeFunctionGradients() const | ||
| { | ||
| return mGetShapeFunctionGradients(); | ||
| } | ||
|
|
||
| } // namespace Kratos | ||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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.