From: Roy Stogner <roystgnr@ic...>  20070216 14:05:41

On Fri, 16 Feb 2007, Ingo Schmidt wrote: > a) Ex 13 is for Reynolds number equal to one (?), i.e. Ex 18 is slightly > more general? (Because I miss some material specific data in ex 13) Yes, that's right. > To handle the problem with that first and second time derivatives in > FEMSystem I'm thinking about extending that time_evolvingswitch > including the order of the time derivative. The time_stepping scheme > should than be able due to that switch to chose the right approximation > order for the different variables. Additional memory for the nodal > values of velocity and acceleration (required for some time stepping > schemes, e.g. Newmark) have to be allocated, e.g. with system.add_vector(). Here's my best current idea: Ignoring the *_constraint functions for now, how you describe a timedependent problem to FEMSystem is basically: M(u) du/dt = F(u) Where the time_evolving switch tells the system that a du/dt term is involved, and the mass_residual term handles the effects of the mass matrix M(u). Most applications, ex18 included, don't redefine mass_residual because they don't need to; the residual is usually just (du/dt,phi_i)_L2. But what if we add a "time_evolving_second_order" flag and "second_mass_residual" method (or, you know, something with better names) as well? Then for problems with second order time derivatives, you turn both flags on, and if you've got a more complicated left hand side than (d2u/dt2 + du/dt) you override the residual functions as well, and that should completely describe: K(u) d2u/dt2 + M(u) du/dt = F(u) There's one concern that comes up for nonconstant K and M  I'll illustrate it for Euler integration since I don't think I understand Newmark yet: When solving M du/dt = F(u), Euler normally does: (M u^(n+1)  M u^(n)) = Dt * F(u^(theta)), where u^(theta) is an interpolant somewhere in between u^(n) and u^(n+1) depending on the specific method. But when M isn't constant, currently you get: (M(u^(n+1)) u^(n+1)  M(u^(n)) u^(n)) on the left hand side, which depending on your formulation probably isn't what you want, and may not even be consistent. We had a student here who had to resort to a few ugly hacks to get a nonlinear stabilizationbased M working, since he needed: (M(u^(theta)) u^(n+1)  M(u^(theta)) u^(n)) If you might have the same problem then maybe it's time to make something like an "elem_theta_solution" variable available in the API. But again, hopefully with a better name. > Beyond that specific insertion concerning FEMSystem I've got a really > general corresponding to the member element_time_derivative(). > Isn't the name too specific to CFD problems or the like? The name is probably bad in general, since it really corresponds to the parts of the equation that aren't attached to any time derivatives! I was just looking for a name that seemed to accurately distinguish it from element_constraint, since my time integrators need to treat those differently. > Thirdly, the TransientSystem<TransientNonlinearImplicitSystem> following > ex13 will be the best alternative to the FEMSystem for me. But is there > a chance to include the SNESSetUpdate > (http://wwwunix.mcs.anl.gov/petsc/petscas/snapshots/petsccurrent/docs/manualpages/SNES/SNESSetUpdate.html) > Option into libmesh at the PetscNonlinearSolver.solve() ? That's a > really general but also important option necessary for residuals > depending on the "current" solution (e.g. nonlinear constitutive laws). You'll have to ask Ben about that  I think he's currently the only one who's done any work with PETSc's nonlinear solvers... and it can't have been too much work, if he's still got that handrolled Newton loop in ex13! He's probably slacking off, finishing his dissertation and graduating next month and such.  Roy 