peprpendicular-flap solid-openfoam case needs tuning

[MakisH]

In #507, while comparing all the solver options we have, it became apparent that the solid-openfoam (solidDisplacementFoam) gives very different results than everything else. Still, it is a nice technical test that only requires OpenFOAM. The solid-solids4foam case is mostly fine already.

I have tried the following:

• Varying the mesh size: I get completely different results
• Varying the fvSolution/stressAnalysis/nCorrectors (suggestion by @philipcardiff): I get a bit different results, but not enough.

Some pictures that might help someone pick it up later on:

For now, I just give up, and document that this is an outlier at the moment.

One case use the plot-all-displacements.sh script to generate these pictures (see #507).

[MakisH]

Even with nCorr 10000 (mesh 6x15):

[MakisH]

For solid-solids4foam, increating nCorr from 1000 to 10k given exactly the same results as before (no improvement):

[philipcardiff]

I suspect another reason for the difference with openfoam-openfoam is that the 1st-order Euler ddt and d2dt schemes are used (solids4foam uses the 2nd-order backward scheme). You could try switching openfoam-openfoam to backward.

[MakisH]

The OpenFOAM documentation states that:

Maybe the d2dt2 backward is a nice solids4foam-only feature? 😄

with:

d2dt2Schemes
{
    default         Euler;
}

ddtSchemes
{
    default         backward;
}

and with also varying the nCorr, I get the following:

At least the ddtSchemes alone does not seem to help.

[philipcardiff]

Ah, OK, yes, we added a d2dt2 backward option in solids4foam.

In that case, using a smaller deltaT would improve the OpenFOAM predictions, if that is an option (probably also the solids4foam predictions).

[MakisH]

Decreasing the deltaT is indeed not really an option, as we wanted to have this consistent among the possible case combinations.

I compared the solid-openfoam and solid-solids4foam directories in a diff tool, and one other difference I found was the tolerances we use:

# solid-openfoam/system/fvSolution
solvers
{
    "(D|T)"
    {
        solver          GAMG;
        tolerance       1e-06;       # Made this 1e-09
        relTol          0.9;         # Made this 0.1
        smoother        GaussSeidel;
        nCellsInCoarsestLevel 20;
    }
}

# solid-solids4foam/system/fvSolution
solvers
{
    "D|DD"
    {
        solver          PCG;
        preconditioner  FDIC;
        tolerance       1e-09;
        relTol          0.1;
    }
}

I changed these tolerances, but it did not seem to improve the results.

Since we have the solid-solids4foam case, it is fine for me if this issue stays open for longer. Maybe someone else stumbles upon it and gets another idea, while it still helps as documentation.

[precice-bot]

This issue has been mentioned on preCICE Forum on Discourse. There might be relevant details there:

https://precice.discourse.group/t/the-perpendicular-flap-model-using-openfoam/1695/25

[philipcardiff]

Coming back to this with something concrete, though it addresses solid-solids4foam rather than the solid-openfoam case this issue is actually about — see the caveat at the end.

TL;DR

Switching solid-solids4foam to the high-order solid solver released in solids4foam v2.4 brings the flap-tip displacement essentially onto the CalculiX / deal.II / FEniCS / DUNE curves. The oscillation period goes from ~11% too short to within 0.3% of the others.

Numbers

All runs: fluid-openfoam (pimpleFoam) + solid-solids4foam, OpenFOAM v2606, preCICE v3.4.1, unchanged precice-config.xml and unchanged 6×15 solid mesh. "Consensus" is the mean of the stored openfoam-{calculix,dealii,fenics,dune}-v2404 watchpoints. relL2 is over the full 5 s of x-displacement; lag is by cross-correlation against the consensus over t ≥ 2 s; "amp ratio" is peak-to-peak relative to the consensus.

configuration	relL2	period [s]	lag [s]	amp ratio
consensus (CalculiX/deal.II/FEniCS/DUNE)	–	1.625	+0.00	1.000
the four of them individually	0.6–2.0 %			0.999–1.014
solids4foam v2.0 / OF-v2012 (stored reference)	47.7 %	1.480	−0.27	0.935
solids4foam v2.4 / OF-v2606, original case config	59.9 %	1.440	−0.37	0.903
solids4foam v2.4 / OF-v2606, high-order (p=3) — proposed	8.6 %	1.620	−0.03	1.076

Reading of this:

• The current solid-solids4foam case is not "mostly fine" — it is the second-largest outlier after solid-openfoam. Its period is ~9% too short, so it drifts steadily out of phase (−0.27 s by the end of the window). That is why it looks acceptable over the first period and diverges afterwards.
• This is a discretisation error, not a version or solver-technology effect. Re-running the original case config under solids4foam v2.4 / OF-v2606 gives 59.9% with period 1.440 s, i.e. no better than the stored v2.0 reference; the difference between those two is the version change alone.
• It is also not the PETSc/SNES solver doing the work. I checked this by keeping the proposed configuration and disabling only the high-order residual: the result is not better than the original case, it is worse (period 0.868 s, amplitude roughly half). The improvement comes from the high-order discretisation itself.

So the standard 2nd-order cell-centred FV discretisation is over-stiffening the flap on this deliberately coarse mesh, and the high-order (p=3) reconstruction removes most of that.

Runtime is unchanged in practice: ~55 s for the solid participant, ~2.0 coupling iterations per time window, same as before.

What is left, and why

It is a good match rather than a converged one: the amplitude is ~8% above the consensus, and after removing phase and amplitude there is a residual waveform-shape difference of ~11% (versus 1.7–4.5% within the FEM group).

My suspicion is that the remainder comes from the relatively low number of interface faces onto which the fluid forces are mapped, and a mesh study would show this. Although not shown here, I previously ran a mesh study with the low-order solver and it does converge towards the other results as the mesh is refined. I would expect the same for solidDisplacementFoam, though it would also need a smaller time step, since it uses the 1st-order Euler scheme.

For the purposes of this tutorial, the current results look acceptable as a demonstration.

Case changes

Three files in solid-solids4foam, no mesh and no preCICE config change:

• constant/solidProperties — solutionAlgorithm PETScSNES, a highOrderCoeffs block (highOrderResidual true, highOrderJacobian false, polynomialOrder 3, faceStencilExtraCells 15, radially-symmetric-exponential weighting with k=6), and an alpha momentum stabilisation.
• system/fvSolution — D solved with solver petsc driving a matrix-free newtonls SNES (snes_mf/snes_mf_operator), lgmres + bjacobi/lu, lagged preconditioner.
• system/fvSchemes — explicit per-term laplacian/snGrad/interpolate entries replaced by default entries, as the high-order path introduces terms the old explicit list does not name.

What this would require

1. solids4foam v2.4 or newer, built with PETSc. Hard requirement — the case will not start without it, since fvSolution selects solver petsc. The README currently says this case needs solids4foam v2.0 on up to OpenFOAM v2012 / OpenFOAM 9; that needs rewriting, and the foam-extend note effectively goes away for this case.
2. System tests: tests/dockerfiles/ubuntu_2404/Dockerfile pins solids4foam/solids4foam:v2.3-openfoam-v2412-latest, which would need bumping to a v2.4 image including PETSc.
3. Reference results need regenerating: perpendicular-flap/reference-results/fluid-openfoam_solid-solids4foam.tar.gz (used by the fieldcompare test via tests/tests.yaml), reference-results/watchpoints/openfoam-solids4foam-v2404.log, the su2-/nutils-/fake-solids4foam watchpoints, and the website plots from plot-all-displacements.sh.
4. Minor: solids4foam writes constant/solidProperties.withDefaultValues and constant/mechanicalProperties.withDefaultValues into the case, and clean_openfoam in tools/cleaning-tools.sh does not remove them, so they survive clean.sh. Worth adding.

The caveat

This does not close #515. That issue is about solid-openfoam (solidDisplacementFoam), which is a much bigger outlier (peak-to-peak 0.0026 vs 0.1620 — essentially dead) and has no high-order option, so none of the above transfers to it directly. What it does support is the diagnosis from earlier in this thread: on this mesh the solid discretisation is the dominant error, with the first-order Euler d2dt2 in solid-openfoam on top of that.

Happy to open a PR for the solid-solids4foam change if you think it is worth it, @MakisH — given the reference-result and Docker-image churn, I would rather ask first. It may also belong in a separate issue, since it is a different case from the one tracked here.

Environment

• OpenFOAM v2606, OpenFOAM-preCICE adapter, preCICE v3.4.1
• solids4foam development at v2.4-4-g9edc0db3 (four commits past the v2.4 tag), built with PETSc
• Serial runs; the proposed configuration reproduced bit-identically on a repeat

[MakisH]

Thanks a lot for looking closer into this, @philipcardiff!

Yes, a PR would be very welcome. I could update the reference results and any testing-related aspects in the same PR (we can definitely change the base Docker image). I started with #893, but happy to close it.

Note that part of the differences you now observe might be due to changed behavior in OpenFOAM v2606: see precice/openfoam-adapter#411