From: Nasser Mohieddin Abukhdeir <nasser.abukhdeir@mc...>  20080519 18:55:11

Hello: A while ago I mentioned something about the time stepping algorithm employed by AdaptiveTimeSolver as being nonideal for my PDE system. I think I have implemented a very simple adaptive time solver based on the local truncation error: core_time_solver>du() / calculate_norm(_system, *_system.solution) where it simply increases the time step up to a default of 2 (max_growth=2). The only advantage of this approach for me application is that the default algorithm is expensive (three nonlinear steps for one effective time step) and, while locally it gives a good estimate of the error, is too optimistic with its adaptive time stepping (my PDEs are highly nonlinear). If this makes any sense for inclusion into LibMesh's AdaptiveTimeSolver I can cleanup and provide the code for review.  Nasser Mohieddin Abukhdeir Graduate Student (A.D. Rey Research Group) McGill University  Department of Chemical Engineering http://webpages.mcgill.ca/students/nabukh/web/ http://people.mcgill.ca/alejandro.rey/ 
From: Nasser Mohieddin Abukhdeir <nasser.abukhdeir@mc...>  20080519 18:55:11

Hello: A while ago I mentioned something about the time stepping algorithm employed by AdaptiveTimeSolver as being nonideal for my PDE system. I think I have implemented a very simple adaptive time solver based on the local truncation error: core_time_solver>du() / calculate_norm(_system, *_system.solution) where it simply increases the time step up to a default of 2 (max_growth=2). The only advantage of this approach for me application is that the default algorithm is expensive (three nonlinear steps for one effective time step) and, while locally it gives a good estimate of the error, is too optimistic with its adaptive time stepping (my PDEs are highly nonlinear). If this makes any sense for inclusion into LibMesh's AdaptiveTimeSolver I can cleanup and provide the code for review.  Nasser Mohieddin Abukhdeir Graduate Student (A.D. Rey Research Group) McGill University  Department of Chemical Engineering http://webpages.mcgill.ca/students/nabukh/web/ http://people.mcgill.ca/alejandro.rey/ 
From: John Peterson <jwpeterson@gm...>  20080519 19:01:49

On Mon, May 19, 2008 at 1:55 PM, Nasser Mohieddin Abukhdeir <nasser.abukhdeir@...> wrote: > Hello: > A while ago I mentioned something about the time stepping algorithm > employed by AdaptiveTimeSolver as being nonideal for my PDE system. I > think I have implemented a very simple adaptive time solver based on the > local truncation error: > > core_time_solver>du() / calculate_norm(_system, *_system.solution) As it's currently implemented, du() is the difference between two successive timesteps, not a measure of the truncation error. While this may be working for your particular application, I don't know of a theoretical justification for this particular method of adaptive timestep selection. J > where it simply increases the time step up to a default of 2 > (max_growth=2). The only advantage of this approach for me application > is that the default algorithm is expensive (three nonlinear steps for > one effective time step) and, while locally it gives a good estimate of > the error, is too optimistic with its adaptive time stepping (my PDEs > are highly nonlinear). > If this makes any sense for inclusion into LibMesh's > AdaptiveTimeSolver I can cleanup and provide the code for review. > >  > Nasser Mohieddin Abukhdeir > Graduate Student (A.D. Rey Research Group) > McGill University  Department of Chemical Engineering > http://webpages.mcgill.ca/students/nabukh/web/ > http://people.mcgill.ca/alejandro.rey/ > > >  > This SF.net email is sponsored by: Microsoft > Defy all challenges. Microsoft(R) Visual Studio 2008. > http://clk.atdmt.com/MRT/go/vse0120000070mrt/direct/01/ > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > 
From: Roy Stogner <roy@st...>  20080519 19:10:43

On Mon, 19 May 2008, Nasser Mohieddin Abukhdeir wrote: > A while ago I mentioned something about the time stepping algorithm > employed by AdaptiveTimeSolver as being nonideal for my PDE system. I > think I have implemented a very simple adaptive time solver based on the > local truncation error: > > core_time_solver>du() / calculate_norm(_system, *_system.solution) This isn't the local truncation error; it will be nonzero even if your solver and system are chosen in such a way that the time integration error is zero. > where it simply increases the time step up to a default of 2 > (max_growth=2). The only advantage of this approach for me application > is that the default algorithm is expensive (three nonlinear steps for > one effective time step) That's three nonlinear steps for two effective time steps, but yes, that's a reasonable complaint. > and, while locally it gives a good estimate of the error, is too > optimistic with its adaptive time stepping (my PDEs are highly > nonlinear). Also a reasonable complaint; I occasionally see the time stepper having to backtrack on an extremely nonlinear system. > If this makes any sense for inclusion into LibMesh's > AdaptiveTimeSolver I can cleanup and provide the code for review. That would be great, thanks! We've currently got the global_tolerance flag for use in choosing between a couple different refinement heuristics; another boolean would probably be the only interface change we'd need to add for your method, right? We might want to refactor that into an enum, a strategy pattern, subclasses, or something like that eventually, but that's what the untested() macro is for... If you've got time to play around with time solver adaptivity, one thing on our todo list is predictorcorrector schemes. Pick an explicit and an implicit scheme of the same order, use the difference between them to estimate your truncation error, then use the implicit scheme for your solution at the next time step to keep stability. That can actually end up being faster than just the implicit scheme alone, since the explicit result should give a better initial guess for your solver.  Roy 
From: Nasser Mohieddin Abukhdeir <nasser.abukhdeir@mc...>  20080519 20:59:58

Roy Stogner wrote: > > On Mon, 19 May 2008, Nasser Mohieddin Abukhdeir wrote: > >> A while ago I mentioned something about the time stepping algorithm >> employed by AdaptiveTimeSolver as being nonideal for my PDE system. I >> think I have implemented a very simple adaptive time solver based on the >> local truncation error: >> >> core_time_solver>du() / calculate_norm(_system, *_system.solution) > > This isn't the local truncation error; it will be nonzero even if > your solver and system are chosen in such a way that the time > integration error is zero. > John pointed this out too, all I can say is "oops," what meant to do was take the error from the last iteration of nonlinear solver which I think makes sense as an estimate of the LTE. I'll go back to the drawing board for this taking into account your following suggestions (I am currently using an enum and default to the original AdaptiveTimeSolver method). > If this makes any sense for inclusion into LibMesh's >> AdaptiveTimeSolver I can cleanup and provide the code for review. > > That would be great, thanks! We've currently got the global_tolerance > flag for use in choosing between a couple different refinement > heuristics; another boolean would probably be the only interface > change we'd need to add for your method, right? We might want to > refactor that into an enum, a strategy pattern, subclasses, or > something like that eventually, but that's what the untested() macro > is for... > > If you've got time to play around with time solver adaptivity, one > thing on our todo list is predictorcorrector schemes. Pick an > explicit and an implicit scheme of the same order, use the difference > between them to estimate your truncation error, then use the implicit > scheme for your solution at the next time step to keep stability. > That can actually end up being faster than just the implicit scheme > alone, since the explicit result should give a better initial guess > for your solver. >  > Roy would this be as simple as just adding an additional UnsteadySolver member to AdaptiveTimeSolver, so for example core_time_solver could be EulerSolver with theta=1 and (let's call it) aux_time_solver would then be another EulerSolver with theta=0. We then just march them in parallel: 1) run the explicit scheme using aux_time_solver 2) use the solution from aux_time_solver to seed the implicit scheme (core_time_solver) and then solve 3) estimate the LTE from the difference of the two solutions and calculate the new timestep 4) synchronize the current solution of the explicit solver (to that of the implicit solution)  Nasser Mohieddin Abukhdeir Graduate Student (A.D. Rey Research Group) McGill University  Department of Chemical Engineering http://webpages.mcgill.ca/students/nabukh/web/ http://people.mcgill.ca/alejandro.rey/ 
From: Roy Stogner <roy@st...>  20080519 21:36:48

On Mon, 19 May 2008, Nasser Mohieddin Abukhdeir wrote: > John pointed this out too, all I can say is "oops," what meant to do was take > the error from the last iteration of nonlinear solver which I think makes > sense as an estimate of the LTE. I'm not sure of that, either. The error on your last algebraic nonlinear solver iteration should depend on those solver tolerances, which don't necessarily have anything to do with your time discretization error. > would this be as simple as just adding an additional UnsteadySolver member to > AdaptiveTimeSolver, so for example core_time_solver could be EulerSolver with > theta=1 and (let's call it) aux_time_solver would then be another EulerSolver > with theta=0. At this point we've differed enough from the original AdaptiveTimeSolver design that I'd refactor the class tree. Maybe make AdaptiveTimeSolver a branch class, with one leaf for the current behavior (with what name?) and another for the new PredictorCorrectorTimeSolver. Then you could rename core_time_solver as corrector_time_solver and call aux_time_solver predictor_time_solver, too. > We then just march them in parallel: > > 1) run the explicit scheme using aux_time_solver Right. I'd want the default diff_solver settings for that to be appropriate for linear problems, too; for forward Euler we're just solving a mass matrix. > 2) use the solution from aux_time_solver to seed the implicit scheme > (core_time_solver) and then solve Right. And since the solution is stored with the system, reusing it should be as simple as "doing nothing". ;) > 3) estimate the LTE from the difference of the two solutions and calculate > the new timestep Right. > 4) synchronize the current solution of the explicit solver (to that of the > implicit solution) Again, "doing nothing" (other than the advance_timestep() you'd do anyway). ;) The only catch that comes to mind is that the time solvers are also responsible for initializing old DiffSystem solution vectors. That's no trouble for forward vs backward Euler, but if you wanted to write an AdamsBashforth predictor for CrankNicholson, you'd have to be careful that the former gets used in initialization since it needs to save more old vectors. That's probably a safe rule of thumb: "use the predictor for initialization"; unless someone tries something weird, the corrector solver is always going to need the same or less old information.  Roy 
From: Nasser Mohieddin Abukhdeir <nasser.abukhdeir@mc...>  20080520 15:26:19

Sounds like a plan, you'll have to figure out what to call the current implementation, I am bad with names. I'm working on a LibMeshbased implementation of one of my liquid crystal models (see my web page if you're interested). Once I have the basics of this worked out, I'll get started on implementing these changes you mentioned in the last email. > > >> would this be as simple as just adding an additional UnsteadySolver >> member to AdaptiveTimeSolver, so for example core_time_solver could >> be EulerSolver with theta=1 and (let's call it) aux_time_solver would >> then be another EulerSolver with theta=0. > > At this point we've differed enough from the original > AdaptiveTimeSolver design that I'd refactor the class tree. Maybe > make AdaptiveTimeSolver a branch class, with one leaf for the current > behavior (with what name?) and another for the new > PredictorCorrectorTimeSolver. Then you could rename core_time_solver > as corrector_time_solver and call aux_time_solver > predictor_time_solver, too.  Nasser Mohieddin Abukhdeir Graduate Student (A.D. Rey Research Group) McGill University  Department of Chemical Engineering http://webpages.mcgill.ca/students/nabukh/web/ http://people.mcgill.ca/alejandro.rey/ 