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

Extend porous navier stokes #11236

Merged
merged 20 commits into from
Oct 19, 2023
Merged

Extend porous navier stokes #11236

merged 20 commits into from
Oct 19, 2023

Conversation

jgonzalezusua
Copy link
Member

@jgonzalezusua jgonzalezusua commented Jun 7, 2023

📝 Description
This PR is aimed to include the porous version of Navier-Stokes via VMS for both projections, ASGS and OSS.

Please mark the PR with appropriate tags:

  • Applications
  • Feature

🆕 Changelog

  • Renovated alternative_qs_vms_dem_coupled.cpp, alternative_d_vms_dem_coupled.cpp, qs_vms_dem_coupled.cpp and d_vms_dem_coupled.cpp
  • Added support for quadratic elements
  • Included tests for quadratic elements to test second derivatives of the shape functions
  • Added porosity term in cpp tests of the alternative_qs_vms_dem_coupled.cpp, alternative_d_vms_dem_coupled.cpp, qs_vms_dem_coupled.cpp and d_vms_dem_coupled.cpp elements.

@@ -229,7 +229,7 @@ class FluidElement : public Element
* @param rRightHandSideVector Local finite element residual vector (output)
* @param rCurrentProcessInfo Current ProcessInfo values (input)
*/
void CalculateLocalVelocityContribution(
virtual void CalculateLocalVelocityContribution(
Copy link
Member

Choose a reason for hiding this comment

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

Why is this required to be virtual?

Copy link
Member Author

Choose a reason for hiding this comment

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

This function is derived in the elements I have renovated. I need to initialize the second derivatives of the shape functions.

Copy link
Member

Choose a reason for hiding this comment

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

The point is that by doing so you're introducing the vtable overhead. Considering that you know at compile time which elements are high order as we've a dimension and number of nodes template argument, I think we can achieve this in a more generic and efficient way by using an if constexpr together with an overload of the method that calculates the kinematics. Note that this approach would make possible to keep a unique implementation for all the elements deriving from the FluidElement.

Copy link
Member Author

Choose a reason for hiding this comment

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

Ok, but we shall consider to implement GetShapeSecondDerivatives in FluidElement and modify 'CalculateLocalVelocityContribution' using an if constexpr.

Copy link
Member Author

Choose a reason for hiding this comment

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

Do you agree with that, @rubenzorrilla?

Copy link
Member Author

@jgonzalezusua jgonzalezusua Jun 8, 2023

Choose a reason for hiding this comment

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

@rubenzorrilla, I am still confused. If I remove GetShapeSecondDerivatives, where do I implement the transform of the second shape derivatives?

Copy link
Member

Choose a reason for hiding this comment

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

In a new overload of the CalculateGeometryData.

Copy link
Member Author

@jgonzalezusua jgonzalezusua Jun 9, 2023

Choose a reason for hiding this comment

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

@rubenzorrilla, this means that we must implement GetShapeSecondDerivatives in the fluid_element and overload the UpdateIntegrationPointData to call the second derivatives from data.

Copy link
Member

Choose a reason for hiding this comment

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

No, you do another version of the CalculateGeometryData (method overload) that does what we're currently doing plus the second derivatives.

Copy link
Member

Choose a reason for hiding this comment

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

It seems that this is not solved yet.

@@ -270,7 +270,7 @@ class BDF2TurbulentScheme : public Scheme<TSparseSpace, TDenseSpace>
* @param Dx Newton-Raphson iteration solution
* @param b Newton-Raphson right hand side (unused)
*/
void Update(
virtual void Update(
Copy link
Member

Choose a reason for hiding this comment

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

Why is the virtual required?

Copy link
Member Author

@jgonzalezusua jgonzalezusua Jun 7, 2023

Choose a reason for hiding this comment

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

Because this function is derived in bdf2_turbulent_schemeDEMCoupled.h, which will be merge to master after this PR.

Copy link
Member

Choose a reason for hiding this comment

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

But this is already virtual in the base class... I mean, it makes no sense to add it also here.

Copy link
Member Author

Choose a reason for hiding this comment

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

It is true

Comment on lines 216 to 226
void InitializeNonLinIteration(
ModelPart &rModelPart,
TSystemMatrixType &A,
TSystemVectorType &Dx,
TSystemVectorType &b) override
{

BaseType::InitializeNonLinIteration(rModelPart, A, Dx, b);

}

Copy link
Member

Choose a reason for hiding this comment

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

This is doing nothing. Please remove it and leave current base class implementation.

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, it was for a testing implementation, I forgot to remove it.

@@ -0,0 +1,113 @@
// | / |
Copy link
Member

Choose a reason for hiding this comment

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

I'd add this test to current test_alternative_qs_vms_dem_coupled_element.cpp. Note that you're testing the very same implementation but with different template instantiation.

@@ -0,0 +1,138 @@
// | / |
Copy link
Member

Choose a reason for hiding this comment

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

I'd add this test to current test_alternative_qs_vms_dem_coupled_element.cpp under the very same reasoning I posted above.

Copy link
Member Author

Choose a reason for hiding this comment

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

I thought that it would be better to split these tests in three in order to ease the detection of errors. However, I can add them in only one.

Copy link
Member

Choose a reason for hiding this comment

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

The tests are already split by the macro.

Copy link
Member Author

@jgonzalezusua jgonzalezusua Jun 7, 2023

Choose a reason for hiding this comment

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

Ok, so you mean to write them in the same .cpp file. I thought you meant to join them in only one test. To me it is ok

Copy link
Member

Choose a reason for hiding this comment

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

Exactly this.

@RiccardoRossi
Copy link
Member

is this done?

@jgonzalezusua
Copy link
Member Author

is this done?

@RiccardoRossi, in theory it is done, but we have decided to move the second derivatives transform to KratosCore. We are waiting for the approval of #11253.

@jgonzalezusua
Copy link
Member Author

It is very strange since in my local machine all tests run OK. But in github, the results in fluid cpp tests are different... How is it possible?

@rubenzorrilla
Copy link
Member

It is very strange since in my local machine all tests run OK. But in github, the results in fluid cpp tests are different... How is it possible?

It could be a memory error. Are you running the tests in FullDebug?

@jgonzalezusua
Copy link
Member Author

It is very strange since in my local machine all tests run OK. But in github, the results in fluid cpp tests are different... How is it possible?

It could be a memory error. Are you running the tests in FullDebug?

Yes, I am running in FullDebug. That's why is so strange.

@jgonzalezusua jgonzalezusua requested a review from a team as a code owner October 13, 2023 16:29
@jgonzalezusua
Copy link
Member Author

Now, tests should work. I have checked the tests with qs_vms.cpp and they give similar results. The difference stems from the calculation of element size. In our case, we are using AverageElementSize instead of MinimumElementSize.

@jgonzalezusua
Copy link
Member Author

@rubenzorrilla tests are runnning. We are ready to merge.

Copy link
Member

@rubenzorrilla rubenzorrilla left a comment

Choose a reason for hiding this comment

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

The issue with the virtual seems to be still pending. Aside of this everything looks good to me.

@jgonzalezusua
Copy link
Member Author

@rubenzorrilla, I have removed the virtual from fluid_element.h

Comment on lines +184 to +191
if(Dim == 2){
if (NumNodes == 9 || NumNodes == 6 || NumNodes == 4)
mInterpolationOrder = 2;
}
else if(Dim == 3){
if (NumNodes == 10 || NumNodes == 27)
mInterpolationOrder = 2;
}
Copy link
Member

Choose a reason for hiding this comment

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

As these are template parameters, in this case you can use if constexpr so you evaluate the if at compile time.

(I'd do this in a future PR considering that this is already taking a while...)

Copy link
Member Author

Choose a reason for hiding this comment

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

Yeah, you are right

Copy link
Member

@rubenzorrilla rubenzorrilla left a comment

Choose a reason for hiding this comment

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

No more changes in base elements. Thanks.

@jgonzalezusua jgonzalezusua merged commit ea7f22a into master Oct 19, 2023
@jgonzalezusua jgonzalezusua deleted the extend_porous_Navier_Stokes branch October 19, 2023 08:00
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants