Open YaqiWang opened 6 years ago
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
[]
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.
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 alphadt or (1-alpha)*dt depending on which stage you are in.
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?
The idea is based on @jwpeterson 's thesis Crank_Nicolson.pdf
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...
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.
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).
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...
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.
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...
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.
Just to summarize my findings regarding the issues with combining LStableDirk2
and VariableTimeIntegrationAux
so far. There are at least two problems:
VariableTimeIntegrationAux::getIntegralValue()
corresponds to the "full" timestep size, but in this context the "stage" timestep should be used while accumulating the integral._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...
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
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...
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.
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...
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.
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.
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:
Solution of
u
ist*t/2
. No spatial discretization error should be seen. Thus the exactintu
is 0.5.Description
Backward Euler: error is
dt/2
, first order; Forward Euler: error isdt/2
, first order; BDF2: error is about3/4*dt^2
, second order; Crank-Nicolson: error is zero; Explicit-midpoint: converges to wrong value 1; DIRK: error isdt/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.