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
------------------------------------------------------------------------------
Everyone hates slow websites. So do we.
Make your web apps faster with AppDynamics
Download AppDynamics Lite for free today:
http://p.sf.net/sfu/appdyn_d2d_feb
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
