Incorrect x_dot (velocity) field with some steppers

[ikalash]

There is an issue discovered recently where the x_dot (velocity) fields computed in Albany/written to Exodus file are incorrect. There seem to be 2 instances of this problem:

1.) The entire x_dot field is 0, when it should be non-zero.
2.) The x_dot field appears correct except at the boundary where a transient Dirichlet BC is applied.

Issue 1.) can be reproduced by running https://github.com/sandialabs/LCM/blob/main/tests/TransientHeat2D/input_be.yaml and https://github.com/sandialabs/LCM/blob/main/tests/TransientHeat2D/input_rk.yaml for the TransientHeat2D problem, where the DBCs applied are all 0 except at the bottom of the domain, where a 1*t time-dependent DBC is prescribed. The x_dot should therefore be 1.0 for all times at the bottom boundary. The correct behavior is observed when using a Backward Euler stepper, but not when using an RK Implicit Midpoint stepper, as shown below.

Issue 2.) can be reproduced by running the LCM/Cuboid problem https://github.com/sandialabs/LCM/blob/main/tests/LCM/Cuboid/cuboid_newmark_tempus.yaml . With the current mesh utilized, the z-velocity field (x_dot) jumps from 0 (top figure below) to 2 (bottom figure below) between time-steps. Since again we have a 1*t time-dependent DBC prescribed, this time at the top of the domain, again x_dot should be 1.0 for all times at the top boundary. This seems to happen independent of which 2nd-order-in-time stepper is used (all variants of Newmark in Tempus, Trapezoidal Rule in Piro). A "fix" is to specify an initial condition that is consistent with the DBC being imposed. However, the DBC value should trump the initial condition and will in general not be constant or known, so this "fix" is not viable in general.

The same issues are present in Albany main, not just LCM. After a lot of detective work, I don't see how x_dot could be modified on the Albany side; it appears to be modified on the Trilinos side, even for issue 2.) (which is a bit surprising since the solvers don't know anything about the DBCs/DBC boundaries and the problem happens for all steppers). For 1.), it seems likely that Tempus is not calculating x_dot correctly for certain steppers.

[ikalash]

I have looked into this a bit more for the Newmark scheme.

For the d-form of Newmark, we have:

dPred = dt*v + dt*dt/2*(1-2*beta)*a + d
vPred = v + dt*(1-gamma)*a 
solve for d (*)
a = 1/beta/dt/dt*(d – dPred) (**)
v = vPred + dt*gamma*a  (***)

What should happen is d, v and a should be solved for only at the unconstrained dofs. The problem is Piro and Tempus have no way to know which dofs are fixed vs. free so they will operate on all the dofs. Consider t = dt, one step into our scheme, and suppose () gave a solution d that satisfies the DBCs. If the DBC is 1t, we will have d[fixed] = dt. Suppose dPred[fixed] = 0 and vPred[fixed] = 0. Then we have, from () and (*) above:

a[fixed] = 1/beta/dt/dt*dt = 1/beta/dt 
v[fixed] = dt*gamma*a[fixed] = dt*gamma/beta/dt = gamma/beta = 0.5/0.25 = 2.0 

This is where the 2.0 in Albany is coming from for xDot at the top boundary.

I think this can be fixed by setting the xdot values to the prescribed ones in Albany. This could be done in the Albany ME, or perhaps elsewhere. One would still have to do some bookkeeping and try to get the relevant info about which DOFs are free vs. constrained to the right place. Also, one would need a way to estimate the correct value of xdot at the constrained DOFs from a BC on x (can it be done using AD for an analytically given x?).

For the other steppers where xdot is always 0 , I think the problem is in Tempus - xdot is not being updated.

[ikalash]

Another update on this issue.

I have talked to @ccober6 about 1.) above. Basically the reason this is happening is because xdot is not implemented for some of the time-steppers. There is no funding to fix this right now on the Tempus side. If it becomes necessary to have xdot for all time-steppers, we can ping @ccober6.

Regarding 2.): it seems feasible to either calculate xdot (and xdotdot) in the using symbolic differentiation here https://github.com/sandialabs/LCM/blob/main/src/evaluators/bc/PHAL_ExprEvalSDBC_Def.hpp or require the user to input an analytic form for xdot (and xdotdot) in their input file when running a dynamic problem. Unfortunately this would not extend to the regular time-dependent DBC where you enter data rather than an analytic function. To handle that case, finite differences would need to be introduced, which is apparently how the DBC values for xdot (and xdotdot) are set in SIERRA.