Hi,
I was wondering if it’s possible to adapt Channelflow for purely 2D configurations? It doesn’t seem that Nz=0 is supported.
Thanks for your help,
Vilda
Hi, Vilda. There is no simple way to eliminate the z coordinate or use a grid of length Nz=1 in channelflow. You can however simulate a flow with no z variation by using Nz=6, a z-constant initial condition, and an Lz value small enough to suppress instability in z variation in the velocity field (but not so small as to cause CFL problems), or make its growth slow. People find this works well in practice.
From a mathematical perspective it seems like Nz=1 ought to work. But comptationally the data layout for the FFTs is really complicated, and Nz=1 turns out not to be a special case of the general structure. Just for example, even in 1d, a real-to-complex Fourier transform in FFTW requires at least two gridpoints, so the two real numbers that tranform to a single complex number occupy the same amount of memory. There are further complications due to the three-dimensionality of the data set and the 2/3 padding for dealiased FFTs. It turns out Nz=6 is the smallest Nz that works. (Maybe Nz=2 or 4 would work if you turn off dealiasing. I don’t know; I’ve never tried.)
Just start with a z-constant field, Nz=6, and Lz set so that dz = Lz/Nz is smaller than the corresponding dx=Lx/Nx by a factor of two or four. Every now and then during our time integration compute a norm of u - u’, where u is u’ shifted in z by exactly dz. If that norm gets above some small number, say 1e-10, replace u with its mean over z and restart the time integration. You won’t have to do this every time step, probably more like every 10 or 100 outer time units.
Dear John,
Thanks for your quick reply.
I’m glad to hear there is a solution that works well in practice. I previously tried using Nz=4 and a small Lz but I ran into problems very quickly, without realising Nz=6 would be the minimum for a 3D geometry. I assume it will be fine not to use de-aliasing in z direction in this case? Thanks very much for the practical advice, I will try it and see how it goes.
Best wishes,
Vilda Markeviciute
Hi again,
I’ve tried implementing the two-dimensionality in your suggested way which seems to work, but I’m now experiencing issues when running the code in parallel on my faculty’s HPC system. In particular, while silencing the unwanted z modes, I’m using the following construct:
Complex zero = 0.0 + 0.0*I;
for (int i=0; i<fields[0].Nd(); ++i){
for (int my=0; my<fields[0].My(); ++my){
for (lint mx=fields[0].mxlocmin(); mx<fields[0].mxlocmin()+fields[0].Mxloc(); ++mx){
for (lint mz=fields[0].mzlocmin(); mz<fields[0].mzlocmin()+fields[0].Mzloc(); ++mz){
if (mz>0){fields[0].cmplx(mx,my,mz,i) = zero;}}}}}
which produces the following error:
channel: /nfs/software/spack/linux-sles12-x86_64/gcc-7.3.0/channelflow-develop-sloyevvcbaztaiv66q2oslnxqdwg6ybb/include/channelflow/flowfield.h:390: int chflow::FlowField::complex_flatten(int, int, int, int) const: Assertion `mz >= mzlocmin_ && mz < Mzloc_ + mzlocmin_’ failed.
I tried my best using the local modes but it seems that something has gone wrong. Is there something obvious that my code is missing?
Kind regards,
Vilda Markeviciute
Hi, Vilda. If you put your code blocks in triple backticks, they’ll format nicely. (https://discourse.julialang.org/t/psa-how-to-quote-code-with-backticks/7530)
Complex zero = 0.0 + 0.0*I;
for (int i=0; i<fields[0].Nd(); ++i){
for (int my=0; my<fields[0].My(); ++my){
for (lint mx=fields[0].mxlocmin(); mx<fields[0].mxlocmin()+fields[0].Mxloc(); ++mx){
for (lint mz=fields[0].mzlocmin(); mz<fields[0].mzlocmin()+fields[0].Mzloc(); ++mz){
if (mz>0){fields[0].cmplx(mx,my,mz,i) = zero;}}}}}
produces error
channel: /nfs/software/spack/linux-sles12-x86_64/gcc-7.3.0/channelflow-develop-sloyevvcbaztaiv66q2oslnxqdwg6ybb/include/channelflow/flowfield.h:390: int chflow::FlowField::complex_flatten(int, int, int, int) const: Assertion `mz >= mzlocmin_ && mz < Mzloc_ + mzlocmin_’ failed.
I’m stumped as to why your loop triggers the assertion error. It certainly looks like your loop starts at mz=mzlocmin and ends with mz < mzlocmin + Mzloc. I even checked the inlined functions FlowField::mzlocmin() and FlowField::Mzloc() and verified they’re correctly returning the FlowField private member variables mzlocmin_ and Mzloc_, which are being used in the assertion. So I really don’t get how it’s being triggered. @florian.reetz? @mirko.farano?
A minor suggestion: pull the function calls out of the loop bounds specifications as follows
int Nd = fields[0].Nd();
int My = fields[0].My();
int Mxloc = fields[0].Mxloc()
int mxlocmin=fields[0].mxlocmin();
int Mzloc = fields[0].Mzloc()
int mzlocmin=fields[0].mzlocmin();
for (int i=0; i<Nd; ++i) {
for (int my=0; my<My; ++my) {
for (lint mx=mxlocmin; mx<mxlocmin+Mxloc; ++mx) {
for (lint mz=mzlocmin; mz<mzlocmin+Mzloc; ++mz) {
if (mz>0){fields[0].cmplx(mx,my,mz,i) = zero;}}}}}
Those are inline functions, and compilers might be capable of hoisting them out of the loop (maybe not, as u.cmplx is a non-const function). But even if so it makes the code more readable.
HI Jhon,
Yes, it looks very weird also to me. From the loop you are referring I don’t see any mistake.
You will need to debug it (if you can reproduce the error with very few number of core such that you can attach a debugger for each proc). Otherwise as first attempt an std::cout will do the job (just output the value of mz and the rank to see where it is failing).
I hope you can easily figure it out!
Best,
Mirko