Crank-Nicolson totally wrong if any PDE in the set missing a time derivative term

[lindsayad]

Bug Description

If a PDE does not have a time derivative (say the mass equation for incompressible Navier-Stokes), then you will end up solving, for that sub-vector of dofs, the equation: residual + residual_old = 0. In my personal experience this yields solutions that look like this in time:

Steps to Reproduce

Solve a system of equations in which at least one equation does not have a time derivative. Your solution is in for a literally bumpy ride

Impact

Makes users distrust MOOSE

[lindsayad]

An additional note: we multiply the Wikipedia form of Crank-Nicolson by 2. I think this likely leads to inconsistent stabilization for SUPG and PSPG in INSFE because du/dt essentially becomes 2*du/dt ... as the mesh is refined, the sum of the strong momentum residuals will likely approach du/dt in this case as opposed to 0.

[YaqiWang]

This does not make sense though. When lacking time derivative term, old non-time residual will be numerically equal to zero, the full residual is the time residual (0 here) plus the current non-time residual plus old non-time residual, so driving that to zero by solvers is equivalent with driving the current non-time residual to zero. I do not see problems here. I admit that I'd like to see the u-dot not having that factor 2 there for a long time but implementation wise I do not see problem.

[YaqiWang]

More thought on this I do see that if the old non-time residual is not equal to zero, the system will tend to converge the current non-time residual to the opposite and having the oscillatory behavior, but that is the numerical effect. CN is known to have this kind of oscillatory solution error even with time derivatives.

[lindsayad]

You don't see a problem with having your new residual pre-determined to be -old_residual?

[YaqiWang]

Not a problem if old residual is zero. Plus I do not think it is correct to switch the code based on if time derivative is there or not.

[lindsayad]

So the user should know to setup an initial condition that magically yields a zero residual to their PDE? Why even use MOOSE then?

[YaqiWang]

Yes, for a variable without time derivative, it should not have initial condition at all. You will need to solve a system equation first to get the starting point correctly. Users do not have the freedom to set the initial condition.

[YaqiWang]

On the other hand, if this is too difficult to set up, we can add an artificial time derivative term for this variable, just to make sure the transient of this variable goes faster than the rest. I am not hundred percent sure if this is a good idea though.

[YaqiWang]

Actually I may take switch the code based on if time derivative is there or not back. We may clear the old residual for primal variables without time derivatives, which can be done fairly easily either in CrankNicolson::init() or in both there and CrankNicolson::postStep().

[lindsayad]

Conceptually time integrators should be applied to variables/equations and their associated dofs as appropriate. For example if solving a dynamic thermal mechanics problem, NewmarkBeta may be a good choice for the displacement variables but totally inappropriate for the transient heat conduction problem. But right now there is only a single time integrator and its methodology gets applied to all variables in all systems, appropriate or not.

[lindsayad]

Users do not have the freedom to set the initial condition.

I don't know what you're trying to say here. Are you arguing that this is how it should be or are you saying that's how it is? Because it is certainly not how it is. ICs are important whether you are solving a steady or transient problem or whether a given equation has a time derivative term or not. The ICs set for velocity and pressure are going to determine the RHS and initial _residual_old for the mass equation dofs in the problem above. The most natural user setup for a flow channel problem is Dirichlet inlet, no slip walls, and natural outlet with implicitly zero ICs. This instantly gives a non-zero _residual_old for the mass equation dofs because the divergence-free condition is violated.

A user of our INS module almost published work basically showing the figure I pasted above with the bcs and ics I outlined. It made MOOSE look really bad, and I can't blame the user here. We can and should do better.

[YaqiWang]

The initial condition I was referring is not ICs in MOOSE. The name in MOOSE is misleading, it would be better named as initialization. If you follow my suggestion on clearing the residual corresponding to variables without time derivatives, things will be ok. I am a long critic of MOOSE time integrators. But I do not think the code here is wrong. We can make it more defensive. Newmark in tensor mechanics module is fine because it is done through a non-time kernel.

[lindsayad]

But I do not think the code here is wrong.

Perhaps you can say it is not wrong, but at a minimum it is insufficient for handling problems that have an equation without a time derivative term. And if used unwittingly by a user, it gives results that can't be described any other way than wrong.

Newmark in tensor mechanics module is fine because it is done through a non-time kernel.

This is also wrong.

// parent class
template <bool is_ad>
using InertialForceParent = typename std::conditional<is_ad, ADTimeKernel, TimeKernel>::type;

template <bool is_ad>
class InertialForceTempl : public InertialForceParent<is_ad>

[YaqiWang]

Yes, unwittingly is wrong because he did not understand the equation well and failed to check the results and made sense out of it.

Yes, I was referring InertialForce. It was not like this before. Let me check when it is changed.

[lindsayad]

Another second order implicit time integrator in our arsenal gives results that are correct for the user who simply tries plug-and-play with our different time integrators. BDF2 and implicit euler give these results

which look a lot better than the ekg pulse that Crank-Nicolson gave.

[YaqiWang]

The change in InertialForce is introduced here 86651914d08a1d63123b5eaef3960743995b355f, before it is derived from Kernel...

IMO, time integration different from the time integrator added by moose is allowed. We have aux kernels can do time integration as well. In neutronics, we typically do the integration for the delayed neutron precursors differently. But some caveats need to be taken like showed in this issue #12347, which has not been resolved.

[lindsayad]

I just think we should make it impossible to generate the figure in my title post. An up-front error would have been better than letting the simulation run. An even better solution is to handle variables without time derivatives in a different way; your suggestion in https://github.com/idaholab/moose/issues/19228#issuecomment-953394282 is a very viable suggestion for that different handling.

IMO, time integration different from the time integrator added by moose is allowed

Yes, you are right. I would still prefer a first-class capability, however, where different TimeIntegrators (or no special residual treatment from TimeIntegrators in this case for the pressure variable) can be applied to different variables/equations.

[YaqiWang]

I agree. The code needs to be more defensive. We could as well warn users that the initial residual for the variables without time derivatives is expected to be zero but it does not, which means the value of those variables is insistent with other solutions.

[GiudGiud]

Likely related, when using the new TimeDerivativeAux with no time derivatives (transient but solving steady states), the results are also really weird.

• 1 process: some noise derivatives (1e-9) in the final time step, uneven spatially
• multiple processes: inconsistent noise, 0s sometimes. So the _dot fields are not great in the absence of time derivative kernels.