Incidentals—Part 1

Yes, I have been working hard, very hard, and have been managing a responsibility, a very difficult and demanding and emotionally very draining a responsibility in a singular capacity, and yes, I have been having problems with people—their irrationalities. And, the irrational scripts they follow. [Mind you, the reference is to scripts and not to scriptures.]


But, just check this out this one, for instance (and, the people I have in mind in the above section wouldn’t do that, I am sure; they never do pursue links from such posts of mine, especially if they are just¬†Indians—they are just self-confident, that’s all): [^].

But, leaving them aside—and I find it very, very easy do that at least in a moment like this—here is a suggestion: For tomorrow and the day after, and may be for a week or so, watch out the physics (esp. astronomy-related) Twitter-feeds, news-feeds, even blogs [complete with words like “kicking” people and all [Indians with ability to speak in English regard it as “unparliamentary,” together with words like “bloody”].


For obvious reasons, no “A Song I Like” section for this posts The news I am reporting about is exciting enough, all by itself.

Best,

–Ajit

Update on 21:35 IST the same day:

A couple of related posts are these: [^][^]

A Song I Like:

(Hindi) “meghaa chhaaye aadhee raat…”
Lyrics: Neeraj
Music: S. D. Burman
Singer: Lata Mangeshkar

 

 

Advertisements

Machine “Learning”—An Entertainment [Industry] Edition

Yes, “Machine ‘Learning’,” too, has been one of my “research” interests for some time by now. … Machine learning, esp. ANN (Artificial Neural Networks), esp. Deep Learning. …

Yesterday, I wrote a comment about it at iMechanica. Though it was made in a certain technical context, today I thought that the comment could, perhaps, make sense to many of my general readers, too, if I supply a bit of context to it. So, let me report it here (after a bit of editing). But before coming to my comment, let me first give you the context in which it was made:


Context for my iMechanica comment:

It all began with a fellow iMechanician, one Mingchuan Wang, writing a post of the title “Is machine learning a research priority now in mechanics?” at iMechanica [^]. Biswajit Banerjee responded by pointing out that

“Machine learning includes a large set of techniques that can be summarized as curve fitting in high dimensional spaces. [snip] The usefulness of the new techniques [in machine learning] should not be underestimated.” [Emphasis mine.]

Then Biswajit had pointed out an arXiv paper [^] in which machine learning was reported as having produced some good DFT-like simulations for quantum mechanical simulations, too.

A word about DFT for those who (still) don’t know about it:

DFT, i.e. Density Functional Theory, is “formally exact description of a many-body quantum system through the density alone. In practice, approximations are necessary” [^]. DFT thus is a computational technique; it is used for simulating the electronic structure in quantum mechanical systems involving several hundreds of electrons (i.e. hundreds of atoms). Here is the obligatory link to the Wiki [^], though a better introduction perhaps appears here [(.PDF) ^]. Here is a StackExchange on its limitations [^].

Trivia: Kohn and Sham received a Physics Nobel for inventing DFT. It was a very, very rare instance of a Physics Nobel being awarded for an invention—not a discovery. But the Nobel committee, once again, turned out to have put old Nobel’s money in the right place. Even if the work itself was only an invention, it did directly led to a lot of discoveries in condensed matter physics! That was because DFT was fast—it was fast enough that it could bring the physics of the larger quantum systems within the scope of (any) study at all!

And now, it seems, Machine Learning has advanced enough to be able to produce results that are similar to DFT, but without using any QM theory at all! The computer does have to “learn” its “art” (i.e. “skill”), but it does so from the results of previous DFT-based simulations, not from the theory at the base of DFT. But once the computer does that—“learning”—and the paper shows that it is possible for computer to do that—it is able to compute very similar-looking simulations much, much faster than even the rather fast technique of DFT itself.

OK. Context over. Now here in the next section is my yesterday’s comment at iMechanica. (Also note that the previous exchange on this thread at iMechanica had occurred almost a year ago.) Since it has been edited quite a bit, I will not format it using a quotation block.


[An edited version of my comment begins]

A very late comment, but still, just because something struck me only this late… May as well share it….

I think that, as Biswajit points out, it’s a question of matching a technique to an application area where it is likely to be of “good enough” a fit.

I mean to say, consider fluid dynamics, and contrast it to QM.

In (C)FD, the nonlinearity present in the advective term is a major headache. As far as I can gather, this nonlinearity has all but been “proved” as the basic cause behind the phenomenon of turbulence. If so, using machine learning in CFD would be, by the simple-minded “analysis”, a basically hopeless endeavour. The very idea of using a potential presupposes differential linearity. Therefore, machine learning may be thought as viable in computational Quantum Mechanics (viz. DFT), but not in the more mundane, classical mechanical, CFD.

But then, consider the role of the BCs and the ICs in any simulation. It is true that if you don’t handle nonlinearities right, then as the simulation time progresses, errors are soon enough going to multiply (sort of), and lead to a blowup—or at least a dramatic departure from a realistic simulation.

But then, also notice that there still is some small but nonzero interval of time which has to pass before a really bad amplification of the errors actually begins to occur. Now what if a new “BC-IC” gets imposed right within that time-interval—the one which does show “good enough” an accuracy? In this case, you can expect the simulation to remain “sufficiently” realistic-looking for a long, very long time!

Something like that seems to have been the line of thought implicit in the results reported by this paper: [(.PDF) ^].

Machine learning seems to work even in CFD, because in an interactive session, a new “modified BC-IC” is every now and then is manually being introduced by none other than the end-user himself! And, the location of the modification is precisely the region from where the flow in the rest of the domain would get most dominantly affected during the subsequent, small, time evolution.

It’s somewhat like an electron rushing through a cloud chamber. By the uncertainty principle, the electron “path” sure begins to get hazy immediately after it is “measured” (i.e. absorbed and re-emitted) by a vapor molecule at a definite point in space. The uncertainty in the position grows quite rapidly. However, what actually happens in a cloud chamber is that, before this cone of haziness becomes too big, comes along another vapor molecule, and “zaps” i.e. “measures” the electron back on to a classical position. … After a rapid succession of such going-hazy-getting-zapped process, the end result turns out to be a very, very classical-looking (line-like) path—as if the electron always were only a particle, never a wave.

Conclusion? Be realistic about how smart the “dumb” “curve-fitting” involved in machine learning can at all get. Yet, at the same time, also remain open to all the application areas where it can be made it work—even including those areas where, “intuitively”, you wouldn’t expect it to have any chance to work!

[An edited version of my comment is over. Original here at iMechanica [^]]


 

“Boy, we seem to have covered a lot of STEM territory here… Mechanics, DFT, QM, CFD, nonlinearity. … But where is either the entertainment or the industry you had promised us in the title?”

You might be saying that….

Well, the CFD paper I cited above was about the entertainment industry. It was, in particular, about the computer games industry. Go check out SoHyeon Jeong’s Web site for more cool videos and graphics [^], all using machine learning.


And, here is another instance connected with entertainment, even though now I am going to make it (mostly) explanation-free.

Check out the following piece of art—a watercolor landscape of a monsoon-time but placid sea-side, in fact. Let me just say that a certain famous artist produced it; in any case, the style is plain unmistakable. … Can you name the artist simply by looking at it? See the picture below:

A sea beach in the monsoons. Watercolor.

If you are unable to name the artist, then check out this story here [^], and a previous story here [^].


A Song I Like:

And finally, to those who have always loved Beatles’ songs…

Here is one song which, I am sure, most of you had never heard before. In any case, it came to be distributed only recently. When and where was it recorded? For both the song and its recording details, check out this site: [^]. Here is another story about it: [^]. And, if you liked what you read (and heard), here is some more stuff of the same kind [^].


Endgame:

I am of the Opinion that 99% of the “modern” “artists” and “music composers” ought to be replaced by computers/robots/machines. Whaddya think?

[Credits: “Endgame” used to be the way Mukul Sharma would end his weekly Mindsport column in the yesteryears’ Sunday Times of India. (The column perhaps also used to appear in The Illustrated Weekly of India before ToI began running it; at least I have a vague recollection of something of that sort, though can’t be quite sure. … I would be a school-boy back then, when the Weekly perhaps ran it.)]

 

Micro-level water-resources engineering—8: Measure that water evaporation! Right now!!

It’s past the middle of May—the hottest time of the year in India.

The day-time is still lengthening. And it will continue doing so well up to the summer solstice in the late June, though once monsoon arrives some time in the first half of June, the solar flux in this part of the world would get reduced due to the cloud cover, and so, any further lengthening of the day would not matter.

In the place where I these days live, the day-time temperature easily goes up to 42–44 deg. C. This high a temperature is, that way, not at all unusual for most parts of Maharashtra; sometimes Pune, which is supposed to be a city of a pretty temperate climate (mainly because of the nearby Sahyaadris), also registers the max. temperatures in the early 40s. But what makes the region where I currently live worse than Pune are these two factors: (i) the minimum temperature too stays as high as 30–32 deg. C here whereas in Pune it could easily be falling to 27–26 deg. C even during May, and (ii) the fall of the temperatures at night-time proceeds very gradually here. On a hot day, it can easily be as high as 38 deg C. even after the sunset, and even 36–37 deg. C right by the time it’s the mid-night; the drop below 35 deg. C occurs only for the 3–4 hours in the early morning, between 4 to 7 AM. In comparison, Pune is way cooler. The max. temperatures Pune registers may be similar, but the evening- and the night-time temperatures fall down much more rapidly there.

There is a lesson for the media here. Media obsesses over the max. temperature (and its record, etc.). That’s because the journos mostly are BAs. (LOL!) But anyone who has studied physics and calculus knows that it’s the integral of temperature with respect to time that really matters, because it is this quantity which scales with the total thermal energy transferred to a body. So, the usual experience common people report is correct. Despite similar max. temperatures, this place is hotter, much hotter than Pune.


And, speaking of my own personal constitution, I can handle a cold weather way better than I can handle—if at all I can handle—a hot weather. [Yes, in short, I’ve been in a bad shape for the past month or more. Lethargic. Lackadaisical. Enervated. You get the idea.]


But why is it that the temperature does not matter as much as the thermal energy does?

Consider a body, say a cube of metal. Think of some hypothetical apparatus that keeps this body at the same cool temperature at all times, say, at 20 deg. C.  Here, choose the target temperature to be lower than the minimum temperature in the day. Assume that the atmospheric temperature at two different places varies between the same limits, say, 42 to 30 deg. C. Since the target temperature is lower than the minimum ambient temperature, you would have to take heat out of the cube at all times.

The question is, at which of the two places the apparatus has to work harder. To answer that question, you have to calculate the total thermal energy that has be drained out of the cube over a single day. To answer this second question, you would need the data of not just the lower and upper limits of the temperature but also how it varies with time between two limits.


The humidity too is lower here as compared to in Pune (and, of course, in Mumbai). So, it feels comparatively much more drier. It only adds to the real feel of a real hot weather.

One does not realize it, but the existence of a prolonged high temperature makes the atmosphere here imperceptibly slowly but also absolutely insurmountably, dehydrating.

Unlike in Mumbai, one does not notice much perspiration here, and that’s because the air is so dry that any perspiration that does occur also dries up very fast. Shirts getting drenched by perspiration is not a very common sight here. Overall, desiccating would be the right word to describe this kind of an air.

So, yes, it’s bad, but you can always take precautions. Make sure to drink a couple of glasses of cool water (better still, fresh lemonade) before you step out—whether you are thirsty or not. And take an onion with you when you go out; if you begin to feel too much of heat, you can always crush the onion with hand and apply the juice onto the top of your head. [Addendum: A colleague just informed me that it’s even better to actually cut the onion and keep its cut portion touching to your body, say inside your shirt. He has spent summers in eastern Maharashtra, where temperatures can reach 47 deg. C. … Oh well!]

Also, eat a lot more onions than you normally do.

And, once you return home, make sure not to drink water immediately. Wait for 5–10 minutes. Otherwise, the body goes into a shock, and the ensuing transient spikes in your biological metabolism can, at times, even trigger the sun-stroke—which can even be fatal. A simple precaution helps avoid it.

For the same reason, take care to sit down in the shade of a tree for a few minutes before you eat that slice of water-melon. Water-melon is nothing but more than 95% water, thrown with a little sugar, some fiber, and a good measure of minerals. All in all, good for your body because even if the perspiration is imperceptible in the hot and dry regions, it is still occurring, and with it, the body is being drained of the necessary electrolytes and minerals. … Lemonades and water-melons supply the electrolytes and the minerals. People do take care not to drink lemonade in the Sun, but they don’t always take the same precaution for water-melon. Yet, precisely because a water-melon has so much water, you should take care not to expose your body to a shock. [And, oh, BTW, just in case you didn’t know already, the doctor-recommended alternative to Electral powder is: your humble lemonade! Works exactly equivalently!!]


Also, the very low levels of humidity also imply that in places like this, the desert-cooler is effective, very effective. The city shops are full of them. Some of these air-coolers sport a very bare-bones design. Nothing fancy like the Symphony Diet cooler (which I did buy last year in Pune!). The air-coolers locally made here can be as simple as just an open tray at the bottom to hold the water, a cube made of a coarse wire-mesh which is padded with the khus/wood sheathings curtain, and a robust fan operating [[very] noisily]. But it works wonderfully. And these local-made air-coolers also are very inexpensive. You can get one for just Rs. 2,500 or 3,000. I mean the ones which have a capacity to keep at least 3–4 people cool.(Branded coolers like the one I bought in Pune—and it does work even in Pune—often go above Rs. 10,000. [I bought that cooler last year because I didn’t have a job, thanks to the Mechanical Engineering Professors in the Savitribai Phule Pune University.])


That way, I also try to think of the better things this kind of an air brings. How the table salt stays so smoothly flowing, how the instant coffee powder or Bournvita never turns into a glue, how an opened packet of potato chips stays so crisp for days, how washed clothes dry up in no time…

Which, incidentally, brings me to the topic of this post.


The middle—or the second half—of May also is the most ideal time to conduct evaporation experiments.

If you are looking for a summer project, here is one: to determine the evaporation rate in your locality.

Take a couple of transparent plastic jars of uniform cross section. The evaporation rate is not very highly sensitive to the cross-sectional area, but it does help to take a vessel or a jar of sizeable diameter.

Affix a mm scale on the outside of each jar, say using cello-tape. Fill the plastic jars to some level almost to the full.

Keep one jar out in the open (exposed to the Sun), and another one, inside your home, in the shade. For the jar kept outside, make sure that birds don’t come and drink the water, thereby messing up with your measurements. For this purpose, you may surround the jar with an enclosure having a coarse mesh. The mesh must be coarse; else it will reduce the solar flux. The “reduction in the solar flux” is just a fancy [mechanical [thermal] engineering] term for saying that the mesh, if too fine, might cast too significant a shadow.

Take measurements of the heights of the water daily at a fixed time of the day, say at 6:00 PM. Conduct the experiment for a week or 10 days.

Then, plot a graph of the daily water level vs. the time elapsed, for each jar.

Realize, the rate of evaporation is measured in terms of the fall in the height, and not in terms of the volume of water lost. That’s because once the exposed area is bigger than some limit, the evaporation rate (the loss in height) is more or less independent of the cross-sectional area.

Now figure out:

Does the evaporation rate stay the same every day? If there is any significant departure from a straight-line graph, how do you explain it? Was there a measurement error? Was there an unusually strong wind on a certain day? a cloud cover?

Repeat the experiment next winter (around the new year), and determine the rate of evaporation at that time.

Later on, also make some calculations. If you are building a check-dam or a farm-pond, how much would be the evaporation loss over the five months from January to May-end? Is the height of your water storage system enough to make it practically useful? economically viable?


A Song I Like:

(Hindi) “mausam aayegaa, jaayegaa, pyaar sadaa muskuraayegaa…”
Music: Manas Mukherjee
Singers: Manna Dey and Asha Bhosale
Lyrics: Vithalbhai Patel

CFD Code Snippets—2: Transients in the Couette flow, Crank-Nicolson

Here is some Python code for modeling transients in the pure shear-driven Couette flow between two infinite horizontal flat plates.

Initially, both the plates are stationary. Then, at time t = 0, the upper plate is suddenly set into motion with a constant speed u = ue. Due to the no-slip boundary condition between the fluid and the upper plate, the topmost layer of the fluid, too, is set into motion with the same speed, i.e., u = ue. Due to viscosity (i.e. friction internal to the fluid) this topmost layer tugs on to the next layer to move in the same direction, too. (That’s sort of like how when you pull a book from the middle of the stack, the books below it tend to come out, too.) Thus, the next layer also begins to move. Then, this second layer, in turn, nudges the third layer to move. …

Yes, it is turtles—but not all the way down.

The bottom-most fluid layer feels the tug not only of the layer immediately above it, but also the cent-per-cent tug of the no-slip boundary condition of the solid plate just below it. And the lower plate remains stationary throughout.

To cut a long story short, if the situation were steady-state, i.e., if the upper plate were always moving (at a constant speed ue) and the bottom plate stationary, then the velocity profile developed would be in the shape of a straight-line.

However, if both plates are initially stationary, and then if only the upper plate is suddenly moved, then the velocity profile would be initially like the Greek letter `\Gamma‘. With the passage of time, the profile would first lose the sharp 90^{\circ} angle, and so become a bit soft at that corner. Then, the curvature at the corner would become more and more diffuse. Then, the curved profile would become less and less curvy, and finally, it would reach the straight-line of the steady-state.

The code below shows how.

The Couette flow is important because it is one of the very, very (, very, very, very) rare cases in which the Navier-Stokes equations can (at all) be solved “analytically.”

The reason is, for this kind of a Couette flow—the one driven purely by the shear force of the upper plate (i.e., when no pressure gradients are applied in the horizontal direction)—the Navier-Stokes equations reduce to nothing but (our dear old friend): diffusion equation!

We did model the transient diffusion equation the last time, here [^]. However, the discretization approach used back then was FTCS—which is an explicit approach. The one now being used is: the Crank-Nicolson method—which is an implicit approach. FTCS is only conditionally stable, whereas CN is (at least for the diffusion equation, and at least within all reasonable limits) unconditionally stable. (Check out the Flow Science’s take on the explicit vs. implicit issue, here [^].) That makes it interesting.

The code once again comes with more comments on more issues than you want to read and/or know about.

The particular test case currently hard-coded is the one from John D. Anderson, Jr.’s celebrated text [^]. You can compare the output values produced by this program with the listing given in Anderson’s text (Table 9.1, p. 430.)

The Python code here even displays a MatPlotLib plot which you can directly compare with the fig. 9.4 in Anderson’s text (p. 431).

Anyway, here is the code, without further ado:

#################################################################
# Transients in the pure shear-driven Couette flow between two
# infinite flat plates, modeled using the Crank-Nicolson method.
# Here we generally follow the treatment and terminology as
# used in the text by John D. Anderson, Jr. However, please
# note, unlike in the book, here we use the _zero_-based
# index, _not_ 1-based. 
# Code written by and copyright (c)Ajit R. Jadhav.
# Version: 03-Feb-2016.
#################################################################

# Importing of modules and routines
# We use numpy arrays for efficiency
import numpy
# We use a neat plotting library
from matplotlib import pyplot as plt
from TDMASolver import TDMASolver

#################################################################
# Program Data. All units are in SI.
#
# Gap between the parallel plates, in m
D = 0.1
# Number of cells along the y-axis
nc = 20  # 4
# Number of nodes along the y-axis
N = nc + 1
# Compute the height of each cell
dy = D / nc
# Material properties
# The values below are for water
# Density, in kg / m^3
rho = 998.2
# viscosity, in Pa-s (or kg/(m.s))
mu = 8.90e-4
# speed of the moving (upper) plate, in m/s
ue = 0.05
# Compute the effective Reynolds' number
Re = rho * ue * D / mu
print("Calculated Re: %lf" % Re)

#################################################################
# Validation
# 1. The program output has been validated against manually
# worked out solutions for the first two time-steps, for these
# values of parameters: N=5, Re = 100.0, E = 1.0. For this set
# of values, and given the default truncation that happens in the
# Python print function, the convergence to the steady-state
# solution occurs at the 30th time-step.
#
# 2. Also validated with the case given in John D. Anderson's
# text. The values of parameters are: N=21, Re = 5000, nt =1.
# The program output matches the solution listing given on
# p. 430 in the book.
#################################################################

# For validation purposes, we will ignore the above-calculated
# Re, and instead set the variable with a simple value.
# Anderson's book uses 5000---quite close to the program data
# assumed above.
Re = 5000.0

# CN being implicit, has unconditional stability. However, for
# ease of comparisons, we define a parameter E anyway.
# $E \equiv \dfrac{(\Delta t)}{(Re)(\Delta y)^2})$
# Anderson's book uses E = 1.0
E = 1.0

# Number of time-steps to take
# Anderson's book shows graph up to 240 time-steps
nt = 241  # 36

print("Assumed Re: %lf, E: %lf, nt: %d" % (Re, E, nt))

# Calculate the time-step duration
dt = E * Re * dy * dy
print("Physical time elapsed in marching one computational time-step: %lf seconds." % dt)

#################################################################
# Create the numpy arrays, initialized to zero, for holding
# the primary variable values (here, u)

# The main array, holds the values after time-marching through one step
u = numpy.zeros(N)
# A temporary array, holds the values at the previous time-step
un = numpy.zeros(N)

#################################################################
#  Setting up of the initial condition:
#
u[N - 1] = 1.0
print("Initial u: ", u)
#################################################################

# Create the y-axis vector, for plotting purposes
yAxis = numpy.linspace(0.0, 1.0, N)
plt.title("Transients in the pure shear-driven Couette flow,\nusing the Crank-Nicolson method")
plt.xlabel("Dimensionless horizontal velocity (u/ue)")
plt.ylabel("Dimensionless distance (y/D)")
plt.grid(True)
plt.plot(u, yAxis, 'r')

# Array containting the indices of time-steps at which the
# solution should be plotted:
ShowPlotInstants = [0, 2, 12, 36, 60, 240]

#################################################################
# Using the CN (Cranck-Nicolson) scheme, for the Couette Flow
#
# Compute A and B; they form a tri-diagonal matrix as below:
# B A
# A B A
#   A B A
A = -E / 2.0
B = 1.0 + E

print("A: %lf, B: %lf" % (A, B))

# a, b, c, d are the 1D arrays that define a tri-diagonal
# matrix. The tri-diagonal matrix in general looks like:
# a0|b0 c0 0 |  |x0| |d0|
#   |a1 b1 c1|  |x1|=|d1|
#   |0  a2 b2|c2|x2| |d2|
# a0 and c2 are not included in the tri-diagonal system, but
# the TDMASolver function assumes that these memory
# locations exist in the a and c arrays. Further, the
# algorithm employs the array indices in such a way that
# the first used element of the array a is a[1], not a[0],
# and similarly, the last used element of the array c is
# c[2], not c[1].
#
# For N-grid points, for the transient Couette flow problem,
# we get a system of N-2 X N-2 size, with the _global_ and
# _zero_-based _index_ going over [1, N-2]. Thus, if the
# number of grid points N = 5, we get [1,2,3].
a = numpy.zeros(N - 2)
b = numpy.zeros(N - 2)
c = numpy.zeros(N - 2)
d = numpy.zeros(N - 2)

# Initialize the a, b and c arrays
a[0:N - 2] = A
b[0:N - 2] = B
c[0:N - 2] = A

#################################################################
# Here begins the actual time-stepping...
nShowPlotIndex = 0
for n in range(nt):
    # We first save the known values of u at time n, to un,
    un = u.copy()
    # and apply the Dirichlet boundary conditions
    # to the solution at the time-step n
    un[0] = 0.0
    un[N - 1] = 1.0

    # The right hand-side vector on the CN is calculated here
    # Anderson calls the variable K
    # Note carefully the indices here. The indices used for Kj
    # are straight-forward from Anderson's book, except that here
    # they are redone as 0-based.
    # But the real change is on the index used on d. It is one
    # less, because here we don't at all store the first and the
    # last rows of the total system in it.
    for j in range(1, N - 1):
        Kj = (1.0 - E) * un[j] + E / 2.0 * (un[j + 1] + un[j - 1])
        d[j - 1] = Kj
    # For the same reason spelt above, note the indices: They are
    # respectively 1 less and 2 less.
    d[0] = d[0] - a[0] * un[0]
    d[N - 3] = d[N - 3] - c[N - 3] * un[N - 1]
    # print( d )

    # Realize, the solution array returned here has values only for
    # the interior nodes. Its elements do not include the two
    # boundary nodes.
    soln = TDMASolver(a, b, c, d)
    for j in range(1, N - 1):
        u[j] = soln[j - 1]

    # We skip plotting for all the time-steps except for those given
    # in the ShowPlotInstants array. The values are for the graph in
    # Anderson's text on p. 431
    if n == ShowPlotInstants[nShowPlotIndex]:
        nShowPlotIndex = nShowPlotIndex + 1
        # Uncomment the following line for erasing out any earlier plots
        # plt.clf()
        # Uncomment the following line if you want to print T to console
        # print( T )
        plt.plot(u, yAxis, 'k')
        plt.pause(0.1)
        print("\nVelocity profile after time-step no. %d (or physical seconds: %lf)" % (n, float(n) * dt))
        print(u)

plt.plot(u, yAxis, 'b')
plt.show()

print("done")
#################################################################
#################################################################

 

The above code calls on a TDMASolver() function written in a Python module of the same name. Here is the listing for that.

Initially, I tried writing my own TDMA (i.e. Thomas algorithm) code, but found that it was getting too time-consuming. So, I did the smart thing. I searched on the ‘net, found Ofan Tau’s code, here [^], and copy-pasted it!

One funny thing here was that Ofan in turn refers to some code on a Wikibooks page (here [^]), but while Ofan’s code is correct, the code on the Wikibooks page (at least the Python version of it) is not.

Incidentally, the best resource on TDMA aka Thomas’ algorithm that I found was one .PDF document by Prof. William T. Lee of the University of Limerick, Ireland, here [(.PDF) ^]. Keep it handy while going through this code.

#################################################################
# TDMASolver.py
#
# Solves a Tri-Diagonal Matrix system.
#
# a, b, c, d are the 1D arrays that define a tri-diagonal
# matrix. The tri-diagonal matrix in general looks like:
# a0|b0 c0 0 |  |x0| |d0|
#   |a1 b1 c1|  |x1|=|d1|
#   |0  a2 b2|c2|x2| |d2|
# a0 and c2 are not included in the tri-diagonal system, but
# the TDMASolver function assumes that these memory
# locations exist in the a and c arrays. Further, the
# algorithm employs the array indices in such a way that
# the first used element of the array a is a[1], not a[0],
# and similarly, the last used element of the array c is
# c[2], not c[1].
#
# The code here was taken from Ofan Tau's blog post:
# http://ofan666.blogspot.in/2012/02/tridiagonal-matrix-algorithm-solver-in.html
# Ofan's blog post makes reference to:
# https://en.wikibooks.org/wiki/Algorithm_Implementation/Linear_Algebra/Tridiagonal_matrix_algorithm
# It so happens that the _Python_ code at the wikibooks URL is
# _wrong_, but Ofan's code is correct (though it does not check
# for stability or even just diagonal dominance.)
# I have mostly kept Ofan's code as is; just added one or two
# comments.
#
# BTW, a wonderfully illustrated (and as far as I can tell,
# correct) reference for understanding TDMA is a .PDF document
# by Prof. William T. Lee (of Industrial Mathematics unit,
# University of Limerick, Ireland):
# http://www3.ul.ie/wlee/ms6021_thomas.pdf
#
#  -- Ajit R. Jadhav
# Version 03-February-2016
#################################################################

import numpy


#################################################################
# This helper function was written by Ajit R. Jadhav
# It uses the a, b and c vectors to fill a full-size matrix.
# It allows for matrix visualization. To be used mainly for
# debugging purposes.
def FillM(DM, a, b, c):
    N = len(a)
    DM[0][0] = b[0]
    DM[0][1] = c[0]
    for i in range(N - 1):
        DM[i][i - 1] = a[i]
        DM[i][i] = b[i]
        DM[i][i + 1] = c[i]
    DM[N - 1][N - 2] = a[N - 1]
    DM[N - 1][N - 1] = b[N - 1]
    return DM


#################################################################
# This function was shamelessly lifted from Ofan's blog-post:
# http://ofan666.blogspot.in/2012/02/tridiagonal-matrix-algorithm-solver-in.html,
# accessed 02-February-2016
#################################################################
def TDMASolver(a, b, c, d):
    '''
    TDMA solver, a b c d can be NumPy array type or Python list type.
    refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

    # a, b, c, d are the 1D arrays that define a tri-diagonal
    # matrix. The triadiagonal matrix in general looks like:
    # a0|b0 c0 0 |  |x0| |d0|
    #   |a1 b1 c1|  |x1|=|d1|
    #   |0  a2 b2|c2|x2| |d2|
    # a0 and c2 are not included in the triadiagonal system, but
    # the TDMASolver function assumes that these memory
    # locations exist in the a and c arrays. Further, the
    # algorithm employs the array indices in such a way that
    # the first used element of the array a is a[1], not a[0],
    # and similarly, the last used element of the array c is
    # c[2], not c[1].

    # One test case is this:
    # a = [0, -1.0, -1.0, -1.0]
    # b = [4.0,  4.0,  4.0,  4.0]
    # c = [-1.0, -1.0, -1.0,  0.0]
    # d = [5.0,  5.0, 10.0, 23.0]
    # result: [2, 3, 5, 7]
    '''
    # number of equations
    nf = len(a)
    # print( "In TDMASolver.py: Size of the system: %d" % nf)

    # Make a copy the arrays
    ac, bc, cc, dc = map(numpy.array, (a, b, c, d))

    ##############################################
    # Make the matrix strictly upper triangular
    for i in range(1, nf):
        mc = ac[i] / bc[i - 1]
        bc[i] = bc[i] - mc * cc[i - 1]
        dc[i] = dc[i] - mc * dc[i - 1]

    ##############################################
    # The back-substitution pass:
    # Note the meaning of the -1 index in Python.
    xc = ac
    xc[-1] = dc[-1] / bc[-1]

    for j in range(nf - 2, -1, -1):
        xc[j] = (dc[j] - cc[j] * xc[j + 1]) / bc[j]

    # DM = numpy.zeros((nf,nf))
    # DM = FillM(DM, a, b, c)
    # print(DM)
    # del DM

    # delete variables from memory
    del bc, cc, dc

    return xc

Some things to try:

  1. Change E to 5.0, 10.0, and then even to 4000.0, and see if you get a solution like the one given in Anderson’s text (fig. 9.5, p. 934)
  2. Remove the hard-coded value of Reynolds’ number, and start using the calculated values for the same. Then, play with the physical parameters: use material properties for honey, glycerine, different SAE grades of motor oils, etc., and see if you get reasonable-looking “settling-in” time periods to approach steady-state or not.
  3. Also change the separation between the flat plates.
  4. What would happen if pressure gradients are applied to the fluid? In what way will the governing differential equation change? The velocity profile? Could you handle it using FDM?
  5. And, oh, yes! Do continue thinking about the SPPU (Savitribai Phule Pune University), and its past, present and future “authorities.”

[As usual, this post could change—improve or deprove—without any prior [or posterior] notice.]

[E&OE]