From: Manav Bhatia <bhatiamanav@gm...>  20130514 15:25:26

Hi, I am working on higherorder simulation using the hierarchich function on quad8s. My errorconvergence plots give me theoretical convergence for first and second order p, but the error stagnates for p > 2. I am speculating that this might be due to the lowerorder geometry, assuming that the x and ycoordinates are interpolated using Lagrange functions associated with nodes of quad8. I am considering adding higher order quads to get x and y interpolation using higherorder Lagrange functions. Any thoughts on how easy/difficult this might be? Thanks, Manav 
From: John Peterson <jwpeterson@gm...>  20130514 17:00:02

On Tue, May 14, 2013 at 9:25 AM, Manav Bhatia <bhatiamanav@...> wrote: > Hi, > > I am working on higherorder simulation using the hierarchich function > on quad8s. > > My errorconvergence plots give me theoretical convergence for first > and second order p, but the error stagnates for p > 2. I am speculating > that this might be due to the lowerorder geometry, assuming that the x > and ycoordinates are interpolated using Lagrange functions associated with > nodes of quad8. What does your domain look like? Unit square? If so I don't think there can be errors due to using even bilinear element maps, the Jacobians are linear in this case. What exact solution are you using for testing the convergence? > I am considering adding higher order quads to get x and > y interpolation using higherorder Lagrange functions. Any thoughts on how > easy/difficult this might be? You are talking about QUAD16? It shouldn't be too hard but we would need to come to a consesus on node numbering for this element.  John 
From: Manav Bhatia <bhatiamanav@gm...>  20130514 17:54:36

On Tue, May 14, 2013 at 12:59 PM, John Peterson <jwpeterson@...>wrote: > On Tue, May 14, 2013 at 9:25 AM, Manav Bhatia <bhatiamanav@...> > wrote: > > Hi, > > > > I am working on higherorder simulation using the hierarchich > function > > on quad8s. > > > > My errorconvergence plots give me theoretical convergence for first > > and second order p, but the error stagnates for p > 2. I am speculating > > that this might be due to the lowerorder geometry, assuming that the x > > and ycoordinates are interpolated using Lagrange functions associated > with > > nodes of quad8. > > What does your domain look like? Unit square? If so I don't think > there can be errors due to using even bilinear element maps, the > Jacobians are linear in this case. > > This is the Gaussian bump problem from the higherorder CFD workshop ( http://dept.ku.edu/~cfdku/hiocfd/case_c1.1.html). You are correct that away from the bump the boundary is straight, so linear elements should be fine. I am looking at the entropy error, since the entropy is supposed to stay constant. The bumpboundary, infact, is adding to the entropyerror. I am able to drop down to 10^7 in the error L2 norm, and then it stagnates. And I have a feeling that this is due to the loworder geometry. > What exact solution are you using for testing the convergence? > > > I am considering adding higher order quads to get x and > > y interpolation using higherorder Lagrange functions. Any thoughts on > how > > easy/difficult this might be? > > You are talking about QUAD16? It shouldn't be too hard but we would > need to come to a consesus on node numbering for this element. > > Yes, I was considering this, and also the 25noded quad. However, I am still considering if the effort might be worth it: Meaning that I may be able to get theoretical order of convergence for this simpler benchmark problem, but I don't know if any practical problem will benefit from a 16 or 25 noded quad. I am not sure if any mesh generator would give me elements with these many nodes. Manav 
From: Roy Stogner <roystgnr@ic...>  20130514 18:11:45

On Tue, 14 May 2013, Manav Bhatia wrote: > This is the Gaussian bump problem from the higherorder CFD workshop ( > http://dept.ku.edu/~cfdku/hiocfd/case_c1.1.html). You are correct that away > from the bump the boundary is straight, so linear elements should be fine. > I am looking at the entropy error, since the entropy is supposed to stay > constant. The bumpboundary, infact, is adding to the entropyerror. I am > able to drop down to 10^7 in the error L2 norm, and then it stagnates. And > I have a feeling that this is due to the loworder geometry. Hmm... the loworder geometry should definitely be hurting your convergence rate here, but it shouldn't be giving you a rate of *zero*. Usually when you start to converge but then stagnate it means you're hitting some previouslynegligible epsilon. I've never seen anything but a penalty BC parameter cause stagnation at a *relative* error as high as 1e7 though. Usually the second limitation I run into is insufficiently tight solver tolerances, and the third is FP roundoff error. From our point of view I don't want to discourage you from the QUAD16 idea, but I don't think it's likely to be the cause of your problem, and it would probably be easy for you to check your solver parameters or to recompile with quad precision before going to the effort of adding a new geometric element to libMesh.  Roy 
From: John Peterson <jwpeterson@gm...>  20130514 18:12:30

On Tue, May 14, 2013 at 11:54 AM, Manav Bhatia <bhatiamanav@...> wrote: > > > On Tue, May 14, 2013 at 12:59 PM, John Peterson <jwpeterson@...> > wrote: >> >> On Tue, May 14, 2013 at 9:25 AM, Manav Bhatia <bhatiamanav@...> >> wrote: >> > Hi, >> > >> > I am working on higherorder simulation using the hierarchich >> > function >> > on quad8s. >> > >> > My errorconvergence plots give me theoretical convergence for first >> > and second order p, but the error stagnates for p > 2. I am speculating >> > that this might be due to the lowerorder geometry, assuming that the x >> > and ycoordinates are interpolated using Lagrange functions associated >> > with >> > nodes of quad8. >> >> What does your domain look like? Unit square? If so I don't think >> there can be errors due to using even bilinear element maps, the >> Jacobians are linear in this case. >> > > This is the Gaussian bump problem from the higherorder CFD workshop > (http://dept.ku.edu/~cfdku/hiocfd/case_c1.1.html). You are correct that away > from the bump the boundary is straight, so linear elements should be fine. I > am looking at the entropy error, since the entropy is supposed to stay > constant. The bumpboundary, infact, is adding to the entropyerror. I am > able to drop down to 10^7 in the error L2 norm, and then it stagnates. And > I have a feeling that this is due to the loworder geometry. Oh yeah, I have heard of this spurious entropy production problem on curved geometries. Are you really using Quad8's? Can you use Quad9's with the hierarchics? Can you modify the test problem slightly so the bump is a quadratic function and verify convergence with higher p's? > Yes, I was considering this, and also the 25noded quad. However, I am still > considering if the effort might be worth it: Meaning that I may be able to > get theoretical order of convergence for this simpler benchmark problem, but > I don't know if any practical problem will benefit from a 16 or 25 noded > quad. I am not sure if any mesh generator would give me elements with these > many nodes. You will probably have trouble viewing the solutions in Paraview too...  John 
From: Manav Bhatia <bhatiamanav@gm...>  20130514 18:37:59

> > > > > This is the Gaussian bump problem from the higherorder CFD workshop > > (http://dept.ku.edu/~cfdku/hiocfd/case_c1.1.html). You are correct that > away > > from the bump the boundary is straight, so linear elements should be > fine. I > > am looking at the entropy error, since the entropy is supposed to stay > > constant. The bumpboundary, infact, is adding to the entropyerror. I am > > able to drop down to 10^7 in the error L2 norm, and then it stagnates. > And > > I have a feeling that this is due to the loworder geometry. > > Oh yeah, I have heard of this spurious entropy production problem on > curved geometries. > > Are you really using Quad8's? Can you use Quad9's with the hierarchics? > > Can you modify the test problem slightly so the bump is a quadratic > function and verify convergence with higher p's? > > This is a great suggestion, John. I will try converting the geometry to second order. I have not tried Quad9, but what would that be any better than Quad8? > > Yes, I was considering this, and also the 25noded quad. However, I am > still > > considering if the effort might be worth it: Meaning that I may be able > to > > get theoretical order of convergence for this simpler benchmark problem, > but > > I don't know if any practical problem will benefit from a 16 or 25 noded > > quad. I am not sure if any mesh generator would give me elements with > these > > many nodes. > > You will probably have trouble viewing the solutions in Paraview too... > As it is, I am having trouble viewing results from higherorder elements in Paraview. I can output only the nodal data for viewing in Paraview, and so the higherorder information of the element solution gets lost. I haven't yet had the chance to look into improving this behavior. Roy: Just to clarify, I do get a reduction in error for higher order elemets after 10^7, but the rate of convergence is considerably lower than the p+1/2 theoretical order. I will try compiling with quad precision to see if that does something. I have played around with the solver parameters, and have also used direct solvers to flush out any potential problems form linear solvers. 
From: John Peterson <jwpeterson@gm...>  20130515 01:01:50

On May 14, 2013, at 12:37 PM, Manav Bhatia <bhatiamanav@...> wrote: >> > >> > This is the Gaussian bump problem from the higherorder CFD workshop >> > (http://dept.ku.edu/~cfdku/hiocfd/case_c1.1.html). You are correct that away >> > from the bump the boundary is straight, so linear elements should be fine. I >> > am looking at the entropy error, since the entropy is supposed to stay >> > constant. The bumpboundary, infact, is adding to the entropyerror. I am >> > able to drop down to 10^7 in the error L2 norm, and then it stagnates. And >> > I have a feeling that this is due to the loworder geometry. >> >> Oh yeah, I have heard of this spurious entropy production problem on >> curved geometries. >> >> Are you really using Quad8's? Can you use Quad9's with the hierarchics? >> >> Can you modify the test problem slightly so the bump is a quadratic >> function and verify convergence with higher p's? > > This is a great suggestion, John. I will try converting the geometry to second order. > > I have not tried Quad9, but what would that be any better than Quad8? Not much other than it would have the x^2 y^2 term in its map. One other thought: do you have some O(h) stabilization terms present? Then if you aren't refining the grid they don't get smaller... >> > Yes, I was considering this, and also the 25noded quad. However, I am still >> > considering if the effort might be worth it: Meaning that I may be able to >> > get theoretical order of convergence for this simpler benchmark problem, but >> > I don't know if any practical problem will benefit from a 16 or 25 noded >> > quad. I am not sure if any mesh generator would give me elements with these >> > many nodes. >> >> You will probably have trouble viewing the solutions in Paraview too... > As it is, I am having trouble viewing results from higherorder elements in Paraview. I can output only the nodal data for viewing in Paraview, and so the higherorder information of the element solution gets lost. I haven't yet had the chance to look into improving this behavior. > > > Roy: Just to clarify, I do get a reduction in error for higher order elemets after 10^7, but the rate of convergence is considerably lower than the p+1/2 theoretical order. I will try compiling with quad precision to see if that does something. I have played around with the solver parameters, and have also used direct solvers to flush out any potential problems form linear solvers. > > > 
From: Manav Bhatia <bhatiamanav@gm...>  20130515 12:36:49

> > > > One other thought: do you have some O(h) stabilization terms present? > Then if you aren't refining the grid they don't get smaller... > > There is a GLS stabilization term with the typical tau matrix. I have experimented with different tau definitions, but get the same behavior. Essentially, almost all definitions in 1D have the form (sum_{i=1,n_shape_funcs}  a dNi/dx  ) ^1 where N is the shape function and 'a' is the velocity. The numeric value of this expression reduces as the polynomial order is increased. So, I think that the influence of the polynomial order is included in the stabilization matrix. Manav 
From: John Peterson <jwpeterson@gm...>  20130515 17:18:12

On Wed, May 15, 2013 at 6:36 AM, Manav Bhatia <bhatiamanav@...> wrote: >> >> >> One other thought: do you have some O(h) stabilization terms present? >> Then if you aren't refining the grid they don't get smaller... >> > > > There is a GLS stabilization term with the typical tau matrix. I have > experimented with different tau definitions, but get the same behavior. > > Essentially, almost all definitions in 1D have the form > > (sum_{i=1,n_shape_funcs}  a dNi/dx  ) ^1 > > where N is the shape function and 'a' is the velocity. The numeric value of > this expression reduces as the polynomial order is increased. Really? The lowest order hierarchics are the Lagrange basis functions, and dNi/dx ~ 1/h is certainly true for those.  John 
From: Manav Bhatia <bhatiamanav@gm...>  20130515 17:58:48

On Wed, May 15, 2013 at 1:17 PM, John Peterson <jwpeterson@...> wrote: > On Wed, May 15, 2013 at 6:36 AM, Manav Bhatia <bhatiamanav@...> > wrote: > >> > >> > >> One other thought: do you have some O(h) stabilization terms present? > >> Then if you aren't refining the grid they don't get smaller... > >> > > > > > > There is a GLS stabilization term with the typical tau matrix. I have > > experimented with different tau definitions, but get the same behavior. > > > > Essentially, almost all definitions in 1D have the form > > > > (sum_{i=1,n_shape_funcs}  a dNi/dx  ) ^1 > > > > where N is the shape function and 'a' is the velocity. The numeric value > of > > this expression reduces as the polynomial order is increased. > > Really? The lowest order hierarchics are the Lagrange basis > functions, and dNi/dx ~ 1/h is certainly true for those. > Hi John, I have attached a pdf which plots the expression for a 2noded and a 3noded element using Lagrange interpolation functions. The element has unit length. The 2noded value is constant at 2 for the entire domain of the element, while the 3noded element shows a variation in the element domain. It equals 2 at xi=0, but rises to about 8 on either end. Each quadrature point uses the inverse of this function value. This is my basis for stating that the the influence of higher polynomial order is included, although I can't say if this is constraining the order of convergence. Manav 
From: lorenzo alessio botti <lorenzoalessiobotti@gm...>  20130515 19:54:39

I'm pretty sure that the geometry approximation is responsible of the convergence degradation. I've experienced this problem with the ringleb flow (2D Euler) using inviscid wall bcs. If you do p refinement instead of h refinement you should see that the entropy error actually stagnates and you get almost zero convergence. This is due to the error in the geometry representation. Imposing the bcs on the approximated geometry according to the exact solution of the ringleb flow allows to recover optimal convergence rates. The geometry approximation error is neglected. Unfortunately I don't think that the bump has an exact solution, does it? Otherwise you could try this way. Lorenzo On May 15, 2013 7:58 PM, "Manav Bhatia" <bhatiamanav@...> wrote: > On Wed, May 15, 2013 at 1:17 PM, John Peterson <jwpeterson@...> > wrote: > > > On Wed, May 15, 2013 at 6:36 AM, Manav Bhatia <bhatiamanav@...> > > wrote: > > >> > > >> > > >> One other thought: do you have some O(h) stabilization terms present? > > >> Then if you aren't refining the grid they don't get smaller... > > >> > > > > > > > > > There is a GLS stabilization term with the typical tau matrix. I have > > > experimented with different tau definitions, but get the same behavior. > > > > > > Essentially, almost all definitions in 1D have the form > > > > > > (sum_{i=1,n_shape_funcs}  a dNi/dx  ) ^1 > > > > > > where N is the shape function and 'a' is the velocity. The numeric > value > > of > > > this expression reduces as the polynomial order is increased. > > > > Really? The lowest order hierarchics are the Lagrange basis > > functions, and dNi/dx ~ 1/h is certainly true for those. > > > > Hi John, > > I have attached a pdf which plots the expression for a 2noded and a > 3noded element using Lagrange interpolation functions. The element has > unit length. The 2noded value is constant at 2 for the entire domain of > the element, while the 3noded element shows a variation in the element > domain. It equals 2 at xi=0, but rises to about 8 on either end. Each > quadrature point uses the inverse of this function value. > > This is my basis for stating that the the influence of higher > polynomial order is included, although I can't say if this is constraining > the order of convergence. > > Manav > > >  > AlienVault Unified Security Management (USM) platform delivers complete > security visibility with the essential security capabilities. Easily and > efficiently configure, manage, and operate all of your security controls > from a single console and one unified framework. Download a free trial. > http://p.sf.net/sfu/alienvault_d2d > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers > > 
From: Manav Bhatia <bhatiamanav@gm...>  20130515 20:18:12

On Wed, May 15, 2013 at 3:54 PM, lorenzo alessio botti < lorenzoalessiobotti@...> wrote: > I'm pretty sure that the geometry approximation is responsible of the > convergence degradation. I've experienced this problem with the ringleb > flow (2D Euler) using inviscid wall bcs. If you do p refinement instead of > h refinement you should see that the entropy error actually stagnates and > you get almost zero convergence. This is due to the error in the geometry > representation. Imposing the bcs on the approximated geometry according to > the exact solution of the ringleb flow allows to recover optimal > convergence rates. The geometry approximation error is neglected. > Unfortunately I don't think that the bump has an exact solution, does it? > Otherwise you could try this way. > Hi Lorenzo, Thanks for sharing your experience in this area. The bump does not have an exact solution, so this approach will be difficult. I did put together the ringleb case, so I will experiment with your suggestion. Manav 
From: John Peterson <jwpeterson@gm...>  20130515 20:13:50

On Wed, May 15, 2013 at 11:58 AM, Manav Bhatia <bhatiamanav@...> wrote: > > I have attached a pdf which plots the expression for a 2noded and a > 3noded element using Lagrange interpolation functions. The element has unit > length. The 2noded value is constant at 2 for the entire domain of the > element, while the 3noded element shows a variation in the element domain. > It equals 2 at xi=0, but rises to about 8 on either end. Each quadrature > point uses the inverse of this function value. (I doubt stabilization is actually causing the pconvergence problems you're seeing, but just to clarify what I was getting at...) You plotted dNi/d(\xi) for these reference elements, where by "\xi" I mean the reference domain coordinate. But to compute the dNi/dx which is actually used in your formula, you have to multiply dNi/d(\xi) by the inverse jacobian, e.g. dNi/dx = dNi/d(\xi) * d(\xi)/dx For example: in 1D, a linear Lagrange element has Jacobian dx/d(\xi) = h/2, inverse d(\xi)/dx = 2/h, and therefore dNi/dx = (2/h) * dNi/d(\xi) which is O(1/h).  John 
From: Manav Bhatia <bhatiamanav@gm...>  20130515 20:27:48

On Wed, May 15, 2013 at 4:13 PM, John Peterson <jwpeterson@...> wrote: > On Wed, May 15, 2013 at 11:58 AM, Manav Bhatia <bhatiamanav@...> > wrote: > > > > I have attached a pdf which plots the expression for a 2noded and a > > 3noded element using Lagrange interpolation functions. The element has > unit > > length. The 2noded value is constant at 2 for the entire domain of the > > element, while the 3noded element shows a variation in the element > domain. > > It equals 2 at xi=0, but rises to about 8 on either end. Each quadrature > > point uses the inverse of this function value. > > (I doubt stabilization is actually causing the pconvergence problems > you're seeing, but just to clarify what I was getting at...) > > > You plotted dNi/d(\xi) for these reference elements, where by "\xi" I > mean the reference domain coordinate. > > But to compute the dNi/dx which is actually used in your formula, you > have to multiply dNi/d(\xi) by the inverse jacobian, e.g. > > dNi/dx = dNi/d(\xi) * d(\xi)/dx > > For example: in 1D, a linear Lagrange element has Jacobian dx/d(\xi) = > h/2, inverse d(\xi)/dx = 2/h, and therefore > > dNi/dx = (2/h) * dNi/d(\xi) > > which is O(1/h). > > Hi John, It is dNi/dx that I plotted. Actually, for this 1D case, the dNi/dxi and dNi/dx are simply scaled variants of each other by a constant (the factor 2/h that you mentioned). I do follow your argument for 1/h, but varying polynomial order also varies the tau value (reducing values of tau with higher p). How would that reflect on identifying the horderdependence? Or would it stay consistent with h^1 for tau, no matter what p is used? Manav 
From: John Peterson <jwpeterson@gm...>  20130515 21:20:33

On Wed, May 15, 2013 at 2:27 PM, Manav Bhatia <bhatiamanav@...> wrote: > > It is dNi/dx that I plotted. Actually, for this 1D case, the dNi/dxi > and dNi/dx are simply scaled variants of each other by a constant (the > factor 2/h that you mentioned). The point being that the inverse of this "constant" gets smaller if you refine the grid, and stays the same if you don't. > I do follow your argument for 1/h, but varying polynomial order also > varies the tau value (reducing values of tau with higher p). How would that > reflect on identifying the horderdependence? Or would it stay consistent > with h^1 for tau, no matter what p is used? That I don't know. Since you are using a consistent scheme like GLS, your stabilization parameter is proportional the strong residual as well. Perhaps one can argue that the strong residual gets smaller at some rate like h^p, so that the overall stabilization term goes like h^{p+1}, and then you'd be fine.  John 