Breaking-dam-2d crashes when using *MODAL DYNAMIC in CalculiX

Official adapters and tutorials
, ,
1

Hi everybody,

I’ve recently started working on modal superposition in an FSI setting, so I tried switching CalculiX in the breaking-dam-2d tutorial (OpenFOAM + CalculiX) from direct time integration to modal dynamics. The coupled simulation crashes. To make sure the problem isn’t in my own setup, I reproduced it on the unmodified tutorial, changing only the analysis step.

System and versions:

OS: Ubuntu 24.04.4 LTS
preCICE: 3.4.1
OpenFOAM: v2512
OpenFOAM adapter: 1.3.1
CalculiX: 2.20
CalculiX adapter: 2.20.2-1

What I changed: In the CalculiX input file, I replaced

*DYNAMIC, ALPHA=0.0, DIRECT

with

*MODAL DYNAMIC DIRECT

Everything else is the same as in the tutorial: mesh, material (E = 40 MPa, ν = 0, ρ = 3000), boundary conditions, *CLOAD, and the 2e-3 s coupling window.

Because the adapter rejects any input file containing *FREQUENCY, I computed the first 10 eigenmodes in a separate, uncoupled CalculiX run with

*FREQUENCY, STORAGE=YES
10

and copied the resulting .eig file to the modal dynamic job name so CalculiX picks it up.

The unmodified tutorial with *DYNAMIC runs without problems using the same preCICE configuration, so the problem seems to be on the structural side or at the interface.

What happens: Before the water reaches the flap at about t = 0.3 s, the loads are small and the watch-point displacement in x rises smoothly to about 1 mm. After impact, the x displacement reverses and swings to -2.32 mm at t = 0.362 s and then to +3.40 mm at t = 0.386 s. The y component of the coupling force has a small high/low alternation between neighbouring time windows from the start of the simulation. After impact this alternation slowly gets larger, and at t = 0.368 s both force components start changing sign every window. Within three windows the y force reaches about ±6 and the x force about ±0.25. The displacement stays smooth the whole time, even while the force changes sign. At t = 0.386 s, IQN-ILS fails to converge within 30 iterations (relative displacement residual 2.37, limit 1e-3). Two iterations into the next window, the simulation crashes.

The Fluid participant crashes first. alpha.water, which should stay between 0 and 1, jumps from about 1.01 to Min = -926.9, Max = 25.85 in a single MULES correction, followed by a floating-point exception in Foam::GAMGSolver::scale. This is the last time step before the crash (t = 0.39):

---[precice] Mapping "Force" for t=0.39 from "Fluid-Mesh" to "Solid-Mesh"
---[precice] Mapping "Displacement" for t=0.39 from "Solid-Mesh" to "Fluid-Mesh"
---[precice] it 2 (min: 1, max: 30), time-window 195, t 0.388 (max: 1), Dt 0.002, max-dt 0.002
Courant Number mean: 0.049499 max: 0.718678
Interface Courant Number mean: 0.00211176 max: 0.52678
Time = 0.39

PIMPLE: iteration 1
smoothSolver:  Solving for cellDisplacementx, Initial residual = 0.997277, Final residual = 1.25549e-17, No Iterations 3
smoothSolver:  Solving for cellDisplacementy, Initial residual = 0.998763, Final residual = 1.67704e-18, No Iterations 3
GAMGPCG:  Solving for pcorr, Initial residual = 1, Final residual = 2.83574e-06, No Iterations 8
GAMGPCG:  Solving for pcorr, Initial residual = 0.192751, Final residual = 5.85416e-06, No Iterations 5
time step continuity errors : sum local = 0.000100864, global = 1.27581e-05, cumulative = 4.31089e-05
smoothSolver:  Solving for alpha.water, Initial residual = 0.0463995, Final residual = 6.47399e-09, No Iterations 25
Phase-1 volume fraction = 0.195618  Min(alpha.water) = 5.60407e-39  Max(alpha.water) = 1.00992
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.212111  Min(alpha.water) = -926.892  Max(alpha.water) = 25.8514
[stack trace]
=============
#1  Foam::sigFpe::sigHandler(int) in /usr/lib/openfoam/openfoam2512/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so
#2  ? in /lib/x86_64-linux-gnu/libc.so.6
#3  Foam::GAMGSolver::scale(Foam::Field<double>&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const in /usr/lib/openfoam/openfoam2512/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so
#4  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const in /usr/lib/openfoam/openfoam2512/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so
#5  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in /usr/lib/openfoam/openfoam2512/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so
#6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) in /usr/lib/openfoam/openfoam2512/platforms/linux64GccDPInt32Opt/lib/libfiniteVolume.so
#7  ? in /usr/lib/openfoam/openfoam2512/platforms/linux64GccDPInt32Opt/lib/libfiniteVolume.so
#8  Foam::fvMesh::solve(Foam::fvMatrix<double>&, Foam::dictionary const&) const in /usr/lib/openfoam/openfoam2512/platforms/linux64GccDPInt32Opt/lib/libfiniteVolume.so
#9  ? in /usr/lib/openfoam/openfoam2512/platforms/linux64GccDPInt32Opt/bin/interFoam
#10  ? in /lib/x86_64-linux-gnu/libc.so.6
#11  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#12  ? in /usr/lib/openfoam/openfoam2512/platforms/linux64GccDPInt32Opt/bin/interFoam
=============
Floating point exception (core dumped)

What I have already tried: I first suspected an added-mass instability, but the added mass is not large enough to explain the divergence. I also tried using more eigenmodes, different time step sizes, and adding Rayleigh damping. None of these fixed the problem; the damping only slightly delayed the crash.

Attached files:

precice-config.xml (2.5 KB)

precice-Solid-convergence.log (46.0 KB)

precice-Solid-iterations.log (10.5 KB)

precice-Solid-watchpoint-Tip.log (22.5 KB)

fluid-openfoam.log (3.3 MB)

precice-Fluid-iterations.log (5.9 KB)

My questions:

  1. Has *MODAL DYNAMIC been tested with the CalculiX adapter?
  2. Is there a known problem with using it with implicit coupling?
  3. Is there anything else I should try?

Thanks in advance for any help :slight_smile:

Bianca

2

Hi @Bianca_Gritti,

we don’t have much experience with this feature, and we don’t have a tutorial either. Any observations or further input here would be very useful to improve the support, and rarely (if at all) someone asks for it, so I see potential here.

You probably already saw the documentation section of the CalculiX adapter, which is mostly there as a hint, rather than a tested feature:

Since the breaking-dam-2d tutorial is also known to have issues when running in parallel (see Breaking dam 2D fails in parallel · Issue #710 · precice/tutorials · GitHub), I would recommend starting with the perpendicular-flap. Does that work fine for you?

When the feature was developed, it was tested with the perpendicular-flap tutorial and with implicit coupling, but it has not seen much activity since then. One of the related PRs:

3

Hi,

thanks for the suggestions! I’ve been working on it, starting from the perpendicular flap (which work in air but not in water). I think I found the cause of my problem.

1. Checkpoint bug dyna_precice.c

In *MODAL DYNAMIC, the load is interpolated linearly over each time window, starting from aanew (the modal load at the start of the window). At an implicit-coupling checkpoint the adapter restores the modal coordinates (bj, bjp) but not aanew. So every repeated iteration ramps the load from the previous iteration’s force instead of F(tn). The coupling operator then changes from iteration to iteration and IQN-ILS struggles.

I first ran the perpendicular flap with run.sh -modal, unchanged: it runs and the result is within 1-2 mm of the direct *DYNAMIC solution. But some windows hit the max number of iterations without converging and it needs about 43% more coupling iterations (1436 vs 1005).

I then tried the same tutorial with water (ρ = 1000, ν = 1e-6, inflow scaled to 0.316 m/s): it crashes after 3 windows, with no convergence at all.

Saving aanew at Precice_WriteIterationCheckpointModal and restoring it at Precice_ReadIterationCheckpointModal should fix the convergence (2-4 iterations per window from the start).

2. Linear load instability with added mass
The interesting part is that the previous adapter version effectively applied an almost constant end-of-window load once converged. With the fix we get now a the linear interpolation which in water it develops a sign-flipping force oscillation that grows by about 3× per window and then crashes after 16 windows so the first fix was still not working for my case.
To solve it I tried applying the end-of-window load as constant over the window (as the direct solver does) and now the coupling is stable for the perpendicular flap case in water. This case now matches the direct *DYNAMIC run within about 2% (peak tip displacement 348 vs 344 mm, frequency 0.362 vs 0.363 Hz). The drawback is that we end up with a first order and that we will have a small energy error that grows with ω*Δt, so for impulsively excited high modes we need a smaller Δt.

perpendicular_flap_comparison
perpendicular_flap_comparison1320×900 167 KB

perpendicular_flap_water_comparison
perpendicular_flap_water_comparison1320×540 46.7 KB

I’ve attached the comparison plots, the first two show the results in an air like fluid with all the adapters, the last shows the same flap in water (with the new adapter working in wet conditions). I’m now working on the dam break+ flap scenario with modal dynamics. For the end-of-window load I currently use an environment variable (CCX_MODAL_END_LOAD=1).

Here are the changes I’ve made: Comparing develop...fix-modal-checkpoint-load · Biancagritti/calculix-adapter · GitHub

I’m happy to talk about any possible doubts or inconsistencies in my reasoning.

Thanks,
Bianca