I am keeping my New Year’s…

I am keeping my NYR [^], made last year.

How about you?


No, really. I AM keeping my NYR. Here’s how.


December is meant for making resolutions. (It doesn’t matter whether it’s the 1st or the 31st; the month is [the?] December; that’s what matters.)

Done.

January is meant for making a time-table. … But it must be something on which you can execute. I have been actively engaged doing that. … You could see that, couldn’t you? … And, what’s more, you could’ve bet about it at any time in the past, too, couldn’t you?

Since execution can only follow, and not precede, planning, it must be February before execution proper itself can begin. As far as I am concerned, I will make sure of that. [And you know me. You know that I always deliver on all my promises, don’t you?]

March is known for madness. To be avoided, of course.

April is known for foolishness. To be avoided, as far as possible, but, hey, as “friends” of this blog, you know, it’s nothing to be afraid of!

May, in this part of the world, is far too hot for any one to handle it right, OK? The work-efficiency naturally drops down. This fact must be factored into any and every good piece of planning, I say! (Recall the British Governors and Other officers of the Bombay Presidency shifting their offices to Matheran/Mahabaleshwar during summer? [Anyone ever cared to measure the efficiency of this measure on their part? I mean, on work?])

Now, yes, June does bring in the [very welcome] monsoons, finally! But then, monsoon also is universally known to be the most romantic of all seasons. [It leaves a certain something of a feeling which ordinarily would require you to down a Sundowner or so. [I am trying to be honest, here!]… And then, even Kalidas would seem to agree. Remember (Sanskrit) “aashaaDhasya pratham…”? Naturally, the month is not very conducive to work, is it?]


OK.


This is [just] January, and my time-table is all done up and ready. Or, at least, it’s [at least] half-way through. …

I will really, really begin work in the second half of the year.

Bye until then.


A Song I Don’t Ever Recall Liking Back Then [When Things Mattered Far More Routinely in Far More Respects than They Do Today]

[Not too sure I like it today either. But there were certain happy isolated instances related to my more recent career which are associated with it. I had registered, but hadn’t known this fact, until recently.

But then, recently, I happened suddenly to “re-hear” the phrase (Hindi) “yeh kaunsaa…”, complete with the piece of the “sax” which follows it…

Then, the world had become [in a [comparatively] recent past] a slightly better place to live in.

So, I’d decided, not quite fully certain but still being inclined to this possibility, that I might actually like this song. … But I still don’t fully, you know… But I still do fully want to run it, you know…

Anyway, just listen to it…]

(Hindi) “chocolate, lime juice, ice-cream…” [No, it really is a Hindi song. Just listen to it further…]
Singer: Lata Mangeshkar [A peculiarity of this song is that precisely when [an aged] Lata sounds [a bit] heavy [of course due to the age not to mention the pressures of the day-to-day work and every one’s normal inability to hit the sweet spot every time!], the directors of the movie and the music together focus your attention on a rather cheerfully smiling and dancing Madhuri. [No, never been one of my favorite actresses, but then, that’s an entirely different story altogether.]]
Music: Ramlaxman
Lyrics: Dev Kohli [?]

[PS: And, coming to the video of this song, did you notice that the hero drives a Maruti Gypsy?

I mean, ask any NRI in USA, and they he will tell you that it was because this was an early 90’s movie; the fruits of the [half-/quarter-/oct-something-/etc.] economic liberalization had still not been had by the general public; the liberalization they [I mean these NRIs] had brought about.

If these [I mean the economic freedoms] were to be brought about , they could easily point out, with good amount of references to Hindi movies of the recent years, that the presence on Indian roads of the [government-subsidized-diesel-driven] SUVs could easily have been seen in the same movie!!!

Hmmm…  Point[s] taken.]


How about your NYR?


[A bit of an editing is still due, I am sure… TBD, when I get the time to do so…]

Advertisements

Yes I know it!

Note: A long update was posted on 12th December 2017, 11:35 IST.


This post is spurred by my browsing of certain twitter feeds of certain pop-sci. writers.

The URL being highlighted—and it would be, say, “negligible,” but for the reputation of the Web domain name on which it appears—is this: [^].


I want to remind you that I know the answers to all the essential quantum mysteries.

Not only that, I also want to remind you that I can discuss about them, in person.

It’s just that my circumstances—past, and present (though I don’t know about future)—which compel me to say, definitely, that I am not available for writing it down for you (i.e. for the layman) whether here or elsewhere, as of now. Neither am I available for discussions on Skype, or via video conferencing, or with whatever “remoting” mode you have in mind. Uh… Yes… WhatsApp? Include it, too. Or something—anything—like that. Whether such requests come from some millionaire Indian in USA (and there are tons of them out there), or otherwise. Nope. A flat no is the answer for all such requests. They are out of question, bounds… At least for now.

… Things may change in future, but at least for the time being, the discussions would have to be with those who already have studied (the non-relativistic) quantum physics as it is taught in universities, up to graduate (PhD) level.

And, you have to have discussions in person. That’s the firm condition being set (for the gain of their knowledge 🙂 ).


Just wanted to remind you, that’s all!


Update on 12th December 2017, 11:35 AM IST:

I have moved the update to a new post.

 


A Song I Like:

(Western, Instrumental) “Berlin Melody”
Credits: Billy Vaughn

[The same 45 RPM thingie [as in here [^], and here [^]] . … I was always unsure whether I liked this one better or the “Come September” one. … Guess, after the n-th thought, that it was this one. There is an odd-even thing about it. For odd ‘n” I think this one is better. For even ‘n’, I think the “Come September” is better.

… And then, there also are a few more musical goodies which came my way during that vacation, and I will make sure that they find their way to you too….

Actually, it’s not the simple odd-even thing. The maths here is more complicated than just the binary logic. It’s an n-ary logic. And, I am “equally” divided among them all. (4+ decades later, I still remain divided.)… (But perhaps the “best” of them was a Marathi one, though it clearly showed a best sort of a learning coming from also the Western music. I will share it the next time.)]


[As usual, may be, another revision [?]… Is it due? Yes, one was due. Have edited streamlined the main post, and then, also added a long update on 12th December 2017, as noted above.]

 

 

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

 

 

Fluxes, scalars, vectors, tensors…. and, running in circles about them!

0. This post is written for those who know something about Thermal Engineering (i.e., fluid dynamics, heat transfer, and transport phenomena) say up to the UG level at least. [A knowledge of Design Engineering, in particular, the tensors as they appear in solid mechanics, would be helpful to have but not necessary. After all, contrary to what many UGC and AICTE-approved (Full) Professors of Mechanical Engineering teaching ME (Mech – Design Engineering) courses in SPPU and other Indian universities believe, tensors not only appear also in fluid mechanics, but, in fact, the fluids phenomena make it (only so slightly) easier to understand this concept. [But all these cartoons characters, even if they don’t know even this plain and simple a fact, can always be fully relied (by anyone) about raising objections about my Metallurgy background, when it comes to my own approval, at any time! [Indians!!]]]

In this post, I write a bit about the following question:

Why is the flux \vec{J} of a scalar \phi a vector quantity, and not a mere number (which is aka a “scalar,” in certain contexts)? Why is it not a tensor—whatever the hell the term means, physically?

And, what is the best way to define a flux vector anyway?


1.

One easy answer is that if the flux is a vector, then we can establish a flux-gradient relationship. Such relationships happen to appear as statements of physical laws in all the disciplines wherever the idea of a continuum was found useful. So the scope of the applicability of the flux-gradient relationships is very vast.

The reason to define the flux as a vector, then, becomes: because the gradient of a scalar field is a vector field, that’s why.

But this answer only tells us about one of the end-purposes of the concept, viz., how it can be used. And then the answer provided is: for the formulation of a physical law. But this answer tells us nothing by way of the very meaning of the concept of flux itself.


2.

Another easy answer is that if it is a vector quantity, then it simplifies the maths involved. Instead of remembering having to take the right \theta and then multiplying the relevant scalar quantity by the \cos of this \theta, we can more succinctly write:

q = \vec{J} \cdot \vec{S} (Eq. 1)

where q is the quantity of \phi, an intensive scalar property of the fluid flowing across a given finite surface, \vec{S}, and \vec{J} is the flux of \Phi, the extensive quantity corresponding to the intensive quantity \phi.

However, apart from being a mere convenience of notation—a useful shorthand—this answer once again touches only on the end-purpose, viz., the fact that the idea of flux can be used to calculate the amount q of the transported property \Phi.

There also is another problem with this, second, answer.

Notice that in Eq. 1, \vec{J} has not been defined independently of the “dotting” operation.

If you have an equation in which the very quantity to be defined itself has an operator acting on it on one side of an equation, and then, if a suitable anti- or inverse-operator is available, then you can apply the inverse operator on both sides of the equation, and thereby “free-up” the quantity to be defined itself. This way, the quantity to be defined becomes available all by itself, and so, its definition in terms of certain hierarchically preceding other quantities also becomes straight-forward.

OK, the description looks more complex than it is, so let me illustrate it with a concrete example.

Suppose you want to define some vector \vec{T}, but the only basic equation available to you is:

\vec{R} = \int \text{d} x \vec{T}, (Eq. 2)

assuming that \vec{T} is a function of position x.

In Eq. 2, first, the integral operator must operate on \vec{T}(x) so as to produce some other quantity, here, \vec{R}. Thus, Eq. 2 can be taken as a definition for \vec{R}, but not for \vec{T}.

However, fortunately, a suitable inverse operator is available here; the inverse of integration is differentiation. So, what we do is to apply this inverse operator on both sides. On the right hand-side, it acts to let \vec{T} be free of any operator, to give you:

\dfrac{\text{d}\vec{R}}{\text{d}x} = \vec{T} (Eq. 3)

It is the Eq. 3 which can now be used as a definition of \vec{T}.

In principle, you don’t have to go to Eq. 3. In principle, you could perhaps venture to use a bit of notation abuse (the way the good folks in the calculus of variations and integral transforms always did), and say that the Eq. 2 itself is fully acceptable as a definition of \vec{T}. IMO, despite the appeal to “principles”, it still is an abuse of notation. However, I can see that the argument does have at least some point about it.

But the real trouble with using Eq. 1 (reproduced below)

q = \vec{J} \cdot \vec{S} (Eq. 1)

as a definition for \vec{J} is that no suitable inverse operator exists when it comes to the dot operator.


3.

Let’s try another way to attempt defining the flux vector, and see what it leads to. This approach goes via the following equation:

\vec{J} \equiv \dfrac{q}{|\vec{S}|} \hat{n} (Eq. 4)

where \hat{n} is the unit normal to the surface \vec{S}, defined thus:

\hat{n} \equiv \dfrac{\vec{S}}{|\vec{S}|} (Eq. 5)

Then, as the crucial next step, we introduce one more equation for q, one that is independent of \vec{J}. For phenomena involving fluid flows, this extra equation is quite simple to find:

q = \phi \rho \dfrac{\Omega_{\text{traced}}}{\Delta t} (Eq. 6)

where \phi is the mass-density of \Phi (the scalar field whose flux we want to define), \rho is the volume-density of mass itself, and \Omega_{\text{traced}} is the volume that is imaginarily traced by that specific portion of fluid which has imaginarily flowed across the surface \vec{S} in an arbitrary but small interval of time \Delta t. Notice that \Phi is the extensive scalar property being transported via the fluid flow across the given surface, whereas \phi is the corresponding intensive quantity.

Now express \Omega_{\text{traced}} in terms of the imagined maximum normal distance from the plane \vec{S} up to which the forward moving front is found extended after \Delta t. Thus,

\Omega_{\text{traced}} = \xi |\vec{S}| (Eq. 7)

where \xi is the traced distance (measured in a direction normal to \vec{S}). Now, using the geometric property for the area of parallelograms, we have that:

\xi = \delta \cos\theta (Eq. 8)

where \delta is the traced distance in the direction of the flow, and \theta is the angle between the unit normal to the plane \hat{n} and the flow velocity vector \vec{U}. Using vector notation, Eq. 8 can be expressed as:

\xi = \vec{\delta} \cdot \hat{n} (Eq. 9)

Now, by definition of \vec{U}:

\vec{\delta} = \vec{U} \Delta t, (Eq. 10)

Substituting Eq. 10 into Eq. 9, we get:

\xi = \vec{U} \Delta t \cdot \hat{n} (Eq. 11)

Substituting Eq. 11 into Eq. 7, we get:

\Omega_{\text{traced}} = \vec{U} \Delta t \cdot \hat{n} |\vec{S}| (Eq. 12)

Substituting Eq. 12 into Eq. 6, we get:

q = \phi \rho \dfrac{\vec{U} \Delta t \cdot \hat{n} |\vec{S}|}{\Delta t} (Eq. 13)

Cancelling out the \Delta t, Eq. 13 becomes:

q = \phi \rho \vec{U} \cdot \hat{n} |\vec{S}| (Eq. 14)

Having got an expression for q that is independent of \vec{J}, we can now use it in order to define \vec{J}. Thus, substituting Eq. 14 into Eq. 4:

\vec{J} \equiv \dfrac{q}{|\vec{S}|} \hat{n} = \dfrac{\phi \rho \vec{U} \cdot \hat{n} |\vec{S}|}{|\vec{S}|} \hat{n} (Eq. 16)

Cancelling out the two |\vec{S}|s (because it’s a scalar—you can always divide any term by a scalar (or even  by a complex number) but not by a vector), we finally get:

\vec{J} \equiv \phi \rho \vec{U} \cdot \hat{n} \hat{n} (Eq. 17)


4. Comments on Eq. 17

In Eq. 17, there is this curious sequence: \hat{n} \hat{n}.

It’s a sequence of two vectors, but the vectors apparently are not connected by any of the operators that are taught in the Engineering Maths courses on vector algebra and calculus—there is neither the dot (\cdot) operator nor the cross \times operator appearing in between the two \hat{n}s.

But, for the time being, let’s not get too much perturbed by the weird-looking sequence. For the time being, you can mentally insert parentheses like these:

\vec{J} \equiv \left[ \left( \phi \rho \vec{U} \right) \cdot \left( \hat{n} \right) \right] \hat{n} (Eq. 18)

and see that each of the two terms within the parentheses is a vector, and that these two vectors are connected by a dot operator so that the terms within the square brackets all evaluate to a scalar. According to Eq. 18, the scalar magnitude of the flux vector is:

|\vec{J}| = \left( \phi \rho \vec{U}\right) \cdot \left( \hat{n} \right) (Eq. 19)

and its direction is given by: \hat{n} (the second one, i.e., the one which appears in Eq. 18 but not in Eq. 19).


5.

We explained away our difficulty about Eq. 17 by inserting parentheses at suitable places. But this procedure of inserting mere parentheses looks, by itself, conceptually very attractive, doesn’t it?

If by not changing any of the quantities or the order in which they appear, and if by just inserting parentheses, an equation somehow begins to make perfect sense (i.e., if it seems to acquire a good physical meaning), then we have to wonder:

Since it is possible to insert parentheses in Eq. 17 in some other way, in some other places—to group the quantities in some other way—what physical meaning would such an alternative grouping have?

That’s a delectable possibility, potentially opening new vistas of physico-mathematical reasonings for us. So, let’s pursue it a bit.

What if the parentheses were to be inserted the following way?:

\vec{J} \equiv \left( \hat{n} \hat{n} \right) \cdot \left( \phi \rho \vec{U} \right) (Eq. 20)

On the right hand-side, the terms in the second set of parentheses evaluate to a vector, as usual. However, the terms in the first set of parentheses are special.

The fact of the matter is, there is an implicit operator connecting the two vectors, and if it is made explicit, Eq. 20 would rather be written as:

\vec{J} \equiv \left( \hat{n} \otimes \hat{n} \right) \cdot \left( \phi \rho \vec{U} \right) (Eq. 21)

The \otimes operator, as it so happens, is a binary operator that operates on two vectors (which in general need not necessarily be one and the same vector as is the case here, and whose order with respect to the operator does matter). It produces a new mathematical object called the tensor.

The general form of Eq. 21 is like the following:

\vec{V} = \vec{\vec{T}} \cdot \vec{U} (Eq. 22)

where we have put two arrows on the top of the tensor, to bring out the idea that it has something to do with two vectors (in a certain order). Eq. 22 may be read as the following: Begin with an input vector \vec{U}. When it is multiplied by the tensor \vec{\vec{T}}, we get another vector, the output vector: \vec{V}. The tensor quantity \vec{\vec{T}} is thus a mapping between an arbitrary input vector and its uniquely corresponding output vector. It also may be thought of as a unary operator which accepts a vector on its right hand-side as an input, and transforms it into the corresponding output vector.


6. “Where am I?…”

Now is the time to take a pause and ponder about a few things. Let me begin doing that, by raising a few questions for you:

Q. 6.1:

What kind of a bargain have we ended up with? We wanted to show how the flux of a scalar field \Phi must be a vector. However, in the process, we seem to have adopted an approach which says that the only way the flux—a vector—can at all be defined is in reference to a tensor—a more advanced concept.

Instead of simplifying things, we seem to have ended up complicating the matters. … Have we? really? …Can we keep the physical essentials of the approach all the same and yet, in our definition of the flux vector, don’t have to make a reference to the tensor concept? exactly how?

(Hint: Look at the above development very carefully once again!)

Q. 6.2:

In Eq. 20, we put the parentheses in this way:

\vec{J} \equiv \left( \hat{n} \hat{n} \right) \cdot \left( \phi \rho \vec{U} \right) (Eq. 20, reproduced)

What would happen if we were to group the same quantities, but alter the order of the operands for the dot operator?  After all, the dot product is commutative, right? So, we could have easily written Eq. 20 rather as:

\vec{J} \equiv \left( \phi \rho \vec{U} \right) \cdot \left( \hat{n} \hat{n} \right) (Eq. 21)

What could be the reason why in writing Eq. 20, we might have made the choice we did?

Q. 6.3:

We wanted to define the flux vector for all fluid-mechanical flow phenomena. But in Eq. 21, reproduced below, what we ended up having was the following:

\vec{J} \equiv \left( \phi \rho \vec{U} \right) \cdot \left( \hat{n} \otimes \hat{n} \right) (Eq. 21, reproduced)

Now, from our knowledge of fluid dynamics, we know that Eq. 21 seemingly stands only for one kind of a flux, namely, the convective flux. But what about the diffusive flux? (To know the difference between the two, consult any good book/course-notes on CFD using FVM, e.g. Jayathi Murthy’s notes at Purdue, or Versteeg and Malasekara’s text.)

Q. 6.4:

Try to pursue this line of thought a bit:

Start with Eq. 1 again:

q = \vec{J} \cdot \vec{S} (Eq. 1, reproduced)

Express \vec{S} as a product of its magnitude and direction:

q = \vec{J} \cdot |\vec{S}| \hat{n} (Eq. 23)

Divide both sides of Eq. 23 by |\vec{S}|:

\dfrac{q}{|\vec{S}|} = \vec{J} \cdot \hat{n} (Eq. 24)

“Multiply” both sides of Eq. 24 by \hat{n}:

\dfrac{q} {|\vec{S}|} \hat{n} = \vec{J} \cdot \hat{n} \hat{n} (Eq. 25)

We seem to have ended up with a tensor once again! (and more rapidly than in the development in section 4. above).

Now, looking at what kind of a change the left hand-side of Eq. 24 undergoes when we “multiply” it by a vector (which is: \hat{n}), can you guess something about what the “multiplication” on the right hand-side by \hat{n} might mean? Here is a hint:

To multiply a scalar by a vector is meaningless, really speaking. First, you need to have a vector space, and then, you are allowed to take any arbitrary vector from that space, and scale it up (without changing its direction) by multiplying it with a number that acts as a scalar. The result at least looks the same as “multiplying” a scalar by a vector.

What then might be happening on the right hand side?

Q.6.5:

Recall your knowledge (i) that vectors can be expressed as single-column or single-row matrices, and (ii) how matrices can be algebraically manipulated, esp. the rules for their multiplications.

Try to put the above developments using an explicit matrix notation.

In particular, pay particular attention to the matrix-algebraic notation for the dot product between a row- or column-vector and a square matrix, and the effect it has on your answer to question Q.6.2. above. [Hint: Try to use the transpose operator if you reach what looks like a dead-end.]

Q.6.6.

Suppose I introduce the following definitions: All single-column matrices are “primary” vectors (whatever the hell it may mean), and all single-row matrices are “dual” vectors (once again, whatever the hell it may mean).

Given these definitions, you can see that any primary vector can be turned into its corresponding dual vector simply by applying the transpose operator to it. Taking the logic to full generality, the entirety of a given primary vector-space can then be transformed into a certain corresponding vector space, called the dual space.

Now, using these definitions, and in reference to the definition of the flux vector via a tensor (Eq. 21), but with the equation now re-cast into the language of matrices, try to identify the physical meaning the concept of “dual” space. [If you fail to, I will sure provide a hint.]

As a part of this exercise, you will also be able to figure out which of the two \hat{n}s forms the “primary” vector space and which \hat{n} forms the dual space, if the tensor product \hat{n}\otimes\hat{n} itself appears (i) before the dot operator or (ii) after the dot operator, in the definition of the flux vector. Knowing the physical meaning for the concept of the dual space of a given vector space, you can then see what the physical meaning of the tensor product of the unit normal vectors (\hat{n}s) is, here.

Over to you. [And also to the UGC/AICTE-Approved Full Professors of Mechanical Engineering in SPPU and in other similar Indian universities. [Indians!!]]

A Song I Like:

[TBD, after I make sure all LaTeX entries have come out right, which may very well be tomorrow or the day after…]

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]

 

CFD Code Snippets—1: 1D transient conduction, FTCS, Dirichlet BCs

This semester, I am teaching an elective UG course on CFD. For illustrating ideas concerning the title problem, I began by borrowing the Python code written by Prof. Lorena Barba (general URL: [^]; the step concerning diffusion: [^] )), and then altered it a bit. I am sharing my modified version here.

In particular, here is a summary of my changes: (i) The Dirichlet BCs are now applied (so that time-stepping can go on for an indefinitely long period), (ii) graph plotting is done right within the time-marching loop, (iii) three different initial conditions can be used, etc. Variable names have been changed, and in fact, even comments have been expanded/corrected. (For instance, I use the term “thermal diffusivity” instead of“viscosity,” and use the actual material property value for it.)


This code can be used for the laboratory component (aka “Term-Work”) for the final-year elective course on CFD, in the S. P. Pune University.

This is the same University that denies me approval as a Professor, now because they think that I don’t have the required amount of work-experience (i.e. 10 years of both industrial, research and teaching experience put together).

The thing is, I didn’t always work for companies with annual turnover of 15 crores or more.

According to the PR guy of my employers, this university is now asking for my previous employers (e.g. startups in India/abroad) to certify in writing on their original letterheads not only the fact that they had employed me, but also to identify whether they had had 15 crores or more of annual turnover during the times when they had employed me or not. The rule is new, but applies retrospectively, i.e., it is supposed to go back even to the late 1990s or early ’00s times, when I worked in the software startups, sometimes with barely 1 or 2 employees too in their initial stages.

Most of these companies are by now wound up (e.g. QueriSoft grew to 500+ people whereas I was may be 13th/14th employee in Pune; nirWANa Inc. and Imagine Tech. where I was the fourth employee, and the second employee etc.—with nirWANa growing to US $1 Million in VC funding (coming from Kanwal Rekhi) but only after I left (at which point it had grown to 60+ employees)). I cannot have their certificates either because they can no longer be contacted, or even if the VC funding was US $ 1 Million, the sales turnover was not, etc.

Thus, the quantum of my total experience, according to the certifiable madness of the “authorities” at the S. P. University of Pune, suddenly dwindles down.

This rule, thus, allows the S. P. University of Pune to happily debar me from becoming a Professor.

It also allows employers to employ people like me not following the UGC scale, but on a consolidated, lower-level, salary.

I do not know precisely who among the two parties wishes to emphasize the rule more, to short-sell me. However, here, I am inclined to believe that it is the former. After all, it has been the government-run COEP and the (now S.P.) University of Pune who have happily supported people like Dr. Ashok Ghatol and Dr. (!) G. K. Kharate.

If professors who have worked nowhere else but in the private engineering colleges affiliated to the university rules are to be applied the same rule, then many of them would disqualify. This contradiction apparently does not matter to the University “authorities.” The rule is retrospective in action, but only for the new applicants, as per what our PR gentleman told me.

Anyway, I am going to share my code, but let me do so with a forthright moral denunciation of those office-bearers of the S. P. Pune University (and/or Maharashtra government) who are (or have been) responsible for framing such idiotic, harassing, anti-life and overall morally putrid rules. (For all I know, there could even be a casteist angle to it all; after all, it’s S. P. Pune University!)


Anyway, here is the Python code. It’s been written in v. 2.7 of Python. You will also need numpy and matplotlib installed. Feel free to use any IDE; I use PyCharm Community Edition.

#################################################################
# 1D Transient Conduction, using FTCS, with Dirichlet BCs.
# Code written by and copyright (c)Ajit R. Jadhav.
# Based on the code by Prof. Lorena Barba (12 Steps to NS).
# Shared with an explicit identification of how morally putrid 
# the S. P. Pune University is.
# Version: 27-Jan-2016.

#################################################################
# Importing Routines
# Use numpy arrays for efficiency
import numpy
# Use a neat plotting library
from matplotlib import pyplot  as plt

#################################################################
#  Program Data. All units are in SI
#
# Domain size, in meters
domLen = 0.1
# Number of cells
nc = 20
# Number of nodes
nx = nc + 1
# Compute the width of each cell
dx = domLen / nc
# Thermal diffusivity. This value is for Alumninium, in m^2/s
alpha = 8.418e-5
# Physical time over which the simulation is to be run
totalTime = 10.0
# Temperatures in degree Celcius, to use in constructing
# the initial and boundary conditions
tempLevel_1 = 0.0
tempLevel_2 = 100.0

# Parrot out the main input data
print("Length of the domain, in m: %lf" % domLen)
print("Number of cells: %d. Number of nodes: %d" % (nc, nx))
print("Therefore, the cell-width, in m: %lf" % dx)
print("Thermal diffusivity, in m^2/s: %lf" % alpha)
print("Duration (in physical seconds) for the simulation: %lf" % totalTime)
print("Temperature levels (deg. C): Level 1: %lf, Level 2: %lf" % (tempLevel_1, tempLevel_2))
# Main program data over
#################################################################

# Compute the duration of each time-step, and the number of time-steps required
# Stability criterion: r = \alpha \dfrac{\Delta t}{(\Delta x)^2}
# r <= 0.5 for stability
r = .35
print("Assumed stability criterion value: %lf" % r)
# Time-interval between steps, in seconds
dt = r * dx ** 2 / alpha
print("Duration of each time-step in s, after applying the stability criterion: %lf" % dt)
# Compute the total no. of time-steps (an integer) to fit within the given time-interval
nt = int(totalTime / dt)
print("Number of time-steps required to step through: %d" % nt)
# No. of time-steps to skip for plotting
stepsToSkipInPlotting = 10

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

# The main array, holds the values after time-marching through one step
T = numpy.zeros(nx)
# A temporary array, holds the values at the previous time-step
Tn = numpy.zeros(nx)
print(T)
#################################################################
#  Setting up of the initial condition:
#
#  Uncomment the required IC and keep others commented

#########################
# IC1: A single pulse in the middle
# TL = tempLevel_1
# TR = tempLevel_1
# T[ nx/2 ] = tempLevel_2

#########################
# IC2: A step function in the middle
# TL = tempLevel_2
# TR = tempLevel_2
# T[ nx/3.0 : 2.0*nx/3.0  ] = tempLevel_1

#########################
# IC3: Left half at tempLevel_1; right half at tempLevel_2
TL = tempLevel_1
TR = tempLevel_2
T[0: nx / 2] = TL
T[nx / 2 + 1: nx] = TR

#
# Setting up of the initial condition is over
#################################################################

# Create the x-axis vector, for plotting purposes
# About the linspace function
#   numpy.linspace(start, stop, num, endpoint=True, retstep=False, dtype=None)
#   Returns evenly spaced numbers over a specified interval.
xAxis = numpy.linspace(0, domLen, nx)
plt.title("1D Transient Conduction, using FTCS")
plt.xlabel("Distance, in m")
plt.ylabel("Temperature, in deg C")
plt.plot(xAxis, T, 'r')

print("Initial condition is:")
print(T)

#################################################################
#  Time-stepping, using FTCS
#
# Multiplying factor to be used at each time-step
MF = alpha * dt / dx ** 2

# Python range(nSize) returns an array of size nSize
# Thus, range(3) returns [0,1,2]
for n in range(nt):
    # We save the values of T at time n, to Tn,
    Tn = T.copy()
    # and apply the Dirichlet boundary conditions
    # to the solution at the time-step n
    Tn[0] = TL
    Tn[nx - 1] = TR

    # Then, we advance time to (n+1) and set T with it
    # A note on the range() function, and the index values.
    # -- The index to the right-most node in the arrays
    #    T and Tn is: nx-1. If we specify range( nx-1 ),
    #    then i will at most become _equal_ to nx-2
    #    i.e., one less than what we specified.
    for i in range(1, nx - 1):
        # This is the FTCS algorithm
        T[i] = Tn[i] + MF * (Tn[i - 1] - 2.0 * Tn[i] + Tn[i + 1])

    # We skip `stepsToSkipInPlotting' number of steps, in between two plots
    if n and not (n % stepsToSkipInPlotting):
        # 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(xAxis, T, 'k')
        plt.pause(0.1)

# The solution after nt number of time-steps
print("The solution after %d steps i.e. %lf physical seconds is: " % (nt, nt * dt))
print(T)
plt.plot(xAxis, T, 'b')
plt.show()

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

Some things to try:

  1. Try different values of the parameter <code>r</code>. For instance, set <code>r = 0.55</code> and see the instability develop.
  2. Try mesh refinement.
  3. Try changing material properties.
  4. Try different initial conditions. For instance, try triangular profile, or a sine-wave profile.
  5. Try setting a higher temperature, and keep the `r’ parameter very close to 0.5. See whether it is the initial condition of a triangle or of a sharp pulse which first develops instability. Try to understand the fact that the von Neumann stability analysis is only linear, and therefore provides only an initial estimate for the stability criterion.
  6. Compare in detail with the analytical solutions. For instance, the analytical solution for the initial triangular profile is worked out in Kreyszig.
  7. If inclined towards programming, try using fft routines to work out analytical solutions, have the program save both the FTCS and the FFT-based solutions to output files, and then write short scripts to plot the errors as functions of, e.g., the shape of the initial profile, space (location), time-step, mesh-refinement, `r’ parameter, etc.

[May be one more editing pass, sometime later this week or so.]

[E&OE]

The Mechanical-vs-Metallurgy “Branch-Jumping” Issue—Part III: A Couple of Stunning Papers

0. The earlier posts in this series may be found here [^] and here [^].

In this post, I am going to talk about two stunning papers.

1. The First Paper:

Before coming to the paper itself, let me discuss its background a bit.

The paper has a topical relevance to the current news; its subject matter in a way refers to those most unfortunate accidents on the Mumbai-Pune Expressway.

In the recent times, there has been a spate of one particular kind of accidents on the Expressway: a tire of a car (or of a tempo or even a container truck) suddenly bursts open, the car gets out of control, crosses over the median and lands in a lane on the opposite side of the traffic, and ends up having a head-on collision.

People have been discussing the matters, building awareness campaigns, trying to analyze causes. For instance, things like: over-speeding as the cause of the tires bursting open; the absence of a tall enough separating wall on either sides; and human factors such as loss of concentration, drunk-driving, etc.

It’s a fact that on the Expressway most every car exceeds the speed limit (of 80 km/hr). You can see many amateur YouTube videos of passenger buses and ordinary passenger cars registering 120–140 km/hr as a routine. BTW, while watching these videos, notice that in India we drive on the left-hand side, and so, the fastest traffic is (usually!) to be found in the rightmost lane.

When a car unfortunately develops a burst tire and lands up on the other side of the road, not only is it the fastest moving car in the forward direction but also the cars in whose path it lands up after crossing over the median, also happen to be the fastest moving ones on the opposite side. Thus, collisions can possibly occur at, say, 100 + 120 = 220 km/hr (even after discounting the reduction in speed during the crossover, say, from 140 to 100 km/hr), making the situation far more grave. Contrast this situation with that of a car moving at 80 km/hr in the left-most lane. Even if the driver loses control and the car veers on the left-hand side, it still collides on to the protective railing (having quite a bit of capacity to absorb the impact energy) at only, say, 60–70 km/hr. The difference of speeds in a head-on collision is, thus, three times over. Little wonder the percentage of fatalities in this kind of accidents has been so high.

Then, there are issues like the lights of the cars on the opposite side shining directly in your eyes, because there are no shades, no dividing walls that are high enough.

Though I am not at all an expert here, another issue I am toying with is the development of a peculiarly wavy kind of surface due to the combined effects of rolling and vibrations,  especially as caused by heavy weight container trucks travelling on a road. A similar effect can be seen also on the railway tracks, where a wavy profile gets developed. However, in the case of the rails, due to the homogeneity of the material, it becomes a somewhat simpler problem to analyse within a mechanics framework. Even then, let me tell you, it’s a very complex problem to analyse. IIT Kanpur, for instance, has been conducing research sponsored by railways on this problem for many years, and one could say that they are still in a somewhat early modelling stage. Now, when it comes to a steel-reinforced cement-concrete road, it becomes, from the applied mechanics viewpoint, a much more complex problem to model. When the road is made of concrete, I am guessing, the peaks of the wavy profile could perhaps be a bit a too sharp, thereby introducing enhanced high frequency components in the car-road interactions, vibrations and what effectively are small-scale but high-speed impacts. Now, if the tires are worn out and are conventional (i.e. if they carry tubes and are not of the tubeless variety), and if a driver has been over-speeding for a long time thereby overheating them, especially on this kind of a more sharply wavy surface, there could be increased chances of the tube of the tire bursting open. (The tubeless tires have a tendency to more gracefully handle an event of a similar kind; the running out of the air or nitrogen is less intense and not so catastrophic.)

So, people have been talking of things like that, and, as you can see, I too have joined them, wondering aloud about what possibly could be the causes of such accidents.

Against this backdrop, here is a paper I found in a recent Web search:

N. P. Dharmadhikari, D. C. Meshram, S. D. Kulkarni, S. M. Hambarde, A. P. Rao, S. S. Pimplikar, A. G. Kharat, and P. T. Patil (2010) “Geopathic stress: a study to understand its nature using Light Interference Technique,” Current Science, vol. 98, no. 5, pp. 695–697

I must say that the paper is simply stunning.

Check it for yourself, and if you think otherwise, do drop me a line. I am sure, though, that you, too, would find it simply stunning. Especially noticeable is the level of rigor in both experimentation and the kind of data that have been gathered therein. Do have a look.

I gather that “Current Science” is “the [sic] leading interdisciplinary science journal from India […] started in 1932 by the then stalwarts of Indian science such as CV Raman, Birbal Sahni, Meghnad Saha, Martin Foster and S.S. Bhatnagar.” I also gather that it is indexed in “Web of Science, Current Contents, Geobase, Chemical Abstracts, IndMed and Scopus” [^].

BTW, the fact that the journal is indexed by so many services, esp. Scopus, goes to satisfy the requirements clarified by Professor Dr. T. P. Singh (Director, Symbiosis Institute of Technology). Recently, in an informal chat, he had mentioned to me that he cannot hire me because I do not have two journal papers already published, and so, I cannot become a PhD guide at Symbiosis from day 1. He then added that it is now his requirement to hire only those people as Full Professors who have two papers published in Scopus/similar indexed journals already, so that they can become PhD guides right from day 1 (his emphasis) and then, he continued in the same conversation, that during the last semester (when I too had applied at Symbiosis) he had hired two people as Full Professors but without PhD degrees so they cannot act as PhD guides from day 1 (here, the emphasis is mine). BTW, at the time that he hired those two people, I, too, had applied for the same post (having PhD already in hand), but was not called for interview at all. … Yes, sometimes Prof. Singh somehow has a charming way of reminding me of the late Jaspal Bhatti, an engineer.

Anyway, coming back to this paper, please note, engineering college professors from JSPM, MIT and Sinhagad have participated in it. There are as many as 7 authors to this paper. (Isn’t this fact stunning by itself?) They all, I suppose, have been approved by the UGC committees for recruitment of professors. In fact, this research began with a PhD granted by the University of Pune, in the first place.

BTW, if you are interested, here is another article which carries more or less the same level of rigor but it is definitely much more readable: [^].

2. The Second Paper:

Another paper, too, has a bit of a topical importance, at least as far as this blog is concerned, in the light of the recent discussion on the Fourier theory and all.

The paper in question is this:

G. K. Kharate, V. H. Patil and N. L. Bhale (2007) “Selection of mother wavelet for image compression on basis of nature of image,” Journal of Multimedia, vol. 2, no. 6, pp. 44–51

This paper, too, is stunning.

Why do I find it stunning?

In this paper, the authors experimentally study the extent to which different types of images can be compressed, using wavelets transforms. In particular, the extent of compression that is possible for a natural image (a woman’s portrait) is compared and contrasted with certain artificial images consisting of simple regular patterns such as horizontal/vertical/oblique lines. The conclusions reached by the authors consist of the following:

“The results demonstrate that for line-based image percentage of zeros is more for db1 as compared to other wavelets and more energy is retained. It shows that the loss of information is less hence the quality is better. Also for natural image percentage of zeros is more for db4 as compared to other wavelets that indicates more compression and the loss of information is less.

If the import of the conclusion is not immediately clear, here are a few passages from the main text of the paper that may help you in ascertaining the nature of the paper:

“It is observed that for periodic functions, Fourier analysis is ideal. However, wavelet transforms are not restricted to only the periodic function, but for any function, provided it is admissible.”

I was left wondering: How about the computational complexity of FFT vis-a-vis wavelets—particularly because the paper deals with compression efficiency? Reading further on, I then saw this:

“We suggest the minimum two steps for the wavelet decomposition. In first step, decompose the image by using any mother wavelet. On the basis of energy content of approximation and detail coefficients, the image is [first] classified as line image or natural image. Then select the proper mother wavelet for further decompositions as mentioned above.” [Emphasis in bold is mine.]

How about entropy? Well, this aspect is not as clearly indicated.

The Journal of Multimedia is, of course, an international journal published from the USA; it is indexed in many services, including Scopus [^].

Anyway, the first author is the current Dean of the Faculty of Engineering at the University of Pune. I mean to say: he was first awarded a doctorate in engineering by the University of Pune (under the guidance of A. A. Ghatol); he was selected by a UGC committee for a Full Professor’s position; he was also democratically elected by the professors coming from all the engineering colleges affiliated to the university of Pune; he now sits on the UGC interview committees that decide whether someone (e.g. me) is to be at all interviewed or not.

Do note that it’s not his the then PhD student (i.e. V. H. Patil), but Kharate himself who is the first author of this paper. Not so stunning. Not any more.

But, yes, the paper, overall, is stunning. Its conclusion simply has a way to force you to do a double take. Its very appearance is against what you expect intuitively—you never expect a paper of this kind to appear in a typical journal covering topics such as multimedia, image or signal processing, or wavelets analysis.

As you might know, the wavelets transform is similar to the Fourier transform though it has a categorically less computational cost involved with it as compared to FFT. Also, take a moment to think about what kind of connections could at all be possible between computational complexity, information entropy and the compression efficiency.

Now, reconsider the cited text and the conclusion of this paper.

Hopefully, now you can see why I say that you must read this journal paper, too.

* * * * *    * * * * *   * * * * *   

BTW, I have finally been rejected in the interviewing process at the JPBTI Phaneesh Murthy’s company, iGATE. The reason communicated to me (over the phone) by Asad Kadri this late evening had to do with the fact that I don’t have documentary evidence for having actually done jobs at three places: (i) Mukand (in 1983), and (ii) Thermax (in 1983–84); at both these places I was a trainee engineer; and then, (iii) QueriSoft (in 1995–96 i.e. about 15+ years ago), a company which is closed by now. BTW, I had arranged for a former colleague to informally verify the fact of my having worked with Querisoft, and anyway did have a xerox copy of the appointment letter for it. However, iGATE apparently thinks that their end-client, Rio Tinto, would find it unacceptable that all of my employment cannot be verified. In any case, they didn’t reply back to that email written by that past colleague of mine (who later on worked as a director in a large American MNC, and has now taken voluntary retirement from the IT field).

On my part, I find both Rio Tinto and iGATE to be very funny companies—tinted, as if they were, with a very dogged and out-of-the-world sort of a sense of humour. (They took three months, two interviews with the VP, three interviews with HR, and many telephonic calls and emails exchanges, in order to reach the conclusion that they cannot appoint me because I do not carry experience certificates for companies where I worked two decades ago and more.)

So, my job search continues, and something on the academia side might materialize—something that allows me at least, say, survival. I will keep you posted.

… But, don’t forget only to see those two papers, na!

* * * * *    * * * * *   * * * * *  

A Song I Like:
[For the above-mentioned reasons, the inclusion this time round of this section is, once again, purely on the whim.]

(Hindi) “machalti hui hawaa main chham chham…” (“hamare sang sang chalein Ganga ki leheren”)
Singers: Kishore Kumar [an unlikely choice, esp. given early ’60s], and Lata Mangeshkar
Music: Chitragupt
Lyrics: Majrooh Sultanpuri [?]

[E&OE]