From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-17 21:01:13
|
The problem is as follows: - You have some parallel, distributed data on an old mesh - You refine the mesh, so you have new elements on your processor - The old_dof_indices for the new elements on the processor are not necessarily contained in the old vectors. The ugly hack is a multi-stage approach as follows: localize() the old vectors on each processor to a serial vector containing *all* the vector values project() these serial vectors copy only the local components into the locations of the new parallel vector. This is exactly the problem I mentioned in the cryptic midnight email I sent out the other day. John, did you test the old implementation in parallel with debugging? Anyway, I'll do a cvs update and see how I can help. -Ben -----Original Message----- From: lib...@li... [mailto:lib...@li...] On Behalf Of Roy Stogner Sent: Tuesday, May 17, 2005 3:45 PM To: lib...@li... Subject: Re: [Libmesh-devel] Unstable code: incoming On Tue, 17 May 2005, Roy Stogner wrote: > How do the add_vector-allocated vectors stand up in parallel > computations? Apparantly the answer is: badly. John's code is now working in serial, and failing an assert() in parallel. If nobody can find the problem soon, I'll revert the code tomorrow afternoon. --- Roy ------------------------------------------------------- This SF.Net email is sponsored by Oracle Space Sweepstakes Want to be the first software developer in space? Enter now for the Oracle Space Sweepstakes! http://ads.osdn.com/?ad_id=7412&alloc_id=16344&op=click _______________________________________________ Libmesh-devel mailing list Lib...@li... https://lists.sourceforge.net/lists/listinfo/libmesh-devel |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-17 21:44:34
|
The reason I ask is because I think the bug should exist in the code from, say, last Wednesday. If so, there is no point in Roy undoing his changes. We can just focus on fixing it now. -Ben -----Original Message----- From: John Peterson [mailto:pet...@cf...] Sent: Tuesday, May 17, 2005 4:05 PM To: KIRK, BENJAMIN (JSC-EG) (NASA) Cc: 'Roy Stogner'; lib...@li... Subject: RE: [Libmesh-devel] Unstable code: incoming KIRK, BENJAMIN (JSC-EG) (NASA) writes: > The problem is as follows: > > - You have some parallel, distributed data on an old mesh > - You refine the mesh, so you have new elements on your processor > - The old_dof_indices for the new elements on the processor are > not necessarily contained in the old vectors. > > The ugly hack is a multi-stage approach as follows: > > localize() the old vectors on each processor to a serial vector > containing *all* the vector values > > project() these serial vectors > > copy only the local components into the locations of the new parallel > vector. > > This is exactly the problem I mentioned in the cryptic midnight email I sent > out the other day. John, did you test the old implementation in parallel > with debugging? I sure didn't. But I'd be happy to revert to whatever date and give it a go. Actually you can reproduce this behavior in example 10 (for sure) and possibly others if you want to see for yourself... > Anyway, I'll do a cvs update and see how I can help. -John |
From: Roy S. <roy...@ic...> - 2005-05-18 03:25:20
|
On Tue, 17 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > The reason I ask is because I think the bug should exist in the code from, > say, last Wednesday. If so, there is no point in Roy undoing his changes. > We can just focus on fixing it now. You're making the kind assumption that there's only that one bug in the code. :-( Even assuming your code doesn't hit the "non-Lagrange coarsening" code path, the new System::reinit isn't behaving as I hoped. For one example, I'd assumed that if refine_and_coarsen_elements() has been called, the following coarsen_elements() and refine_elements() calls would be no-ops. Instead, it looks like it's possible for the make_{coarsening,smoothing}_compatible calls in coarsen_elements() to decide to add an Elem::REFINE flag, which triggers the refine_elements()! I think I can get everything working in the serial case tomorrow, but anyone who doesn't want to futz with my bugs might want to revert to Monday's system.C and mesh_refinement.C in the meantime... --- Roy |
From: Roy S. <roy...@ic...> - 2005-05-18 14:00:26
|
On Tue, 17 May 2005, Roy Stogner wrote: > I think I can get everything working in the serial case tomorrow, but > anyone who doesn't want to futz with my bugs might want to revert to > Monday's system.C and mesh_refinement.C in the meantime... Yeah, I've triggered some subtle new bugs - apparantly of the difficult to debug "scribble on the wrong memory now and so crash the program much later" variety. I'd rather try to fix them today than revert the changes immediately. Anyone can avoid them by running: cvs update -D '5/17/05 05:00 CST' src/solvers/system.C src/mesh/mesh_refinement.C If anyone else wants to revert those files immediately instead, don't let my delusions of competence stop you. --- Roy |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-17 21:46:32
|
Yes, and MeshBase::prepare_for_use() repartitions the mesh. If you comment out the partition call then the existing code will probably work, albeit with a horrible load imbalance. -Ben -----Original Message----- From: Roy Stogner [mailto:roy...@ic...] Sent: Tuesday, May 17, 2005 4:25 PM To: KIRK, BENJAMIN (JSC-EG) (NASA) Cc: lib...@li... Subject: RE: [Libmesh-devel] Unstable code: incoming On Tue, 17 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > The problem is as follows: > > - You have some parallel, distributed data on an old mesh > - You refine the mesh, so you have new elements on your processor > - The old_dof_indices for the new elements on the processor are not > necessarily contained in the old vectors. That seems to be an obvious consequence of repartitioning the mesh, but just refining? Unless repartitioning is done, shouldn't new fine elements end up on the same processor as their parent? Their parent's DoFs are the only ones needed for projection. --- Roy |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-18 14:25:03
|
Is there an easy way to reproduce that bug from the examples? I'll see if I can help... -Ben -----Original Message----- From: lib...@li... [mailto:lib...@li...] On Behalf Of Roy Stogner Sent: Wednesday, May 18, 2005 9:00 AM To: lib...@li... Subject: RE: [Libmesh-devel] Unstable code: incoming On Tue, 17 May 2005, Roy Stogner wrote: > I think I can get everything working in the serial case tomorrow, but > anyone who doesn't want to futz with my bugs might want to revert to > Monday's system.C and mesh_refinement.C in the meantime... Yeah, I've triggered some subtle new bugs - apparantly of the difficult to debug "scribble on the wrong memory now and so crash the program much later" variety. I'd rather try to fix them today than revert the changes immediately. Anyone can avoid them by running: cvs update -D '5/17/05 05:00 CST' src/solvers/system.C src/mesh/mesh_refinement.C If anyone else wants to revert those files immediately instead, don't let my delusions of competence stop you. --- Roy ------------------------------------------------------- This SF.Net email is sponsored by Oracle Space Sweepstakes Want to be the first software developer in space? Enter now for the Oracle Space Sweepstakes! http://ads.osdn.com/?ad_id=7412&alloc_id=16344&op=click _______________________________________________ Libmesh-devel mailing list Lib...@li... https://lists.sourceforge.net/lists/listinfo/libmesh-devel |
From: Roy S. <roy...@ic...> - 2005-05-18 14:42:59
|
On Wed, 18 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > Is there an easy way to reproduce that bug from the examples? I'll > see if I can help... On my system, running ex10 in debug mode from the CVS head crashes after step 2. I don't know how reproduceable this is, though, because the crash location doesn't appear to be related to any recent changes - it appears that something elsewhere in the code has scribbled over a FE object's memory. If that's the case, then the crash might occur elsewhere or not at all depending on your linker... It also appears that commenting out one of the two Mesh::contract() calls (from either refine_and_coarsen_elements or System::reinit) prevents the bug on my system. That suggests that either something's wrong with Mesh::contract() (whose second call should be read only if there's been no mesh coarsening since the first call) or the bug is even earlier and has scribbled over the mesh's or geometric elements' memory too. --- Roy |
From: Roy S. <roy...@ic...> - 2005-05-18 14:52:17
|
On Wed, 18 May 2005, Roy Stogner wrote: > It also appears that commenting out one of the two Mesh::contract() > calls (from either refine_and_coarsen_elements or System::reinit) > prevents the bug on my system. Whoops - only commenting out both calls prevents it. There may be something wrong with the Elem::active() changes I checked in a while ago, in that case! Other than wasting memory and keeping around children for project_vector to play with, an uncontracted mesh is supposed to behave the same as a contracted one. --- Roy |
From: Roy S. <roy...@ic...> - 2005-05-18 16:15:16
|
On Wed, 18 May 2005, Roy Stogner wrote: > There may be something wrong with the Elem::active() changes I checked > in a while ago, in that case! Other than wasting memory and keeping > around children for project_vector to play with, an uncontracted mesh > is supposed to behave the same as a contracted one. Well, I've found one place where an uncontracted mesh fails: DoFMap::create_dof_constraints was happily calling FE::compute_constraints() on every element. The tests inside compute_constraints() skip over ancestor elements, but when acting on subactive elements could end up with disaster. With that fixed (already committed to CVS), adding a prepare_for_use() call to the end of Mesh::contract() appears to make the bug go away. I'm not going to commit that second change until I fully understand it, though - I want to make sure it's actually fixing a bug and not just hiding one again. --- Roy |
From: John P. <pet...@cf...> - 2005-05-18 16:19:56
|
Roy Stogner writes: > On Wed, 18 May 2005, Roy Stogner wrote: > > > There may be something wrong with the Elem::active() changes I checked > > in a while ago, in that case! Other than wasting memory and keeping > > around children for project_vector to play with, an uncontracted mesh > > is supposed to behave the same as a contracted one. > > Well, I've found one place where an uncontracted mesh fails: > DoFMap::create_dof_constraints was happily calling > FE::compute_constraints() on every element. The tests inside > compute_constraints() skip over ancestor elements, but when acting on > subactive elements could end up with disaster. > > With that fixed (already committed to CVS), adding a prepare_for_use() > call to the end of Mesh::contract() appears to make the bug go away. > I'm not going to commit that second change until I fully understand > it, though - I want to make sure it's actually fixing a bug and not > just hiding one again. It sounds like you know what's going on, but I have no idea what you are talking about half the time :) The beauty of open source? Just don't get hit by a bus any time soon!! -John |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-18 16:25:21
|
In that case I bet there is an issue with the neighbor information getting out of sync? An element can try to access its neighbor, which has previously been deleted in contract(). Of course, dereferencing the neighbor pointer is invalid. Just a thought. The call to prepare_for_use() will update the neighbor information and prevent that. BTW, WTF is an ancestor? :-) -Ben -----Original Message----- From: Roy Stogner [mailto:roy...@ic...] Sent: Wednesday, May 18, 2005 11:15 AM To: KIRK, BENJAMIN (JSC-EG) (NASA) Cc: lib...@li... Subject: RE: [Libmesh-devel] Unstable code: incoming On Wed, 18 May 2005, Roy Stogner wrote: > There may be something wrong with the Elem::active() changes I checked > in a while ago, in that case! Other than wasting memory and keeping > around children for project_vector to play with, an uncontracted mesh > is supposed to behave the same as a contracted one. Well, I've found one place where an uncontracted mesh fails: DoFMap::create_dof_constraints was happily calling FE::compute_constraints() on every element. The tests inside compute_constraints() skip over ancestor elements, but when acting on subactive elements could end up with disaster. With that fixed (already committed to CVS), adding a prepare_for_use() call to the end of Mesh::contract() appears to make the bug go away. I'm not going to commit that second change until I fully understand it, though - I want to make sure it's actually fixing a bug and not just hiding one again. --- Roy |
From: Roy S. <roy...@ic...> - 2005-05-18 17:00:55
|
On Wed, 18 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > BTW, WTF is an ancestor? :-) It's the reason John doesn't understand what I'm talking about: I'm not any smarter, I just make up my own vocabulary, like those pairs of twins that invent their own languages and then have to be painstakingly retaught to speak to normal people. The distinction of an element as "active/inactive" didn't make as much sense once I introduced "subactive" elements, elements which were technically contained in the mesh but which should be ignored by almost all code except project_vector. I called them "subactive" because they were children of active (or of other subactive) elements and thus "below" the active layer. So when I needed a word for "inactive but not subactive"; "superactive" had all the wrong connotations, and "ancestor" was the only thing I could think of that fit an existing ("parent"/"child") metaphor: an "ancestor" element is an element whose children are either active or other "ancestors". > In that case I bet there is an issue with the neighbor information getting > out of sync? An element can try to access its neighbor, which has > previously been deleted in contract(). Of course, dereferencing the > neighbor pointer is invalid. Just a thought. The call to prepare_for_use() > will update the neighbor information and prevent that. That turned out not to be the case. For backwards compatibility's sake, only subactive elements can have subactive neighbors. Active elements always have active neighbors, even if that neighbor is coarser and there's a subactive element that would be a closer fit; similarly ancestor elements still always have active or ancestor neighbors. The only place we could hit a deleted neighbor is from a subactive element, and AFAIK there's no code that looks for subactive neighbors except in find_neighbors itself. Or, to ramble less: even prepare_for_use() is overkill; for some reason just adding renumber_nodes_and_elements() to the end of contract() seems to fix this bug. --- Roy |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-18 17:12:31
|
That is with the fix you made in the DofMap for recomputing the constraints, right? Can you send me a patch for that particular file? I think the issue is probably that contract() does not set the elem->id(). renumber_nodes_and_elements() ensures that 0 <= elem->id() < mesh.n_elem() so you can do things like make an error_vector which is mesh.n_elem() long. I've changed contract() locally so that it properly sets the id, but I want to check it first. -Ben -----Original Message----- From: Roy Stogner [mailto:roy...@ic...] Sent: Wednesday, May 18, 2005 12:01 PM To: KIRK, BENJAMIN (JSC-EG) (NASA) Cc: lib...@li... Subject: RE: [Libmesh-devel] Unstable code: incoming On Wed, 18 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > BTW, WTF is an ancestor? :-) It's the reason John doesn't understand what I'm talking about: I'm not any smarter, I just make up my own vocabulary, like those pairs of twins that invent their own languages and then have to be painstakingly retaught to speak to normal people. The distinction of an element as "active/inactive" didn't make as much sense once I introduced "subactive" elements, elements which were technically contained in the mesh but which should be ignored by almost all code except project_vector. I called them "subactive" because they were children of active (or of other subactive) elements and thus "below" the active layer. So when I needed a word for "inactive but not subactive"; "superactive" had all the wrong connotations, and "ancestor" was the only thing I could think of that fit an existing ("parent"/"child") metaphor: an "ancestor" element is an element whose children are either active or other "ancestors". > In that case I bet there is an issue with the neighbor information > getting out of sync? An element can try to access its neighbor, which > has previously been deleted in contract(). Of course, dereferencing > the neighbor pointer is invalid. Just a thought. The call to > prepare_for_use() will update the neighbor information and prevent > that. That turned out not to be the case. For backwards compatibility's sake, only subactive elements can have subactive neighbors. Active elements always have active neighbors, even if that neighbor is coarser and there's a subactive element that would be a closer fit; similarly ancestor elements still always have active or ancestor neighbors. The only place we could hit a deleted neighbor is from a subactive element, and AFAIK there's no code that looks for subactive neighbors except in find_neighbors itself. Or, to ramble less: even prepare_for_use() is overkill; for some reason just adding renumber_nodes_and_elements() to the end of contract() seems to fix this bug. --- Roy |
From: Roy S. <roy...@ic...> - 2005-05-18 17:29:42
|
On Wed, 18 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > That is with the fix you made in the DofMap for recomputing the constraints, > right? Can you send me a patch for that particular file? I've already committed that to CVS; there was no reason to call compute_constraints on an inactive element even if doing so wasn't a bug. Update src/base/dof_constraints.C, or just change the iterator initializations to use mesh.active_* on lines 71 and 72. > I think the issue is probably that contract() does not set the elem->id(). > renumber_nodes_and_elements() ensures that > 0 <= elem->id() < mesh.n_elem() so you can do things like make an > error_vector which is mesh.n_elem() long. I've changed contract() locally > so that it properly sets the id, but I want to check it first. That sounds like it could be the problem - all the failing programs I've tested have been using such an error estimator. Let me know if it works! Oh, BTW, there's one more bug outstanding: currently my new assert(!this->eliminate_unrefined_patches()) in refine_and_coarsen_elements() is failing on timestep 27 of example 10; you may want to comment those asserts out. I'm still trying to figure out why that's happening; I thought my changes to eliminate_unrefined_patches() were sufficient to prevent unrefined patches from being created to begin with. This might be a result of non-updated neighbor pointers. BTW, why aren't we updating those on the fly? It seems like overkill to loop over every single mesh element in find_neighbors to fix just the pointers around a couple newly refined or coarsened cells. --- Roy |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-18 17:52:36
|
Check out out.gmv.020 from ex10... It seems that somehow the refinement flags are getting changed in a way they should not... The mesh starts as uniformly refined and should coarsen immediately everywhere except around the initial pulse. Now there are refined regions that persist, they seem to be coarsened one layer at a time on the perimeter. My guess is something is being overlooked in make_coarsening_compatible(). Updating the neighbors for all the cells is serious overkill. There is a TODO at the beginning of find_neighbors() that shouldn't be there, provided elem->nullify_neighbors() gets called in the right spot. Last time I tried that (a *long* time ago) I had some trouble, but I haven't messed with it recently. -Ben -----Original Message----- From: Roy Stogner [mailto:roy...@ic...] Sent: Wednesday, May 18, 2005 12:30 PM To: KIRK, BENJAMIN (JSC-EG) (NASA) Cc: lib...@li... Subject: RE: [Libmesh-devel] Unstable code: incoming On Wed, 18 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > That is with the fix you made in the DofMap for recomputing the > constraints, right? Can you send me a patch for that particular file? I've already committed that to CVS; there was no reason to call compute_constraints on an inactive element even if doing so wasn't a bug. Update src/base/dof_constraints.C, or just change the iterator initializations to use mesh.active_* on lines 71 and 72. > I think the issue is probably that contract() does not set the > elem->id(). > renumber_nodes_and_elements() ensures that > 0 <= elem->id() < mesh.n_elem() so you can do things like make an > error_vector which is mesh.n_elem() long. I've changed contract() locally > so that it properly sets the id, but I want to check it first. That sounds like it could be the problem - all the failing programs I've tested have been using such an error estimator. Let me know if it works! Oh, BTW, there's one more bug outstanding: currently my new assert(!this->eliminate_unrefined_patches()) in refine_and_coarsen_elements() is failing on timestep 27 of example 10; you may want to comment those asserts out. I'm still trying to figure out why that's happening; I thought my changes to eliminate_unrefined_patches() were sufficient to prevent unrefined patches from being created to begin with. This might be a result of non-updated neighbor pointers. BTW, why aren't we updating those on the fly? It seems like overkill to loop over every single mesh element in find_neighbors to fix just the pointers around a couple newly refined or coarsened cells. --- Roy |
From: Roy S. <roy...@ic...> - 2005-05-18 18:31:34
|
On Wed, 18 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > Check out out.gmv.020 from ex10... It seems that somehow the refinement > flags are getting changed in a way they should not... The mesh starts as > uniformly refined and should coarsen immediately everywhere except around > the initial pulse. Now there are refined regions that persist, they seem to > be coarsened one layer at a time on the perimeter. My guess is something is > being overlooked in make_coarsening_compatible(). Something was overlooked in my new eliminate_unrefined_patches(). I've committed the fix to CVS, and it appears to work. Note that even if eliminate_unrefined_patches() is working now, there's one behavior change: instead of just eliminating unrefined patches, now AMR will refuse to create them in the first place. I seem to recall some "flickering" in one of John's experiments, of a cell that was coarsened every odd timestep by the error indicator and refined every even timestep by eliminate_unrefined_patches(); hopefully the change will prevent situations like that. > Updating the neighbors for all the cells is serious overkill. There > is a TODO at the beginning of find_neighbors() that shouldn't be > there, provided elem->nullify_neighbors() gets called in the right > spot. I don't think it's where elem->nullify_neighbors() gets called, but rather what it does. If it set neighbor pointers to the correct values rather than to NULL, that would handle coarsening at least. > Last time I tried that (a *long* time ago) I had some > trouble, but I haven't messed with it recently. With the luck I'm having this week, I'm not going to mess with a troublesome change just for efficiency's sake, not when the inefficient code is still O(n). In fact, I think I may add an additional find_neighbors() call (wrapped in #ifdef DEBUG) to the middle of refine_and_coarsen_elements() so redundant calls of eliminate_unrefined_patches() will continue to be no-ops and I can leave in those assert()s. Fixing nullify_neighbors to set pointers for a newly coarsened cell itself should be easy, but making sure every element's neighbor's descendants all get set correctly sounds ugly, and writing another function to handle neighbors of newly refined cells correctly would be even uglier. --- Roy |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-19 15:05:13
|
Most everything I've tested so far is good, but there is still a problem in ex10. It now dies at step 43 with Solving time step 43, time=1.1000... Refining the mesh... ex10: src/mesh/mesh_refinement.C:198: bool MeshRefinement::refine_and_coarsen_elements(bool): Assertion `!this->eliminate_unrefined_patches()' failed. Thoughts? -----Original Message----- From: Roy Stogner [mailto:roy...@ic...] Sent: Wednesday, May 18, 2005 1:31 PM To: KIRK, BENJAMIN (JSC-EG) (NASA) Cc: lib...@li... Subject: RE: [Libmesh-devel] Unstable code: incoming On Wed, 18 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > Check out out.gmv.020 from ex10... It seems that somehow the > refinement flags are getting changed in a way they should not... The > mesh starts as uniformly refined and should coarsen immediately > everywhere except around the initial pulse. Now there are refined > regions that persist, they seem to be coarsened one layer at a time on > the perimeter. My guess is something is being overlooked in > make_coarsening_compatible(). Something was overlooked in my new eliminate_unrefined_patches(). I've committed the fix to CVS, and it appears to work. Note that even if eliminate_unrefined_patches() is working now, there's one behavior change: instead of just eliminating unrefined patches, now AMR will refuse to create them in the first place. I seem to recall some "flickering" in one of John's experiments, of a cell that was coarsened every odd timestep by the error indicator and refined every even timestep by eliminate_unrefined_patches(); hopefully the change will prevent situations like that. > Updating the neighbors for all the cells is serious overkill. There > is a TODO at the beginning of find_neighbors() that shouldn't be > there, provided elem->nullify_neighbors() gets called in the right > spot. I don't think it's where elem->nullify_neighbors() gets called, but rather what it does. If it set neighbor pointers to the correct values rather than to NULL, that would handle coarsening at least. > Last time I tried that (a *long* time ago) I had some trouble, but I > haven't messed with it recently. With the luck I'm having this week, I'm not going to mess with a troublesome change just for efficiency's sake, not when the inefficient code is still O(n). In fact, I think I may add an additional find_neighbors() call (wrapped in #ifdef DEBUG) to the middle of refine_and_coarsen_elements() so redundant calls of eliminate_unrefined_patches() will continue to be no-ops and I can leave in those assert()s. Fixing nullify_neighbors to set pointers for a newly coarsened cell itself should be easy, but making sure every element's neighbor's descendants all get set correctly sounds ugly, and writing another function to handle neighbors of newly refined cells correctly would be even uglier. --- Roy |
From: Roy S. <roy...@ic...> - 2005-05-19 19:43:00
|
On Thu, 19 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > Most everything I've tested so far is good, but there is still a problem in > ex10. It now dies at step 43 with > > Solving time step 43, time=1.1000... > Refining the mesh... > ex10: src/mesh/mesh_refinement.C:198: bool > MeshRefinement::refine_and_coarsen_elements(bool): Assertion > `!this->eliminate_unrefined_patches()' failed. > > Thoughts? Hopefully the patch I committed this morning fixes that. I decided to remove the extra eliminate_unrefined_patches() assertion rather than insert the extra find_neighbors() necessary to ensure it would work. --- Roy |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-19 18:14:22
|
Multiple systems do not seem to work with the current coarsening/refinement scheme. No single system should be responsible for coarsening/refining the mesh, rather the EquationSystems should handle that (or the user directly?) and then perform restriction/projection on a system-by-system basis. -Ben |
From: Roy S. <roy...@ic...> - 2005-05-19 20:01:41
|
On Thu, 19 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > Multiple systems do not seem to work with the current coarsening/refinement > scheme. Wow. I didn't think about that... The shame goes on... > No single system should be responsible for coarsening/refining the > mesh, rather the EquationSystems should handle that (or the user directly?) > and then perform restriction/projection on a system-by-system basis. I'd rather not make the user handle the coarsening and refining directly, except that the current functions which do so should be taken into account for backwards compatibility. The less boilerplate users have to write, the better, and the easier it will be to make underlying but backwards compatibile changes in the future. I'm not certain what the correct implementation is. We can't just have EquationSystems::reinit call System::reinit for every system. The current System::reinit behavior is three-pass: one attempted restrict-and-prolong projection in case the user has already done the coarsening/refining, then a coarsen+restrict, then a refine+project. That three pass behavior needs to be done one pass at a time, all systems in each pass. It ought to be done in EquationSystems::reinit but what do we call the functions that get called? Perhaps reinit() for the first, which will have to be called for every EquationSystems::reinit - even if we end backwards compatibility and assume the user has only flagged rather than altered the mesh, the re_update() (and possibly other functions in derived classes) will have to be called every time. How about "coarsen_and_restrict" for the second (which will only be called if the EquationSystems::reinit finds cells to coarsen) and "refine_and_project" for the third (which will only be called if the EquationSystems::reinit finds cells to refine)? --- Roy |
From: Roy S. <roy...@ic...> - 2005-05-20 05:39:32
|
On Thu, 19 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > Multiple systems do not seem to work with the current coarsening/refinement > scheme. No single system should be responsible for coarsening/refining the > mesh, rather the EquationSystems should handle that (or the user directly?) > and then perform restriction/projection on a system-by-system basis. Making the user handle it directly would break existing user level code. I'm still trying to keep backwards compatibility for the moment; a big API change (the subclassable System) every now and again may be necessary, but small code-breaking PETSc-esque changes ought to be avoided when possible. I'd like to put out something like an "0.5.0" with a backwards compatible (if occasionally deprecated) API, so there'll be a fixed release with C1 elements, non-Lagrange AMR, OS X support, triangle/tetmesh, etc. Ironically, one reason I wanted such a thing soon was to be able to point people to it if the CVS API was ever unstable while we were designing and implementing the System changes. After all, if you give CVS write access to jerks like me, who knows what's liable to be broken... :-( Speaking of which, I've commited System and EquationSystems changes that I hope fix the multiple systems problem. Let me know if it's insufficient (or if there are any more regressions...) --- Roy |
From: KIRK, B. (JSC-E. (NASA) <ben...@na...> - 2005-05-19 20:20:26
|
See if this makes any more sense given the current situation... http://sourceforge.net/mailarchive/forum.php?thread_id=6762215&forum_id=2462 8 I'm not at all opposed to changing the user interface slightly from what it has been to date. It seems we will need to add some functions now and change the user interface without much concern about backwards compatibility, and then in the future the details can be hidden within a smart solver class. Given all the steps that have to performed, I still think the best interface might be something like mesh_refinement.flag_elements_by_error_fraction(...); equation_systems.restrict(); mesh_refinement.refine_and_coarsen_elements(); equation_systems.project(); or mesh_refinement.flag_elements_by_error_fraction(...); mesh_refinement.coarsen_elements(); equation_systems.restrict(); mesh_refinement.refine_elements(); equation_systems.project(); When the user knows they are working only with Lagrange elements, for example, they can omit the coarsen/restrict calls in the second option. -Ben -----Original Message----- From: Roy Stogner [mailto:roy...@ic...] Sent: Thursday, May 19, 2005 3:02 PM To: KIRK, BENJAMIN (JSC-EG) (NASA) Cc: lib...@li... Subject: RE: [Libmesh-devel] Unstable code: incoming On Thu, 19 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > Multiple systems do not seem to work with the current > coarsening/refinement scheme. Wow. I didn't think about that... The shame goes on... > No single system should be responsible for coarsening/refining the > mesh, rather the EquationSystems should handle that (or the user > directly?) and then perform restriction/projection on a > system-by-system basis. I'd rather not make the user handle the coarsening and refining directly, except that the current functions which do so should be taken into account for backwards compatibility. The less boilerplate users have to write, the better, and the easier it will be to make underlying but backwards compatibile changes in the future. I'm not certain what the correct implementation is. We can't just have EquationSystems::reinit call System::reinit for every system. The current System::reinit behavior is three-pass: one attempted restrict-and-prolong projection in case the user has already done the coarsening/refining, then a coarsen+restrict, then a refine+project. That three pass behavior needs to be done one pass at a time, all systems in each pass. It ought to be done in EquationSystems::reinit but what do we call the functions that get called? Perhaps reinit() for the first, which will have to be called for every EquationSystems::reinit - even if we end backwards compatibility and assume the user has only flagged rather than altered the mesh, the re_update() (and possibly other functions in derived classes) will have to be called every time. How about "coarsen_and_restrict" for the second (which will only be called if the EquationSystems::reinit finds cells to coarsen) and "refine_and_project" for the third (which will only be called if the EquationSystems::reinit finds cells to refine)? --- Roy |
From: John P. <pet...@cf...> - 2005-05-19 20:26:16
|
KIRK, BENJAMIN (JSC-EG) (NASA) writes: > See if this makes any more sense given the current situation... > http://sourceforge.net/mailarchive/forum.php?thread_id=6762215&forum_id=2462 > 8 > > I'm not at all opposed to changing the user interface slightly from what it > has been to date. Yeah, no worries... we are still at version 0.x :) OK, sorry just had to inject my 2$CURRENCY_UNIT. -John |
From: John P. <pet...@cf...> - 2005-05-17 21:04:57
|
KIRK, BENJAMIN (JSC-EG) (NASA) writes: > The problem is as follows: > > - You have some parallel, distributed data on an old mesh > - You refine the mesh, so you have new elements on your processor > - The old_dof_indices for the new elements on the processor are > not necessarily contained in the old vectors. > > The ugly hack is a multi-stage approach as follows: > > localize() the old vectors on each processor to a serial vector > containing *all* the vector values > > project() these serial vectors > > copy only the local components into the locations of the new parallel > vector. > > This is exactly the problem I mentioned in the cryptic midnight email I sent > out the other day. John, did you test the old implementation in parallel > with debugging? I sure didn't. But I'd be happy to revert to whatever date and give it a go. Actually you can reproduce this behavior in example 10 (for sure) and possibly others if you want to see for yourself... > Anyway, I'll do a cvs update and see how I can help. -John |
From: Roy S. <roy...@ic...> - 2005-05-17 21:25:07
|
On Tue, 17 May 2005, KIRK, BENJAMIN (JSC-EG) (NASA) wrote: > The problem is as follows: > > - You have some parallel, distributed data on an old mesh > - You refine the mesh, so you have new elements on your processor > - The old_dof_indices for the new elements on the processor are > not necessarily contained in the old vectors. That seems to be an obvious consequence of repartitioning the mesh, but just refining? Unless repartitioning is done, shouldn't new fine elements end up on the same processor as their parent? Their parent's DoFs are the only ones needed for projection. --- Roy |