Cook's membrane

[philipcardiff]

Hi @iBatistic,

Thanks for pushing the hyper-elastic Cook's membrane case.

It seems to work for me now with OF.com when using the SNES solver; see the case and scripts I pushed at https://github.com/solids4foam/solid-benchmarks/blob/main/hyperElasticity/cooksMembrane, and the predictions at https://github.com/solids4foam/solid-benchmarks/blob/main/hyperElasticity/cooksMembrane/times.mac-studio-m1.hex.lu.txt:

# Mesh Time (in s) Max-Memory (in kB) Disp
1 1 78176 0.0140949
2 2 78224 0.0139142
3 3 80176 0.01409
4 6 90208 0.0141904
5 18 122432 0.01425
6 71 243792 0.014282
7 337 812592 0.0142993

Is 0.0142993 close to the correct answer?

Note that you need to use the latest feature-petsc-snes branch of solids4foam; I was using OF-v2312 but it should also work with OF-v2012 or anything in between.

I have not yet checked the standard segregated solver in OF.com but I am fairly sure that the problem is the gradient scheme, as I had some issues with this when using the SNES solver (and hence created leastSquaresS4f).

By the way, in terms of a fair comparison of calculation time and memory in the paper, it might be better if we compare SNES-JFNK and SNES-segregated (Picard iterations). These will give exactly the same answer (apart from iteration error). Maybe somewhere in the paper, we can also comment on or compare it with the standard segregated approach (the gradient discretisation is slightly different at the boundaries).

It seems like a good idea for us to create a different GitHub issue to discuss each case for the paper. Once they all run and give reasonable answer, then it should be straightforward to prepare the results section of the paper.

[iBatistic]

Is 0.0142993 close to the correct answer?

Yes, it is fairly close. Below are the results from the deal.II:

It seems to work for me now with OF.com when using the SNES solver

That is great.

By the way, in terms of a fair comparison of calculation time and memory in the paper, it might be better if we compare SNES-JFNK and SNES-segregated (Picard iterations). These will give exactly the same answer (apart from iteration error). Maybe somewhere in the paper, we can also comment on or compare it with the standard segregated approach (the gradient discretisation is slightly different at the boundaries).

Sure, I agree.

It seems like a good idea for us to create a different GitHub issue to discuss each case for the paper. Once they all run and give reasonable answer, then it should be straightforward to prepare the results section of the paper.

Good idea!

p.s. I like the way you made those scripts for mesh study. Having one setup file is a neat idea.

[philipcardiff]

OK, great. It seems that we may be converging to a slightly different answer than deal.ii, but it is close.

By the way, I got the standard segregated solver working in OF.com: see the case here: https://github.com/solids4foam/solid-benchmarks/tree/main/hyperElasticity/cooksMembrane/base/hex.seg

Notice this part of the Allrun script:

# Use new version of extendedLeastSquares fo OF.com
# Once we are confident this works, we should update solids4FoamScripts
if [[ "${WM_PROJECT_VERSION}" == *"v"* ]]
then
    echo "OpenFOAM.com specific: using 'extendedLeastSquares 0' for the gradSchem in in system/fvSchemes"
    sed -i "s/ leastSquares;/ extendedLeastSquares 0;/g" "${CASE_DIR}"/system/fvSchemes
    sed -i "s/ pointCellsLeastSquares;/ extendedLeastSquares 0;/g" "${CASE_DIR}"/system/fvSchemes
fi

You will see I pushed changes to extendedLeastSquares in solids4foam for OF.com/OF.org and I made a fix for the mesh motion. These results seem to be (almost) the same as the FE41 results. You will see in my comments that maybe we should now default to extendedLeastSquares for all versions, or maybe we can rename this scheme to be solids4foam specific and make all solvers uses that. We can decide this later.

[iBatistic]

By the way, I got the standard segregated solver working in OF.com: see the case here: https://github.com/solids4foam/solid-benchmarks/tree/main/hyperElasticity/cooksMembrane/base/hex.seg

Yes, it is working. However, for some reason, I'm experiencing the same crash with v2012, while v2312 works like a charm.

You will see in my comments that maybe we should now default to extendedLeastSquares for all versions, or maybe we can rename this scheme to be solids4foam specific and make all solvers uses that. We can decide this later.

Nice, it is extremely useful to have the same name

[philipcardiff]

Hmnn, I am not sure about OFv2012. For the paper, let's stick with OFv2312.

[iBatistic]

Philp, how do you set the maximum number of iterations for snes.seg? Currently, the solver stops at 9999 iterations, and adjusting the snes_max_it parameter does not seem to affect this limit.

[philipcardiff]

-snes_max_it sets the maximum number of outer iterations and -ksp_max_it sets the maximum iterations for the linear (Krylov Subspace Methods = KSP) solver.

Which one is reaching 9999 iterations?

[iBatistic]

Outer iterations are reaching max value; -snes_max_it is working for numbers below 9999. Any number above is ignored and 10000 is set. I will let you know when I push the code so you can try. Hope that this limit is not hardcoded somewhere in PETSc code.

[iBatistic]

In the PR I opened, you can run hex.snes.seg to replicate this. The last two cases with the finest mesh are reaching this limit.

In the code, I have added:

• [ ] hex.segregated
• [ ] hex.snes.seg
• [ ] tet.segregated
• [ ] tet.snes
• [ ] tet.snes.seg

Tet cases can also be run as poly mesh (polyDualMesh).

I was thinking to also adding those base cases to hyperelastic and plastic Cooks membrane cases, but first I will finish this (you can leave PR open).

p.s. I slightly changed plot scripts by using load USER_VARIABLE_NAMES. By doing this, I managed to have different names of plots (name depends od inputs in USER_VARIABLE_NAMES).

[philipcardiff]

Outer iterations are reaching max value; -snes_max_it is working for numbers below 9999. Any number above is ignored and 10000 is set. I will let you know when I push the code so you can try. Hope that this limit is not hardcoded somewhere in PETSc code.

I can reproduce this behaviour; here is the error output I get:

9998 SNES Function norm 3.598369866926e-01 
9999 SNES Function norm 3.595468568952e-01 
--> FOAM Warning : 
    From void foamPetscSnesHelper::checkConvergence(SNES snes) const
    in file solidModels/solidModel/foamPetscSnesHelper/foamPetscSnesHelper.C at line 524
    PETSc SNES solver return error check disabled
The SNES nonlinear solver did not converge.
 PETSc SNES convergence error code: -2
 PETSc SNES convergence reason: DIVERGED_FUNCTION_COUNT

I Googled for DIVERGED_FUNCTION_COUNT, and I found that by default, PETSc only allows the user-defined function (our "residual" function) to be called 10,000 times. But we can increase this number with -snes_max_funcs. So I think we can added -snes_max_funcs 1000000 to the petscOptions file and it should work. I will check.

Wow, convergence is slow in this case. I will see if changing the relative tolerance of the linear solver will help (changing -ksp_rtol 0.99 to -ksp_rtol 0.1). And I found a MOOSE page about HYPRE saying that a good value for -pc_hypre_boomeramg_strong_threshold is 0.25` 2-D and 0.7 in 3-D. I will check if this helps.

[philipcardiff]

Here is the MOOSE page about HYPRE Boomerang: https://mooseframework.inl.gov/releases/moose/2021-09-15/application_development/hypre.html

[philipcardiff]

Convergence seems better with -ksp_rtol 0.1 and -pc_hypre_boomeramg_strong_threshold 0.25. From a quick check, the -ksp_rtol is the more important change.

I will push my updates.

[iBatistic]

Yep, it works fine now with those settings

[iBatistic]

Note for the cases I've tested with the AllrunMeshStudy script. I will update it as I continue running the cases.		Hookean	Neo-Hookean	Neo-Hookean-MisesPlastic
hex.snes	✔	✔	✘ 5 -8 SNES diverged
hex.snes.seg	✔	✔	✔
hex.seg	✔	✔	✔
tet.snes	✘**	✘ * **
tet.snes.seg	✘*	✘ *	✘ **
tet.seg	✔	✔	✔

Notes: nonLinearGeometryUpdatedLagrangian and nonLinearGeometryTotalLagrangian are crashing at the beginning. nonLinearGeometryTotalLagrangianTotalDisplacement is working fine.

* diverging at the begining. Setting the Rhie-Chow parameter to a higher value can help initiate the simulation. However, one cell in the upper-right corner shows an unphysical solution for stress

** Some simulations finished and some diverged. Related to issue *. It is interesting that always meshes 1,6,7 are failing.

@philipcardiff, below you can see the problem with the solution at the right upper corner cell. Based on this, I think that the problem is not related to the material type, solver, or SNES. It may be due to gradient calculation in the corner cells, where there is one internal face and two boundary faces.

[philipcardiff]

Great, thanks, Ivan!

FYI, I also noticed that our current JFNK approach does not work well with plasticity; I suspect that the "impK" preconditioning matrix is not a good approximation of the problem. Maybe that is a useful point to make in the paper. There may be some smart way we can scale the Laplacian with plasticity (be inspired by the definition of the true Jacobian) but this can be future work.

By the way, you will notice that I made some changes to Section 4 of the paper; I realised that we had no place to present the results of the different cases (e.g. stress/displacement fields, etc.). Now, I have arranged Section 4 as

• 4.1: Testing the accuracy and order of accuracy
  • Describe all cases and present the predictions (and comparisons with any reference). Analysis of time/memory and JFNK parameters are left to subsequent sections
  • 4.2: Time and memory requirements
  • 4.3: Effect of mesh type => I think I will merge this with 4.2
  • 4.4: Effect of the precondition choice => LU vs ILU vs HYPRE
  • 4.5: Effect of Rhie-Chow stabilisation
  • 4.6: Effect of globalisation strategies => to be removed or merged
  • 4.7: Parallelisation

[iBatistic]

FYI, I also noticed that our current JFNK approach does not work well with plasticity; I suspect that the "impK" preconditioning matrix is not a good approximation of the problem. Maybe that is a useful point to make in the paper. There may be some smart way we can scale the Laplacian with plasticity (be inspired by the definition of the true Jacobian) but this can be future work.

I reached the same conclusion. I spent some time adjusting the settings, experimenting with line search damping and various tolerance levels, but nothing seems to work once the region affected by plastic deformation begins to expand. As you said, we can leave this for later.

By the way, you will notice that I made some changes to Section 4 of the paper; I realised that we had no place to present the results of the different cases (e.g. stress/displacement fields, etc.). Now, I have arranged Section 4 as

I think it looks even better now.

I've made a few changes to the Ivan branch. The updates are minor, but you might find them useful. I've added bibliography entries for the references I believe you plan to include, as well as pictures of the geometry for all cases. You can use them or feel free to leave them out, as you already have some of these images from previous articles.

[philipcardiff]

@iBatistic : what page in the Zienkiewicz book is the linear elastic Cook's membrane described? I had a quick look but could not find it (it is a big book :D).

I found it described at https://www.researchgate.net/publication/268613538_On_the_main_properties_of_the_primal-mixed_finite_element_formulation (they set E=1 and tau=F=1), where they cite the original Cook paper (which I cannot get access to).

[iBatistic]

@iBatistic : what page in the Zienkiewicz book is the linear elastic Cook's membrane described? I had a quick look but could not find it (it is a big book :D).

I had the same problem (he doesn’t refer to it as Cook's membrane). The results can be found in section 11.5.4, on page 297.

In the book, they don't provide the dimensions for 𝐸 or specify the force amount. However, in this link, it shows that both the load and E are given in MPa, with the dimensions in mm. We have the same in our case. I saw your comment in the article mentioning that the disp are large—do you want us to make some adjustments?

For the comparison of tip displacement convergence (Figure 14), would you like three separate graphs including the results from the literature?

[iBatistic]

Philip, do you want diagrams like this?

We can reduce the load to 6250 Pa for the linear elastic case, and we will get 0.032 mm displacement. This allows us to use the results from Zienkiewicz and Taylor, but they will be more physical, i.e., smaller than the hyperelastic results.

[philipcardiff]

These graphs look ideal. They show that our solver converges to the same results as other solvers and give some idea of the order of convergence. You can add our solver to the legend, e.g., "Present work" or similar.

0.032 m is more physical, so yep, let's go with that for the small strain case.

For the second graph, our method seems less mesh-sensitive than Pelteret and McBride but converges to a slightly different solution. If needed, we can ask Dylan to run this case in Abaqus, too.

Feel free to add these figures to the paper.

It looks like we have almost finished Section "4.1 Testing the Accuracy and Order of Accuracy". In addition to Cook's membrane, I am waiting on the fine mesh for Hron-Turek to finish, and I also need to run the axial turbine blade; the more simple turbine blade case works fine, but the more complex one is not converging (or is converging very slowly). I will try a bit more to get the complex one to converge; if not, I will go with the less complex one.

Once they are done, I think the remainder of Section 4 will be relatively easy to complete. We're almost there!

[iBatistic]

Feel free to add these figures to the paper.

Sure, I'm waiting for one simulation to finish, and then I will make PR.

[iBatistic]

I came across article (incompressilbe neoHookean) where they actually used the midpoint to plot the graph xD.

For the elastoplastic case, in all the papers I’ve looked at, the solution is consistently just below 7 mm simplas Simo and Armero Cesar et al. Glaser and Armero

From what I understand, they are using the same material model that we are. Specifically, they refer to the Simo and Armero model, which is the same model described in the Computational Inelasticity book.

I tried using the useUndeformedArea option, and observed a 4.17% lower displacement with a 48x48 mesh. If this ratio holds even with finer meshes, we should obtain a displacement of around 7 mm, which aligns very closely with the results reported in the literature.

@philipcardiff , do you think this could be a similar case to what we encountered with the hyperelastic model? Today, I checked other factors—such as load increments, RC scaling, pressure equation, tolerances, etc.—but none of them affected the results as significantly as this.

[iBatistic]

update: With the useUndeformedArea option enabled, the finest mesh shows a displacement of 6.92 mm. César et al. presented their results in tabular form, reporting values of 6.92 mm and 6.97 mm depending on the element used. @philipcardiff, do you think that's the cause? According to the results, it seems that they all apply load in this manner

[philipcardiff]

update: With the useUndeformedArea option enabled, the finest mesh shows a displacement of 6.92 mm. César et al. presented their results in tabular form, reporting values of 6.92 mm and 6.97 mm depending on the element used. @philipcardiff, do you think that's the cause? According to the results, it seems that they all apply load in this manner

Ok, great. Yep, it sounds like that is the reason, and now we are getting good agreement (within the range of values reported). This will be sufficient for the paper. Thanks!

[iBatistic]

@philipcardiff, just a quick question. Below are the convergence results for the TLTD and UL solvers in case of elastoplastic case. it looks like the UL solver has a more "favorable curve". Should we continue using the TLTD solver for the elastoplastic case (we used TLTD in other two cases), or would it be better to switch to the UL solver given its more appealing convergence?

p.s. ignore the offset between results as they were generated without useUndeformedArea option. If we decide to go with the UL solver, we’ll need to add the code for the useUndeformedArea feature in the UL solver.

[philipcardiff]

I don't mind. Whatever is easiest.

For UL, we may need to store the original area vectors (since we will lose them). It may be easier to just use TL... I don't mind.