User Tools

Site Tools


This is an old revision of the document!


The FlowField class represents vector-valued fields of the form

{\bf u}({\bf x}) &= \sum_{k_x,k_y,k_z,j}  \hat{u}_{k_x k_y k_z j}\bar{T}_{k_y}(y) \; e^{2 \pi i (k_x x/L_x + k_z z/L_z)} {\bf e}_j

and also scalar and tensor fields with appropriate changes in the dimensionality of the coefficients. The barred T function is a Chebyshev polynomial scaled to fit the domain y ∈ [a,b]. 1) The spatial domain of a FlowField is Ω = [0,Lx] x [a,b] x [0,Lz], with periodicity in x and z.

In channelflow programming, fields such as velocity, pressure, stress tensors, vorticity, etc. are stored as variables of type FlowField. The main functionality of the FlowField class is

  • algebraic and differential operations, +/-, +=, ∇, ∇2, norms, inner products, etc.
  • transforming back and forth between spectral coefficients  \hat{u}_{k_x k_y k_z j} and gridpoint values u_j (x_{n_x}, y_{n_y}, z_{n_z})
  • serving as input to DNS algorithms, which map velocity fields forward in time: u(x,t) → u(x, t+Δt)
  • setting and accessing scpetral coefficients and gridpoint values
  • reading and writing to disk

Constructors / Initialization

FlowFields are initialized with gridsize and cellsize parameters, read from disk, or assigned from computations. Examples:

   FlowField f;                                   // null value, 0-d field on 0x0x0 grid 
   FlowField u(Nx, Ny, Nz, Nd, Lx, Lz, a, b);     // Nd-dim field on Nx x Ny x Nz grid, [0,Lx]x[a,b]x[0,Lz]
   FlowField g(Nx, Ny, Nz, Nd, 2, Lx, Lz, a, b);  // Nd-dim 2-tensor  
   FlowField h("h");                              // read from file "h.ff"
   FlowField omega = curl(u);

Algebraic and differential operators

Assume f,g,h etc. are FlowField variables with compatible cell and grid sizes. Examples of possible operations

   f += g;                  // f = f + g
   f = curl(g); 
   f = lapl(g);
   f = div(g);
   f = diff(g, j, n);       // f_i  = d^n g_i /dx_j
   f = grad(g);             // f_ij = dg_i / dx_j
   f = cross(g,h);      
   f *= 2.7;                // f = 2.7*f
   Real c = L2IP(f,g);      // L2 inner product of f,g
   Real n = L2Norm(u);
   Real D = dissipation(u);
   Real E = energy(u); 
   Real I = wallshear(u);

The latter functions are defined as

 $ \begin{align*}
L2IP(f,g) &= \frac{1}{L_x L_y L_z} \int_{\Omega} {\bf f} \cdot {\bf g} \,\, d{\bf x} \\
L2Norm(u) &= \left(\frac{1}{L_x L_y L_z} \int_{\Omega} \|{\bf u}\|^2\,\, d{\bf x} \right)^{1/2}\\
     E(u) &= \frac{1}{2 L_x L_y L_z} \int_{\Omega} \|{\bf u}\|^2\,\, d{\bf x}  \\
     D(u) &= \frac{1}{L_x L_y L_z} \int_{\Omega} \|\nabla \times {\bf u}\|^2 \,\, d{\bf x} \\
     I(u) &= \frac{1}{L_x L_z} \int_{y=a,b} \frac{\partial u}{\partial y} \, dx dz
\end{align*} $

Transforms and data access

FlowField transforms are a complicated subject –there are transforms in x,y, and z and implicit symmetries in complex spectral coefficients due to the real-valuedness of the field, for instance. This section outlines the bare essentials of transforms and data access methods. For further details see the Channelflow User Guide.

The FlowField class has a large data array that contains the spectral coefficients of the expansion listed at the top of this page. Most operations on FlowFields are calculated in terms of these spectral coefficients. But sometimes we need to know the value of the field at gridpoints. Rather than directly evaluating the expansion sum for given values of (x,y,z) (which would be very slow), we use fast Fourier transforms to transform the array of spectral coefficients into an array of gridpoint values. The main FlowField class transform methods are

   u.makePhysical(); // transform spectral coeffs to gridpt values
   u.makeSpectral(); // transform gridpt values to spectral coeffs

Because the transforms change the meaning of the FlowField's internal data array, you make sure the FlowField is in the proper state before trying to access either its spectral coefficients or its gridpoint values.

For example, to print out the entire set of gridpoint values of a FlowField, you would want to make the field Physical first.

   char s = ' ';
   for (int i=0; i<u.Nd(); ++i)
     for (int nx=0; nx<u.Nx(); ++nx)
        for (int ny=0; ny<u.Ny(); ++ny)
           for (int nz=0; nz<u.Nz(); ++nz)
              cout << nx <<s<< ny <<s<< nz <<s<< u(nx,ny,nz,i) << endl;

To print out its spectral coefficients, you need to make it Spectral first

   for (int i=0; i<u.Nd(); ++i)
     for (int mx=0; mx<u.Mx(); ++mx) {
        int kx = u.kx(mx);
        for (int my=0; my<u.My(); ++my)
           for (int mz=0; mz<u.Mz(); ++mz) {
              int kz =;
              cout << kx <<s<< my <<s<< kz <<s<< u.cmplx(mx,my,mz,i) << endl;

Note that the loop variables for mx,mz are not the same as the wavenumbers kx,kz. That's because the Fourier transforms leave the data in a peculiar order. Channelflow tries to ease the pain of this difference by providing functions int kx = u.kx(mx) and int mx = that translate between data ordering mx and Fourier wavenumbers kx, and similarly for mz,kz. 2). Note also that the data access method for spectral coefficients is the complex-valued u.cmplx(mx,my,mz,i), compared to the real-valued gridpoint access method u(nx,ny,nz,i) and that the bounds of the indexing variables are different.

If you really want to loop in kx,kz order (at the slight cost in efficiency), do this

   for (int i=0; i<u.Nd(); ++i)
     for (int kx=u.kxmin(); kx<u.kxmax(); ++kx) {
        int mx =;
        for (int my=0; my<u.My(); ++my)
           for (int kz=u.kzmin(); kz<u.kzmax(); ++kz) {
              int mz =;
              cout << kx <<s<< my <<s<< kz <<s<< u.cmplx(mx,my,mz,i) << endl;

But in general it is better to use built-in FlowField operations such as curl and diff than to loop over the data arrays, if you can.

We could eliminate the issue entirely, but at the cost of run-time efficiency
docs/classes/flowfield.1234830157.txt.gz · Last modified: 2009/02/16 16:22 by wikiadmin