From: <ti...@ce...> - 2006-11-03 15:50:25
Attachments:
test.cpp
|
Dear all, The attached example code does the following: * Initialize a system with one component. * Set the initial value to the constant value of 1. * Compute the min and max value: Both are 1 as expected. * Refine the grid uniformly. * Compute the min and max values again. In the last item, the max value is always 1, but the min value is 0 on more than one processor. What am I doing wrong? Best Regards, Tim |
From: John P. <pet...@cf...> - 2006-11-03 16:13:57
|
Tim Kr=F6ger writes: > Dear all, >=20 > The attached example code does the following: > * Initialize a system with one component. > * Set the initial value to the constant value of 1. > * Compute the min and max value: Both are 1 as expected. > * Refine the grid uniformly. > * Compute the min and max values again. >=20 > In the last item, the max value is always 1, but the min value is 0 = on=20 > more than one processor. What am I doing wrong=3F You are looping over non-active elements in computeMinMax. Can you try= it looping over active elements only=3F -J > Best Regards, >=20 > Tim#include <iostream> > #include <fstream> > #include <math.h> >=20 > #include "libmesh.h" > #include "mesh.h" > #include "mesh=5Fgeneration.h" > #include "equation=5Fsystems.h" > #include "fe.h" > #include "quadrature=5Fgauss.h" > #include "dof=5Fmap.h" > #include "sparse=5Fmatrix.h" > #include "numeric=5Fvector.h" > #include "dense=5Fmatrix.h" > #include "dense=5Fvector.h" > #include "linear=5Fimplicit=5Fsystem.h" > #include "transient=5Fsystem.h" > #include "linear=5Fsolver.h" > #include "mesh=5Frefinement.h" > #include "elem.h" >=20 > void assemble (EquationSystems& es, > =09=09const std::string& system=5Fname) > { > } >=20 > void init (EquationSystems& es, > =09 const std::string& system=5Fname) > { > const Mesh& mesh =3D es.get=5Fmesh(); > const unsigned int dim =3D mesh.mesh=5Fdimension(); > TransientLinearImplicitSystem & system =3D es.get=5Fsystem<Transie= ntLinearImplicitSystem> ("a"); > const unsigned int var =3D system.variable=5Fnumber ("u"); > FEType feType =3D system.variable=5Ftype(var); > AutoPtr<FEBase> fe (FEBase::build(dim, feType)); > const DofMap & dof=5Fmap =3D system.get=5Fdof=5Fmap(); > std::vector<unsigned int> dof=5Findices; > MeshBase::const=5Felement=5Fiterator el =3D mesh.active=5F= local=5Felements=5Fbegin(); > const MeshBase::const=5Felement=5Fiterator end=5Fel =3D mesh.activ= e=5Flocal=5Felements=5Fend();=20 > for ( ; el !=3D end=5Fel; ++el) > { =20 > const Elem* elem =3D *el; > dof=5Fmap.dof=5Findices (elem, dof=5Findices, var); > const unsigned int n=5Fdofs =3D dof=5Findices.size(); > for (unsigned int l=3D0; l<n=5Fdofs; l++) > { > if(dof=5Findices[l]>=3Dsystem.solution->first=5Flocal=5Fin= dex() && > dof=5Findices[l]<system.solution->last=5Flocal=5Findex(= )) > { > system.solution->set(dof=5Findices[l], 1.0); > } > } > } > system.update(); > } >=20 > void computeMinMax(EquationSystems& es) > { > const Mesh& mesh =3D es.get=5Fmesh(); > TransientLinearImplicitSystem & system =3D es.get=5Fsystem<Transie= ntLinearImplicitSystem> ("a"); > const DofMap& dof=5Fmap =3D system.get=5Fdof=5Fmap(); > std::vector<unsigned int> dof=5Findices; >=20 > Real minVal =3D 100000.0; > Real maxVal =3D -100000.0; >=20 > MeshBase::const=5Felement=5Fiterator el =3D mesh.local=5F= elements=5Fbegin(); > const MeshBase::const=5Felement=5Fiterator end=5Fel =3D mesh.local= =5Felements=5Fend();=20 > for ( ; el !=3D end=5Fel; ++el) > { =20 > Elem* elem =3D *el; > dof=5Fmap.dof=5Findices (elem, dof=5Findices); > const unsigned int n=5Fdofs =3D dof=5Findices.size();=20 > for (unsigned int l=3D0; l<n=5Fdofs; l++) > { > =09 const Real value =3D (*system.current=5Flocal=5Fsolution)(dof=5F= indices[l]); > =09 if(maxVal<value) > =09 { > =09 maxVal =3D value; > =09 } > =09 if(minVal>value) > =09 { > =09 minVal =3D value; > =09 } > =09} > } >=20 > { > std::vector<double> buf(2,0.0); > buf[0] =3D -minVal; > buf[1] =3D maxVal; > std::vector<double> buf2(buf); > MPI=5FAllreduce(&buf2[0], > =09=09 &buf[0], > =09=09 buf.size(), > =09=09 MPI=5FDOUBLE, > =09=09 MPI=5FMAX, > =09=09 libMesh::COMM=5FWORLD); > maxVal =3D buf[1]; > minVal =3D -buf[0]; > } >=20 > std::cout << " Min/Max: " << minVal << ", " << maxVal << std::end= l; > } >=20 > int main(int argc, char** argv) > { > libMesh::init (argc, argv); > { =20 > Mesh mesh (3); > MeshTools::Generation::build=5Fcube(mesh,2,2,2,0.0,1.0,0.0,1.0,0= .0,1.0,HEX8); > EquationSystems equation=5Fsystems (mesh); >=20 > TransientLinearImplicitSystem & system =3D=20 > equation=5Fsystems.add=5Fsystem<TransientLinearImplicitSystem>= ("a"); > system.add=5Fvariable ("u", FIRST, LAGRANGE); > system.attach=5Fassemble=5Ffunction (assemble); > system.attach=5Finit=5Ffunction (init); > equation=5Fsystems.init (); > =20 > MeshRefinement meshRefinement(mesh); >=20 > computeMinMax(equation=5Fsystems); >=20 > meshRefinement.uniformly=5Frefine(); > equation=5Fsystems.reinit(); >=20 > computeMinMax(equation=5Fsystems); > } > return libMesh::close (); > } |
From: Roy S. <roy...@ic...> - 2006-11-03 16:27:54
|
On Fri, 3 Nov 2006, John Peterson wrote: > Tim Kr=F6ger writes: > > Dear all, > > > > The attached example code does the following: > > * Initialize a system with one component. > > * Set the initial value to the constant value of 1. > > * Compute the min and max value: Both are 1 as expected. > > * Refine the grid uniformly. > > * Compute the min and max values again. > > > > In the last item, the max value is always 1, but the min value is 0 o= n > > more than one processor. What am I doing wrong? > > You are looping over non-active elements in computeMinMax. Can you try > it looping over active elements only? Looping between active_local_elements_* fixes it for me. Let me guess: the partition between processors becomes uneven, so that parent elements' children are split between processors - and when a parent element is asked for the solution value at a node that lies solely on another processor, it returns zero? That's the only scenario I can think of, but it seems like it shouldn't happen on an 8 hex mesh; it should be practically impossible for the partitioning to happen in such a way that the node being asked for isn't at least a ghost node. Maybe this is technically a bug in user code, but it still looks like something that should be fixed in the library; it would be nice if at least there was an assert failing somewhere. --- Roy |
From: <ti...@ce...> - 2006-11-06 15:25:01
Attachments:
test.cpp
|
Dear Roy and John, On Fri, 3 Nov 2006, Roy Stogner wrote: > Looping between active_local_elements_* fixes it for me. For me also. However, only in that special situation. Please find attached a new version of the test code. The changes are as follows: * Loop over active cells only, as suggested; * Use HEX27 cells instead of HEX8; * Refine only one cell instead of all; * Repeat refinement 10 times. The error again occurs, namely after the fourth refinement step (if run on two processors). What did I do wrong this time? Best Regards, Tim |
From: Roy S. <roy...@ic...> - 2006-11-06 15:50:23
|
On Mon, 6 Nov 2006, Tim Kr=F6ger wrote: > On Fri, 3 Nov 2006, Roy Stogner wrote: > >> Looping between active_local_elements_* fixes it for me. > > For me also. However, only in that special situation. Please find att= ached=20 > a new version of the test code. The changes are as follows: > > * Loop over active cells only, as suggested; > * Use HEX27 cells instead of HEX8; > * Refine only one cell instead of all; > * Repeat refinement 10 times. > > The error again occurs, namely after the fourth refinement step (if run= on=20 > two processors). > > What did I do wrong this time? This time you didn't do anything wrong that I can see at first glance - and it's not that complicated a program. I'm playing around with your test case, but for debugging in parallel I'm only familiar with stone-age methods; if anyone else is handy with more than std::cerr, now's a good time to speak up. This kind of bug would wreak havoc on a time-dependent parallel adaptive code. ... which actually makes me wonder why nobody's seen it before. My own time-dependent parallel adaptive code mostly uses Hermite and Clough elements, so perhaps there's a problem with the Lagrange code path in System::project_vector... but John and others do this sort of thing with the Lagrange elements all the time. You'd think someone else would have run into this bug already - it even affects the 0.5.0 release, not just CVS! I'm going to forward this to Ben in case he's not following the lists right now; he's busier than hell these days but I think he's the one who wrote most of the parallel code originally, and he uses it enough that he definitely won't want to miss this. --- Roy |
From: Kirk, B. \(JSC-EG\) <ben...@na...> - 2006-11-06 17:20:12
|
I've been loosely following this, but you now have my full attention... =20 -----Original Message----- From: Roy Stogner [mailto:roy...@ic...]=20 Sent: Monday, November 06, 2006 9:50 AM To: Tim Kr=F6ger Cc: John Peterson; lib...@li...; Kirk, Benjamin = (JSC-EG) Subject: Re: [Libmesh-users] Solution projection On Mon, 6 Nov 2006, Tim Kr=F6ger wrote: > On Fri, 3 Nov 2006, Roy Stogner wrote: > >> Looping between active_local_elements_* fixes it for me. > > For me also. However, only in that special situation. Please find=20 > attached a new version of the test code. The changes are as follows: > > * Loop over active cells only, as suggested; > * Use HEX27 cells instead of HEX8; > * Refine only one cell instead of all; > * Repeat refinement 10 times. > > The error again occurs, namely after the fourth refinement step (if=20 > run on two processors). > > What did I do wrong this time? This time you didn't do anything wrong that I can see at first glance - and it's not that complicated a program. I'm playing around with your test case, but for debugging in parallel = I'm only familiar with stone-age methods; if anyone else is handy with = more than std::cerr, now's a good time to speak up. This kind of bug = would wreak havoc on a time-dependent parallel adaptive code. ... which actually makes me wonder why nobody's seen it before. My own = time-dependent parallel adaptive code mostly uses Hermite and Clough = elements, so perhaps there's a problem with the Lagrange code path in = System::project_vector... but John and others do this sort of thing with = the Lagrange elements all the time. You'd think someone else would have = run into this bug already - it even affects the 0.5.0 release, not just = CVS! I'm going to forward this to Ben in case he's not following the lists = right now; he's busier than hell these days but I think he's the one who = wrote most of the parallel code originally, and he uses it enough that = he definitely won't want to miss this. --- Roy |
From: Kirk, B. \(JSC-EG\) <ben...@na...> - 2006-11-06 19:15:40
Attachments:
patch
|
OK, so this was a tricky one.... Check out the attached patch. Here is a summary of what I think is = going on: When solution vectors are projected to a new mesh the old solution = vector is used. It must therefore contain valid values for all local = degrees of freedom. For an implicit system these DOF constraints are = usually built-in to the linear problem so that the solution vector = contains the proper values for the constrained nodes. In your example, there is no linear solve going on, hence there is no = mechanism for ensuring the proper constrained DOF values to make it back = in to the system.solution vector. So, after a few steps, the = system.solution vector becomes tainted with 0s. =20 I'm still trying to track down why this only appears in parallel, but = that is proving a little difficult... =20 See if this fix works for you. =20 -Ben -----Original Message----- From: Roy Stogner [mailto:roy...@ic...]=20 Sent: Monday, November 06, 2006 9:50 AM To: Tim Kr=F6ger Cc: John Peterson; lib...@li...; Kirk, Benjamin = (JSC-EG) Subject: Re: [Libmesh-users] Solution projection On Mon, 6 Nov 2006, Tim Kr=F6ger wrote: > On Fri, 3 Nov 2006, Roy Stogner wrote: > >> Looping between active_local_elements_* fixes it for me. > > For me also. However, only in that special situation. Please find=20 > attached a new version of the test code. The changes are as follows: > > * Loop over active cells only, as suggested; > * Use HEX27 cells instead of HEX8; > * Refine only one cell instead of all; > * Repeat refinement 10 times. > > The error again occurs, namely after the fourth refinement step (if=20 > run on two processors). > > What did I do wrong this time? This time you didn't do anything wrong that I can see at first glance - and it's not that complicated a program. I'm playing around with your test case, but for debugging in parallel = I'm only familiar with stone-age methods; if anyone else is handy with = more than std::cerr, now's a good time to speak up. This kind of bug = would wreak havoc on a time-dependent parallel adaptive code. ... which actually makes me wonder why nobody's seen it before. My own = time-dependent parallel adaptive code mostly uses Hermite and Clough = elements, so perhaps there's a problem with the Lagrange code path in = System::project_vector... but John and others do this sort of thing with = the Lagrange elements all the time. You'd think someone else would have = run into this bug already - it even affects the 0.5.0 release, not just = CVS! I'm going to forward this to Ben in case he's not following the lists = right now; he's busier than hell these days but I think he's the one who = wrote most of the parallel code originally, and he uses it enough that = he definitely won't want to miss this. --- Roy |
From: Roy S. <roy...@ic...> - 2006-11-06 21:10:14
|
On Mon, 6 Nov 2006, Kirk, Benjamin (JSC-EG) wrote: > Check out the attached patch. Here is a summary of what I think is > going on: > > When solution vectors are projected to a new mesh the old solution > vector is used. It must therefore contain valid values for all > local degrees of freedom. For an implicit system these DOF > constraints are usually built-in to the linear problem so that the > solution vector contains the proper values for the constrained > nodes. > > In your example, there is no linear solve going on, hence there is > no mechanism for ensuring the proper constrained DOF values to make > it back in to the system.solution vector. So, after a few steps, > the system.solution vector becomes tainted with 0s. > > I'm still trying to track down why this only appears in parallel, > but that is proving a little difficult... You're debugging with std::cerr too, aren't you? ;-) There's a related bug (well, what I would consider a bug, anyway) with adaptive coarsening: when coarsening an element causes nodes of it's more refined neighbors to become hanging nodes, there's no mechanism to immediately adjust those hanging node DoF values to be properly constrained. For most libMesh codes, that's not a big problem; the next linear solve will re-constrain everything nicely. I've had to add a call to enforce_constraints_exactly() to the start of my NewtonSolver class, though. It was only enforcing constraints on the deltas between iterates and not on the Newton iterates themselves, and having perpetually misconstrained hanging nodes was causing me all sorts of problems. What I never figured out, though, was why the misconstrained dofs seemed to only be giving me trouble in parallel... Anyway, depending on an application's formulation, it's possible for systems that aren't using my framework to have trouble with misconstrained inputs too. Should we just throw a call to enforce_constraints_exactly() at the end of System::project_vector() and kill two or three birds with one stone? --- Roy |
From: <ti...@ce...> - 2006-11-07 09:58:37
|
Dear Roy and Ben and all, Thank you very much for your help. In fact, Ben's patch helps. However, since in my production code, I don't compute the minimum and maximum (at least not so often), I have written the following function, which I call after each refinement operation: void postpareRefinement(EquationSystems& es) { const Mesh& mesh = es.get_mesh(); std::vector<unsigned int> dof_indices; for(unsigned int i=0; i<es.n_systems(); i++) { System& system = es.get_system(i); const DofMap& dof_map = system.get_dof_map(); MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); for ( ; el != end_el; ++el) { Elem* elem = *el; dof_map.dof_indices (elem, dof_indices); const unsigned int n_dofs = dof_indices.size(); for (unsigned int l=0; l<n_dofs; l++) { system.solution->set(dof_indices[l], (*system.current_local_solution)(dof_indices[l])); } } system.solution->close(); } } This works, but I wonder whether there might be something simpler. Note that in my production code, I have several systems in one EquationSystem object, and I typically solve one of them, then refine, solve again etc. I won't touch the other systems before the first system is solved satisfactorily. Best Regards, Tim |
From: <ti...@ce...> - 2006-11-07 15:33:28
Attachments:
test.cpp
|
Dear Roy and Ben and all, Another problem is that I don't know how to generalize Ben's patch to the situation that I have additional vectors in my systems. This is demonstrated with the attached test code; it's still the same (with Ben's patch integrated), but this time I also compute the minimum and maximum of old_local_solution. This is still filled with some zeros (after the fourth refinement step on two processors). This is not surprising, since I don't know how to adopt Ben's patch to the old_local_solution. In my real application, I do the following: 1.) Set *old_local_solution = *current_local_solution; 2.) Solve system. 3.) Refine grid. 4.) Solve system again. Step 4 needs to know about the old_local_solution vector (since it's a transient system) and gets wrong values. Any ideas? Best Regards, Tim |
From: Roy S. <roy...@ic...> - 2006-11-07 19:43:13
|
On Mon, 6 Nov 2006, Kirk, Benjamin (JSC-EG) wrote: > Check out the attached patch. Here is a summary of what I think is > going on: > > When solution vectors are projected to a new mesh the old solution > vector is used. It must therefore contain valid values for all > local degrees of freedom. For an implicit system these DOF > constraints are usually built-in to the linear problem so that the > solution vector contains the proper values for the constrained > nodes. The more I think about it, the more I think this can't be the problem. System::project_vector, if it's working as designed, should never produce a misconstrained degree of freedom when refining a finite element based on nested function spaces. The LAGRANGE codepath uses interpolation at each node; the generic codepath uses L2 projection, and either way the solution on hanging nodes should be exact (except for floating point error) unless there's simultaneously coarsening being done on a neighboring element. Anyway, I'm committing a DofMap::max_constraint_error() function to CVS in case anyone else wants to play around with it to test these theories. I haven't used it on that test case program yet, but now that I know what I'm looking for I think I've tracked down a bug in one of my application codes that may have the same root cause... --- Roy |
From: Roy S. <roy...@ic...> - 2006-11-06 15:54:38
|
On Mon, 6 Nov 2006, Tim Kr=F6ger wrote: > * Loop over active cells only, as suggested; > * Use HEX27 cells instead of HEX8; > * Refine only one cell instead of all; > * Repeat refinement 10 times. Here's an oddity: I can't seem to find a 2D test case that exhibits this bug, and in 3D the HEX20 seems immune too. Is it possible that there's something wrong with the HEX27 element specifically? --- Roy |
From: Roy S. <roy...@ic...> - 2006-11-06 16:06:20
|
On Mon, 6 Nov 2006, Roy Stogner wrote: > On Mon, 6 Nov 2006, Tim Kr=F6ger wrote: > >> * Loop over active cells only, as suggested; >> * Use HEX27 cells instead of HEX8; >> * Refine only one cell instead of all; >> * Repeat refinement 10 times. > > Here's an oddity: I can't seem to find a 2D test case that exhibits > this bug, and in 3D the HEX20 seems immune too. Is it possible that > there's something wrong with the HEX27 element specifically? And oddity 2: I don't have any problems with second-order Lagrange FE on HEX27. --- Roy |
From: Roy S. <roy...@ic...> - 2006-11-09 08:32:43
|
On Mon, 6 Nov 2006, Tim Kr=F6ger wrote: > * Loop over active cells only, as suggested; > * Use HEX27 cells instead of HEX8; > * Refine only one cell instead of all; > * Repeat refinement 10 times. > > The error again occurs, namely after the fourth refinement step (if run= on=20 > two processors). Okay, it looks like there's at least two bugs at work here. For the first bug, I've committed a patch to CVS - basically what seems to be happening is that we would get hanging node degrees of freedom that were owned by processors which didn't express those degrees of freedom on any elements. So when System::project_vector() time came around, the hanging degrees of freedom would be correctly calculated only on processors which would never contribute those calculations to the new global vector. For the second bug... I haven't found it yet. Filling the _send_list (e.g. in DofMap::sort_send_list) seems to be a successful workaround, though (in combination with a fix for the first bug), so the problem has to be missing _send_list entries somewhere. I've never worked on that section of the code, though, so I'm going to cross my fingers and hope Ben fixes it. ;-) Everyone using CVS should keep an eye out to make sure your apps don't break after the next update. I've dome some testing on this patch before committing it, of course, but it's a tricky change, and if there's one thing this mess has taught me it's that it's easy for uncommon cases to get missed by test code... --- Roy |
From: <ti...@ce...> - 2006-11-09 09:40:35
|
Dear Roy, On Thu, 9 Nov 2006, Roy Stogner wrote: > Okay, it looks like there's at least two bugs at work here. > > For the first bug, I've committed a patch to CVS - basically what > seems to be happening is that we would get hanging node degrees of > freedom that were owned by processors which didn't express those > degrees of freedom on any elements. So when System::project_vector() > time came around, the hanging degrees of freedom would be correctly > calculated only on processors which would never contribute those > calculations to the new global vector. > > For the second bug... I haven't found it yet. Filling the _send_list > (e.g. in DofMap::sort_send_list) seems to be a successful workaround, > though (in combination with a fix for the first bug), so the problem > has to be missing _send_list entries somewhere. I've never worked on > that section of the code, though, so I'm going to cross my fingers and > hope Ben fixes it. ;-) What do I have to fill the _send_list with? Could you give an example for that? Since your fix of the first bug alone does not seem to cure the problem (in particular, the same test program as before does even stranger things than before), it would be nice to have at least a workaround so that I can go on with my work. I'm attaching the test program again and the output it currently produces on two processors -- just to make sure that we're talking about the same thing. Note that Ben's change is already integrated. Best Regards, Tim |
From: Roy S. <roy...@ic...> - 2006-11-11 02:29:22
|
On Thu, 9 Nov 2006, Tim Kr=F6ger wrote: > What do I have to fill the _send_list with? Could you give an example = for=20 > that? Since your fix of the first bug alone does not seem to cure the=20 > problem (in particular, the same test program as before does even stran= ger=20 > things than before), it would be nice to have at least a workaround so = that I=20 > can go on with my work. Okay; I've rewritten the _send_list construction in CVS to be a little more liberal about the DoFs it adds; you should be able to run your code without any workarounds now. Let me know if you run into any more problems. While playing around with parallel adaptive corner cases, I found another nasty bug - the caching trick I tried to do in FE::reinit should have used a per-object cache, and for some incomprehensible reason I made the cache a static variable in the reinit function instead. The bug was only triggered in my application code when using a particular adaptivity scheme in parallel, but there's no telling what effects it could have on a different app. Anyway, anyone who's using CVS should probably either fix that caching by hand or do a CVS update; the risks of using undertested new System::project_vector code are probably outweighed by the risks of using broken FE::reinit code. --- Roy |
From: Roy S. <roy...@ic...> - 2006-11-09 14:21:32
|
On Thu, 9 Nov 2006, Tim Kr=F6ger wrote: > What do I have to fill the _send_list with? Could you give an example = for=20 > that? Every degree of freedom index; what I did was to comment out the body of sort_send_list() and replace it with:o _send_list.resize(n_dofs()); for (unsigned int i=3D0; i!=3D n_dofs(); ++i) _send_list[i] =3D i; > Since your fix of the first bug alone does not seem to cure the=20 > problem (in particular, the same test program as before does even stran= ger=20 > things than before), it would be nice to have at least a workaround so = that I=20 > can go on with my work. I wish I could offer you one, but it looks like my "fix" was broken. I'll let you know if I can find the problem this morning, but if not I'll have to revert the change before I break too many other cvs junkies' libMesh installations. --- Roy |
From: Roy S. <roy...@ic...> - 2006-11-09 15:30:25
|
On Thu, 9 Nov 2006, Roy Stogner wrote: > I wish I could offer you one, but it looks like my "fix" was broken. > I'll let you know if I can find the problem this morning, but if not > I'll have to revert the change before I break too many other cvs > junkies' libMesh installations. Or maybe I'll split things right down the middle, fixing the bug in my patch but simultaneously discovering that it doesn't yet work correctly for vectors other than the primary solution... --- Roy |
From: Roy S. <roy...@ic...> - 2006-11-09 18:47:30
|
On Thu, 9 Nov 2006, Roy Stogner wrote: > On Thu, 9 Nov 2006, Roy Stogner wrote: > >> I wish I could offer you one, but it looks like my "fix" was broken. >> I'll let you know if I can find the problem this morning, but if not >> I'll have to revert the change before I break too many other cvs >> junkies' libMesh installations. > > Or maybe I'll split things right down the middle, fixing the bug in my > patch but simultaneously discovering that it doesn't yet work > correctly for vectors other than the primary solution... With the DofMap send_list filled, the code I just committed to CVS works on your testcase (even if you don't use Ben's solution.set workaround in the user code). Give it a shot. --- Roy |
From: <ti...@ce...> - 2006-11-10 07:29:29
|
Dear Roy, On Thu, 9 Nov 2006, Roy Stogner wrote: > With the DofMap send_list filled, the code I just committed to CVS works on > your testcase (even if you don't use Ben's solution.set workaround in > the user code). Give it a shot. Thank you very much for your help. Everything seems to work now. Please let me know when the bug is completely fixed and the workaround can thus be removed. By the way, how much does the workaround slow down the code? Best Regards, Tim |
From: Roy S. <roy...@ic...> - 2006-11-10 07:46:10
|
On Fri, 10 Nov 2006, Tim Kr=F6ger wrote: > On Thu, 9 Nov 2006, Roy Stogner wrote: > >> With the DofMap send_list filled, the code I just committed to CVS wor= ks on >> your testcase (even if you don't use Ben's solution.set workaround in >> the user code). Give it a shot. > > Thank you very much for your help. Everything seems to work now. Pleas= e let=20 > me know when the bug is completely fixed and the workaround can thus be= =20 > removed. By the way, how much does the workaround slow down the code? Hard to say. The cost of the full send_list will depend on your application and system, since it doesn't add any FLOPS but it may add a lot of extraneous bandwidth. At the same time as committing the half-fix for your problem, though, I added a call to enforce_constraints_exactly() to the end of project_vector(), and that does add a few FLOPS. I'm still not sure whether we should leave that in, since it shouldn't be necessary for your fix. On the one hand, if adaptively coarsening your mesh can leave you with unexpected discontinuities at hanging nodes - well, unexpected discontinuities can be a very very bad thing. On the other hand, for many users those discontinuities will be patched up by the next linear solve, which will proceed fine since your basic time-stepping schemes only require the previous timestep's solution to be in L2; so why force everybody's code to do spend even a couple percent of its runtime doing potentially unnecessary work? After the time I spent hunting this bug down, though, I'm leaning toward "fix up everyone's hanging node constraints whether it's necessary or not". It's easier to optimize inefficient correct code than it is to fix efficient incorrect code. --- Roy |