Trapezoidal rule in Euler implicit solver

[alxbilger]

Hi,

I would like to discuss the trapezoidal rule option in EulerImplicitSolver. I tried to derive the maths, but I cannot find the same result than in the code.

Here is my maths:

Here is the code:

mFact = 1 + tr * h * d_rayleighMass.getValue();
bFact = -tr * h;
kFact = -tr * h * (h + d_rayleighStiffness.getValue());

with $h=1/2$

If we remove the Rayleigh damping:

mFact = 1;
bFact = -tr * h;
kFact = -tr * h * h;

There is a discrepancy between the code and the maths for the factor of $K$. The code should be kFact = -tr * tr * h * h; according to the maths.

The code was introduced in c1acea6

[hugtalbot]

Hi @alxbilger

After "derouling" all equations, I definitely agree with you. In fact, there is more than just the factor K to be fixed.
Do you think about a scene illustrating the scheme? We could run a pendulum case which should be less numerically dissipated with the Trapezoidal rule.

See my associated PR #5169