## libmesh-devel

 Re: [Libmesh-devel] Re=0 From: Roy Stogner - 2006-06-06 02:28:23 ```On Sun, 4 Jun 2006, vikramvgarg@... wrote: > Quoting Roy Stogner : > >> On Sun, 4 Jun 2006, vikramvgarg@... wrote: >> >>> Ex 13 wont work with Re=0 unless the pressure terms are non-zero. >> >> Explain what you mean by "pressure terms are non-zero"? Do you mean >> that pinning a pressure node fixes things? Or do you mean that you've >> got Re attached to a pressure term? - that wouldn't be correct. > > I mean that if I solve just vector Laplacian u = 0 by turning the > appropriate K matrix entries to zero the I get the u = v = 0 solution as the > other day. If the Kpu Kpv etc terms are non-zero the solver works fine. Okay, then - the problem is just that you're turning off too much. The Reynolds number is a coefficient of the convection term and the corresponding Newton terms, but it doesn't scale the pressure term or the continuity equation. >> In the meantime, I'd be curious to know more about how this bug >> behaves. What happens if you set Re=0.0001, for example? Do the >> solutions at 0.001 and 0.0001 look essentially the same, like they >> should? > Re = 0.001 > Works fine > Re = 0.0001 > Works fine, solution similar to Re = 0.001 > Re = 0.00001 > [0]PETSC ERROR: Detected zero pivot in LU factorization! > Re = 0.00001 with zero pivot tolerance = 1.e-30 > Works fine, correct solution > Re = 0 > Will not work even with zero pivot tolerance = 1.e-50 > > It seems to be a solver issue. Okay, this just sounds like the unpinned pressure mode is screwing up your preconditioner. I didn't realize that was a Reynolds number dependent problem, but it is an expected problem. We probably ought to explicitly pin a pressure node in example 13, rather than passing the buck on to the linear solver. >>> The same is true for my steady state and transient codes. > > About matching those two, well it turns out that the velocity is extremely > transient initially and then quickly settles to a steady state, so by using a > very small timestep (0.001) for the first 25 steps and then switching to dt = > .1 I was able to get the correct steady state solution. What happens if dt is too large initially? This is odd; if you use theta=1, then as dt -> infinity you should be solving the steady state problem. An implicit time-stepping scheme shouldn't be so sensitive to dt. --- Roy ```
 Re: [Libmesh-devel] Re=0 From: Roy Stogner - 2006-06-06 02:59:20 ```On Sun, 4 Jun 2006, vikramvgarg@... wrote: > Ex 13 wont work with Re=0 unless the pressure terms are non-zero. Explain what you mean by "pressure terms are non-zero"? Do you mean that pinning a pressure node fixes things? Or do you mean that you've got Re attached to a pressure term? - that wouldn't be correct. And I assume that by "won't work" you mean the same thing we were seeing with your code the other day - that the solution becomes u = v = 0? I do wonder if there's something incorrect in the ex13 Navier-Stokes formulation. I'm currently refactoring ex13 to use the "element_residual()" type implementation which we've been talking about putting in libMesh for years, and the bilinear form I saw didn't quite look right. I'm used to solving for deltau and getting rid of the pressure variable entirely with div-free elements, though, so I'd just assumed that the u^new solution and velocity-pressure formulation had confused me. On the other hand, it's possible that you're just not turning off the right set of terms when you go to Re=0. You're getting rid of both the convection terms and their associated Newton terms, right? I'm going to copy this to libmesh-devel, in case Ben or John or someone wants to try switching ex13 to Stokes flow and see if your bug can be duplicated. In the meantime, I'd be curious to know more about how this bug behaves. What happens if you set Re=0.0001, for example? Do the solutions at 0.001 and 0.0001 look essentially the same, like they should? > The same is true for my steady state and transient codes. Well, you're basically building straight off of ex13 with more interesting boundary conditions, right? So that's not surprising. It was a good idea to go back to ex13 itself to check that. --- Roy ```
 Re: [Libmesh-devel] Re=0 From: Roy Stogner - 2006-06-06 05:35:43 ```I'm not sure what's going haywire with the mailing list, but I sent this stuff like 36 hours ago! Hmmm... maybe I should say, "I sent this stuff on Sunday!", otherwise my complaint won't make sense when it's finally received 36 hours from now... --- Roy On Sun, 4 Jun 2006, Roy Stogner wrote: > On Sun, 4 Jun 2006, vikramvgarg@... wrote: > >> Quoting Roy Stogner : >> >>> On Sun, 4 Jun 2006, vikramvgarg@... wrote: >>> >>>> Ex 13 wont work with Re=0 unless the pressure terms are non-zero. >>> >>> Explain what you mean by "pressure terms are non-zero"? Do you mean >>> that pinning a pressure node fixes things? Or do you mean that you've >>> got Re attached to a pressure term? - that wouldn't be correct. >> >> I mean that if I solve just vector Laplacian u = 0 by turning the >> appropriate K matrix entries to zero the I get the u = v = 0 solution as the >> other day. If the Kpu Kpv etc terms are non-zero the solver works fine. > > Okay, then - the problem is just that you're turning off too much. > The Reynolds number is a coefficient of the convection term and the > corresponding Newton terms, but it doesn't scale the pressure term or > the continuity equation. > >>> In the meantime, I'd be curious to know more about how this bug >>> behaves. What happens if you set Re=0.0001, for example? Do the >>> solutions at 0.001 and 0.0001 look essentially the same, like they >>> should? > >> Re = 0.001 >> Works fine >> Re = 0.0001 >> Works fine, solution similar to Re = 0.001 >> Re = 0.00001 >> [0]PETSC ERROR: Detected zero pivot in LU factorization! >> Re = 0.00001 with zero pivot tolerance = 1.e-30 >> Works fine, correct solution >> Re = 0 >> Will not work even with zero pivot tolerance = 1.e-50 >> >> It seems to be a solver issue. > > Okay, this just sounds like the unpinned pressure mode is screwing up > your preconditioner. I didn't realize that was a Reynolds number > dependent problem, but it is an expected problem. > > We probably ought to explicitly pin a pressure node in example 13, > rather than passing the buck on to the linear solver. > >>>> The same is true for my steady state and transient codes. >> >> About matching those two, well it turns out that the velocity is extremely >> transient initially and then quickly settles to a steady state, so by using a >> very small timestep (0.001) for the first 25 steps and then switching to dt = >> .1 I was able to get the correct steady state solution. > > What happens if dt is too large initially? This is odd; if you use > theta=1, then as dt -> infinity you should be solving the steady state > problem. An implicit time-stepping scheme shouldn't be so sensitive > to dt. > --- > Roy > > > _______________________________________________ > Libmesh-devel mailing list > Libmesh-devel@... > https://lists.sourceforge.net/lists/listinfo/libmesh-devel > ```
 Re: [Libmesh-devel] Re=0 From: John Peterson - 2006-06-06 05:54:22 ```Unknown, but they recently upgraded the mailing list management software. -J Roy Stogner writes: > > I'm not sure what's going haywire with the mailing list, but I sent > this stuff like 36 hours ago! > > Hmmm... maybe I should say, "I sent this stuff on Sunday!", otherwise > my complaint won't make sense when it's finally received 36 hours from > now... > --- > Roy > > On Sun, 4 Jun 2006, Roy Stogner wrote: > > > On Sun, 4 Jun 2006, vikramvgarg@... wrote: > > > >> Quoting Roy Stogner : > >> > >>> On Sun, 4 Jun 2006, vikramvgarg@... wrote: > >>> > >>>> Ex 13 wont work with Re=0 unless the pressure terms are non-zero. > >>> > >>> Explain what you mean by "pressure terms are non-zero"? Do you mean > >>> that pinning a pressure node fixes things? Or do you mean that you've > >>> got Re attached to a pressure term? - that wouldn't be correct. > >> > >> I mean that if I solve just vector Laplacian u = 0 by turning the > >> appropriate K matrix entries to zero the I get the u = v = 0 solution as the > >> other day. If the Kpu Kpv etc terms are non-zero the solver works fine. > > > > Okay, then - the problem is just that you're turning off too much. > > The Reynolds number is a coefficient of the convection term and the > > corresponding Newton terms, but it doesn't scale the pressure term or > > the continuity equation. > > > >>> In the meantime, I'd be curious to know more about how this bug > >>> behaves. What happens if you set Re=0.0001, for example? Do the > >>> solutions at 0.001 and 0.0001 look essentially the same, like they > >>> should? > > > >> Re = 0.001 > >> Works fine > >> Re = 0.0001 > >> Works fine, solution similar to Re = 0.001 > >> Re = 0.00001 > >> [0]PETSC ERROR: Detected zero pivot in LU factorization! > >> Re = 0.00001 with zero pivot tolerance = 1.e-30 > >> Works fine, correct solution > >> Re = 0 > >> Will not work even with zero pivot tolerance = 1.e-50 > >> > >> It seems to be a solver issue. > > > > Okay, this just sounds like the unpinned pressure mode is screwing up > > your preconditioner. I didn't realize that was a Reynolds number > > dependent problem, but it is an expected problem. > > > > We probably ought to explicitly pin a pressure node in example 13, > > rather than passing the buck on to the linear solver. > > > >>>> The same is true for my steady state and transient codes. > >> > >> About matching those two, well it turns out that the velocity is extremely > >> transient initially and then quickly settles to a steady state, so by using a > >> very small timestep (0.001) for the first 25 steps and then switching to dt = > >> .1 I was able to get the correct steady state solution. > > > > What happens if dt is too large initially? This is odd; if you use > > theta=1, then as dt -> infinity you should be solving the steady state > > problem. An implicit time-stepping scheme shouldn't be so sensitive > > to dt. > > --- > > Roy ```