Fix issue in warm start.

[tupek2]

For linear problems with body or traction loads, the warm start was not including the initial residual contributions. Then we were solving the linear system again with the traction loads appropriately included in the nonlinear solver after. This was doubling the cost of some solves.

This PR also simplifies the warm start logic by assuming its only the displacement boundary conditions that may lead to element inversion. Parameter changes are simply directly included in the new residual evaluation. Always the previous stiffness is used.

An option to disable the warm start was added. This is recommended for linear problems, and perhaps even for problems with relatively small deformations.

[samuelpmishLLNL]

It looks like our documentation of this "warm start" part of the code isn't documented very well in the sourcefiles or the theory reference, so I'll try to write down my understanding of what's going on here for future reference.

---

We assume that at some point in the past, we were given $p_n$ (parameters), $U_n$ (essential b.c. values) and time $t_n$, and found $u_n$ satisfying

$$ r(u_n, p_n, U_n, t_n) = 0 $$

Now, we want to solve

$$ r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) = 0 $$

for $u_{n+1}$, given new values of parameters, essential b.c.s and time. The problem is that if we use $u_n$ as the initial guess for this new solve, most nonlinear solver algorithms will start off by linearizing at (or near) the initial guess. But, if the essential boundary conditions change by an amount on the order of the mesh size, then it's possible to invert elements and make that linearization point inadmissible (either because it converges slowly or that the inverted elements crash the program).

So, we need a better initial guess. This "warm start" generates a guess by linear extrapolation from the previous known solution:

$$ 0 = r(u_{n+1}, p_{n+1}, U_{n+1}, t_{n+1}) \approx \cancel{r(u_n, p_n, U_n, t_n)} + \frac{dr_n}{du} \Delta u + \frac{dr_n}{dp} \Delta p + \frac{dr_n}{dU} \Delta U + \frac{dr_n}{dt} \Delta t $$

Move all the known terms to the rhs and solve for $\Delta u$,

$$ \Delta u = - \bigg( \frac{dr_n}{du} \bigg)^{-1} \bigg(\frac{dr_n}{dp} \Delta p + \frac{dr_n}{dU} \Delta U + \frac{dr_n}{dt} \Delta t \bigg) $$

This assumed that the old residual was converged ($r_n = 0$), which means that the change in this PR:

$$ \Delta u = - \bigg(\frac{dr_n}{du} \bigg)^{-1} \bigg(r_n + \frac{dr_n}{dp} \Delta p + \frac{dr_n}{dU} \Delta U + \frac{dr_n}{dt} \Delta t \bigg) $$

is almost certainly more stable, since there were small errors which were previously being accumulated over time (although for reasonable solver tolerances, this was likely not a noticeable issue at the timescales we've been considering).

However, if the initial conditions were not close to equilibrium then the error would be significant.

[tupek2]

Yes, this looks correct. I'll add that an additional property of the warm start I'd like to see is:

For linear problems, the warm start should solve the system exactly.

In the case of a single static linear mechanics solve, there is no previous solution, so the initial residual will certainly be significant.

[tupek2]

For the case of dynamics, we also should probably be including the inertia/mass term. Similarly for contact, contact stiffness from the previously converged (stable) solution should be included in the warm start calculation,

[tupek2]

I don't think we are computing dr/dt * dt? So, this pr might be wrong still for time-varying body/Neumann loads.
I think we may want to evaluate

r(u^n, p^n, t+dt) instead of r(u^n, p^n, t) + dr/dt * dt. Similarly, we could change to
r(u^n, p^{n+1}, t+dt) instead of r^n + dr/dt * dt + dr/dp * dp.

[samuelpmishLLNL]

an additional property of the warm start I'd like to see is [that] for linear problems, the warm start should solve the system exactly.

That sounds like a reasonable expectation, although it potentially conflicts slightly with Kenny's "always do a nonlinear solve" request.

For the case of dynamics, we also should probably be including the inertia/mass term.

It's also worth noting that for the case of dynamics, there is already a natural choice of predictor from the time integrator itself.

instead of r(u^n, p^n, t^n) + dr/dt * dt, we just use r(u^n, p^n, t^{n+1}). We could change it to be
r(u^n, p^{n+1}, t^{n+1}) and get rid of the parameter sensitivities as well. This would be assuming that transient parameter changes dont cause element inversion. I slightly prefer this change, but I'm not set on it.

We could try that too. The reason for linearizing at the previous step is twofold:

• it's safe, since you know it's a "converged" state without inverted elements
• you don't need to recalculate the residual or jacobian (although it looks like the implementation may not be taking advantage of this aspect)

I don't think we are computing dr/dt * dt?

We're not right now, but we should be in general (I'm okay with omitting this term for now).

[tupek2]

The Jacobian for sure needs to always be evaluated at u^n, p^n, t^n (and I agree we should be using a previous one, if it exists). The residual has some more flexibility without strictly hurting the order of accuracy of the warm start, e.g., evaluating at u^n, p^{n+1}, t^{n+1}.

[btalamini]

• I also agree with Sam's math
• Agree that the question of using the warm start in dynamics is unresolved. I've tried both the Newmark predictors and a warm start in the past, but I don't have any hard data on relative effectiveness.
• However, I do think that if the warm start is activated for dynamics, it should include the inertia
• As for Kenny's expectation, I think it's just that the configuration be updated if there's any change to the inputs. This does not necessarily require a nonlinear solve, though it does with our current implementation. See next point.
• As for warm start being used as the main solver for linear problems: we need to be a little careful. Linear in which variables? The warm start doesn't (and shouldn't) update internal variables. We need to make sure that happens, which we currently do in the nonlinear solver step.

[samuelpmishLLNL]

A thought from yesterday: the correctness issue fixed in this PR was related to the assumption that the initial configuration was in equilibrium. This was a bad assumption because it puts the burden on the user to figure that out (and even worse, it wasn't even documented).

There is a related assumption for the dynamics problem, that the initial velocity values should be consistent with the time rate of the prescribed essential displacements. Is there something similar we should do to handle that case as well (not necessarily in this PR)? Maybe that case is easier and the initial velocities are just overwritten by du/dt @ t == 0

[tupek2]

I'm not sure how bad the velocity and displacements being out of sync is. I'm guessing the velocities are essentially ignored on the dirichlet nodes, but maybe not. I suppose a material model with viscosity might be looking at the wrong velocity in cases like that? Regarding what to do for dynamics, I think we need to sketch out the math a bit more. I think the warm start may still be appropriate to use there, as its really just a linear extrapolation, which may still be valid for dynamics. But I need to write it out.