I write regarding an issue that popped up in my usage of both the old and the new channelflow. My modified files are attached to this post. The focus of my research is the Poiseuille flow.
The first changes I made were to ensure that I could work in an averaged moving frame (Navier Stokes are Galilean invariant and hence I could simplify my work while obtaining the same physics). This was somehow not achievable by simply imposing a bulkv constraint, while requesting a parabolic base flow (or alternatively, picking out appropriate boundary values) and so I had to create my own flags in simulateflow to achieve it. I write this because on inspecting the Ubase.asc files, after imposing the above conditions, I saw that the base profile was somehow either obeying the bulkv constraint or the parabolic profile, but never both. My modifications took care of that and the Ubase was thereafter found to be as expected.
For my purposes, the objects of interest were not the actual velocities of the flow, but the deviation of the velocity from the expected base parabolic flow. Hereon, I shall refer to them as velocity values.
The actual issue popped up when I tried to output values of velocity fields along certain directions to what were called the history files. These are simply files that contain velocity field values along the z axis for some fixed (x, y) and then along the y axis for some fixed point (x, z) (value of x doesn’t have to be the same in both). These velocity field values are stored for each time at {1, 2, 3 … T1}.
It was noticed that these history files have a peculiar property of having 0 deviation from the parabolic flow at some collection of (x, y, z) values for all times. The troubling part of this observation was that this is reproduced by CF1.5 and CF2 for any initially chosen Re or random initial velocity field. As I am merely outputting values calculated by the program and not doing any calculations myself, I take this to mean that CF1.5 and CF2 are telling me that no matter the initial physical variables or velocity profiles picked, the same points observe 0 deviation from the parabolic profile for all times. Unless I am somehow always picking out flows that somehow have the same fixed points, that is in contradiction to observed reality.
It is possible that my changes somehow caused the above paradoxical behavior, even though as far as I am aware, I am merely outputting values that CF1.5 and CF2 produce. Hence, I am attaching all my changes and data from the most recent of runs to this post. As I rely on the results produced by channelflow for my research, it is imperative I can rely upon its output. Thus, I would be grateful for a look into the matter whenever possible.
Best Link to files
Can you be more specific --mathematically more specific-- about the modifications you made to the base flow conditions? E.g. specify the base flow you are trying to achieve (is it U(y) = (1 - y^2) - 2/3 = 1/3 - y^2 ?). Then specify what base flow conditions (wall speeds, bulk velocity constraint, or pressure constraint) you tried to use to achieve that base flow, and how you modified the channelflow baseflow flags in order achieve it.
And can you be more specific about what (x,y,z) point has zero (deviation from laminar) velocity for all time? Just specify a single point where the anomaly happens, what channelflow code you’re using to extract that data from the velocity field. I realize I could probably figure this out from analyzing the source code you sent, but it would save me a lot of time if you just told me.
What do you mean by CF1.5 and CF2?
Lastly (and I don’t mean to sound grumpy or discourage you in any way at all), it would really be helpful if you cleaned up the source code of extraneous commented-out code and used conventional indentation so that it wasn’t so hard to read. Start with simulateflow.cpp, modify it as minimally as possible to reproduce the proble, and mark your modifications with comments // MODIFIED HERE or something like that. That would really help.
Also, quote the complete command-line command you used to start the simulation that gives erroneous results, including how you constructed the initial velocity field.
looking forward to helping further, with further information!
Ok, I’ve dug into your code a bit. I can see you did mark you r modifications with comments, and by the lower/upper wall speed of -0.666666 I can see you are trying to get U(y) = 1/3 - y^2.
I think the likely issue here is that you’re using the gridpoint evaluation functions u(nx,ny,nz,i) on the flow fields during the timestepping loops when these flowfields are surely stored as spectral coefficients instead of gridpoint values. The internal data representation of a flowfield is either girdpoint values of the velocity or spectral coefficients of the velocity field. At any given time a Flowfield is physical (gridpoint) or spectal (coefficient). Flowfields transform back and forth between these representations using FFT and IFFT. It’s an error to ask for a gridpoint value when the Flowfield is in spectral representation (and vice versa). But checks for such errors happen only when channelflow is compiled in debugging mode.
Long story short: if you want to get a gridpoint value inside the timestepping loop, you need to do it like this
If you want to output multiple gridpoint values, put the for loop between the calls to makePhysical() and makeSpectral(). They’re relatively expensive.
This post is a follow up to the suggestions of Prof Gibson. To make matters easy to follow for everyone, I will first address the points raised already. The following contains the various files: link
CF1.5 and CF2 are acronyms for Channelflow_1.5 and Channelflow_2.
The professor rightly found that the base profile I was trying to achieve is the U(y) = 1/3 - y^2.
The commands used to test the code was the following:
randomfield -sd 232 -m 0.30 -Lx 5 -Lz 20 -Nx 32 -Ny 33 -Nz 128 u0
simulateflow -mc bulkv -R 600 -T 109 -dtmin 0.0005 -dt 0.001 -Ubulk 0 -theta 0.4188790204786391 u0.nc
The professor’s advice on structuring the for loops inside a umakePhysical() and umakeSpectral() pair does remove the occurrence of fixed point of velocity (no 0s are fixed anymore in history file). This holds only for CF2. (Case for CF1.5 is dealt separately, but in summary, there the problem is not resolved at all)
In CF2, this however still results in a problem, namely that files being produced are not producing consistent results. There are 3 cases:
Using default CF2: The energy_no_mods.asc shows behavior that is expected.
Using default CF2, but with the addition to allow printing of history files only: Even though this step should not change the development of the flow, as all the addition does is print out values, compared to case 1, it does produce different values completely in the energy file.
Using default CF2, but with addition to allow the desired base profile only: Quantitative difference is expected in the files in this case, due to the fact that we are in a completely different frame of reference this time. However, surprisingly, just changing to this frame apparently caused such a drastic difference that the code ascertained that the energy had blown up to infinity after only 3 steps.
Using default CF2, with changes from point 2 and 3 together: the energy blows up even faster than case 3 and the code basically doesn’t run (produces Nan at the very first step)
Making the suggested changes in couette.cpp of CF1.5 seems to have had no effect at all and still results in fixed points of velocity in the same gridpoint (0 6 1) at all times.
As an addendum to the problems observed above, in CF1.5, employing the Professor’s method of nesting the for loops in couette.cpp in umakePhysical and umakeSpectral produces a history file with different data than the one produced by code not nesting the loop in the two commands. The only exception is the fixed point, which doesn’t change at all. The fixed gridpoint is (0 6 1)