[mn14cat]

Hello follow foamers,

I'm a PhD student at the university of leeds and I am currently working on a large strain poroelastic solver for analysis of cartilage tissue. I have successfully compiled my own large strain poroelastic solver developed in foam-extend 4.0.

However, I am currently running into floating point exception errors when running my case. The issue is I do not receive any further information other than floating point exception.

I have done some reading on the issue on other related threads and a lot of the post show the stack trace which gives some indication as to the cause of the error.

Is there any functionality in foam extend 4 which would provide me with the stack trace to give me some hope of being able to solve the floating point exception.

For any one who is interested in my research or my solver I am more than happy to provide the necessary files and further explanations of my project

Thank you in advance

[bigphil]

Originally Posted by mn14cat

Hello follow foamers,

I'm a PhD student at the university of leeds and I am currently working on a large strain poroelastic solver for analysis of cartilage tissue. I have successfully compiled my own large strain poroelastic solver developed in foam-extend 4.0.

However, I am currently running into floating point exception errors when running my case. The issue is I do not receive any further information other than floating point exception.

I have done some reading on the issue on other related threads and a lot of the post show the stack trace which gives some indication as to the cause of the error.

Is there any functionality in foam extend 4 which would provide me with the stack trace to give me some hope of being able to solve the floating point exception.

For any one who is interested in my research or my solver I am more than happy to provide the necessary files and further explanations of my project

Thank you in advance

Hi Craig,

The easiest way to debug is to add Info statements to the code, e.g.:

    Info<< "I am here!" << endl;

to find where the error comes from. Or if you want to be fancier:

    Info<< "File: " << __FILE__ << " at line: " << __LINE__ << endl;

If that does not help then you might find it useful to use a debugger with your code, e.g. as described here https://openfoamwiki.net/index.php/HowTo_debugging.
You can run the solver with the gdb debugger and then use the backtrace command to show exactly where the error comes from.

Philip

[mn14cat]

Hi Philip,

First can I say its an honour to receive a response from yourself, I've used a number of your papers related to solid mechanics with openfoam as a bible for the past two years whilst trying to develop my large strain poroelastic solver.

I have spent the better part of two years first learning openfoam from a complete novice to being able to modify/compile and debug a lot of openfoams fairly cryptic error messages haha.

My project is based around the poroelastic response of soft rubber like materials such as articular cartilage(Poisson ratio of around 0.45). So far I have couple together the pore pressure fluid equation to the total lagrangian large strain solid solver and I am able to run a full solution but for a very small parameter space not really representing material properties of cartilage.

However when I push my case model to more realistic values for the fluid equation I run into the floating point errors and as discussed in my previous post I receive no other information as to the cause of the floating point excepetion.

After doing some reading (I have previously read the web page you have posted) I understand that it is typically caused by dividing by zero (illegal mathematical operation). I was hoping to be able to get the print stack trace is this not possible using the pre-compiled functionality in openfoam ?

I have tried your suggestion of using the backtrace command with the gdb debugger and have now received so more information as to the casue:

Program received signal SIGFPE, Arithmetic exception.
0x00007ffff45818fd in Foam::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned int) const ()
from /home/mn14cat/foam/foam-extend-4.0/lib/linux64GccDPOpt/libfoam.so

(gdb) backtrace
#0 0x00007ffff45818fd in Foam::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned int) const ()
from /home/mn14cat/foam/foam-extend-4.0/lib/linux64GccDPOpt/libfoam.so
#1 0x00007ffff45a1065 in Foam::GAMGSolver::solveCoarsestLevel(Foam::Field<d ouble>&, Foam::Field<double> const&) const ()
from /home/mn14cat/foam/foam-extend-4.0/lib/linux64GccDPOpt/libfoam.so
#2 0x00007ffff45a1ddd in Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMa trix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned int) const ()
from /home/mn14cat/foam/foam-extend-4.0/lib/linux64GccDPOpt/libfoam.so
#3 0x00007ffff45a3acc in Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned int) const ()
from /home/mn14cat/foam/foam-extend-4.0/lib/linux64GccDPOpt/libfoam.so
#4 0x00007ffff712702b in Foam::fvMatrix<double>::solve(Foam::dictionary const&) ()
from /home/mn14cat/foam/foam-extend-4.0/lib/linux64GccDPOpt/libfiniteVolume.so
#5 0x000055555559358f in main ()

I can somewhat make out that the issue is coming from the PCG and GAMG solvers in my case but after that I'm not sure what the error message is telling me, If its not to much trouble could you help me understand exactly what the output is telling me ??

I'm more than happy to share case and solver files if you have the time to look through

[bigphil]

Hi Craig,

No problem. Always nice to discuss solver development!

OK, so the floating point error is coming from the linear solver. This is good and bad: good because it means it is not directly from some custom piece of code you wrote, but bad because this sort of error is very general and can be caused by many things (mesh, boundary conditions, material properties).

Basically there is something about your case that the solid solver is struggling to converge. It might be also that your code is all perfect but the current iterative solution algorithm is not robust enough for this type of case (from experience with poor-elasticity sometimes the convergence can be quite slow).

Have you tried under-relaxation for p and D (or U or whatever you call displacement)?

Also, are you able to describe your case (or a representative version of it) here (geometry/mesh, boundary conditions, and material law and properties)?

Philip

[mn14cat]

Good Morning Philip,

Apologies for the delay in response. Thank you for the response, so far I haven't tried under-relaxation for p and U but I ave seen this functionality in a number of the pre-compiled solid mechanics tutorials so I will have ago at that.

In terms of my case:

I am trying to replicate a confined compression test of cartilage, my geometry is a disk shaped of 20 mm diameter and 10 mm thickness. I am displacing the top surface by 4 mm over 1 sec to give 40% compressive strain. The tangential surfaces (sides) are confined laterally and are displaced as a function of the top surface.

I have implemented groovyBC to give me a time varying displacement BC on the top and side surfaces for U and for pressure the sides and BOTTOM surfaces are zerogradient and the TOP surface is fixedValue to represent it being open.

I am implementing a Saint Vennant Kirchhoff material law for the solid and material properties of a Youngs of 1.452 MPa and a Poissons ratio of 0.452

Here is the link for access to my case/solver files:

https://leeds365-my.sharepoint.com/:...-Byobl4ZEZ7L_A

Please let me know if you have any issue accessing the files.

Thanks in advance

[bigphil]

Do you set up the case as axisymmetric or with symmetries?

[mn14cat]

Hi Philip,

I've set the case as axis symmetric about y axis .

[bigphil]

I have generally found axisymmetric to be a bit fiddly in OpenFOAM; I am not sure if it is still the case but I think the case may need to be set up to symmetrically straddle the XY plane (min Z = max Z) with the X axis as the axis. You can find a couple of examples in solids4foam if you search for the "solidWedge" boundary condition.

In addition, I found that the mesh errors in the cells at the axis tend to be quite large due to the lack on non-orthogonal or skew correction on wedge patches and so I implemented the solidWedge boundary conditions to fix this. This may affect convergence.

Also, have a look at the under-relaxation values used for p and D in the poor-elastic tutorials in solids4foam.

Finally, St Venant Kirchhoff is a simple generalisation of linear elasticity to finite strains however its predictions are not very physical for truly larges. Laws like neo-Hookean tend to be more physical (although it obviously depends on the material to be simulated).

[mn14cat]

Thanks for your suggestions, I agree with you with respect to the St Venant Kirchhoff, however the choice of this material law was to offer the simplest route to developing a large strain poroelastic solver.

would you suggest implementing a change in the mesh to have a central block and then have 4 surround blocks and use face matching to generate a disk space geometry similar to an o-block geometry?

I really appreciate you help, I will let you know how get on if that's okay with you.

[bigphil]

Originally Posted by mn14cat

Thanks for your suggestions, I agree with you with respect to the St Venant Kirchhoff, however the choice of this material law was to offer the simplest route to developing a large strain poroelastic solver.

Sure. Understood.

Originally Posted by mn14cat

would you suggest implementing a change in the mesh to have a central block and then have 4 surround blocks and use face matching to generate a disk space geometry similar to an o-block geometry?

It is worth a try. There is no need to simulate the full geometry; you could use two symmetry planes and simulate one quarter of it (or even one eighth if you add another symmetry plane). The symmetry planes will avoid problems with rigid body translations.