Re: [Libmesh-users] transient system of second order

[Roy S. <roy...@ic...>]

On Thu, 21 Feb 2013, Manav Bhatia wrote:

> What is throwing me off is that the class documentation of Euler2_Solver says: 
> 
>  * Euler solves u' = f(theta*u_new + (1-theta)*u_old),
>  * Euler2 solves u' = theta*f(u_new) + (1-theta)*f(u_old)
> 
> which seems to imply that no mass matrix is accounted for.

This is a flaw in the documentation; I'll fix it ASAP.

> Also, while going through the code, I came across the following in euler_solver.C (lines 124-129)
> 
>   // We do a trick here to avoid using a non-1
>   // elem_solution_derivative:
>   context.elem_jacobian *= -1.0;
>   jacobian_computed = _system.mass_residual(jacobian_computed, context) &&
>     jacobian_computed;
>   context.elem_jacobian *= -1.0;
> 
> 
> What is the need to multiply the elem_jacobian by -1.0 both before
> and after the function call? Wouldn't the call to mass_residual
> overwrite the elem_jacobian anyways? 

The call to mass_residual (or anything else) is supposed to *augment*
the Jacobian, rather than writing to a separate jacobian that would
then just be added to the original.  This was my bright idea for
minimizing L1 cache use and is one of the "premature optimization"
examples I was talking about.

A problem with that is, if mass_residual is computing a jacobian
w.r.t. elem_solution, but elem_solution is storing something whose
derivative w.r.t. the real solution isn't 1.0, then the user has to
remember to multiply those jacobian terms by elem_solution_derivative
to take the chain rule into account.  The way we do scaling with
deltat and 1.0 bits was intended to try and avoid this complication
for most codes.
---
Roy