From: Roy Stogner <roystgnr@ic...>  20130103 17:21:58

To do a decent stabilized problem with nonlinearities on the time derivative term, our mass_residual() functions ought to have access to the current element solution u(t) at the simulationtime being evaluated, a fixed element solution u_f which is constant from t_n to t_{n+1}, and some measure of the rate of change du/dt at t. We currently pack the rate of change (as deltat * du/dt) into elem_solution (a woefully confusing decision originally made for backwards compatibility reasons), we pack a u_f into elem_fixed_solution, and nonlinear mass_residual calls can reuse elem_fixed_solution as u(t) so long as they're not doing any time integration more complicated than the singleevaluation theta method. We currently can't add more accurate time integrators (or even use a simple trapezoidal rule!) on problems with nonlinear mass matrices without the results being inaccurate. I'd like to change this, by: Adding an elem_solution_rate (or elem_solution_derivative? opinions welcome) to FEMContext which holds du/dt. Making elem_solution consistently hold u(t). The problem is that this will thoroughly break the results of anyone with a FEMSystembased code complicated enough to require overriding mass_residual(). I've talked to the few developers who I personally know to be in such a situation, and so far there's unanimous support for "just break it and we'll update", but if any other users disagree I'd like to know about it. Adding a "transition release" with forwards and backwards compatibility for this change would be a pain in the neck, but it can be done if someone really wants it.  Roy 