[GeoMechanicsApplication] Extract a static utility function for the calculation of the Damping Matrix (D) #12368
Conversation
| const Matrix& rMassMatrix, | ||
| const Matrix& rStiffnessMatrix) | ||
| { | ||
| return rMassMatrix * RayleighAlpha + rStiffnessMatrix * RayleighBeta; |
There was a problem hiding this comment.
Correct, but the formula is usually written as: alphaM + betaK
Please do the same here.
Is there a need for checking the sizes of M and K before doing this operation, or do we accept that will be o.k. as those are generated for the same element and therefore have the same and correct size.
There was a problem hiding this comment.
Thank you. Indeed ublas provides both options. Changed it. Yes the same size. In all three cases the matrices are generated just before CalculateDampingMatrix function.
| noalias(rDampingMatrix) = GeoEquationOfMotionUtilities::CalculateDampingMatrix( | ||
| rCurrentProcessInfo[RAYLEIGH_ALPHA], rCurrentProcessInfo[RAYLEIGH_BETA], mass_matrix, stiffness_matrix); |
There was a problem hiding this comment.
Apparently, the damping for the "structural" elements within geo have an implementation that differs from the continuum elements. I wonder whether we should make the possibility for a material based alpha and beta possible or to avoid that effort as we hope to let the "structural" elements within geo slowly disappear.
| // Compute Damping Matrix | ||
| if (rDampingMatrix.size1() != N_DOF) rDampingMatrix.resize(N_DOF, N_DOF, false); | ||
| noalias(rDampingMatrix) = ZeroMatrix(N_DOF, N_DOF); | ||
|
|
||
| const PropertiesType& rProp = this->GetProperties(); | ||
|
|
||
| if (rProp.Has(RAYLEIGH_ALPHA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_ALPHA] * MassMatrix; | ||
| else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * MassMatrix; | ||
|
|
||
| if (rProp.Has(RAYLEIGH_BETA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_BETA] * StiffnessMatrix; | ||
| else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_BETA] * StiffnessMatrix; | ||
| const PropertiesType& r_prop = this->GetProperties(); | ||
| noalias(rDampingMatrix) = GeoEquationOfMotionUtilities::CalculateDampingMatrix( | ||
| r_prop.Has(RAYLEIGH_ALPHA) ? r_prop[RAYLEIGH_ALPHA] : rCurrentProcessInfo[RAYLEIGH_ALPHA], | ||
| r_prop.Has(RAYLEIGH_BETA) ? r_prop[RAYLEIGH_BETA] : rCurrentProcessInfo[RAYLEIGH_BETA], | ||
| mass_matrix, stiffness_matrix); |
There was a problem hiding this comment.
Here, I think it might be nice to simplify by just assigning rDampingMatrix to GeoEquationOfMotionUtilities::CalculateDampingMatrix without noalias, such that we can remove the resizing and assigning it a zeromatrix.
There was a problem hiding this comment.
Let's start removing resize and noalias ;) Indeed it is not clear a benefit of noalias nowadays.
| @@ -58,4 +58,12 @@ Vector GeoEquationOfMotionUtilities::CalculateDetJsInitialConfiguration(const Ge | |||
| return det_Js_initial_configuration; | |||
| } | |||
|
|
|||
| Matrix GeoEquationOfMotionUtilities::CalculateDampingMatrix(double RayleighAlpha, | |||
There was a problem hiding this comment.
Do we ever want to check the inputs for the damping matrix (e.g. in this case I'd expect alpha and beta to be >0). We can also just assume the numbers to be correct, but since they are direct user inputs in this case, it might be a good idea to check them?
This might also be a broader discussion, since we don't do input validation for the jsons now, we might want to do that in a central place and avoid expensive checks for each time we call this function (which is done for each time step and each element of course)
There was a problem hiding this comment.
Checks should be as early as possible and once, so not here but at element initialization and once per material property/constitutive law.
There was a problem hiding this comment.
Fully agreed with Wijtze-Pieter.
| auto mass_matrix = scalar_matrix(n, n, mass_matrix_value); | ||
|
|
||
| constexpr double stiffness_matrix_value = 20; | ||
| auto stiffness_matrix = scalar_matrix(n, n, stiffness_matrix_value); |
There was a problem hiding this comment.
These can be const
| auto mass_matrix = scalar_matrix(n, n, mass_matrix_value); | |
| constexpr double stiffness_matrix_value = 20; | |
| auto stiffness_matrix = scalar_matrix(n, n, stiffness_matrix_value); | |
| const auto mass_matrix = scalar_matrix(n, n, mass_matrix_value); | |
| constexpr double stiffness_matrix_value = 20; | |
| const auto stiffness_matrix = scalar_matrix(n, n, stiffness_matrix_value); |
| rayleigh_alpha, rayleigh_beta, mass_matrix, stiffness_matrix); | ||
|
|
||
| const double expected_matrix_value = rayleigh_alpha * mass_matrix_value + rayleigh_beta * stiffness_matrix_value; | ||
| auto expected_damping_matrix = scalar_matrix(n, n, expected_matrix_value); |
There was a problem hiding this comment.
Can be const
| auto expected_damping_matrix = scalar_matrix(n, n, expected_matrix_value); | |
| const auto expected_damping_matrix = scalar_matrix(n, n, expected_matrix_value); |
| namespace Kratos::Testing | ||
| { | ||
|
|
||
| KRATOS_TEST_CASE_IN_SUITE(CalculateDampingMatrixGivesCorrectResults, KratosGeoMechanicsFastSuite) |
| #include "custom_utilities/equation_of_motion_utilities.h" | ||
| #include "includes/checks.h" | ||
| #include "testing/testing.h" | ||
| #include <boost/numeric/ublas/assignment.hpp> |
There was a problem hiding this comment.
I don't see assignment being used, could we remove this include?
There was a problem hiding this comment.
Thank you. It is removed.
applications/GeoMechanicsApplication/tests/cpp_tests/test_damping_matrix.cpp
Outdated
Show resolved
Hide resolved
| @@ -327,26 +328,19 @@ void UPwBaseElement<TDim, TNumNodes>::CalculateDampingMatrix(MatrixType& | |||
| const unsigned int N_DOF = this->GetNumberOfDOF(); | |||
|
|
|||
| // Compute Mass Matrix | |||
There was a problem hiding this comment.
Feel free to remove this comment, since it doesn't add anything that the code already tells.
There was a problem hiding this comment.
Thank you. Indeed the comments just repeat the function names. Removed
|
|
||
| this->CalculateMassMatrix(MassMatrix, rCurrentProcessInfo); | ||
| MatrixType mass_matrix(N_DOF, N_DOF); | ||
| this->CalculateMassMatrix(mass_matrix, rCurrentProcessInfo); | ||
|
|
||
| // Compute Stiffness matrix |
There was a problem hiding this comment.
Feel free to also remove this comment.
|
|
||
| this->CalculateMaterialStiffnessMatrix(StiffnessMatrix, rCurrentProcessInfo); | ||
| MatrixType stiffness_matrix(N_DOF, N_DOF); | ||
| this->CalculateMaterialStiffnessMatrix(stiffness_matrix, rCurrentProcessInfo); | ||
|
|
||
| // Compute Damping Matrix |
There was a problem hiding this comment.
And the same thing applies to this comment.
There was a problem hiding this comment.
Thanks for eliminating the duplication; this looks much better. I have the same suggestions about code comments as for the previous file.
There was a problem hiding this comment.
removed the comments
| rDampingMatrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix( | ||
| rCurrentProcessInfo[RAYLEIGH_ALPHA], rCurrentProcessInfo[RAYLEIGH_BETA], mass_matrix, stiffness_matrix); |
There was a problem hiding this comment.
The thing that struck me here was that apparently for "geo structural" elements the element properties need not be checked for the Rayleigh coefficients (contrary to the U-Pw elements). They are always taken from the "current process info".
There was a problem hiding this comment.
Same suggestions about the code comments.
There was a problem hiding this comment.
removed the comments
| damping_matrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix( | ||
| rayleigh_alpha, rayleigh_beta, mass_matrix, stiffness_matrix); | ||
|
|
||
| const double expected_matrix_value = rayleigh_alpha * mass_matrix_value + rayleigh_beta * stiffness_matrix_value; |
There was a problem hiding this comment.
Excellent unit test! The only thing I'm a bit hesitant about is repeating the implementation of the damping matrix computation here. Instead, we could also express the expected matrix like you did in the unit test above this one (see expected_mass_matrix).
There was a problem hiding this comment.
Since the expected matrix only has a single value, I would not express the matrix in the same way as expected_mass_matrix. We might want to give this double just the value it should have (which I guess is 15.0) and use expected_damping_matrix = scalar_matrix(n, n, expected_matrix_value) to fill the matrix as we do now. I'm not sure this approach is needed for such a simple calculation though. Not sure what you both think?
There was a problem hiding this comment.
Agreed with Richard. Now the calculated damping matrix is compared with expected damping matrix only. Each time the expected matrix is filled with a specific single value.
馃摑 Description