Program Error: time step size is zero to machine precision.

[Chen Yicheng]

Dear all,

Thank you for reading this discussion. In my IBFE codes, a strange error repeated after running for some time. Does anyone know the source of this bug?
____________________________________________________________________________________
P=00005:
P=00008:Program abort called in file ``../../../IBAMR/ibtk/lib/../src/utilities/HierarchyIntegrator.cpp'' at line 258
P=00008:ERROR MESSAGE:
P=00008:IBHierarchyIntegrator::advanceHierarchy():
P=00008:  at time = 0.0111208: time step size dt = 9.96447e-11 is zero to machine precision.
P=00008:
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode -1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.

______________________________________________________________________________________

This error was not encountered in the beginning of program, but after it has been running for some time.  A larger time step size could delay it, but cannot eliminate the bug and also limited by the "time step size change". Increasing the number of solid grids also only delays bugs.

Thank you so much for your kind help!

Best,

Yicheng

[Wells, David]

Hi Yicheng,

If the time step is this small then something has gone wrong with the fluid solver - in particular, I suspect that the velocity is blowing up (has unphysically large values) at some point in the domain. Can you visualize the velocities and forces in visit? In addition - what kind of solid mechanics models are you using? They may be stiffer than we can reasonably handle with explicit time integration.

Best,

David

---

From: ibamr...@googlegroups.com <ibamr...@googlegroups.com> on behalf of Chen Yicheng <cychi...@gmail.com>
Sent: Monday, February 14, 2022 12:28 AM
To: IBAMR Users <ibamr...@googlegroups.com>
Subject: [ibamr-users] Program Error: time step size is zero to machine precision.

 

--
You received this message because you are subscribed to the Google Groups "IBAMR Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ibamr-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/9e699f64-f7bb-4809-b5eb-e40cd3b39b6dn%40googlegroups.com.

[Chen Yicheng]

Hi David,

Thank you for your help. I am trying to perform a pressure drive flow in a flexible tube, the wall was set to be a thin plate as follow. The normal pulsating pressure was applied to boundary_0 on the left.
Height of tube = 1.0, Length = 25.0, Thickness of wall = 0.01, excepted Re = 100.

Here are the velocities and pressure before program broke (left, U; right, P).

[Wells, David]

Hi Yicheng,

Thanks. It looks like there is a checkerboarding problem with the pressure: the maximum value and the minimum value are both huge and adjacent. How fine is the grid resolution relative to the spacing of those peaks?

Best,

David

---

From: ibamr...@googlegroups.com <ibamr...@googlegroups.com> on behalf of Chen Yicheng <cychi...@gmail.com>

Sent: Monday, February 14, 2022 9:29 PM
To: IBAMR Users <ibamr...@googlegroups.com>
Subject: Re: [ibamr-users] Program Error: time step size is zero to machine precision.

 

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/16ef7c41-3425-445b-9100-8e6ac30475ban%40googlegroups.com.

[Boyce Griffith]

This is probably not a fluid solver problem — I’ll follow up when I get a chance later today.

On Feb 15, 2022, at 11:22 AM, Wells, David <drw...@email.unc.edu> wrote:



To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/MN2PR03MB4879957BCDBAAB60FE7B4416ED349%40MN2PR03MB4879.namprd03.prod.outlook.com.

[Yicheng Chen]

Hi David, 

The minimum Cartesian grid size was 0.0125 while the distance between two peak value was about 0.05. Dose the grid resolution not not applicable?

The maximum and minimum value are adjacent near (5.0, 1.5), where is the segment point of PK1 force. The upper wall was set to be rigid ≤5.0 (using penalty force) and flexible between 5.0 and 10.0. Could this be the cause of sharp changing?

Best,

Yicheng

[Yicheng Chen]

Thank you so much, Professor! 

Best,

Yicheng

[Boyce Griffith]

Sorry for the delay in getting back to you.

In IBAMR, the time step size is determined to satisfy what is called the CFL condition: dt <= C dx / u_max. (For most of the methods implemented in IBAMR, I recommend starting with a CFL number of C = 0.25 or less.) If the flow speed increases, dt may need to be decreased to satisfy this condition.

For FSI models, the CFL condition is not the only stability condition that must be satisfied. In our current FSI time stepping schemes, the structural dynamics are actually treated explicitly, and so there will be a largest stable time step size that is associated with the structure model. We don't currently have a way to dynamically determine this largest stable time step size, however. Usually increases in the flow velocity are actually an indication that the structural dynamics are becoming unstable.

When setting up a new model, I recommend running with the option "error_on_dt_change" set to "TRUE". (In most of the examples in IBAMR, this can be set by a top-level variable, "ERROR_ON_DT_CHANGE".) This will cause the code to abort immediately if there is a time step size change. In this case, I would recommend reducing the time step size and re-running the model. If the time step size change is resulting from a structural instability, eventually you should be able to find a stable time step size.

On the other hand, it also is possible that the flow speed is large enough to require a time step size change to satisfy the CFL condition. In this case, you should be able to set "error_on_dt_change" to "FALSE". However, if the time step size is going to zero, this almost always means that there is an instability in the structural dynamics part of the model.

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/53e4de7e-c3f2-4863-bce5-d44ea07af4e4n%40googlegroups.com.

[Yicheng Chen]

Oh! Thank you so much for your detailed explanation! Can I take it this way? In IBAMR, time step errors (whether Time step size change or time step os zero to machine precision) are due to structural instability, and modification of time step size can improve the robustness of code to predict structural dynamics. In other simulations I've tried, it's true that we can increase or decrease the time step size to keep the simulation stable. 

One surprise error is that smaller or bigger time step size will still bring the 'time-to-0' bug. To figure out what happened, structural deformation are checked based on your analysis. No matter how little drive (inlet) pressure I apply, there is always asymmetric deformation of the structure (thickness = 0.01), then induced time-to-0. However, the previous simulation with larger thickness (0.1) performed symmetrically and worked very well (See the figure below). Is it possible for me to improve my calculation without changing the thickness?

[Boyce Griffith]

On Feb 18, 2022, at 7:37 AM, Yicheng Chen <cychi...@gmail.com> wrote:

Oh! Thank you so much for your detailed explanation! Can I take it this way? In IBAMR, time step errors (whether Time step size change or time step os zero to machine precision) are due to structural instability, and modification of time step size can improve the robustness of code to predict structural dynamics.

It depends. If the flow speed increases enough, that can require a time step size change to satisfy the CFL condition that is required by the fluid solver --- in that case, running at the same time step size would result in a fluid solver instability. However, for FSI models with relatively stiff structures, the dominant stability constraint is usually related to the structural dynamics.

The reason that the time step size can shrink to zero is that if the time step size violates the structural time step size restriction, this will cause energy in the model to blow up, which generates large forces that generates large velocities and, consequently, very small time step sizes to maintain a CFL-type condition. (The model problem I think about is y' = -lambda y. If you discretize via forward Euler and pick a time step size that is too large, then the solution blows up instead of decaying to zero, independent of the initial condition.) One way to monitor whether this kind of blow up is happening is to track the CFL number in the log file. If you see the CFL number increasing by a consistent multiplicative factor from time step to time step (e.g. consistently going up by a factor of 2 or 5 or 10 in each time step), then that usually means that the time step size is too large for the structural model.

In other simulations I've tried, it's true that we can increase or decrease the time step size to keep the simulation stable. 

If things are unstable, you usually cannot improve thing by increasing the time step size.

On the other hand, it is beneficial to find the largest stable time step size, since it impacts the overall runtime for the model --- if I have a model that is stable that I want to run for a long time, I will usually spend some time to see if the time step size is close to optimal.

One surprise error is that smaller or bigger time step size will still bring the 'time-to-0' bug.

Let me know if the explanation above makes sense for how the time step size can shrink to zero.

To figure out what happened, structural deformation are checked based on your analysis. No matter how little drive (inlet) pressure I apply, there is always asymmetric deformation of the structure (thickness = 0.01), then induced time-to-0.

The structural stability conditions are not strongly related to the flow conditions --- they are determined by the stiffness of the material model, the grid spacings, and the form of the regularized delta function. Specifically, if you have a very stiff material but no flow or loading --- and so no induced motion --- you still would expect that any small perturbations to the model (which can result from finite precision effects) will eventually cause the structure to "blow up" if the time step size is too large.

 However, the previous simulation with larger thickness (0.1) performed symmetrically and worked very well (See the figure below). Is it possible for me to improve my calculation without changing the thickness?

I think it should be possible to make the thin structure work. What is the material model and structural discretization?

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/6c0d2e0b-6384-4924-bbc6-66038790c696n%40googlegroups.com.
<微信图片_20220218203533.png>

[Yicheng Chen]

Thank you so much!  

                 It depends. If the flow speed increases enough, that can require a time step size change to satisfy the CFL condition that is required by the fluid solver --- in that case, running at the same time step size would result in a fluid solver instability. However, for FSI models with relatively 

                  stiff structures, the dominant stability constraint is usually related to the structural dynamics.

                  The reason that the time step size can shrink to zero is that if the time step size violates the structural time step size restriction, this will cause energy in the model to blow up, which generates large forces that generates large velocities and, consequently, very small time step sizes 

                  to maintain a CFL-type condition. (The model problem I think about is y' = -lambda y. If you discretize via forward Euler and pick a time step size that is too large, then the solution blows up instead of decaying to zero, independent of the initial condition.) One way to monitor whether 

                  this kind of blow up is happening is to track the CFL number in the log file. If you see the CFL number increasing by a consistent multiplicative factor from time step to time step (e.g. consistently going up by a factor of 2 or 5 or 10 in each time step), then that usually means that the 

                  time step size is too large for the structural model.

In the last 100 step before 'time-to-0', the CFL is almost fixed around 0.00112487, but suddenly increased to 1421.58 in the last step as follow. Is it a 'blow up' in this time step?

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

At beginning of timestep # 48458
Simulation time is 0.121145
IBHierarchyIntegrator::advanceHierarchy(): time interval = [0.121145,0.121148], dt = 2.5e-06
IBHierarchyIntegrator::preprocessIntegrateHierarchy(): performing Lagrangian forward Euler step
IBHierarchyIntegrator::advanceHierarchy(): integrating hierarchy
IBHierarchyIntegrator::integrateHierarchy(): computing Lagrangian force
IBHierarchyIntegrator::integrateHierarchy(): spreading Lagrangian force to the Eulerian grid
IBHierarchyIntegrator::integrateHierarchy(): solving the incompressible Navier-Stokes equations
INSStaggeredHierarchyIntegrator::integrateHierarchy(): stokes solve number of iterations = 1
INSStaggeredHierarchyIntegrator::integrateHierarchy(): stokes solve residual norm        = 11.9153
IBHierarchyIntegrator::integrateHierarchy(): interpolating Eulerian velocity to the Lagrangian mesh
IBHierarchyIntegrator::integrateHierarchy(): performing Lagrangian midpoint-rule step
IBHierarchyIntegrator::postprocessIntegrateHierarchy(): interpolating Eulerian velocity to the Lagrangian mesh
IBHierarchyIntegrator::postprocessIntegrateHierarchy(): CFL number = 0.00112487
IBHierarchyIntegrator::postprocessIntegrateHierarchy(): Eulerian estimate of upper bound on IB point displacement since last regrid = 0.00224975
IBHierarchyIntegrator::advanceHierarchy(): synchronizing updated data
IBHierarchyIntegrator::advanceHierarchy(): resetting time dependent data

At end       of timestep # 48458
Simulation time is 0.121148
+++++++++++++++++++++++++++++++++++++++++++++++++++


+++++++++++++++++++++++++++++++++++++++++++++++++++
At beginning of timestep # 48459
Simulation time is 0.121148
IBHierarchyIntegrator::advanceHierarchy(): time interval = [0.121148,0.12115], dt = 2.5e-06
IBHierarchyIntegrator::preprocessIntegrateHierarchy(): performing Lagrangian forward Euler step
IBHierarchyIntegrator::advanceHierarchy(): integrating hierarchy
IBHierarchyIntegrator::integrateHierarchy(): computing Lagrangian force
IBHierarchyIntegrator::integrateHierarchy(): spreading Lagrangian force to the Eulerian grid
IBHierarchyIntegrator::integrateHierarchy(): solving the incompressible Navier-Stokes equations
INSStaggeredHierarchyIntegrator::integrateHierarchy(): stokes solve number of iterations = 6960
INSStaggeredHierarchyIntegrator::integrateHierarchy(): stokes solve residual norm        = 5.31686e+10
IBHierarchyIntegrator::integrateHierarchy(): interpolating Eulerian velocity to the Lagrangian mesh
IBHierarchyIntegrator::integrateHierarchy(): performing Lagrangian midpoint-rule step
IBHierarchyIntegrator::postprocessIntegrateHierarchy(): interpolating Eulerian velocity to the Lagrangian mesh
IBHierarchyIntegrator::postprocessIntegrateHierarchy(): CFL number = 1421.58
IBHierarchyIntegrator::postprocessIntegrateHierarchy(): Eulerian estimate of upper bound on IB point displacement since last regrid = 1421.58
IBHierarchyIntegrator::advanceHierarchy(): synchronizing updated data
IBHierarchyIntegrator::advanceHierarchy(): resetting time dependent data

At end       of timestep # 48459
Simulation time is 0.12115
+++++++++++++++++++++++++++++++++++++++++++++++++++


+++++++++++++++++++++++++++++++++++++++++++++++++++
At beginning of timestep # 48460
Simulation time is 0.12115
P=00000:Program abort called in file ``../../../IBAMR/ibtk/lib/../src/utilities/HierarchyIntegrator.cpp'' at line 258
P=00000:ERROR MESSAGE:
P=00000:IBHierarchyIntegrator::advanceHierarchy():
P=00000:  at time = 0.12115: time step size dt = 1.75861e-10 is zero to machine precision.
P=00000:

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

                    The structural stability conditions are not strongly related to the flow conditions --- they are determined by the stiffness of the material model, the grid spacings, and the form of the regularized delta function. Specifically, if you have a very stiff material but no flow or loading --- and 

                    so no induced motion --- you still would expect that any small perturbations to the model (which can result from finite precision effects) will eventually cause the structure to "blow up" if the time step size is too large.

                    I think it should be possible to make the thin structure work. What is the material model and structural discretization?

My material is not very rigid and has loadings - a prestress and an outside pressure. Here is the material model: E = 100 dyne/cm2, Nu = 0.49. The prestress was 1000*E to simplify the vibration modes, and the outside pressure was 0.1*E to bring disturbance to the flexible wall. The IB_DELTA_FUNCTION was IB_3, and the structural discretization was QUAD4 in ELEM_TYPE.  

[Boyce Griffith]

Try using a QUAD9 mesh and, if you see the same behavior, see what happens if you increase the bulk modulus.

On Feb 18, 2022, at 11:39 PM, Yicheng Chen <cychi...@gmail.com> wrote:



To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/256e2d2f-654e-4aab-8fdf-134b572deff2n%40googlegroups.com.

[Yicheng Chen]

Great! The change of mesh type eliminates this problem! The 9 noded quadrilaterals can better keep the structural dynamic stable, but I can not find any explaination in libMesh.

Thank you so much, professor!

Best,

Yicheng

[Boyce Griffith]

On Feb 20, 2022, at 12:31 AM, Yicheng Chen <cychi...@gmail.com> wrote:

Great! The change of mesh type eliminates this problem! The 9 noded quadrilaterals can better keep the structural dynamic stable, but I can not find any explaination in libMesh.

Great. I’d suggest using QUAD9 for the thicker case as well since your (numerical) Poisson ratio is close to 0.5.

Note that, in principle, IBFEMethod treats all structures as exactly incompressible. It is still useful in practice to include volumetric energies/stresses in material models to reduce discretization errors that impact volume conservation. See this paper for discussion and tests: https://doi.org/10.1016/j.cma.2020.112978 or https://arxiv.org/abs/1811.06620.

To view this discussion on the web visit https://groups.google.com/d/msgid/ibamr-users/492537fa-4605-4282-a617-12a4da938462n%40googlegroups.com.

[Yicheng Chen]

Thanks so much for providing references!! I'll study them right away. ^_^