Stateful physics error when execute_on nonlinear

[zachmprince]

Bug Description

There was a recent libmesh implementation that detects if there are stateful quantities in a problem. This is really useful to detect cyclic dependencies on auxkernels, postprocessors, materials, etc. However, it brought up an issue when trying to evaluate objects on nonlinear. It seems this is because objects with execute_on = nonlinear are evaluated when the Jacobian is updated, i.e. after nonlinear residual computation but before the first linear residual computation. This causes a mismatch in these two residuals and triggers the aforementioned error, despite there being no cyclic dependencies.

Steps to Reproduce

Compiling in dbg and running

./moose_test-dbg -i tests/auxkernels/element_aux_var/element_aux_var_test.i AuxKernels/constant/execute_on=nonlinear

causes the following error:

 0 Nonlinear |R| = 1.224745e+00
      0 Linear |R| = 1.224745e+00
You requested to reuse the nonlinear residual vector as the base vector for computing the action of the matrix-free Jacobian, but the vectors are not the same. Your physics must have states; either remove the states from your code or make sure that you set_mf_reuse_base(false)

Impact

I am unsure of the implications when trying to reuse the base vector during this type of simulation, but it is quite annoying to have this error when trying to debug models with objects that have execute_on = nonlinear. It can be avoided by adding snesmf_reuse_base = false in the executioner, but this results in a non-optimal solve.

Tagging @YaqiWang @lindsayad

[lindsayad]

I am unsure of the implications when trying to reuse the base vector during this type of simulation, but it is quite annoying to have this error when trying to debug models with objects that have execute_on = nonlinear. It can be avoided by adding snesmf_reuse_base = false in the executioner, but this results in a non-optimal solve.

What do you mean by non-optimal? The error message is actually being really helpful here...the evaluation of your objects on execute_on = nonlinear has changed the state of your function and I guess more relevantly here, the state of your Jacobian/preconditioner. It would be wrong to reuse the nonlinear residual as the base vector for approximating the Jacobian action; you've changed what the Jacobian "is" since the nonlinear residual evaluation and so it's more accurate to recompute the function in order to get an accurate computation of the base vector for approximating the Jacobian action.

Perhaps what you really want and what would eliminate the error message you're seeing is that execute_on = nonlinear would actually execute on nonlinear residual evaluation instead of during the Jacobian evaluation. And then your nonlinear residual and jacobian would actually match up (and then also your base vector would be reusable!)

[zachmprince]

I completely agree, the error message was very helpful and not reusing the base vector is right thing to do with the current implementation. Ideally, yes, the objects would be computed on the nonlinear residual evaluation. To me, there are three options: 1) we change when EXEC_NONLINEAR is called, but I don't know if this goes against what EXEC_NONLINEAR means as whenever the Jacobian is updated. 2) We automatically set snesmf_reuse_base = false when the system detects that there are objects executed on nonlinear. 3) Leave the current implementation, assume that the user knows what they are doing with the Krylov solve, and have an option to disable the check. What would you think is the best course of action?

[zachmprince]

Actually, I don't really like the third option, so disregard that one. Also, thinking about it, there would be cases that base vector can be used even if objects are executed on nonlinear, say if those objects only affect the Jacobian and not the residual. This would be the case where the current EXEC_NONLINEAR is appropriate as well. This is tricky.

[lindsayad]

In terms of the math, there should never be an object that affects the Jacobian and not the residual, right? The Jacobian is just the derivative of the function. The function is king. Now of course in MOOSE, the execution of nonlinear on Jacobian doesn't follow the math, which is sad and unfortunate. Back when we created execute_on = nonlinear, we didn't have a way to tell apart nonlinear and linear residual evaluations, so doing nonlinear on Jacobian was the best we could do. However, we can now tell them apart, so I would love to see execute_on = nonlinear moved to nonlinear residual evaluation.

[lindsayad]

The only problem with moving this to nonlinear residual evaluation, is that presumably the user has selected execute_on = nonlinear for performance reasons, e.g. they don't want it done for every residual evaluation in PJFNK. If we moved nonlinear to mean nonlinear residual evaluation, we would get additional evaluations when we employ a line search that we don't when we execute nonlinear on Jacobian.

[zachmprince]

Okay, that's great! I will start thinking about where and how to move the evaluation. For the line search, fundamentally, do we want to execute a "nonlinear" object on a line search evaluation? It seems like the answer would be yes, so that those objects are updated for next set of linear iterations.

[lindsayad]

In terms of correctness, you would definitely want it evaluated on every line search evaluation. But also in terms of correctness, it's unclear to me why you would want to just evaluate an object on nonlinear, as opposed to both nonlinear and linear. We do this conceptually in mechanical contact (e.g. only update the contact set on nonlinear), but I don't like it. I believe we should either be using a semi-smooth method where the function is independent of whether you're evaluating on nonlinear or linear or maybe a reduced active set method that doesn't have function call-backs during linear evaluation.

Anyway we don't need to talk about contact. Why do you guys just execute objects on nonlinear?

[zachmprince]

Sounds like our reasoning is pretty similar to the one you gave about contact. We have transient simulation that has neutronics coupled with temperature feedback. It is a adiabatic heat up model, so the temperature can be solved in an auxkernel. The transient starts with very low power, so almost no temperature feedback through a lot of the simulation and it would be wasteful to be evaluating that auxkernel and recalculating materials on linear iterations during this part. The compromise was to evaluate this auxkernel on nonlinear so it would still be somewhat strongly coupled when the power increases. @YaqiWang might have a better explanation.

[YaqiWang]

execute_on=nonlinear is like mixed Picard iteration with Newton. Let us say you have a nonlinear problem, with two sources of nonlinearity from two user objects. You can make both user objects on timestep_end, solving a linear problem and two UOs with Picard iterations. If we put the two UOs on nonlinear, we are basically doing Picard iteration, where each Newton iteration is equivalent with Picard. If we put both on linear, we are doing normal Newton solve. If we put one on linear, another on nonlinear we are doing a hybrid one, it is neither Newton or Picard. This can result into a faster solve. The final solution is dictated by the nonlinear residual. Because of this, I think the nonlinear object should not be done in line search.

[lindsayad]

We should have execute_on = nonlinear_residual and execute_on = jacobian. And of course execute_on = picard_end