Some time integrators do not work as expected

[YaqiWang]

Rationale

I am able to create a simple test to show the unexpected behavior of time integrators. The test is small so I will just copy it here:

[Mesh]
  type = GeneratedMesh
  dim  = 2
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  nx   = 2
  ny   = 2
[]

[Variables]
  [./u]
  [../]
[]

[AuxVariables]
  [./v]
  [../]
  [./c]
    initial_condition = 1
  [../]
[]

[Kernels]
  [./ie]
    type = TimeDerivative
    variable = u
  [../]

  [./cv]
    type = CoupledForce
    variable = u
    v = v
  [../]
[]

[AuxKernels]
  [./vt]
    type = VariableTimeIntegrationAux
    variable = v
    variable_to_integrate = c
  [../]
[]

[Postprocessors]
  [./intu]
    type = ElementIntegralVariablePostprocessor
    variable = u
  [../]
[]

[Executioner]
  type = Transient
  start_time = 0.0
  end_time   = 1.0
  dt = 0.1
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-12
[]

[Outputs]
  csv = true
[]

Solution of u is t*t/2. No spatial discretization error should be seen. Thus the exact intu is 0.5.

Description

Backward Euler: error is dt/2, first order;
Forward Euler: error is dt/2, first order;
BDF2: error is about 3/4*dt^2, second order;
Crank-Nicolson: error is zero;
Explicit-midpoint: converges to wrong value 1;
DIRK: error is dt/2, first order;
Explicit-tvd-rk-2: converges to wrong value 1.

Looks like BE, FE, BDF2, CN work as expected, but not the rest. I think these integrators cannot handle stateful aux variables properly. I did not try stateful material properties, but guess we will see similar issue.

Impact

Fix buggy time integrators.

[jwpeterson]

Just wanted to note that a modified version of @YaqiWang's input file that solves the same problem without involving the "stateful" AuxKernel VariableTimeIntegrationAux behaves the same way for both Backward Euler and Crank-Nicolson, and similarly to Crank-Nicolson, LStableDirk2 is also exact for this version of the problem.

So I think the issue is definitely traceable to AuxVariables (and possibly MaterialProperties) not being updated correctly for the "multi-stage" methods (ExplicitMidpoint, LStableDirk2, explicit-tvd-rk-2, etc).

# This is a modification of Yaqi's original problem that simply solves du/dt = t
# without involving stateful AuxVariables/Postprocessors.
[Mesh]
  type = GeneratedMesh
  dim  = 2
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  nx   = 2
  ny   = 2
[]

[Variables]
  [./u]
  [../]
[]

[Functions]
  [./forcing_fn]
    type = ParsedFunction
    value = 't'
  [../]

  [./exact_fn]
    type = ParsedFunction
    value = '0.5*t*t'
  [../]
[]

[Kernels]
  [./ie]
    type = TimeDerivative
    variable = u
  [../]

  [./cv]
    type = BodyForce
    variable = u
    function = forcing_fn
  [../]
[]

[Postprocessors]
  # Estimate spatial norm of error at fixed time, ||e||_{L2}
  [./l2_err]
    type = ElementL2Error
    variable = u
    function = exact_fn
  [../]
  # Estimate spacetime norm ||e||_{L2, \infty}
  [./max_l2_err]
    type = TimeExtremeValue
    value_type = max
    postprocessor = l2_err
  [../]
  # Estimate spacetime norm ||e||_{L2, L1}
  [./cumulative_l2_err]
    type = CumulativeValuePostprocessor
    postprocessor = l2_err
  [../]
[]

[Executioner]
  type = Transient
  start_time = 0.0
  end_time   = 1.0
  # dt = 1.0
  # dt = 0.5
  # dt = 0.25
  # dt = 0.125
  # dt = 0.0625
  dt = 0.03125
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-12

  # Use to set different TimeIntegrators
  [./TimeIntegrator]
    # type = ImplicitEuler
    # type = CrankNicolson
    type = LStableDirk2
  [../]
[]

[Outputs]
  csv = true
[]

[jwpeterson]

I added a print statement at the top of AuxiliarySystem::compute() and around each stage solve of LStableDirk2 to verify whether the AuxKernels are being updated correctly before the second stage starts, and it appears that they are:

Called AuxiliarySystem::compute(TIMESTEP_BEGIN)

Time Step 2, time = 1
                dt = 0.5

Called AuxiliarySystem::compute(LINEAR)
Starting LStableDirk2 first stage solve.
Called AuxiliarySystem::compute(LINEAR)
1st stage
 0 Nonlinear |R| = 1.098350e-01
Called AuxiliarySystem::compute(NONLINEAR)
      0 Linear |R| = 1.098350e-01
Called AuxiliarySystem::compute(LINEAR)
Called AuxiliarySystem::compute(LINEAR)
      1 Linear |R| = 2.886246e-03
Called AuxiliarySystem::compute(LINEAR)
      2 Linear |R| = 1.595813e-11
Called AuxiliarySystem::compute(LINEAR)
Called AuxiliarySystem::compute(LINEAR)
 1 Nonlinear |R| = 1.086367e-10
Called AuxiliarySystem::compute(NONLINEAR)
      0 Linear |R| = 1.086367e-10
Called AuxiliarySystem::compute(LINEAR)
Called AuxiliarySystem::compute(LINEAR)
      1 Linear |R| = 8.779707e-12
Called AuxiliarySystem::compute(LINEAR)
      2 Linear |R| = 2.713448e-19
Called AuxiliarySystem::compute(LINEAR)
Called AuxiliarySystem::compute(LINEAR)
 2 Nonlinear |R| = 7.757919e-18
Finished LStableDirk2 first stage solve.
Starting LStableDirk2 second stage solve.
Called AuxiliarySystem::compute(LINEAR)
2nd stage
 0 Nonlinear |R| = 2.651650e-01
Called AuxiliarySystem::compute(NONLINEAR)
      0 Linear |R| = 2.651650e-01
Called AuxiliarySystem::compute(LINEAR)
Called AuxiliarySystem::compute(LINEAR)
      1 Linear |R| = 6.968015e-03
Called AuxiliarySystem::compute(LINEAR)
      2 Linear |R| = 1.284181e-10
Called AuxiliarySystem::compute(LINEAR)
Called AuxiliarySystem::compute(LINEAR)
 1 Nonlinear |R| = 8.321540e-10
Called AuxiliarySystem::compute(NONLINEAR)
      0 Linear |R| = 8.321540e-10
Called AuxiliarySystem::compute(LINEAR)
Called AuxiliarySystem::compute(LINEAR)
      1 Linear |R| = 5.005130e-11
Called AuxiliarySystem::compute(LINEAR)
      2 Linear |R| = 6.056120e-12
Called AuxiliarySystem::compute(LINEAR)
      3 Linear |R| = 5.572965e-19
Called AuxiliarySystem::compute(LINEAR)
Called AuxiliarySystem::compute(LINEAR)
 2 Nonlinear |R| = 7.277573e-17
Finished LStableDirk2 second stage solve.
 Solve Converged!
Called AuxiliarySystem::compute(TIMESTEP_END)

I am now looking more closely at the VariableTimeIntegrationAux class itself. It's possible that it is using the "wrong" _dt in its getIntegralValue(). As I understand how this AuxKernel works, it needs to integrate from zero to the "current" time, which in a multi-stage TimeIntegrator means the current stage time, not the end of the current timestep, so I'm checking on that now.

[jwpeterson]

I am now looking more closely at the VariableTimeIntegrationAux class

I think this class might be the problem. After putting a print statement in this class, I see:

In VariableTimeIntegrationAux::getIntegralValue(), using _dt= 0.5

every time. This corresponds to the dt I'm setting in the input file, but I don't think it's correct for multistage TimeIntegrators like LStableDirk2 which evaluate the first stage at time t_{alpha} = t_n + alpha * dt (i.e take an initial step of size alpha * dt).

For this AuxKernel, which is accumulating the integral of t * dt numerically, we should probably see the value of the current stage timestep, which is either alpha*dt or (1-alpha)*dt depending on which stage you are in.

[YaqiWang]

I tried to save old time mass matrix in Crank-Nicolson integrator in the above commit and hope to restore the convergence even we have time dependent quantities in time kernel, but failed to make it work with mesh adaptation. Must be missing something. I think the idea should be ok. Anyone can see any obvious issue in my commit?

[YaqiWang]

The idea is based on @jwpeterson 's thesis
Crank_Nicolson.pdf

[jwpeterson]

but failed to make it work with mesh adaptation.

Hmm... this probably isn't surprising. The matrix you saved is now the wrong size if you subsequently adapted the mesh. The matrix would need to be "projected" onto the new mesh, similarly to the way vectors are projected, but projecting matrices isn't really a thing...

[YaqiWang]

I think when mesh is adapted, equation system reinitialization will reinit all the matrix including the one I added. It also project the solution onto the new mesh. So I can call computeJacobian to evaluate the Jacobian for the new mesh. But somehow the newly calculated Jacobian, which should be the same as the Jacobian on the new mesh in the new time step for this particular problem, does not give the time residual when multiplied with udot.

[jwpeterson]

including the one I added.

Oh, that may be right, if you added it to the system with an add_matrix() call.

Before you get too far down this path though, remember that it still won't necessarily give you 2nd-order convergence (Eq. B-2.8) unless the dependence of the mass matrix on time is fairly weak (Eq. B-2.9).

[YaqiWang]

Oh. @jwpeterson do you have any experience on handling the time dependent stabilization parameter in SUPG? I guess, I will just stay with Euler scheme for easy life...

[YaqiWang]

Ok, I will stop working on this. I cannot figure out why assembling time matrix and multiplying it with udot, does not give the same the same time residual with mesh being adapted. Updated my commit to make the issue more clearer.

[jwpeterson]

do you have any experience on handling the time dependent stabilization parameter in SUPG?

You mean with regard to verifying second-order accuracy in time? I have not looked at it in detail in MOOSE-base codes, as first-order convergence or just time-stepping to steady state was usually good enough for my purposes. In traditional libMesh codes or the FEMSystem stuff, it is possible to assemble the different terms at different points in time, and I believe I used that to implement the version of Crank-Nicolson in (B-2.16), but that was back during my dissertation work...

[YaqiWang]

Yes. I will just use Euler or our current Crank-Nicolson. MOOSE Crank-Nicolson is not wrong though just could potentially lose convergence order. Thanks.

[jwpeterson]

Just to summarize my findings regarding the issues with combining LStableDirk2 and VariableTimeIntegrationAux so far. There are at least two problems:

1. The value of "_dt" used in VariableTimeIntegrationAux::getIntegralValue() corresponds to the "full" timestep size, but in this context the "stage" timestep should be used while accumulating the integral.
2. We don't call _fe_problem.advanceState() in LStableDirk2 so the old/new values of the AuxKernel are not updated between stage 1 and stage 2, and therefore the wrong AuxVariable value is used in the 2nd stage of each timestep, leading to the observed first-order global convergence.

Since ExplicitRK2 does call advanceState(), it's possible that it could work if problem 1 above is fixed. I'm testing that out now. If it works, then we may be able to fix LStableDirk2 the same way...

[jwpeterson]

The reason that the ExplicitRK2-derived classes (Heun, Ralston, ExplicitMidpoint) don't work for this test problem (BTW: remember to set implicit=false in the CoupledForce integral when running with an explicit TimeIntegrator!) is slightly different than the implicit case.

In the explicit case, the problem is that the _fe_problem.advanceState(); call causes the integral in VariableTimeIntegrationAux to accumulate (up to) twice as much as it should. In the case of Heun for instance, the first stage is an explicit Euler step of size _dt, followed by the aforementioned advanceState() call, followed by what looks like a trapezoidal integration which uses the solution from the first stage. The VariableTimeIntegrationAux therefore accumulates exactly twice as much as it is supposed to. (In the output below, we want the value returned by computeValue() to always match the value of _t, but it clearly does not.)

In VariableTimeIntegrationAux::computeValue(), returning: 2
In VariableTimeIntegrationAux::getIntegralValue(), using _dt= 0.5
In VariableTimeIntegrationAux::getIntegralValue(), using _t= 1
_fe_problem.timeOld()=0.5

[jwpeterson]

I'm not really sure how to fix all these issues. I think it's clear that VariableTimeIntegrationAux needs to be redesigned somehow.

One idea is to split its functionality into 2 classes, one of which can compute \int_{t_n}^{t_alpha} f(t) dt for any t_alpha within the current step, but doesn't keep a running sum and so doesn't get messed up by advanceState() calls? Need to think about this further...

[YaqiWang]

Looks like the root cause is the independent time integration in aux kernels, which users are responsible for, and the time integrator. I am not even sure if it is possible to have a consistent time integration in aux kernel with the time integrator.

[jwpeterson]

Hi @YaqiWang, yes we were discussing this over here as well, and agree with you. Our current idea is that one probably has to use the full "backup & restore" capability (like what is used in MultiApps) to get the multistage TimeIntegrators to work correctly.

In the meantime, I will be adding a mooseWarning() to the various multistage TimeIntegrators to let people know that they are known not to work with Materials and AuxKernels that accumulate time-dependent "state" values, and they should currently be considered experimental in general...

[YaqiWang]

Sad part is that in Rattlesnake, we often use aux kernels to evaluate delayed neutron precursors since they do not have diffusion or convection. Also we often treat temperature equation without conduction (diffusion) thus use aux kernel to evaluate the temperature. Things can get more complicated in multiphysics. I remembered I once did coupled tensor mechanics and radiation transport, where I was using MOOSE time integrators and using Newmark for tensor mechanics time kernels. We possible just need to disable coupledOlder in AuxKernel and Material to make our life easier.

[YaqiWang]

Maybe we should have an aux kernel like TimeAuxKernel, which evaluate _u_dot[_qp] instead of _u[_qp], then time integrator can use the evaluated _u_dot to update _u. If TimeAuxKernel couples _u, we will have to have local iterations to resolve it. Similarly for stateful material property, we have declarePropertyDot. getMaterialPropertyOld and getMaterialPropertyDot can still exist. Stateful postprocessors can be handled in the same way.