From: Nasser Mohieddin Abukhdeir <nasser.abukhdeir@mc...>  20080428 14:02:19

Hi: Thanks, I ended up figuring that out over the weekend, my background in C is not serving me very well with C++! So far I have run into two major difficulties which I think require me to make some code additions. Note, I am posting this to the dlist in case I am using LibMesh incorrectly, plus, I plan on eventually submitting this code to the LibMesh developers if all goes well: 1) The AdaptiveTimeSolver implementation is not working well for my system of PDEs. I think there are two reasons behind this: a) the method involves solving three nonlinear steps per overall timestep which is extremely expensive for my system (on the order of millions of DOFs) b) my system is "very" nonlinear, so the provided algorithm is too optimistic and I end up having to set the max_growth variable to limit the predicted timestep. The solution, I think, is to implement a more simple onestep algorithm based on a discrete PIcontrol method (see the reference) using the LTE estimate from the nonlinear step. Basically I would end up creating a new AdaptiveTimeSolver class where you can chose the timeadaptive algorithm and I would add the one I just mentioned plus maybe another simple onestep algorithm. 2) Using adaptive meshing is a huge advantage, but I have an issue finding a good initial mesh. I plan on implementing a simple routine that repeats the initial time step, reprojecting the initial solution on the refined mesh. This would get around the problem I have now where my initial coarse mesh has to be notsocoarse to resolve my initial solution (which involves a heaviside function). I'm not sure where this would fit in to the LibMesh class structure. ref "Adaptive Multilevel Solutions of Nonlinear Parabolic PDE Systems," Jens Lang, Springer 2001, pages 3235 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/ John Peterson wrote: > Hi, > > Since your initialization function is a member function of a class you > can't pass a pointer to it for project_solution the same way you would > for a "global" function. An easy fix is to use a global function > instead, but I believe there may be another, more correct way to > achieve what you are trying to do. > > J > > On Wed, Apr 23, 2008 at 3:14 PM, Nasser Mohieddin Abukhdeir > <nasser.abukhdeir@...> wrote: > >> Thanks for the informative answer, I moved over the DiffSystem and I >> think I worked out some incorrect math, but now I am back to square one >> with the initial condition of the transient problem. I added this line >> to the end of my Nem2DSystem:init_data () call (Nem2DSystem is the name >> of my class that inherits FEMSystem ): >> >> this>project_solution(initial_value, >> NULL,(this>get_equation_systems()).parameters); >> >> where the initial value function is defined similar to the >> nonDiffSystem example with an initial condition: >> >> Number Nem2DSystem::initial_value (const Point& p, >> const Parameters& parameters, >> const std::string& sysname, >> const std::string& unkname) >> >> and I get a compile error that I cannot figure out: >> >> nem2dsystem.C: In member function \u2018virtual void >> Nem2DSystem::init_data(): >> nem2dsystem.C:149: error: no matching function for call to >> \u2018Nem2DSystem::project_solution(<unresolved overloaded function >> type>, NULL, Parameters&) >> /home/grad/nasser/nemfem/libmesh0.6.2/include/solvers/system.h:165: >> note: candidates are: void System::project_solution(Number (*)(const >> Point&, const Parameters&, const std::string&, const std::string&), >> Gradient (*)(const Point&, const Parameters&, const std::string&, const >> std::string&), Parameters&) const >> >> I guess I have to study C++ templates a little more, but any ideas? >> >> 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/ >> >> >> >> >> >> Roy Stogner wrote: >> > >> > On Wed, 23 Apr 2008, Nasser Mohieddin Abukhdeir wrote: >> > >> >> *** Solving time step 0, time = 1e07 *** >> >> ERROR: LASPACK Error: >> >> in ILUFactor for Mat: >> >> Factorization produces zero pivot elements. >> >> [0] src/numerics/laspack_linear_solver.C, line 276, compiled Apr 21 2008 >> >> at 14:01:19 >> >> Aborted >> >> >> >> Basically I added my expressions for the residual and Jacobian and >> >> removed all of the boundary condition code because I need to impose >> >> Neuman BCs (which from what I understand is implicitely taken care of by >> >> the FEM algorithm). >> > >> > Homogeneous natural BCs are implicitly taken care of by the FEM >> > formulation, but just what the "natural boundary condition" is depends >> > on your precise formulation; if you've got something complicated, take >> > a look at what boundary integrals pop out of your integration by parts >> > to doublecheck that what you're setting to zero is what you want to >> > be setting to zero. >> > >> >> I initialized t=0 with an acceptable solution to the equations, so I >> >> just expected some smooth sailing... >> > >> > You're solving a transient problem, then? In an allNeumann problem >> > with no reaction terms you can often get a singular matrix doing a >> > steady solve (since the problem itself isn't uniquely defined); but >> > for a transient problem, the mass matrix coming from the time >> > integration should prevent anything like that from happening. >> > >> >> I am also wondering about DiffSystem, should I just scrap my work with >> >> ex13 and start with DiffSystem? >> > >> > There are tradeoffs. DiffSystem makes it easier to add some new >> > features (e.g. adaptive time stepping, more complex nonlinear >> > equations) to your code, but it also means you'll have to make >> > more changes to your code to keep up with the library. >> > >> > When you construct a DiffSystem subclass, you read: >> > "Warning, This code is untested, experimental, or likely to see future >> > API changes" >> > >> > The "untested" and "experimental" adjectives no longer apply, but >> > "future API changes" certainly does. I'm probably going to alter even >> > the most basic class methods to make it work better multithreaded. >> > The changes should be straightforward (methods will take reference >> > parameter(s) instead of using certain member variables), but it's >> > never fun when your code suddenly requires repair because a new >> > library version made backwards incompatible API changes. >> >  >> > Roy >> >>  >> This SF.net email is sponsored by the 2008 JavaOne(SM) Conference >> Don't miss this year's exciting event. There's still time to save $100. >> Use priority code J8TL2D2. >> http://ad.doubleclick.net/clk;198757673;13503038;p?http://java.sun.com/javaone >> _______________________________________________ >> Libmeshusers mailing list >> Libmeshusers@... >> https://lists.sourceforge.net/lists/listinfo/libmeshusers >> >> 