On reducing spurious oscillations in discontinuous Galerkin methods for convection-diffusion equations
This is a modal window.
The media could not be loaded, either because the server or network failed or because the format is not supported.
Formal Metadata
| Title |
| |
| Title of Series | ||
| Number of Parts | 22 | |
| Author | ||
| License | CC Attribution 3.0 Germany: You are free to use, adapt and copy, distribute and transmit the work or content in adapted or unchanged form for any legal purpose as long as the work is attributed to the author in the manner specified by the author or licensor. | |
| Identifiers | 10.5446/57505 (DOI) | |
| Publisher | ||
| Release Date | ||
| Language |
Content Metadata
| Subject Area | ||
| Genre | ||
| Abstract |
|
Leibniz MMS Days 20228 / 22
11
19
00:00
Galerkin-MethodeNichtlineares GleichungssystemMathematical analysisOscillationReduction of orderStandard errorThermal conductivityTerm (mathematics)MereologyObservational studyFocus (optics)Multiplication signResultantNichtlineares GleichungssystemHydraulic jumpEquationGalerkin-MethodeRight angleSteady state (chemistry)Computer animation
01:25
Time domainNeumann boundary conditionTerm (mathematics)ForceLogical constantEquationEquationMechanism designScaling (geometry)Water vaporTerm (mathematics)Time domainPopulation densityConcentricOctahedronConnected spaceStandard errorSet theoryThermal conductivityDistanceOrder of magnitudeState of matterDirection (geometry)Cartesian coordinate systemMaß <Mathematik>Different (Kate Ryan album)Point (geometry)Model theoryCharge carrierKörper <Algebra>PhysicalismScalar fieldSteady state (chemistry)Right angleNumerical analysisFluidPartial differential equationComputer animation
03:45
CuboidComputabilityTime domainPositional notationComplete metric spaceCircleWater vaporModel theoryMass flow rateBoundary value problemConcentricMereologyDirection (geometry)Mass diffusivityRight angleConnected spaceComputer animation
05:19
Order of magnitudeCircleDirection (geometry)MathematicsDiffuser (automotive)Computer animation
05:51
Basis <Mathematik>Galerkin-MethodeMaxima and minimaBasis <Mathematik>OscillationPoint (geometry)TriangleOrder of magnitudeFunctional (mathematics)Standard deviationCategory of beingOrder (biology)PhysicalismAnalytic continuationConcentricUniformer RaumWeightMilitary baseSet theoryContrast (vision)Matching (graph theory)Maxima and minimaIndependence (probability theory)TorusGraph coloringGreatest elementMathematicsComputer animation
08:36
Galerkin-MethodeMaxima and minimaProcess (computing)PolynomialLimit (category theory)ApproximationNichtlineares GleichungssystemBasis <Mathematik>Functional (mathematics)Analytic continuationSlide ruleOscillationCategory of beingDegree (graph theory)PolynomialStandard deviationContrast (vision)TrianglePhysical systemDifferent (Kate Ryan album)Local ringLinearizationNominal numberOpen setStatistical hypothesis testingSet theoryMilitary baseMusical ensembleComputer animation
11:49
Linear mapKeilförmige AnordnungFunction (mathematics)INTEGRALSquare numberFunctional (mathematics)TriangleQuadratic functionLinearizationPoint (geometry)Hydraulic jumpParabolaArithmetic meanSpecial unitary groupObservational studyProcess (computing)Military baseKontraktion <Mathematik>Computer animation
15:03
Hydraulic jumpINTEGRALSquare numberApproximationLogical constantHydraulic jumpArithmetic meanContrast (vision)TriangleDifferent (Kate Ryan album)Square numberMassOscillationComputer animation
16:03
Hydraulic jumpSquare numberINTEGRALMaxima and minimaPairwise comparisonProcess (computing)Galerkin-MethodeOscillationStandard deviationTime domainMaxima and minimaMortality rateApproximationArithmetic meanContrast (vision)MeasurementWater vaporResultantCircleMereologyRight anglePoint (geometry)Galerkin-MethodeReal numberMaxwell's equationsDistribution (mathematics)RectangleMass flow rateWell-formed formulaComputer animation
18:41
Galerkin-MethodeParameter (computer programming)ApproximationResultantParameter (computer programming)Normal (geometry)Arithmetic meanOrder (biology)Hydraulic jumpLimit (category theory)OscillationObservational studyTerm (mathematics)Maxima and minimaMeasurementGalerkin-MethodePosition operatorShooting methodStandard deviationNichtlineares GleichungssystemProcess (computing)Slide ruleMultiplication signMass flow rateConcentricLine (geometry)Degrees of freedom (physics and chemistry)Spectrum (functional analysis)Basis <Mathematik>Conservation lawNumerical analysisAxiom of choiceCircleVelocityKörper <Algebra>Functional (mathematics)SpacetimeCategory of beingMassDegree (graph theory)Matching (graph theory)Connected spaceTheorySet theory2 (number)Food energyGroup representationState of matterLogical constantPoint (geometry)FrustrationInterior (topology)Group actionGreatest elementComputer animation
Transcript: English(auto-generated)
00:00
Okay. Hello, welcome back to the session part two of the CFD and GFD session today. The first speaker in the second part is Dirk Freireich-Muho. He works at the Weierstrass Institute as well, one of my colleagues. And he will talk about DG methods to reduce spurious oscillations and a special method for that. Okay.
00:28
Okay. All right. So hello and welcome everybody. I'm Dirk, as Ulrich said. And I want to present today some results of my PhD so far that I conduct under the supervision of Professor Volker Jones.
00:40
So I present today not only my work, but also joint work together with Volker. And before we start, let me just short mention that even though that I'm in the third year of my PhD, this is my first onsite participation on conference. And I'm quite happy that this is possible and want to thank especially all the organizers for making this possible.
01:06
So as Ulrich said, and also as the title suggests, I'm going to talk about a method how to reduce spurious oscillations in discontinuous Galerkin methods for convection diffusion equations. And I would say without further ado, let's right jump into this talk by explaining or by let me explaining you the terms in the title.
01:24
So let's start with the steady state convection diffusion reaction equation. So this is a basic PDE with which one can model the spread of a scalar quantity inside a flowing medium. So this can be, for instance, the charge carrier density in semiconductor devices. This can be the temperature in the air when the air is moving, or this can be, for instance, a chemical inside a fluid.
01:50
So and as I said, it's a basic model that can be applied to many different fields. I guess since most of you come from a more applied or physical point of view or
02:01
engineering point of view, you might know it even better where it can be applied in your case. And yeah, it models the spread of a scalar quantity inside a fluid medium, as I said, and it models three different mechanisms. Maybe let me explain this using this right picture here. So the first mechanism is diffusion.
02:20
So we know that the in this case, the dye in the water spreads in the water due to concentration differences. So we've got regions where we've got high concentrations and we've got regions with low concentrations. And due to this concentration differences, the dye spreads in the water. This is one mechanism. The second mechanism is convection.
02:41
So if I would stir the water, then the water itself is moving, right? And due to the movement of the water, also the dye in this picture gets transported in the domain. And yeah, this is called convection. So the movement of the underlying medium. And last but not least, but this is more like also a technical issue.
03:03
So we also can model or can have a reaction term. So if, for instance, somewhere in the domain, dye is created or deleted, this is not of particular importance for this talk, but this can be added here as well. This also fits into the framework. So this is, as I said, or this is the steady state convection diffusion reaction equation.
03:23
And in most of the practical applications, the convection is much, much stronger. So orders of magnitude stronger than the diffusion. So this is what is needed in almost all practical applications. And I want to show you what happens with our numerical algorithms in this case.
03:45
So this is the benchmark problem that I want to use in this talk also later. So like, yeah, this is the only problem that we consider and forget about all the notation on top. This is just for the sake of completeness. Let me explain this using this picture. So the domain of interest.
04:02
So it's a 2D model problem. So, yeah, it's a 2D domain. The domain is this rectangular box over here. So this is the domain without the circle. So the circle is not part of the domain. OK, so there's a hole basically. And what happens here is that maybe let's stick again to dye and water.
04:22
So the water flows, comes from the left and flows constantly. So with a constant convection to the right. So just enters the domain on the left hand side. So it enters the domain here and leaves the domain on the right hand side. When it enters the domain, it has a concentration of dye of zero.
04:40
So there's no dye in the water. No dye at all when the water enters the domain. Here at the circle, here's pure dye. That means here the dye is part of the domain and then gets transported due to the convection. Transported within the flow direction to the right hand side.
05:04
And therefore also most of the dye will be behind the circle here. So in the direction of the convection. But due to the diffusion, also some parts of the dye spread to the upper and lower boundary. And so this is how a solution looks like for a rather moderate diffusion coefficient.
05:23
So as I said, here at the circle the dye is added and then is transported within the direction of the convection. So most of the dye is behind the circle here. And a little bit spreads also to the top and bottom edge due to the diffusion.
05:42
And so here now, the diffusion is only one order of magnitude smaller compared to the convection. But if I now decrease the diffusion, you see that the solution changed. So first of all, it gets here a little bit steeper because we've got less diffusion. So less spread to the top and bottom side.
06:01
And if I do this further, until let's say 10 to the minus six. So now we are six orders of magnitude smaller. So the diffusion is six orders of magnitude smaller than the convection. We see the following. So first of all, we've got a steep slope here, which is in principle okay. But unfortunately, we've got here overshoots.
06:24
We've got unphysical values. You've got values of 1.2, which means we've got a concentration of 120%, which does not make sense at all because we cannot have more than 100% of concentration. And also, you can see here that we've got negative concentrations of minus 38%, which also does not make sense from a physical point of view.
06:44
And these over and undershoots are also often called spurious oscillations. So this unphysical values. And my talk is basically about a method how to reduce these oscillations. So this is what I'm trying to explain you. A method how to reduce these over and undershoots, these unphysical values.
07:05
So as also was stated in the title, I am using a discontinuous Galerkin method. So basically one could also use the standard continuous Galerkin method. So I hope most of you know just the standard P1 continuous Galerkin method with continuous basis function.
07:23
This can be used as well. But there are even more spurious oscillations. There are also other methods that might be non-oscillatory, but they are either computationally very expensive or they are also lowest order only. So we are interested in a computationally cheap method also of higher order.
07:42
And here comes discontinuous Galerkin into play. It is extremely flexible due to the design of the method. And I don't want to go too much into detail here, but one important property that we need later in the talk.
08:01
Let me explain this also here with this picture. So you see here basis functions. And you can see that the basis function in this triangle and the basis function in the neighboring triangle are independent of each other. And we see here also the discontinuity. So the basis function are discontinuous as the name sets.
08:21
And this can be also seen in the previous picture. You can see here the solution also, which is also discontinuous. So we haven't got a continuous solution anymore, but it can have jumps. So discontinuities along the edges of the mesh. And so this is in contrast to the continuous Galerkin method where the basis
08:41
function, also the solution, is continuous over the edges of your triangles in this case. And yeah, as I said, the most important property for my method on the next slide is that the basis functions are, let me mention this again, discontinuous.
09:01
So the basis function in this cell are completely independent of the basis function in the neighboring cell. So they are completely independent. And this can be used now. So our idea for reducing the oscillations is to use a post-processing technique, which is, and let me also emphasize this, it's computationally very cheap.
09:22
So it's a linear method that we apply after solving the standard discontinuous or with a standard discontinuous Galerkin method. And so if you want to take away just only one thing of this talk, this is your slide. So this is like the basic principle of what we're trying to do.
09:41
This works as follows. So as I said, it's a post-processing technique. So the zeroth step is just solving with standard discontinuous Galerkin method and maybe having these spurious oscillations over here. And then the post-processing looks as follows. So in the first step, we want to mark all the cells. So we want to identify the regions where the spurious oscillations occur.
10:03
So we want to mark the cells where this happened. And then after we have marked the cells, we just change the solution locally in this region by a polynomial of a lower degree and hopefully reducing or in the best case, of course, removing these spurious oscillations.
10:25
And this is now where we need this discontinuity property because, as I said, we can change the solution in this region without affecting the solution in the other regions. This wouldn't be possible with continuous Galerkin method. And this method is basically designed for discontinuous Galerkin methods.
10:43
And here is where we need this discontinuity property of the basis function also of the solution. So this is basically the idea. And Volker and I published two papers on that where we investigated different post-processing techniques, so different slope limiters.
11:02
And yeah, I think it's not like that important to go into detail with all the methods. But if you're interested, you can have a look at these papers. We investigated, as I said, several methods. And you can look up the details there if you're interested.
11:23
Or that you get at least an idea of how these methods in principle look like. I will just present you roughly two methods so that you have at least an idea how these methods could look like. And as I said, if you're interested in other, please look them up or ask me later.
11:40
And to counteract publication bias, I want to show you a method that works not so good and a method that works quite well. So let's start with one method from the literature. So it's a rather old method from the late 90s of the last century by Cockburn and Shue. And we call it linear reconstruction on triangles.
12:02
So let me again explain this using this picture. So assume that this is your grid. So we've got these four triangles here. So this is the grid. In this triangle here and also in this triangle, the solution is constant zero. So it's a constant zero solution here and here. Here in the back triangle, here at the end, the solution is constant one.
12:24
And here in the middle triangle, we've got this parabolic, this P2 function. This is the solution. So this parabola here or parabolic function is the solution in this middle triangle. So I hope you get the idea of this picture.
12:41
It's a little bit difficult to draw these things, but I hope you understand. And Lin Trierico now works as follows. So we now, we are in this situation on being on that cell and want to check if we want to mark that cell. So this is, as I said, the first step in this post-position technique. We want to mark this cell. And for this, or for Lin Trierico, we just evaluate the function here in the edge midpoint.
13:05
So, yeah, we evaluate the discrete solution just in a single point in the edge midpoint and check if this value is too large, whatever too large means. So it's basically compared to the integral values of the function in the cell and also in the neighbors, but this is more like technical details.
13:22
So we check the solution here in this point coming from the cell. So we are evaluating the solution in the cell at the midpoint and check if this value is too large. Because this hopefully indicates some overshoots in this case. And after we have marked this, we change the solution to a linear or rather affine function.
13:47
So this is why this method is called linear reconstruction on triangles because we change the solution locally to a linear function. This method basically computes three different options and then chooses one.
14:02
So this could be one, or this is one option how the solution is reconstructed. But as you can see, it's an affine function, so a P1 function. We change the solution to a P1 function. So this is one method or basically or roughly the idea of one method. The next method, we call it constant jump norm.
14:22
And yeah, this is like our method or one of our method. Here, the situation is similar. So again, we've got the constant zero solution on this triangle and this triangle. So here and here. And constant one solution here on the triangle at the back.
14:40
And here in this middle triangle, the solution is again a quadratic function. And we now want to again check if we need to mark this cell. And as you can see, we certainly want to because we've got some overshoots here. Okay, as I said, the solution jumps between these triangles.
15:04
So whenever I'm on an edge, I have let's say two values of the solution. One solution on the left triangle and one on the right triangle. But what we can of course compute is exactly the jump. So the difference of the solution in this triangle and the difference in this triangle.
15:24
So we just take the difference of these two solutions. And then basically we integrate this jump along the edge. So we're computing the integral or basically the mean value of the square of the jump. And we want to compare this again to some reference value.
15:42
And we want to check if this jump or this mean jump is too large. Then we might again have detected spurious oscillations. Okay, and after the cell is marked, it gets just replaced by a constant approximation. That is the integral mean. So basically the integral mean to preserve the mass.
16:02
Okay, so we just replace it here with a constant approximation. In contrast to the previous method, which has used a linear or affine approximation. So yeah, these are now two methods. And maybe let's see how they work in our experiments. So soon on the next slide, I'll show you the results again for the Hamka example.
16:24
But we want to somehow evaluate what is a good solution. So which solution is good? And basically we want to measure somehow these spurious oscillations. These over and under shoots. We want to somehow quantify how much are we off or how much unphysical values do we have?
16:46
How many unphysical values? And for this, we've got these two measures called OC max and OC mean. I guess you also don't have to understand this formula. Let me just explain. So OC mean, the idea of OC max.
17:00
It's called OC max because it tries to measure the maximal appearing oscillation. So again, this picture here. Let's say here, this is the largest value that the discrete solution has. And so we take this value and compare it to the maximum value of the exact solution, of the real solution.
17:22
So this is one contribution to OC max. We also do this with minimum. But basically this measures the maximal appearing over and under shoots. So this is what OC max does. But if you compare the left solution and the right solution, they have the same OC max value.
17:43
Because it just takes into account a single point. But certainly this solution looks worse compared to the right solution. Here we've got less oscillations. And for this we've got the second quantity called OC mean, which measures the mean oscillation. So it just computes the mean of all oscillations occurring in the domain.
18:05
And we've got these two measures to have an idea of maximal oscillation and mean oscillations. Now, as I said, we tried out these methods and rate these methods using our measures.
18:21
And compare this to the standard DG method without any post-processing, which I call Galerkin. And just to remind you shortly, we tried this out at the Hamke example. So again, this is the rectangular domain without the circle. So the circle is not part of the domain. The water flows from the left constantly to the right. And so this is how a reference solution would look like.
18:46
So we've got here no over and under shoots. So this is what we want to, this is like a sketch of the reference solution. Because the exact solution in this example is not known. So this is also just an approximation, but a sketch of the solution.
19:00
Okay, and let's come to the results. So as you can see, this is a snapshot from the paper. As I said, we investigated several methods. But since I just talked about the first two methods compared to the standard Galerkin method, I also want to just represent the results of these two methods. If you're interested in these, please ask me later or look them up in the paper.
19:23
And what we can see here is the two measures. So OSC mean and OSC max, compared to the number of degrees of freedom that we have in our mesh. And let's start maybe on the left hand side. So the mean oscillation. The blue line here is the standard Galerkin method.
19:43
So the standard DG method without any post-processing. If we apply Lin-Trier-Reco, so the linear reconstruction, we certainly reduce the mean oscillations. But not as much as if we use cons jump norms. So even though that it's quite an easy method,
20:03
this reduces the mean oscillation significantly compared to Galerkin and also compared to Lin-Trier-Reco. So Lin-Trier-Reco, as I said, also reduces the mean oscillation. But if we have a look at OSC max, it even increases the maximum or maximal occurring oscillations. So we're reducing the mean oscillation with Lin-Trier-Reco,
20:23
but at the price of introducing even larger maximum oscillations. We'll also see a picture on the next slide or second next slide. But yeah, so if you have a look at the maximum oscillations, Lin-Trier-Reco might not be the best choice.
20:42
But cons jump norm also here improves the situation significantly from around like 30% off. We are just now only let's say 10% off, which is well better on the one hand. But on the other hand, we're still not able to remove the oscillations completely, unfortunately.
21:00
This was P1, I forgot to mention, so this is like for piecewise affine basis function. We also have higher degrees. So for instance, P4, and here the situation looks similar. So again, in the mean, Lin-Trier-Reco and cons jump norm reduce the oscillation significantly, while cons jump norm is way better than Lin-Trier-Reco.
21:23
And compared to OSCmax, again, Lin-Trier-Reco introduces even larger maximum oscillations. And cons jump norm here really, really significantly decrease the oscillations also in the maximums norm or measure, let's say measure.
21:40
So here we've got like 60% off, and here we are only like 5% off, 4 to 5% off, which is way better. And so that you have an idea how this looks like, how the solution now looks like. Here you see the standard DG method without pause processing with these over and under shoots.
22:00
And here you see the solution of Lin-Trier-Reco. And I hope, I'm not sure if you can see this, but if we have a look here at these over shoots, we can see here that most of the over shoots are detected and also corrected. But here there are some small peaks, and this is exactly the position where we unfortunately introduce larger maximum oscillations.
22:23
So as also the measures suggested in the mean, we reduce the oscillation, but at the price of introducing even larger maximum oscillations. And the second solution here looks way better. So as I said, we are only 4% off and it detects all the, so like
22:43
the region with the spurious oscillation and change the solution to the constant value to the integral mean. And this is how the solution looks like. And compared to the standard Galerkin method, as I guess you can see, this is just way better as also our measures suggested. And that's basically it. Let me just shortly summarize.
23:03
So we saw that the standard DG method may show spurious or large spurious oscillation for convection dominated convection diffusion equations. And we saw several or at least two slope limiting techniques that are a cheap way to significantly reduce these oscillations. But on the other hand, even though we reduced them, none of them is able to remove them completely.
23:25
And in the future, we want to have a look at maybe a parameter studies because some of these limiters also depend on parameters. And yeah, there we can have a look if we can further tune these methods. And we want to also investigate what happens if we combine several of these methods, if this brings us better results maybe.
23:43
So with this, I'd like to thank you all for your attention and I'm happy to take questions. Thank you, Dirk. We have time for a few questions. Anyone?
24:05
I have a short question. This velocity field for your example, this is not constant. Because they had to go around this circle or something like that because you were. Well, how is this constructed? Like like the convection term.
24:23
So this B is a constant value of just one zero. So it's basically what happens at the circle, by the way. So but then it's not divergent for you to go somewhere around the circle. Yeah.
24:42
Yeah. Okay. So if you. Yeah. But then it's it goes there and then how it goes to be left to right only by diffusion. To the top and bottom goes only due to the diffusion.
25:04
Exactly. To the top and bottom edge. The concentration spreads only due to diffusion. Yeah. Okay. Maybe it's stationary. The flow goes through the obstacle. Time dependent then it goes from left to right and then it's a different story.
25:51
The main advantage of this continuous collecting to go to a higher order, it's easier than to as for other methods. And if you go to a higher order, what about the slope limiter in that case?
26:05
What do you mean? I like the slope limiter work like for P1, but also for P2. Even higher. You go to six, seven order. So we investigated the slope limiters until four. So of order one, two, three and four. And the results are all similar.
26:20
So const jump norm reduces the oscillations. Yeah, pretty much. And lintier reco also not that good. So it introduces larger oscillations. But all these methods work also for higher orders. Okay. Thank you. Does this answer your question? Okay. Anyone else?
26:43
Maybe I have one. Sorry. So how does this post-fix solution affect conservation properties of your scheme? And how well does one understand the origin of this oscillation? I mean, there's related to a spectral representation.
27:02
It's a Gibbs phenomenon, I suppose. So you cannot represent this continuous solution in these space functions. Like fixing it after the fact, is that the best approach? Or does it actually preserve the properties you would like to have?
27:21
Like conservation of mass and so on? This is preserved. This is guaranteed by this method. So by construction, the mass is conserved. So we just lose some order of approximation in these discontinuous regions. So this is basically what we lose, of course.
27:42
Because we just clipped the polynomial degree in these regions. But the mass is preserved. Okay. So if you have any more questions, then I think we're done with this one.
28:01
Okay. Thank you, Dirk. Once more.