Flash and JavaScript are required for this feature.
Download the video from Internet Archive.
Description: This session introduces finite volume methods, comparing to finite difference. After discussing scalar conservation laws, and shockwaves, the session introduces an example of upwinding.
Instructor: Qiqi Wang
The recording quality of this video is the best available from the source.
Session 10: Numerical Metho...
The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high quality educational resources for free. To make a donation or to view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.
QIQI WANG: OK. So today is our lecture on finite volume method. So we talked about finite difference method. And we are going to be talking about finite volume method and finite element method. So I'm going to-- there is a request for me to go over what did I do on the matrix form of the two dimensional finite difference. So I'll go over that. But before I do that, let me show you what is the difference between finite volume method and finite difference method.
Because there are two concepts, they are two conceptually different methods of describing and discretizing partial differential equations. So I have a few demos. I have a finite volume demo. I have a finite difference demo. I have also a finite element demo. So they are really three different conceptually different methods of solving PDEs.
So first I have a finite difference demo. So let me explain what is happening here. We are representing a function of x, a function of space in solving partial differential equations. And this function is going to evolve as a function of time. But in order to [INAUDIBLE] that function of x, we need to discretize it. Because a function is an infinite dimension [INAUDIBLE]. Like that. Huh?
AUDIENCE: [INAUDIBLE].
QIQI WANG: I think I-- I think there is something I can get rid of. I just don't know-- Advanced Settings. Monitor. Yeah. I don't know why. Oh. That went away? Oh OK.
[INTERPOSING VOICES]
AUDIENCE: Don't close that dialogue box.
QIQI WANG: OK. Advanced Settings. Monitor. I think I didn't do anything except for--
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah, just like that. Yeah.
AUDIENCE: [INAUDIBLE].
QIQI WANG: OK. Wow.
[LAUGHTER]
I have no idea. This is weird. Technology. OK.
So I have a function of x. I need to discretize it. Can somebody come and draw me a function? And visualize how this is kind of different, how finite difference discretizes a function. So can somebody come up and draw a function? And maybe you can describe to me before you draw the function of how finite difference discretizes the function, how will you [INAUDIBLE] that function.
How would the-- in a finite difference method, how do we represent, like, obviously a function? Let me see. These blue lines are where the grid is. [INAUDIBLE], Right. OK. So who wants to draw it for me? We just did the project. It's a time to relax a little bit. So [INAUDIBLE] the difference in [INAUDIBLE].
[INTERPOSING VOICES]
QIQI WANG: All right. Good. So when we draw a function, it looks like this. Right? [? If this ?] curve is the function [INAUDIBLE] a little bit, right? On the-- if not the [INAUDIBLE]. In finite difference, we forget about that the function is. We only remember the value of the function at these distinct points. This is finite difference, though.
Do we all have a mental picture of what this does? What finite difference does? And then when you have approximated derivative, right? For example the derivative here, we can approximate it either as the slope of this line or as the slope of this line. Depending on if we do backwards [INAUDIBLE] or [INAUDIBLE]. And in particular we may approximate the derivative here, [INAUDIBLE] taking this slope. If it's backwards [INAUDIBLE] then we are taking this slope [INAUDIBLE]. So there are some approximations [INAUDIBLE]. This is what finite difference does.
OK. So can I invite somebody to do something that is completely unpredictable? Because we are going to be looking at the same thing. But we're going to be using the same grid. I'm also going to ask somebody to draw me a function. But now we have a visual of how finite volume is going to discretize this function. OK. I'm going to type FvDemo. It's the same plot before-- yes, I'm going to ask you to draw the function of the same plot. But can somebody draw--
[INTERPOSING VOICES]
QIQI WANG: And let's see how a finite volume is going to discretize the function. Whoa. OK.
AUDIENCE: [INAUDIBLE].
QIQI WANG: OK. Well, we all see the function in the background, OK? See that? This is how finite volume [INAUDIBLE] this function if we do that. What do we do? What does finite volume do? [INAUDIBLE] the function value anywhere?
Would somebody describe to me what we did? What finite volume did? And we're taking like an average value in between the points, each two grid points. Right? In finite difference, we talked about grid points. Because we are exploring the value of the function at these grid points. In finite volume, we are talking about cells or controlled volumes, right? Each space is in two points, or controlled volume. We are discretizing the function by storing the average value of the function inside each controlled volume, or inside each cell. OK?
So in particular, if we look at this one, we will have the [INAUDIBLE]. But [INAUDIBLE] for the average difference. OK. And as we see, when we use finite volume, we need to evolve these averages. And we know the function value. In finite difference, when we look at the function value, and the function value, of course, will satisfy the PDE. And the PDE in order to solve the PDE, we need to approximate the spatial derivative.
Here, we need to evolve this average. And we are going to look at what the differential equations on this average can [INAUDIBLE], how these averages cannot satisfy the PDE. It satisfies something else. And this is what finite volume does. It is going to derive from the PDE a different set of equations than describe the evolution of these averages. And so this equation set describes the evolution of the averages.
OK. So this is a conceptual difference between finite difference and finite volume. And when we start finite element, I'm also going to-- well, let me not show you today. OK? If you're interested in it, you can just go run it at home. I'm going to give you the three scripts. If interested, you can run it at home and see what you get. [INAUDIBLE] who don't have a tablet, you can go with your mouse on the function.
So I'm going to give you the scripts later on. But right now let me first focus on and go to our last lecture folder which is 09. And I'm going to show you a little bit of what exactly we did. The two dimensional advection question.
OK. So this is the two dimensional advection equation. I think Professor Willcox ran that last class, right? We have this wavelike thing going on there. Of course, this is not as advanced as your project. Because we are looking at the constant cx and cy. We're not looking at a cx and cy responding to the [INAUDIBLE]. So this is actually a simpler script than what you guys wrote in the project.
But one thing I'm doing differently from the project is that I'm using implicit methods in solving this equation. [INAUDIBLE]. So in using the implicit time-integration method, I have to use the matrix form of the finite difference. Right? Does everybody understand why by switching from explicit to implicit I need to use the matrix form of the finite difference? You don't?
AUDIENCE: [INAUDIBLE]. Does the matrix only have a times b plus something else times [INAUDIBLE]?
QIQI WANG: Exactly. The reason I have to use the matrix form is because in an implicit time-integration scheme, the spatial derivative is operated not only on time step n but also on time step n plus 1, which is not known before I do the time step. I have a spatial differential operator operating on some unknown solution. So I cannot compute it explicitly. I have to move that curve to the left hand side of the question and solve it.
So let me make a new page and show you what this is. Page options. Insert. OK. It's more like this. [INAUDIBLE]. OK. So in an explicit scheme which would be partial u-partial t equals to, for example, a simplified version of what you guys did is this, right? This is like the simplified version of what you did. And if I do explicit time step, and if I use forward Euler, then I can do U of n plus one minus Un over delta t. That is an approximation of the time derivative.
OK. I can do matrix form. I can do non-matrix form. But whatever I do, I'm making-- I'm either evaluating distance in matrix form or not in matrix form. But I'm actually computing the differential operator operated on Un, which is known. The fact that Un is known makes it possible for me to evaluate this without [INAUDIBLE] the matrix A. Right? This is explicit.
But what if I need to do this? What if I want to use an implicit method? For example, if I need to do backward Euler, then the difference is that these spatial finite difference operator is operated on not Un but Un plus one. Well, yes. Where a is finite difference, approximation of minus partial u-partial x minus partial u-partial y. Right?
OK. So in this case, I do not know what we have for [INAUDIBLE]. Therefore, I cannot just evaluate-- I can no longer do-- most of you did, which is taking a bunch of indexing with this U and subtract from doing the indexing or doing the circular shifting of U and subtracting it from each other. I cannot do that. Because I don't have the n plus one.
What I need to do instead is to move all the Un plus one terms, is I need to move all these blue terms to the left hand side of the equation and move the green terms to the right hand side of the equation. So this is what I need to do. I need to move these to this side. So what I end up getting is identity over delta t minus A of Un plus one is equal to Un over delta t.
So I [? initially ?] was so busy [INAUDIBLE] I did not [INAUDIBLE] must have [? deflated ?] A, constructed. Right? So that is why, when I am solving-- when I'm using an implicit scheme, I need to have the matrix form of the finite difference operator.
OK so this is how I did it in the finite difference matrix. I am taking-- so I am constructing the same derivative matrix I constructed in the [INAUDIBLE] case. This is-- if you let me-- OK. Let me run this. I get a grid that looks like this. And then when I run this, getting-- If I do a spy on Ax, I'm getting the same matrix as we get in a one-dimensional case. On the y diagonal, I think we have something on the diagonal [INAUDIBLE]. Remember this? And Ay has the same structure. It is almost identical, except for the dimension are a little bit different because the y direction has a different number of grid points as the x direction.
Then I did something interesting. I did a run of Ax with Iy. So Iy is just an identity matrix. So I did a [INAUDIBLE] of Ax in Iy. To show what [? connect a product ?] does, I think I just want to do a spy of the product. This is what gets me. OK. You see, what I'm getting-- let me zoom in. Let me zoom in using my pen. I think it's a little better. So what I'm getting is a bi-diagonal matrix. But the [INAUDIBLE] I'm replacing each element in Ax with that identity matrix. So now this [INAUDIBLE] this is a series of dots. Each dot is one entry of a matrix. There is a bunch of elements over here. And there is a bunch of elements over here.
So when I'm eventually doing is, OK. If I modify something with this matrix, I'm subtracting-- so let me look at this one-- I'm subtracting the value corresponding to this column from the elements corresponding to this column. So I'm subtracting two values that are spaced apart with the spacing exactly equal to the number of x elements. Right?
So in terms of spacing, if I have-- let me insert another one. If I'm organizing U2 1, U2 2, et cetera, U 2nx, et cetera, to Uny1, et cetera, Uny2. So if I'm organizing my array like this and I'm going through-- [INAUDIBLE].
OK. And if I'm organizing my array like this, I'm organizing my array as this one goes to this one goes to this one goes to this one. So I'm storing this first in memory and then this one and then this one. And then you want to two in memory and following this. So I am storing one column the other column then another column, et cetera. Then if my two elements are spaced apart with ny elements in memory, then it is the exact same, the adjacent elements in the x direction. That is how my elements are stored in memory.
So which means, if I'm taking derivative in the x direction, my spacing should be ny. And what happens if I take the derivative in the y direction? What if I do kron of Ix with Ay? Oh. OK. Yes, I mean this.
OK. When I go to something like this, so if you look at each block-- so I'm getting a periodic structure. And each block [INAUDIBLE] is exactly the same as the one dimensional derivative matrix. And once we zoom in over here-- oh, you can't see anymore. But these are individual dots. So when we take the derivative in the y direction, we are taking [INAUDIBLE] adjacent values in the array, which is the adjacent value actually in the y direction. So adjacent values in the x direction are actually spaced over by ny in the memory. So these values in the y-direction are spaced with only one in memory. That is why we need to use the kron in the matrix form.
OK. Questions on that? OK. Now let's go to our finite volume analysis. OK. In finding the volume, we are ultimately going to two dimensional too. Today, we are just going to focus on one dimension finite volume. And again, in finite volume-- OK. Let me go back to our last [INAUDIBLE]. In finite volume, we draw a function. And obviously in the function we are only storing the averages inside each cell.
All right. So let's look at why that is, what kind of equations can finite volume solve? And why does it solve this kind of equations well? Finite volume only works if you have a conservation law. Finite volume is not for all differential equations. It works only when you have something that can be formulated in terms of conservation laws.
And what does conservation law mean? Conservation laws mean the equations that are in this form. OK. You can have a [INAUDIBLE]. But like you can have a [INAUDIBLE] here. All right. So this is conservation laws. It's some kind of a conservative quantity. Can be mass, momentum energy, some molecules, number of [INAUDIBLE], number of people. Something that can be described as conservative.
This F is the flux. The flux is the number of quantities u that goes through a boundary per unit amount of time. OK. That's what flux means. All right. And the specific force, the force is how much quantity, again described in u, is generated, comes out from nowhere, per unit amount of time.
This is a case, for example, if you have some kind of a-- if u is describing the amount of molecules and F is the number of molecules generated in an amount of time if you have a chemical reaction or things like that. All right. And if u is describing momentum, then F is the flux of momentum through the boundary. And S would be something like the [INAUDIBLE] of gravity. Something like that. Right? Is it clear what the conservation law means?
The conservation law really makes even more sense if you try to integrate it. OK? So this is the differential form. Now if you try to integrate the conservation law, it even makes more sense. And by integrate, I don't mean integrating time. I mean integrating space. So let's integrate over any control volume. And control volume is described as an interval from L to R. L is the left hand of the interval. R is the right hand of the interval. So the integral of left hand side has to be equal to the integral of the right hand side of S, both are dx.
So let's do some manipulation to this integral form. So first of all, the integrative time derivative in space. Because time and space are two independent variables, the integral in space [INAUDIBLE] with the derivative in time. So the integral in space of a time derivative is equal to the time derivative of the integral in space. The first term can be written as this. And this is the time derivative of [INAUDIBLE] using that control volume. If U is density, then this is the amount of mass inside the control volume. If U is the density of [INAUDIBLE] like molecules, density is the total number of molecules in that control volume.
OK. The second term. The second term is not a timed derivative. It is a spatial derivative. Now what is the spatial integral of the spatial derivative? The function itself. Exactly. So it is F of U at R minus F of U at L. All right? And let's actually move this to the right hand side. Because it makes more sense. Let's move it to the right hand side and put a minus sign here, put a plus sign here.
OK. So this is the result of the spatial derivative term. And the [INAUDIBLE] that's [INAUDIBLE]. OK. So this is what the integral is going to give us. Now let's try to interpret that integral. The time derivative of the total quantity within the control volume is equal to the flux at the left minus the flux at the right. This is what goes into the control volume. This is what goes on the control volume.
So makes sense? The rate of change of the amount of stuff inside the control volume is equal to the rate of change for the stuff going in minus the rate of stuff coming out plus the rate of stuff generated inside the control volume. So this is the integral form. This is the integral form of the conservation law.
Now why is it important? And how does it relate to finite volume method? Because this is the equation we are solving in finite volume method. This is the equation we're solving in finite difference method. Right? We stick with this kind of equation. Because in finite difference method, the U, we are storing the value of the solution U at individual points. Right? And that is why we need to approximate the distance, approximate the spatial derivatives [? and it ?] kind of a finite difference method. And then we can do [INAUDIBLE] for U [INAUDIBLE].
With finite volume method, we are not interested in the solution at individual grid points. Instead, we are storing the average of the solutions. And the average of the solutions, what is a mathematical description of an average of a solution in a control volume? What is the average when we look at the amount of stuff [INAUDIBLE]?
Yeah. It's the amount of stuff in space. If I ask you to write a mathematical formula of the value of the average as a function of-- in terms of the function within these intervals, what is it?
AUDIENCE: [INAUDIBLE].
QIQI WANG: Exactly. The average would be equal to the integral of this function from 0.6 is 0.7 divided by the width of the integral, which is 0.1. That is what the average means. And now we have derived the integral form of the differential equation. And it describes the evolution of the integral. Do we have the evolution of the average?
Of course. We just divide this thing by R minus L, right? We just divide the whole thing by R minus L. Then this whole term actually-- let me include the time derivative-- this is the d dt of the U average. Right? This whole term would be the time derivative of the average value of the function between L and R. The time derivative of the average value between L and R is just equal to the amount of stuff going in minus the amount of stuff going out divided by L over R. Plus the average amount of stuff generated in the control volume. Right?
So the reason we derive the integral form from the differential form is because this integral form is what is satisfied by the finite volume method, by the volume average. So we use this form for the finite volume method. Now this form has distinct advantage over this form. Why? Why is sometimes this form valid while this form is invalid? Can you guys think of a case?
AUDIENCE: [INAUDIBLE].
QIQI WANG: Yeah. Of course, this has derivatives [INAUDIBLE]. This has a spatial derivative. This does not have a spatial derivative. Right? So in what cases can you not have a derivative?
AUDIENCE: [INAUDIBLE].
QIQI WANG: Yeah. We can have singularities. What kind of singularities have we seen that have no derivative?
AUDIENCE: [INAUDIBLE].
QIQI WANG: Huh.
AUDIENCE: [INAUDIBLE].
QIQI WANG: [INAUDIBLE]. Shocks. Shocks. Right? We saw-- when you say shocks in solving this kind of differential equation. Once we have a shock wave, you have a discontinuity over x. Do you accept you can take the derivative of a shock wave? No. Actually, I'm going to show you. If it applies to you to define a different method to solve a differential equation with a shock wave, sometimes you can still [INAUDIBLE]. And you only get correct solution when you use finite volume method, because finite volume method doesn't care if it has shock waves or not. It solves the integral form of the equation, which is immune to shock waves, and not even taking spatial derivative. As long as I can integrate, I can solve the equation. And the solutions with shock waves have their integral.
All right. So before we go into finite volume method and look at its solution, let's review how the solutions behave. And I think we can do that review best if we just ran the Burgers equation phase again. All right. So Burgers' equation. I'm going to draw the initial condition this time. I haven't drawn it. So I'm just going to draw something that goes up and goes down and goes up again. OK. And let's look at how the solution evolves.
OK. So this is a typical conservation law. [INAUDIBLE] is one of the simplest conservation laws. And remember, the solution moves with a speed proportional to the value itself. So if I have a high value and a positive value, I move toward the right. If I have a negative value, I move towards the left. Right? This is even more dramatic than the case where everybody's [INAUDIBLE].
This is the solution we are [INAUDIBLE]. So if I have something moving towards the right and I have something moving towards the left, and you can very clearly see a shock wave forming over here. And we can know what the solution is really intuitively because we know at each point, the solution moves proportional to the value, which if you look at the plot with [? accuracy ?] in space while at the same time, then the characeristic lines-- which are straight lines in space and time-- the flow-- or I should say the inverse flow because their speed is delta x over delta y. So the delta in the x axis over the delta in the y axis.
AUDIENCE: [INAUDIBLE].
QIQI WANG: [INAUDIBLE] is always delta x over delta y. In this case, delta x is [INAUDIBLE] delta y. Delta t [INAUDIBLE]. Delta x over delta t is kind of the inverse. Does that make sense? I mean, a speed equal zero is like a vertical line. A vertical line has no delta x but lots of delta t. So it's low speed. And the line that is very shallow is the high speed line because it has a lot delta x over small delta t. So we see these are high speed. And these near vertical lines are low speed. And these lines sloping towards the left are negative speed. They are moving towards the left. Right?
OK. The solution behaves like this. The characteristics are proportional to the solution itself. The redder it is, the larger the position is. And the more quantity, the larger the speed is. The bluer the solution is, the more negative the solution is. And the solution moves backwards at a negative velocity. This is how the solution behaves in a Burgers equation.
And we also have this shock wave that moves at an average velocity of [? 6.54 ?]. Right? So we this shock wave moves kind of slowly towards the right. Because if you look at the average left value and right value, it's slightly positive. So the shock waves moves slightly towards the right.
OK. That's how the solution for the Burgers' equation behaves. And we need to have a method to solve the Burgers equation. And we're interested not only in the smooth part of the solution, but also in the discontinuous part of the solution.
OK let me just do another case and ask you how the solution is going to behave. All right. So if I have a function u of x for the Burgers equation, the Burgers equation is partial u-partial t plus U times partial-partial x equals zero. Or in a conservative form, it's partial u-partial t plus partial u over partial x and half of u squared equal to zero. Right? So this is-- this form is [INAUDIBLE] u squared out of the spatial derivative.
After a while, can anybody tell me how the solution is going to look? Yeah?
AUDIENCE: Is this periodic?
QIQI WANG: Let's say they are periodic,
AUDIENCE: [INAUDIBLE].
QIQI WANG: Yeah.
AUDIENCE: The shock waves are formed and [INAUDIBLE].
QIQI WANG: Yes. Exactly. So this part moves at a positive speed. This part moves at a negative speed. So you can really see the solution later on should become something like this. And later on, the shock wave should form. I actually don't know if it's moving towards the left or moving towards the right. But the shock wave is going to form over here.
AUDIENCE: [INAUDIBLE].
QIQI WANG: Uh huh.
AUDIENCE: [INAUDIBLE].
QIQI WANG: Yeah. Good point. This part is a little bit [? difficult ?]. So the positive part is still dominant and the shock wave is going to move towards the right. Good point. OK. So this is how the equation behaves. And let's see how do we solve it using a finite volume method.
OK. So let me write down the equation again. Partial u-partial t plus partial-partial x of half of u squared equal to zero. Right? Or in the more general form, I should be writing partial u-partial t plus partial F of u partial x equal to zero. And to specialize this general conservation law into the Burgers equation, we just need to set F of u equal to half of u squared.
OK. And again, you find the volume. We are solving for the derivative. OK. Now we need to, we are solving the-- we are tracking the evolution of this Uk average. And the Uk average is defined as one over delta x. So delta x in this case is the size of each control volume. OK? Control volume. OK. Times the integral from the left of the k-th control volume to the right of the k-th control volume. OK. And the size of the control volume is of course equal to Rk minus Lk.
OK. The integral is the u of x and t dx. This is my cell average. And in general, this is delta xk, right? Because each control volume can be different. And here we are just going to talk about the case with a uniform, uniform mesh. And the uniform mesh responds to delta Xk or r equal to delta x. And the L of k is equal to k minus one times delta x. And the r of k is equal to k times delta x.
All right. So the L of 1, the left boundary of the first [INAUDIBLE] is at zero. The right boundary of the k-th control volume is always the same as the left boundary of the next control volume. That has to be satisfied for finite volume method. You have to divide the entire space, the entire domain into non-overlapping into non-overlapping control volumes and leave no empty space in between. You have to basically partition the entire [? combination ?] domain into non-overlapping control volumes. And you can't leave anything open.
OK. So now let's take a look at what is the time derivative of this average? To figure out the time derivative of this average, we have to use the integral form of the conservation law. It is equal to one over delta k times d/dt of Lk Rk Uxk dx. And by applying the integral form of the conversion law, we get the flux at the right hand side minus the flux at the left hand side of the control volume divided by-- oh, we don't need to divide anything because we are looking at the derivative of the integral. Right? And plus the integral L to R of a [INAUDIBLE] which is zero in the Burgers equation case. Because in the Burgers equation, we don't have any source then. The source is zero.
OK. Now we have an evolution equation for the cell average. And up to here, we've made no approximations. Right? In finite difference, we are approximating the spatial derivative. But here, we don't need to make any approximations in getting to here. And the flux in this case is half of u squared at the cell boundaries.
Now this is the point we have to start making approximations. Because unlike in finite difference, we have the value of the function x at [? these ?] points. But here we don't. We don't have the solution at either Lk or Rk. Right? We don't have the solution at the control volume boundaries. We only have the cell average solution values.
So we have to approximate, we have to do the finite volume approximation. I have to approximate the flux at rk. I have to approximate it as F of a function. The most simple cases u bar of x Uk minus one and u bar of k. F at-- sorry, this has to be at Lk because I'm on the left. Rk has to be approximated as a flux of Uk and Uk plus one.
So let me draw what is happening here. So I have a bunch of cells. This is, let's say, k minus 1. This is k. k plus one. k plus two. Right? OK.
So this point is the left boundary of k minus one. Here is the right boundary of k minus one. But it is also corresponding to the left boundary of k. This point is the right boundary of k but also corresponds to the left boundary of k plus one. This point is the right boundary of k plus one but is also the left boundary of k plus two. This point is the right boundary of k plus two but also the left boundary of k plus 3.
I need to attach a flux value at each of these boundaries in order to solve this equation. For this equation, I need to not only resolve the cell average in each cell, I need to have the flux at these boundary points. But I don't have the solution at these boundary points. All I have is the cell average inside of the cells.
How do I approximate the flux at these boundaries? It makes no sense to approximate the flux at this point using the cell averages around that cell boundary. At the left hand side of [INAUDIBLE], and this is [? like the cell interface ?]. On the left hand side interface is [INAUDIBLE] right of the interface [INAUDIBLE]. So I need to approximate the flux using the values.
OK. The easiest way to approximate it is what we call first order upwinds. And upwind means the reverse of the direction towards which the solution goes. OK. In this situation, can somebody tell me what is the upwind direction? Yeah?
AUDIENCE: [INAUDIBLE].
QIQI WANG: This direction is the upwind direction? Always? It depends. Right. The upwind direction in the Burgers equation actually means almost all [INAUDIBLE] conservation laws ends on the solution. In this part, the solution was here and now moved here. The solution is [INAUDIBLE] towards the right. It's kind of blown, like blown by the wind, and that comes from the left.
So what is the optimum direction for this part in the red region? It's towards the left, right? What is the optimum direction over here in the blue region? Towards the right, right? The wind is coming from the right. So the wind is going this way, if this were the wind. So upwind direction is looking where the wind is coming from. And in the Burgers equation, the upwind direction is towards minus x if u is positive. It is plus x if u is negative. This is upwind.
It makes a lot of sense to do upwind because that's where the solution comes from, right? So you want to-- I'm going to talk about upwinding, why upwinding makes sense really later when we talk about stability. OK. So upwinding schemes. At this point, I just need to tell you it gives stability. And in the upwinding scheme, I'm just going to say F of k-- so F at the interface-- so I'm going to define F of k plus 1/2, defined as F at the right boundary of the k value, which is also the same as the left boundary of the k plus one cell. It is the flux at the interface between the k and k plus one cell. Right? I call it k plus 1/2, because it's kind of a halfway between the k and k plus 1 cell. Right?
And the first order upwind scheme is to compute F of k plus 1/2 is equal to 2 cases. It is either F of u bar at k, which in the Burgers equation is just u bar of k squared over 2. If the minus direction is the upwind direction, I'm here again at-- I'm going to draw the same cells again. So k minus 1, k, plus one.
OK. I am looking at k plus 1/2, which is the interface between k and k plus 1. I am going to use the average method of k to approximate here. If k is the upwind direction of the cell, which means if the wind travels towards the right, if locally u is greater than zero. And here locally-- I can say locally u is greater than zero if the average between the two cells is greater than zero. And else, I'm going to approximate it using the upwind direction, which is now towards the right. It is this.
Does it make any sense? This is what I'm going to write in the code. OK. This is what I'm going to implement in my code. That is, in the [INAUDIBLE] code [INAUDIBLE]. The F at interface is either the flux at the left value or the right value where [? it is ?]. Any questions on this? No? OK.
Let's do it then. What I'm going to do is I'm going to pull up a skeleton. And I'm going to ask you to fill in the blanks. So I'm going to go to the shared border. Open your computers if you haven't. I'm going to create an order 2014. Today is the 10th, right?
OK. I'm going to start from scratch. I'm going to make a function just called the du dt Burgers. OK. And like our-- so we're assuming PDEs. But as long as we can convert the PDEs into a set of ODEs, we know how to integrate it at this point. I'm sure you guys all can.
So what is it going to be? The input is the u bars. And the output is d u bar dt. And what I need to do is I need to compute from the u bars the du bar dt. So what do I do first? I'm first going to do F bar is equal to u bar squared over two. What am I doing here?
AUDIENCE: [INAUDIBLE].
QIQI WANG: Computing the flux. This is the only thing I need to change if I switch from the Burgers equation to some other equation. OK? This is really defining the flux. If I have a flux that is equal to another form, this is what I need to change. OK. And here, I'm going to have you fill in the blanks. Compute the flux at cell interfaces. So I need to compute F. Let me call this F of [INAUDIBLE] k. So this is F at k plus 1/2 equal to blah, blah, blah. This is what you need to do. OK?
And then what I'm going to write down is this is the finite volume method. This is the same for all finite volume methods. What is it? It is the d u bar dt equal to-- according to our formula-- du bar dt is equal to one over delta x times the F, the flux, at the right boundary of k, which is k plus 1/2, minus the F at the left boundary of k, which is F of k minus 1/2. Right?
So this is F k plus 1/2 minus F k minus 1/2 divided by delta x. And I need to compute F of k minus 1/2. It's basically a [INAUDIBLE] of F k plus 1/2. Let's assume they are periodic boundaries for now. So how should I compute F of k minus 1/2 from F of k plus 1/2?
I'm going to use circshift, right? And which direction should I shift it? So I already have F of k-- I already have F of k plus 1/2 over here. I need to have-- so if I already have F of k up here, I need to shift it towards the right to get-- so that's the same point I have, F of k minus 1/2.
So we assume-- let's assume it's a row vector, so that when I shift it, I shift it to 0 and 1. OK. And let's assume my boundary goes from 0 to 1. OK. So if I have a boundary that goes from 0 to 1, my delta x is equal to the length of the domain divided by the length of u bar, how many control volumes I have. So that's my delta x.
The size of each control volume? I have F of k plus 1/2, F of k minus 1/2. Then I get the du bar dt. OK. Go ahead and just copy what I wrote here. Just rename it to whatever. Just rename it to dudtBurgers underscore your name. You can do it in the same directory. And tell me what you get.
All right. Anybody have anything I did wrong? I just wrote a script. So it's writing [INAUDIBLE]. So this is basically making [INAUDIBLE]. OK. And using all the input [INAUDIBLE]. So [INAUDIBLE]. All right? So I'm [INAUDIBLE].
AUDIENCE: [INAUDIBLE].
QIQI WANG: OK. You can try. All right.
AUDIENCE: [INAUDIBLE].
QIQI WANG: OK.
AUDIENCE: [INAUDIBLE]. The function in--
QIQI WANG: Oh. Oh. OK. Oops. I added just--
AUDIENCE: [INAUDIBLE].
QIQI WANG: OK. So this is-- let's see. [INAUDIBLE]. So I goes from one to the length of u bar minus 1. If u bar is [INAUDIBLE]. You also put I here.
AUDIENCE: [INAUDIBLE].
QIQI WANG: OK.
AUDIENCE: It's be a good idea.
QIQI WANG: Yeah. I think I need to close this so that I don't get a [INAUDIBLE]. OK. It's updated. OK. Let's run the script again. The driver. Run. Ah. OK. This is--
[INTERPOSING VOICES]
QIQI WANG: OK. Update it again.
AUDIENCE: [INAUDIBLE].
[INTERPOSING VOICES]
QIQI WANG: Oh, you must-- OK. So that's my problem. I need to change it so that I [? deleted ?] a column vector. Let me [? delete ?] a column vector and [INAUDIBLE] a column vector. OK. Now let's see. Ah. You're returning a vector of length 63. But the length I should get is 64.
AUDIENCE: Oh, right. So it's returning [INAUDIBLE].
QIQI WANG: Ah. We are assuming the domain is periodic. So at the last control volume, if we need to take towards the right, we should be taking the first value. Right? OK. Let me see, I see to [? copy ?] the network, oh.
AUDIENCE: [INAUDIBLE].
QIQI WANG: Into the [INAUDIBLE] [? space ?]. Let me see. Hm. OK. For the last cell, here I think I [INAUDIBLE] from the [INAUDIBLE], I think that will end after-- so let me see. For r equals from 1 to 10, i n. Yet, if I do i again, it is still 10. It is not 11.
AUDIENCE: [INAUDIBLE].
QIQI WANG: Yes.
AUDIENCE: [INAUDIBLE].
QIQI WANG: OK. [INAUDIBLE] update. [INAUDIBLE] update. OK. Let me do this again. Oh. I need a minus sign here. Right?
[LAUGHTER]
Is that right? OK. I need a minus sign in my-- OK. I need a minus sign in my script because--
AUDIENCE: Yeah.
QIQI WANG: Right? Is that right? Let me see. So it is equal to the right value minus the left one, which I may be--
AUDIENCE: Yes, it would be--
QIQI WANG: Should have been the left value minus the right one. Sorry. I have been deriving it incorrectly. So yes. So for the time derivative of the average is equal to what comes in from the left minus what goes out from the right. So I need a-- yes, I need a minus. So don't change your code. I just put a minus sign here.
All right. Now it looks like it's behaving a little bit weirder at the end. Otherwise it seems to be working.
AUDIENCE: [INAUDIBLE].
QIQI WANG: Let me-- OK Right. So we do see a shock wave developing, right? So let me change the delta t to be a little smaller. And then run it again. Ooh. Right? That's the correct behavior of the Burgers equation. Right? OK. Very nice.
So who else-- I think the-- let's see. Yeah. OK. Let me run it again. All right. Great. That's great. And let me see if I can break it without betraying initial conditions. And let me take it as x is zero equal to draw periodic function. Do I have that? Oh no. I don't have that in here. So I'm just going to copy it.
AUDIENCE: [INAUDIBLE].
QIQI WANG: Draw periodic function.
AUDIENCE: [INAUDIBLE].
[INTERPOSING VOICES]
[INTERPOSING VOICES]
QIQI WANG: Let me do that. OK. Let me do something. Anyway, OK. So I see a lot of you are doing good work. Anybody else? I can try your solution. OK. Ah.
[LAUGHTER]
OK. I think I'm going over time. I had a lot of fun trying your scripts. But let's meet on Wednesday and talk about more of the finite volume method. So if you want to stay, I'm going to try a little bit more of your scripts. [INAUDIBLE]. Oops. OK.
AUDIENCE: [INAUDIBLE].
QIQI WANG: Oops. All right. OK . Anyway so once you have a script, you can always try my driver and see how they work. All right. I'll see you on Wednesday. Yes. If you didn't send me a project here, please put it over here. And the homework [INAUDIBLE] with the homework.
