Mod - 01 lec - 36 foundation of scientific computing - 36

Lecture 36, and we continue our discussion on properties of filters with the help of 1D convection equation. We pick up a basically unstable Euler-CD 2 scheme and we design a transfer function in such a way, that we can modify the numerical amplification factor G with the help of this and solve this problem in a perfectly neutrally stable manner. However, in doing all these, we have already seen that basically, most of the discrete computational method introduced spurious upstream propagating waves, and this is extreme form of a dispersion error, and as we have talked about filters can alter the dissipation and dispersion properties. So, we need to relook at q waves in the context of filters; but, before we do that, we take a look the g contours with and without filter; and we pick up some examples of rotary oscillation of a cylinder in a uniform flow and we find out that, if we perform a filtering in a particular direction, then this vertical structures are actually elongate stretches in that particular direction; and of course, filter does attenuate the solution, so, we cannot avoid that. This can be somewhat mitigated by reducing the frequency, increasing the frequency of filtering, as well as considering the directionality of the filters in terms of azimuthal verses radial filter for this rotary oscillation problem. And we notice, that azimuthal filters applied globally, gives rise to serious issue. Now, in the following, we discuss about a new class of filters that we have designed, where the filters are themselves upwinded in the interior. Once again, this is done using Pade's schemes, and we workout the transfer function real and imaginary part and identify what is alteration of the dispersion properties; and, please note that, in this particular effort, we are trying to alter the dispersion properties in such a way that we can remove q x; so, the motivation for this upwind filter design is to remove q waves by a, is clever design of the filter itself. We also note that filters can be very interestingly used for large a d simulation or detached a d simulation, used in fluid dynamics; because, essentially, filters do band limit the solution and that is exactly what is done in LES. But LES is far too complicated a process, where you actually filter the governing equation and that brings in newer stress terms which have to be modeled; so, that could be completely removed by the present approach of using this implicit filters as a post processing operation. This should conclude our discussion on filters and we can then begin our discussion on Introduction to Finite Element Methods, which is basically, identifying the nature of global verses local method, and what are weighted residual methods? That is what we will briefly introduce and will conclude this lecture with. I would still a discuss little more about filters because the at the one of the assignment, and I think this is one of the most beautiful tools that has emerged over last 15 years, to do many things which were otherwise thought impossible. So, let us say if you are trying to solve a problem, an example, a convection dominated flow, and you take this model equation once again; and we know it is property that does not disperse, it does not attenuate; and you apply an algorithm, let's say, which we know is problematic, namely the central scheme on space and forward in time and going through this motion of defining functions with respect to the Fourier Laplace Transform; we define the numerical amplification factor taken with the help of the CFL number. This is the algorithm that you actually use, and when you actually use the the spectral representation of a in d, you get this; this is what we have noted; that it is an unstable scheme. So, what do you? Do people have probably thought earlier that it is a no-win situation? But we will see today itself, that by applying additional explicit filter, we can work out a stable algorithm and may be a neutrally stable algorithm; that is the one of the ideas. So, basically, if you have a solution obtained at an advanced time level, you want to filter it indicated by this cap quantity. Then, we will define the amplification factor at the end of numerical integration; and filtering is given by what you obtain after filtering, divided by the solution that you had before time integration, right? So, this we split it into two parts; and you can notice that we have artificially introduced this factor u of k hat t n plus 1, so, this is basically the output of your numerical integration and that is what we have called as the transfer function of the filter, right? So, this is what the filter does. It takes here, this denominator and produces the numerator, so that is the transfer function. And then, this is the basic numerical method by which you have reached t n plus 1 from t n, so, this is the usual language with which we speak. When we talk about explicit filters, we need to look at e j of a, so, show you a result which is somewhat different. Then, what you are going to do in your assignment? Your assignment talks about propagation of a packet, but here, we have looked at even a simpler problem taken a domain and we say that we have periodic problem, so, the periodic problem is something like this, you know? A sinusoidal wave propagating, entering through the left and exiting through the right, so that you can use this same algorithm CD2 Euler and we know it is going to be unstable; but what we do is, we apply a second order periodic filter. And for this combination of k h equal to 0.2 pi, if we choose delta t in such a way that n c is 0.1 for this, value of c equal to 0.1. This choice of filter coefficient of 0.48 would make your method neutrally stable, so that, if I identify one period of the wave by this solid red line here, that after a sufficient number of steps when you reach at t equal to 30, it is identified here at the exit plane; and you see, there is absolutely no attenuation here, so, this should give you hope in thinking that filter is a technique by which you can take an unstable method and make it work. So, this is one example that you are seeing here. Now, there is something that we already have seen, and that is about this dashed region here, in this figure. And, in this figure, where we are showing the group velocity, normalized contour for this type of combination; so, this is a fourth stage Runge-Kutta method, used with that compact scheme that we have talked about and this is CD2 Scheme with RK4; and as you see, that replacing Euler by RK4 gives you a neutrally stable region which is shown by this horizontally hatched region. But you can see that region exist over a very small value of k h; if you go little on this side, what will happen is, you will get a damped region, right? That may not add to accuracy of the solution, that is what we talked about. Now, when it comes to these spurious upstream propagating waves which we have called the q waves, this situation is worse when you look at CD2 Scheme as compared to the compact scheme. And this is what we have seen earlier some animations also, what this entails. Today, we are basically going to talk about that what filter can do to alleviate similar such problem. So, if I look at OUCS3 RK4 scheme on the right hand side, you can see that the g contours are shown here in the KHNC plane and you see only this region to the extreme left where you have neutral stability, right? Now, suppose I apply a sixth order filter with a value of alpha as 0.45, then this g counters actually change like this. What is interesting for us to notice is that, of course this neutral region degrades, you get a neutral region, now, which is sandwiched here only; whereas, this 0.99 line that you you could see here, actually terminates up to here; and what is important for us to realize is that, for this higher value of k h, the filter actually attenuates it; significantly, the filter affect that. You would see mostly on this high k h range, even when your n c is small here, we had a neutrally stable region, right? This part, but here, it is all gone; the g has become close to zero, right? So, this is something we have to realize why filter is sort of welcomed so often. Because, most of the numerical problem originates from high wave number region, and if there are problems, applications of filter actually removes it quite drastically. However, you may lose accuracy of the solution; that is another issue that we have to keep in mind. Now, show you an example of a problem; this appears very simple. What you have is basically a cylinder, and this cylinder is performing some kind of a oscillation; it is going back and forth. So, that is what we call by rotary oscillation; it is not just simple rotation, it is doing a rotary oscillation. So, what happens? Half the cycle part of the flow on one surface will be will be opposing the convection direction. So, suppose, let us say, the flow is going from left to right here. Then, when it is going anticlockwise, then the top surface you will have an opposing motion, and that opposing motion actually enhances the chance of flow separation; and in another half of the cycle, you will see a similar thing happening on the lower surface. And when you keep doing this, this cylinder keeps performing this rotary oscillation; you will see a very neat train of vertices come out. These vertices are much much stronger than what you are probably told, if you have taken a course on fluid mechanics; if you keep a stationary cylinder. So, these are extremely strong vertices and capturing this flow is not so easy. I will not go about CFD of that path, but what I am trying to show you here is basically, look at this top figure where we have shown the lift coefficient, the upward force that the cylinder experiences as a function of time. And, if you look at the three curves, that those have been shown here; the solid line is the one that is obtained without using any filter and it seems to go on, go on. Then, at t equal 47, it actually blows up. So, you see what happens is, sometimes you have to really enlarge your window observation to make sure that the things are alright. Suppose, I would have stopped at 40 and say, oh! I could get a good solution and I could present such results. But here, it is clearly seen that, at t equal to 47, all of a sudden the c l nose dives and the solution blows up. Why does... so that is given in this vorticity contour plot. What happens is, this is one of the vortex that has been shared in recent times; but prior to that, there was another vortex that was shared from the lower side, which goes and convex downstream; and there are certain parameters of this problem, is the amplitude of oscillation and the frequency of oscillation. Those are given by this a and f f, so please do not bother too much about exactly, about those parameters, but let me tell you this; flow is at a very low Reynolds number; it is as low as 150, and computing it is one of the toughest challenge that we have ourselves encountered, so, it just simply blows up, right? Now, what we could do is, we could then try to use filter and see what happens, and this is what has happened. Here, we have applied a fourth order filter with a alpha value very close to 0.4, that is what is given here, with this square symbol here, with this dash line, with this solid symbol here; that is, a fourth order filter and goes to the value of alpha is 0.495, and what happens if you do that? The solution goes little further, then again that also blows up. So, this is something that we should keep in mind, that in an actual flow, computations it is not like 1D wave equation that you are solving in your assignment; everything is known, you can work it out. For a real problem, you would not know a priory, what is the order, what is the filter coefficient that you have to choose. So, here is a choice of a fourth order filter with a alpha equal to 0.495, does not suffice. Why that? You can see here, the vortex which was causing that solution breakdown, this vortex keeps on growing un-physically continues to happen; so, here is a case, actually I think it is a sixth order filter; it is a sixth order filter, the filtering direction. This is another issue that I would like to bring to your attention; it is a two dimensional flow problem, right? So, we have a radial direction and we have a azimuthal direction, and the filters that we have talked about, they are all one dimensional filter, so, what we could do is, we could of course, apply it in the azimuthal direction, then we could also apply it in the radial direction. What happens is, the here you are seeing some results, azimuthal filters, but a sixth order filter with a value of 0.495 for alpha does not suffice. At the same time, if you apply a second order filter and take the same value of alpha 0.495, what you notice is that, you can compute indefinitely, how good are the result. Now, that is a very very a valid question. Because, the offending vortex which is causing the numerical breakdown here, in these two cases, attenuate to this; this is a much more weaker vortex, but it also has a very funny attribute to the solution. The vortex has been stretched in the theta direction, because that is the direction along which we are doing filtering. So, if I keep doing azimuthal filter, I may get around the instability problem. But then, the vertigo structures would take some unphysical attributes; like here, you can see this vortex has been detached already and this vortex has also been stretched in the theta direction. So, these are some of the issues that one needs to worry about and the another thing that happens is, there is no need for you to do the filtering at every time step, so, you can do it infrequently, so that this frequency of filtering itself is another parameter at your disposal, that you could make use of. So, what is being seen here is, the solution that we are showing for another case, this is a milder case, so, here we could compute without any problem. You see, the amplitude has been reduced from 5 to 2 and the frequency has also changed to 4 and the unfiltered solution is shown here. Now, if you start applying the filter in this frame, you are seeing that a second order filter with alpha equal to 0.4 has been applied, and the interval of application of filter is about every thousand steps you are doing it now. What happens what happens is that, if I do infrequent filtering, then the basic numerical method has an amplification factor which is given by the first factor on the right hand side g j, and if I would have done it at every step, I would have multiplied by t j. But since I am doing after n times step, so, the resultant resultant amplification factor would be g j of the original method times the transfer function raise to the power 1 over n, so you could actually reduce the intensity of filtering by doing infrequent filtering, right? That is what is being suggested here, and you could apply this procedure for any order filter in any direction and that is what is shown in this last three frames. Suppose I try to solve the problem using a second order filter and take the value of alpha same as 0.4, and if I keep doing the filtering at every time level, this is the solution that you get. And I told you that if you do an azimuthal filter, you see this unphysical attribute, the vortices are stretched in the theta direction and that is what you are getting. So, you are getting so called stable solution, but not necessarily a good correct solution. If you keep filtering every ten steps, that is, this frame that you are seeing here; well, you can see little more features coming out, those unphysical stretching in the theta direction has been prevented somewhat, but it is still not completely gone. Because, you see, these vertices which are all there, they are all absent in this two cases; so, filtering actually removes signal also, you have to keep that in mind; and then of course, when you increase the filtering frequency further, you do it every hundred steps; well, you start recovering back some of those lost signal, but it still, you see there is a missing structures; whereas, if you increase it, increase the frequency ten times, that is, your filtering every thousand time-step; well, you can see more or less the structures there, but they are much more weaker because of the factor that we are applying a filter. Please do understand that alpha equal to 0.4, is a quite strong filter; it is its not a very small quantum of filtering that one uses. Now, having talked about the frequency of filtering, here are some results where we are talking about the direction of filtering. Once again on top, you have the unfiltered solution for the case; just now, we have seen at a different time, we are looking at; and now, instead of applying a second order filter with alpha equal to 0.4, which actually makes the solution totally unphysical here, we decide to take a high order filter a sixth order filter and also take alpha larger at 0.495. And what we need to look at, is what we get if we do a filtering in different direction. For example, this one here, we have used a azimuthal filter, so, basically, the filtering is done only in the theta direction; and you can see, there is some effect some effects of the solutions; there is a wrong speed of propagation of this vortical structure, this strength also changes the speed of convection of the vortical structure changes. However, near field structure are modified less, but they are indeed modified, but somewhat less ok. Now, if you switch over from a azimuthal filter to a radial filter and you start from the cylinder surface all the way in the full domain in the free stream, then you get a solution which actually looks like what you have in the unfiltered solution. So, basically, it tells you that you have to understand the physics of the problem, because what happens? The flow is coming from left to right and then a boundary layer kind of forms and that is separated by this rotary oscillation. However, if you now apply a filter in the theta direction, it kind of performs a mixing operation in the theta direction, right? And that mixing, numerical mixing, actually prevents separation and that is what you have seen here, that you would not have to worry about numerical separation at all; the results are wrong, that is a different issue. Whereas, so this basically counteracts the action of rotary oscillation in the flow here, the physical mechanism of the flow is determined by this rotary oscillation; and if I now filter in the theta direction, I am actually trying to take out the effect of rotary oscillation, whereas, if you do a radial filter, you see that that action of rotary oscillation is not interfered that much; and that is why you can see that you get fairly decent result. And, of course, you notice that this filtering operation itself is done over the whole domain, so you take all the points together; then, you do a filtering. So, it does matter at times, that whether you are doing the filtering from the cylinder to the free stream, that is going from the surface to outwards or from outside to the cylinder. And you can see there would be some differences and why would that difference be? That is because you are solving a tri-diagonal matrix equation, and it does depend on the operation accumulates, those errors round off errors, right? If I go from j equal to 1 to n, then I have one sequence of error accumulation; and if I go from j equal to n to 1, I have a different feature. So, there are some minute differences but which is not seen in these two frames; in a sense, it is good that at least for a problem of this kind, you do not need to worry tremendously about the direction of radial filter. But, for some combination of physical parameters, this might affect the solution also. So, basically, this is a kind of a portrait which will tell you how you are doing as compared to what people may have noticed in an experiment. This an experiment done in France, published in j f m in 2006, nice experimental visualization pictures; this is the cylinder, performing rotary oscillation, and you see nice vortical structures. These are not like your common vortex street, so that is what you can see very clearly, the vortical structures are interleaved; they are all connected together, that is exactly what you actually get if you solve it with lot of care; and care is done in this way, taken this way, this is a quite a fine grid calculation is the middle column, so, do not have any worry of filtering to do. So, that gives you a very good math with the experiment; this is kind of computation one would like to do, where your computation and experiments, they actually look rather very close to each other. Now, on the last column, we are showing some filtered solution obtained with a sixth order filter; and this interesting, because we have applied filter on the first 30 lines, close to the surface of the cylinder, with a value of alpha 0.49; whereas, the radial filter has been applied in the complete domain. So, it yes/ no. They there is not There is not, that is, the whole intension is to show that if you choose the parameters, choose the directions and choose the way you want to do, then you should not be seeing much of a difference in this case. There are actually no differences at all. So, this solution that you have shown in the middle, is to convince you, I mean, if you do not apply a filter, what sort of solution you should get? This stable case, this does not lead to instability; whereas, if I do something what I have shown you earlier, if I would have used second order, I have used azimuthal filter all over the domain, then, things would have gone pretty bad. But, if you do a high order filter and even if you apply azimuthal filter, you applied very close to the surface only, and then that does not cause any problem. And of course, the radial filter is applied over the whole domain and with a filter value that is quite drastic, 0.45 is quite a strong filtering, but even then you can see there is no difference at all between these two sets of results. So, basically, my intension is to, sort of, convince you, that choose a filter correctly and it would work to your advantage and you it would not lead to unphysical result. It is not some how to get so called stable solution, should be our goal; we need accuracy and that is rather important. Now, this is something which we done in very very recent times. This is the only work that has been reported so far. What we do is, what we had done earlier for central schemes. When we are developing compact schemes, we realize that sometimes numerical in stabilities can be cured by upwinding. So here, we try to develop an up wind filter like Yogesh did. And you actually add an additional fourth order dissipation term with a floating constant eta, that should be what we will call as the upwind coefficient, and that should give you some additional degree of freedom, and why we are doing it? I will explain to you, this is important. Because this is not just simply a another filter that you have. Now, what you are seeing is the transfer function of this upwind filter at different nodes and for different values of the upwind coefficient. the this This is eta equal to 0.001, both of them are. So, at top, you are seeing the real path; the real path will tell you how things are the first and the last point, the variables do not change at all, they are they just remain same. Now, if you look at, let us say, j equal to q and n minus 1, you get to this value, that is, this value as the lowest sets of points and as you go inwards, it keeps improving; so, filtering is respected only to the high k h and this imaginary part, actually would add on to the dispersion property; that is what we have talked about, if you recall what we said that g hat is g times t. Now, suppose I have this; of course we will have this, right? g i of the basic numerical method will be complex, so it has a real part, and let us say, these are looking at a particular node, j th node, that is what we are looking at. This I should multiply with t j; now, t j also has a real part and it will also have an imaginary part. Now, what happens to your g hat? Real would, of course, come from here; let me the g j r, and so this is the real part. And so, basically, you can see what has happened here is this; I could call it g j, or and this I will call it as g j i, right? Remember, that with the original method that we had that defined your numerical phase speed through that beta j beta j was minus tan inverse of g j i by g j r, right? So, now what will happen here? You are after filtering; you are going to get this, this will change because the real and imaginary parts have changed and that would be this tan inverse g hat of j i by g hat of j r, so, that is what happens. So, you can see, if you have a transfer function with the imaginary part, then you are going to see that beta j and beta hat j, they are not same; but, you can very clearly see that if t j i is equal to, identically equal to 0, then you will get beta j should be equal to beta j. So, there is no change in dispersion property, because that is what eventually this defines your c n by c n v g n by c. If you look at your overload, you will see that, that is what you get, c n by c is beta j by k delta t. Remember that expression that we have. So, now things do change; so, what happens is, the two waves we can bring in the transfer function to become complex; one is of course this boundary closure and the other one here is by design. We are purposely introducing a upwind filter which will have an imaginary part and which will show you this kind of a behavior; so, you will see that t j i has this kind of a property, and what does it do? That is what we like to see. Now, the previous case was for a eta equal to 0.001 and alpha was a 0.45 here. We are basically seeing the effect of real path for different upwind coefficients and you can see its its got quite a bit of a control for us, because if I take eta equal to 0.01, I can see the filtering operations start from a rather a moderate value of k h of one itself; whereas, if I keep reducing the value of eta, the filtering gets delayed to higher k, right? So, this is what you see at 0.45, alpha equal to 0.45 If you increase alpha to even higher value of 0.485, you will see that this is further delayed the filtering action for lower values of eta. So, basically, it gives you a confidence to really be able to have another parameter by which you can perform. Now, this these are the corresponding imaginary path and you can see the kind of thing is that, larger values of eta actually brings in a very very strong value of the imaginary path, and that can change dispersion property significantly, very very significantly; whereas, smaller and smaller values of eta, the effects would be marginal and restricted to very high values of k h. So, having gotten that lead, we could basically define the benefits of upwind filters. We have seen that problems arise due to numerical instability near the boundary and excessive damping. This, we could rectify by using upwind filer, and we also can allow controlled amount of dissipation in the interior; and the absolute control of this, allows one to really perform what is called as a sub grid scale model, used for large a d simulation, a very specialized branch of computing which has been in use for many decades. But, with the this kind of analysis, now we are seeing a better explanation of the same. So, basically, this sub grids scale models are related to that part of the flow which are not resolving because of the grid size. So, to be not able to resolve that scale is equivalent to adding some having some extra stress, and this is what is attempted in this sub grid scale stress models. You say a fertile area of research, but let us not worry too much about it. But, let us see what we can do, if we look at a central filter verses an upwind filter on the, right. Now, you are seeing basically the numerical amplification contours, and as you can see, that there is not much of a difference between the the two, not much of a difference between the two. Only thing is, with the upwind filter, let us say, this 0.999 value has come slightly below 1; earlier, it was above 1, so, there is a marginal difference that you see. However, what was intended in doing this exercise was is given in this figure. If we used a central filter, then we identified a region of q waves here, shown by this shaded area; so, this is your v g n by c equal to 0 line. Now, what you do is, you keep the same filter coefficient 0.45, but use an upwinding of the filter, so, eta becomes 0.001. And then, what happens? Then you see what happens is that, this zero line actually turns back and then you have a region of m c as shown by this; the green shaded area on the right, for which there is not going to be any q waves at all. So, this is something that was intended, that was the whole motivation behind designing this upwind filter; and it just does show that you have a very narrow window, but you do have where you will not suffer from that extreme dispersion effect of q waves. So, this is one thing that was achieved. And this is, we have seen before, this is another example, a tougher case, you know an aerofoil section which looks like this, is exposed to a flow coming from left to right. So, this is like a normal flat plate type of geometry and you do get very very significant vortical structure originating from the ends of this aerofoil; and there are some comparison of results, you know, very high order filters have been used; sixth and eight order filter in the azimuthal direction with alpha equal to 0.48, and these are the corresponding experimental visualization. It is very difficult to visualize these types of flows. So, the aerofoil is somewhere there, and this is the vortex that is shed from top, and this is another vortex that is shed from bottom. So, there is some kind of asymmetric of shed vortices. Because of the asymmetry of the geometry itself, the top side is rounded; whereas, the bottom side has got a very sharp edge, and this sharp edge actually brings about much more diffused larger vortex; whereas, the top surface has additional secondary separation occurring from on the surface itself, where this is the very well defined. In fact, if we have, those of you who do not belong to aerospace engineering, you may have wondered why the wing shape is like this; it is for this, that if you have a sharp trailing edge, you can actually force the flow to change locally there; whereas, if you have a rounded edge, then the effects are diffused and it can actually cause lot of unsteadiness. Whereas, this kind of sharp change in trailing edge actually helps in a better performance of the wing, and that is why the shape like this; you know, this is very wrongly understood that aircraft flight is more like a swimming of a fish than flying of a bird, because the physics is totally different a when bird flies, it actually flaps its wing; so, it is basically unsteady motion; whereas, a fish glides through the water, and if you look at a cross-section of a fish, that would be more like this. than Then, if you take a section of a bird's wing, that would look like this. Anyway, the whole idea is this; is again a very tough problem to solve, and which one could solve here by by filtering. With both, we have done it; with azimuthal filters as well as upwind filters in vulnerable directions and some of these results are shown here for different time. But, let me just simply close this discussion by looking at what we should be doing. We must develop a strategy that we have to find out an optimum filter so that optimal filter should be a combination of azimuthal filter that is applied close to the wall; and whereas, in the wall number direction, you would have a non-periodic filter. See, azimuthal direction is a special direction because you have a perfect periodicity of the problem, right? So, you do not need to worry about this so called boundary closure. So, what happens is, many a times people are very much attracted to use filter in the azimuthal direction because it is easier; you do not have to do additional problems and then, it is also has the property of retaining its property, uniformly across all nodes. However, we have seen azimuthal filters have this tendency to change the physical nature of the flow. So, what you should do is, basically, apply azimuthal filter which can be a central filter, but apply it only close to the wall; whereas, the radial direction, you can you have to use a non-periodic filter, and an upwinded filter is found to be the best. I do not wish to go through this second point; it is something to do with computing fluid mechanics. So, let us not worry about it; third point also, as I told you, it actually circumvent sometimes to do these fancy turbulence or modeling. So, this is a very very good tool for one to learn and use it imaginatively; and you could also note that whenever you use non-periodic filter, you still have the problem; numerical instability in an inflow and excessive damping in outflow. This has to be investigated analyzed and sorted out up front. So, that can be worked out by using upwind composite filter which actually, even allows you to add controlled dissipation in the interior; upwind filters has better dispersion property. We showed that it can circumvent q wave formation, and this approach of using upwind filter does not require using different equations in different parts of the flow. This is one of the problem with large a d simulation or direct a d simulation people use. So, I think will stop here on this and I would like to go into the new topic which is going to be on finite element method. So, basically, it's our way of looking at this very special branch of computing. A finite element method is all pervasive; people like to know more about finite element method, but, we have kept our focus, in fact, on the scientific computing aspects of it, which is often not very well understood by many. Let us start by discussing about the various things that we know of, that we have focused most of our attention in whole of this course on finite difference methods. One of the problems that we have noted is it cannot, well, we have not noted, because we have not talked about geometry aspect. We did not do anything about grid generation etcetera etcetera; that is a separate sub branch of computing, but finite differences are essentially more difficult to use on complex geometries. However, they have excellent resolution properties and we have seen that finite difference methods, when well-designed and produce results which are extremely accurate, which are otherwise not possible by other method. For example, we have finite element and finite volume method which are extremely attractive in handling complex geometries. Therefore, boundary conditions become easier to apply. Specifically, in finite element method, what you do is, you breakdown the computational domain into smaller subdomains and on each of these subdomains, you approximate the governing equation by either variational methods or by what is called as a weighted residual methods. Calculus of variation is a very classic field; it started off with Euler, and lots of people have contributed. Unfortunately, when it comes to some special branches of computing, like in fluid mechanics, we do not pay much attention to variational methods, because the variational principle does not exist when you are looking at dissipative flows. For example, the governing Navier-Stokes equation does not allow you to develop a variational method, whereas, the complimentary aspect of it is via this weighted residual method, that could be quite easily done. And, we will spend time talking about the next four days that we have at our disposing. Now, what you are trying to do is, in FEM, you take this smaller subdomains, and in this smaller subdomains, you try to fit the solution locally by simpler polynomials, that is, the essential idea. So, what do you have? You have a big problem. You break it into smaller pieces and then, on each of the small pieces you locally fit the solutions by polynomials and try to find out the coefficients in the polynomial that is essentially is done. So, basically, in FEM, it would appear that overall solution is approximated by a local representation which is different from global method. So, basically, let us say if I have, say some unknown u of x, I could write it by, let us say, Fourier series and I could write it as u hat of k, and I could write, sorry, let us say k j, right? And here, of course, j could go from 1 to infinity. In the actual case, this is what you know, what you do in a Fourier Spectral Method. So, this is what is called as a Fourier Spectral Method. So, what you keep doing is, of course, you cannot accept an infinite series, so, you will stop at some finite number of terms, and your intension would be in a problem to be able to define these Fourier coefficients. What you notice in this? Global method, so, this is a kind of a global method. Why do we call it a global method? I will explain a global method is one where, if you make some kind of a change, then that effect is felt globally in the whole domain. So, for example, if I change it from n to n plus 1, I have just added one more term in the Fourier representation. You would be surprised to see that this quantity will change across the whole scale; that is, that is your global method, right? That is what we are saying, that let us say, if we use a spectral method, this is a global method. But we have already seen; we are now, more or less convinced that if we take a Fourier Spectral Method, what do is, we find for c n by c is equal to 1 v g n by c equal to 1, so, spectral methods are very very attractive. However, they have some restrictions. What is the main restriction? The main restriction is, as we have written here, the problem has to be periodic, right? If we are summing by Fourier serious, it has to be a periodic problem. So, it is so happens that many physical problems, this imposition of periodicity lends itself to some sort of a unphysical attribute to the problem. So, not all problems can be... Well, most of the problems cannot be replaced by a periodic extension of the same problem; so, that is one thing. The second thing is of course, in Fourier method, you end up working in the physical domain itself. Suppose, I am looking at the problem in the Cartesian frame, then I stay there, and I am also forced to use equal spacing, and that is a major major constraint. Whereas, in contrast, what can you do in most of the other computing method? We have seen all of them suffer from various sources of error; that is what we have spent the whole course, almost talking about various sources of error; and all these other discrete method suffers from those problems. And what happens is, of course, despite that, we live with them, is simply for the reason that we can actually handle problems which are not restricted by the constraint of this special type of global method. We can talk about non-periodic problem; we can talk about non-uniform spacing where we need more accuracy; we can cluster more number of points where we do not, we cannot. So, this is something that we really need to appreciate, and that is one of the reasons that finite element from finite volume method enjoys its reputation of what it is. Let us briefly discuss about what this weighted residual methods are, so that, let us talk about a problem; It is evolution equation is written like this, l of u could be any problem that you talk about, let us say, you have a space-time dependent problem. You perform all the spatial discretization; at the end of the spatial discretization, you will have a time varying equation of this kind. And let us say this talk about a one dimensional problem for simplicity; x belongs in a domain omega and for all the problems starts at t equal to 0. So you are interested in finding out solution, evolution and the initial condition is given by u 0, defined everywhere; and in addition, you would have the boundary conditions like what you did in your second assignment; you have the time dependent boundary condition at the boundary defined as delta omega. Now, what you need to do in developing a weighted residual method, is to start with defining a trail function, so, the numerical solution we indicate by subscript capital n, basically that defines how many terms we are taking. We are showing it, this trail solution to be composed of two parts. One is here, c j t times u j x, and there is this additional term u b of x, what is this u b of x? It is that is a special function that is used there to specifically satisfy the boundary condition. You see, we are not going to do periodic condition; this will be non-periodic condition, so that, boundary condition would be given to you as has been given in three; like f of b as a function of time, so, this second part of the trail solution actually helps you in satisfying the boundary condition; whereas, at the boundary, this u j's are 0, so that, this exactly satisfies the boundary condition, right? Now, this base, these functions are called the basis function in this series, that summation term which has a, let us say, time dependent coefficient c j of t times a space dependent function. Now, this u js are the basis function and they are such that, you can satisfy the prescribed boundary condition, but that would not satisfy the initial conditions and the governing differential equation. So, that means, what you will have to use equation 4 into the governing equation and try to see if you can derive some relations among all these c js. I think I will continue with this in the next class. This needs a little bit of understanding, all these things more.