#98 fvc::d2dt2 reads its scheme from ddtSchemes, while fvm::d2dt2 reads d2dt2Schemes

[Philip Cardiff]

Summary

The explicit and implicit second time derivative operators look up their
discretisation scheme in different dictionaries:

• fvm::d2dt2 uses mesh.d2dt2Scheme("d2dt2(...)"), i.e. d2dt2Schemes;
• fvc::d2dt2 uses mesh.ddtScheme("d2dt2(...)"), i.e. ddtSchemes.

So a case whose d2dt2Schemes default differs from its ddtSchemes default
silently uses different time schemes for the implicit and explicit forms of the
same term. The explicit form falls back to the ddt default, and nothing warns
the user.

Location

src/finiteVolume/finiteVolume/fvc/fvcD2dt2.C, both overloads (with and
without density):

// fvc::d2dt2(vf)
return fv::d2dt2Scheme<Type>::New
(
    vf.mesh(),
    vf.mesh().ddtScheme("d2dt2(" + vf.name() + ')')   // expected: d2dt2Scheme
).ref().fvcD2dt2(vf);

// fvc::d2dt2(rho, vf)
vf.mesh().ddtScheme("d2dt2(" + rho.name() + ',' + vf.name() + ')')

The corresponding lines in fvm/fvmD2dt2.C call d2dt2Scheme(...).

Versions checked:

• OpenFOAM.com v2412 and v2512: ddtScheme(...). (i have also opened bug reports there)
• OpenFOAM.org 9 and 13: schemes().ddt(...). (i have also opened bug reports there)
• foam-extend 4.1 and 5.0: schemesDict().ddtScheme(...).

Expected behaviour

fvc::d2dt2 should use the scheme in d2dt2Schemes, as fvm::d2dt2 does.
Otherwise the implicit and explicit forms of the same operator can differ.

Impact / how it was found

Solid mechanics codes built on OpenFOAM (e.g. solids4foam) assemble a nonlinear
residual with fvc::d2dt2, for example for a Newton-Krylov (PETSc SNES) solve.
With d2dt2Schemes { default backward; } and ddtSchemes { default Euler; },
the residual's inertia term is silently computed with the Euler scheme, so the
solution is first order in time instead of second.

A quasi-static case with d2dt2Schemes { default steadyState; } and a transient
ddt default gets a spurious inertia term in the explicit residual.

Minimal reproduction

1. In any transient case, set ddtSchemes { default Euler; } and
   d2dt2Schemes { default backward; }.
2. After a few time steps with a time-varying field D, evaluate both operators on
   the same field: fvc::d2dt2(D), and the implicit form
   (fvm::d2dt2(D) & D) (i.e. fvm::d2dt2(D).A()*D - fvm::d2dt2(D).H() or
   equivalent).
3. The two differ, because fvc::d2dt2 has used the Euler scheme and
   fvm::d2dt2 the backward scheme. Setting both defaults to the same scheme
   makes them agree.
4. Adding "d2dt2\(.*\)" backward; to ddtSchemes also makes them agree. That
   confirms fvc::d2dt2 takes its scheme from ddtSchemes.

Suggested fix

Replace ddtScheme("d2dt2(...") with d2dt2Scheme("d2dt2(...") in both
fvc::d2dt2 overloads. On OpenFOAM.org use the schemes().d2dt2(...)
equivalent.

This could change results for cases whose ddtSchemes and d2dt2Schemes
defaults differ, in the explicit form only. Cases that set both the same are
unaffected.

Downstream tracking: https://github.com/solids4foam/solids4foam/issues/502

[Hrvoje Jasak]

Fixed in foam-extend-5.1, git commit 47f1f009ba945eb6b209d57ff06e3cd48b44b0a6

Thank you!

[Philip Cardiff]

Thank you!