From: John P. <pet...@cf...> - 2007-05-04 17:33:49
|
Hi, I've done some numerical experiments with 3D Navier-Stokes and hanging edges. I'm using prisms since they are simpler to visualize. Here's the test mesh used: http://www.cfdlab.ae.utexas.edu/~peterson/hanging_prisms_full_domain.png The boundary condition is a C1-continuous lid velocity which decays to zero at the edge. This avoids pressure singularities and discontinuous boundary conditions. The edges in the center of the box have a level-2 mismatch even though level-1 face mismatch is enforced throughout. On coarse grids, I initially noticed some slightly higher velocities along the level-2 mismatched edge. Here's a y-cutplane of u: http://www.cfdlab.ae.utexas.edu/~peterson/hanging_prisms_u.png We can't really interpret this too closely, since these are quadratic elements plotted as linears. Under uniform refinement, this region of lower-than-average velocity is not present. http://www.cfdlab.ae.utexas.edu/~peterson/hanging_prisms_x2_u.png To get some idea of the relative sizes of the discontinuity, I also plotted as a surface here (again, quadratics plotted as 4 linears!): http://www.cfdlab.ae.utexas.edu/~peterson/discontinuous_u_surface.png Obviously, there aren't any problems when u=const across the interface. The (potentially) more interesting thing appears to be the linear/nonlinear solver convergence when such a mismatched edge is present. Here, I present the relative and absolute Newton step sizes, and the nonlinear and linear system residuals for the first timestep of the once uniformly-refined version of the mesh with the level-2 mismatch hanging edge: |du|/|u| |du| nonlinear |R| linear |R| 1 114.377 2.32242e+08 7.74163e-11 0.376716 41.3106 0.000908376 7.27909e-12 0.236782 27.9334 3.08814e-05 4.52935e-13 0.0958868 11.8114 7.95717e-06 6.28991e-14 0.040872 5.13872 3.45092e-06 2.43645e-14 0.0176168 2.23529 1.50099e-06 1.30614e-14 0.0076324 0.972337 6.52953e-07 5.33396e-15 0.0033142 0.422959 2.841e-07 2.29361e-15 0.00144055 0.183984 1.23752e-07 9.92329e-16 5.40359e-08 4.30634e-16 The poor convergence of Newton does not appear to be due to insufficient solution of the linear subproblems. This is the first timestep, so it is the hardest solve. Timesteps nearer steady state start with smaller initial residuals but still converge poorly. For comparison, I also ran the problem with the opposing central prisms also refined: http://www.cfdlab.ae.utexas.edu/~peterson/matched_prisms_x.png In this case, the convergence of the solver during the first timestep is much better |du|/|u| |du| nonlinear |R| linear |R| 1 158.132 1.54155e+08 8.68055e-11 0.171657 26.5923 0.000375404 6.72722e-12 0.00567181 0.878692 9.36429e-06 3.1895e-13 6.53108e-09 2.82395e-16 So while I don't find any evidence of spurious velocities "blowing up" and trashing the solution, I do think it is highly possible that the level-2 mismatched edge is causing Newton to fail to converge, possibly due to ill-conditioning of the linear systems (?) -John Lorenzo Botti writes: > >.) Were the spurious velocities coming particularly from the > > non-level-1 constrained nodes? Or were they arising from hanging > > (edge or face) nodes? > > It's hard to tell right away... It seems that the problem involves a small > region surrounding a non level-one edge, > but it's tricky to tell for sure just using gmv. The point is that the > problem doesn't seem to be systematic. > What I can tell is that there are non-level-1 constrained nodes that don't > involve spurious velocity vectors, but > without limiting nodal level mismatch we always find a region were the > solution blows up. > I'll investigate this and I'll send you some screenshots... > > >.) (Curiosity) How many more elements does flagging based on level-1 > > at the nodes introduce vs. not maintaining level-1 at the nodes? > > It seems that your mesh size would blow up quickly doing this. > > Not so quickly but actually with my 2Gb of ram i can't go further 80000 > elements in my simulations > (using GMRES/ILU with a restart of 400). > The resulting refined mesh is much more uniform than without > level-1-mismatch control, with an > increase of 10000/15000 elements over 80000. > > > Possibility 1 should be easy enough to test - run a > > DiscontinuityMeasure over your solutions and make sure the > > amount discontinuity it sees is on the order of floating point error. > > I'll try this... > Thanks a lot for your help. > > Lorenzo > > 2007/5/3, John Peterson <pet...@cf...>: > > > > Roy Stogner writes: > > > On Thu, 3 May 2007, John Peterson wrote: > > > > > > > I'm glad you found that this works, but I wonder if it means there is > > > > a problem with the way we are constraining hanging "face" nodes on > > > > Tet10's, and this fix is just masking the bigger problem. > > > > > > It sounds like there's two possible problems on tets here. Either: > > > > > > 1. Our hanging nodes aren't producing conforming function spaces. > > > > > > Or: > > > > > > 2. The function spaces you get by only enforcing level-one > > > restrictions on faces aren't well-suited for Navier-Stokes. > > > > > > Possibility 1 should be easy enough to test - run a > > > DiscontinuityMeasure over your solutions and make sure the > > > amount discontinuity it sees is on the order of floating point error. > > > > > > For possibility 2, we might want to make "maintain_nodal_level_one" an > > > option along with "maintain_level_one". I don't want to make the > > > former behavior standard, since it would be a backwards-incompatible > > > behavior change and it might be too expensive for many users. However > > > making it an option would be a good idea. I'm thinking we deprecate > > > the maintain_level_one method argument and then make both options into > > > private MeshRefinement variables with accessor methods, like > > > refine_fraction and such. > > > > This is a good solution in the meantime, IMO. I'd still like to know > > what's going on with Navier-Stokes ... including why this doesn't show > > up (so far) on hexes. It seems that just having "fewer" > > greater-than-level-1 mismatched nodes in Hex meshes wouldn't be enough > > to prevent the problem, i.e. *any* greater-than-level-1 mismatched > > nodes should seemingly trigger it. > > > > I never got around to setting up a smooth (with smooth BCs) 3D > > Stokes/Navier-Stokes test problem. It seems that would be an important > > first step in testing your hypotheses above. If anyone has a suggestion > > I'd be happy to try it... > > > > -J > > |
From: Roy S. <roy...@ic...> - 2007-05-04 21:23:38
|
On Fri, 4 May 2007, John Peterson wrote: > The poor convergence of Newton does not appear to be due > to insufficient solution of the linear subproblems. This still seems to suggest a bug to me, but in our repeated testing we can only get similarly odd behavior out of the Stokes problem, not out of the Poisson problem on the same mesh. It's possible that there is no bug but that there's some LBB violation here which is making the velocity-pressure system underdefined or extraordinarily ill-conditioned; we'll try looking into that. In the meantime, we've changed the MeshRefinement API somewhat to try and make it easier to avoid meshes which trigger such problems. The method argument "maintain_level_one" has been replaced by accessors to options "face_level_mismatch_limit", "edge_level_mismatch_limit", and "node_level_mismatch_limit", which default to 1, 0 (i.e. "off"), and 0 respectively. Setting edge_level_mismatch_limit to 1 will hopefully prevent these sorts of awkwardly constrained hanging nodes from occuring, but without as much added expense as enforcing a node_level_mismatch_limit would entail. I've committed these options to CVS without giving them much testing, and John and I probably won't give them as much testing as they deserve, so Lorenzo, please give them a try. --- Roy |
From: Roy S. <roy...@ic...> - 2007-05-04 23:59:45
|
On Fri, 4 May 2007, Roy Stogner wrote: > This still seems to suggest a bug to me, but in our repeated testing > we can only get similarly odd behavior out of the Stokes problem, not > out of the Poisson problem on the same mesh. It's possible that there > is no bug but that there's some LBB violation here which is making the > velocity-pressure system underdefined or extraordinarily > ill-conditioned; we'll try looking into that. Interestingly enough, there was a bug (introduced by myself, in fact...) when running with Lagrange elements on meshes without the level-1 rule true across faces. This shouldn't have effected the case where level-1 regularity was true across faces but not edges, however. The bug still exists for non-Lagrange elements on non-level-1 regular meshes, so beware. I'll try to have a general fix posted to CVS this weekend, but it needs more testing. --- Roy |
From: Shun W. <shu...@ui...> - 2007-05-06 04:44:53
|
Hi, I didn't read through the whole post on the issue. But I found some similar problem for Hex8 elements. If you checked the attached the picture. This is a local mesh I get from libMesh AMR. There, both A and B are hanging nodes, while A depends on B. So when I do both dofmap.enforce_constraints_exactly(), A is interpolated from B and C, while B as a hanging node has value 0. Therefore, A gets a wrong value. I am guessing that libMesh level one condition is enforced only on surface, i.e. the level difference of two elements that share surfaces can not exceed one. If this is true, it could be a potential problem for hanging nodes. Is my conjecture right? How can I solve this problem? Thank you very much! -Shun |
From: Roy S. <roy...@ic...> - 2007-05-06 05:03:50
|
On Sat, 5 May 2007, Shun Wang wrote: > If you checked the attached the picture. The picture didn't come through, but I take it you're talking about a situation where the mesh is level-1 regular on faces but there is a level-2 irregularity on edges? > This is a local mesh I get from libMesh AMR. There, both A and B > are hanging nodes, while A depends on B. So when I do both > dofmap.enforce_constraints_exactly(), A is interpolated from B and > C, while B as a hanging node has value 0. Therefore, A gets a > wrong value. Hmm... the bug I found Friday should only have affected meshes which have level-2 irregularities on sides, not just edges. If A is interpolated from B and C, while B is interpolated from (suppose) D and E, then enforce_constraints_exactly() is supposed to recursively build up a constraint matrix which directly interpolates A from C, D, and E. It's of course possible that there are bugs in this process. > I am guessing that libMesh level one condition is enforced only on > surface, i.e. the level difference of two elements that share surfaces can > not exceed one. > Is my conjecture right? How can I solve this problem? There are now options for enforcing level one regularity between edge and/or node neighbors in addition to face neighbors. If level-2 edge irregularities are causing problems for your particular formulation, using either of those options should fix it. However, even if your mesh has *no* regularity, libMesh is still supposed to ensure that your constrained finite element spaces have the appropriate level of continuity. If you're seeing an exception to that rule, it's a bug and we'd appreciate if you can give us more details or a reproducable test case. --- Roy |
From: Roy S. <roy...@ic...> - 2007-03-31 04:09:34
|
On Thu, 29 Mar 2007, Lorenzo Botti wrote: > We're testing our solver on the lid driven cavity problem 3D but we're > experiencing some problems, maybe with hanging nodes constraints...? Well, you're not looping over inactive elements, you haven't forgotten to multiply any terms by JxW, you haven't forgotten the constrain_ call before matrix assembly, and you don't have a nonlinear loop combined with inexact constraint enforcement. That exhausts my list of "things I've seen new users do which break adaptivity". All I can suggest now is that you try to simplify the problem as much as possible; I don't think any of us have time to debug something so complicated unless we can be sure it's a library problem. I only do 3D problems on hexes, so you might try running on a hex mesh to see if we have some bug in AMR on tets. I wouldn't bet on it, though; I believe John uses AMR on tets for some of his double-diffusion code. Oh, and if you haven't already, try the latest CVS head and try running with METHOD=dbg, just to see if there are any bugfixes you need or assert() calls you might trip. And, when you do find the problem, even if it's just a user mistake and not a library bug, let us know. If it's something we can prevent in the future by adding more assert() statements that would be nice, and at the very least it might help future users if we had a more complete "things which can break adaptivity" list. --- Roy |