User Tools

Site Tools


docs:classes:flowfield

This is an old revision of the document!


FlowField

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, differential, and symmetry 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

For a complete description of FlowField functionality, see the header file flowfield.h.

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 f(g);                                // make a copy of g
   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

<code>
   f *= 2.7;                // f = 2.7*f
   f += g;                  // f = f + g
   f -= g;                  // f = f - g

   f = xdiff(g);            // f_i  = d g_i/dx
   f = ydiff(g);            // f_i  = d g_i/dy
   f = zdiff(g);            // f_i  = d g_i/dz
   f = diff(g, j, n);       // f_i  = d^n g_i/dx_j
   f = diff(g, j, n);       // f_i  = d^n g_i/dx_j
   f = grad(g);             // f_ij = dg_i/dx_j
   f = curl(g); 
   f = lapl(g);
   f = div(g);
   f = cross(g,h);      
  
   xdiff(g, dgdx);          // same as dgdx = xdiff(g), but often more efficient
   curl(g, curl_g);         // ditto
   lapl(g, lapl_g);         // ditto
   ...

   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*} $

Note that expressions such as f = g + h or f = 0.5*(g + h) are not allowed on FlowFields, since these would generate temporary FlowField variables g + h and 0.5*(g+h) during expression evaluation. Instead, use sequences such as

  // A sequence that results in f = 0.5*(g + h);
  f = g;
  f += h;
  f *= 0.5;

As C++ objects, FlowFields are huge monsters. It is best to minimize the amount of construction, copying, assignment of FlowFields by reusing temporaries and figuring out the minimal sequence of operations to get the desired result.

Symmetry operations

The symmetry group of FlowFields is represented by the FieldSymmetry class. Briefly, the symmetries of 3D FlowFields can be parameterized as

 $ \begin{align*}
  \sigma &= (s_x, s_y, s_x, a_x, a_z, s)\\
  s_x, s_y, s_z, s &= \pm 1\\
  a_x, a_z &\in [-0.5, 0.5)
\end{align*} $

with the action of σ on a velocity field u as


  \sigma [u, v, w](x,y,z) = s (s_x u, s_y v, s_z w)(s_x x + a_x L_x, s_y y, s_z z + a_z L_z)

A FieldSymmetry can be constructed and applied to a FlowField as follows

 FieldSymmetry sigma(sx, sy, sz, ax, az, s);  // construct sigma = (sx,sy,sz,ax,az,s)
 FlowField sigma_u = sigma(u);                // apply symmetry sigma to u

Or, the symmetric component of a field can be obtained by

 FlowField Pu = u;  
 Pu += sigma(u);    // Pu now equals u + sigma u
 Pu *= 0.5;         // Pu now equals (u + sigma u)/2

For more examples of FlowField and FieldSymmetry usages, see the FieldSymmetry documentation.

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 need to 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 = ' ';
   u.makePhysical();
   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

   u.makeSpectral();
   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 = u.kz(mz);
              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 = u.mx(kx) 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

   u.makeSpectral();
   for (int i=0; i<u.Nd(); ++i)
     for (int kx=u.kxmin(); kx<u.kxmax(); ++kx) {
        int mx = u.mx(kx);
        for (int my=0; my<u.My(); ++my)
           for (int kz=u.kzmin(); kz<u.kzmax(); ++kz) {
              int mz = u.mz(kz);
              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.

You can also perform the $x,z$ and the $y$ transforms independently. For example, if u is representing pure gridpoint values you could do this

  // get a gridpoint value
  Real u_nxnynzi = u(nx,ny,nz,i); 
 
  u.makeSpectral_xz();
 
  // get kx,kz Fourier coefficient at ny-th gridpoint in y
  Complex ukxnykzi = u.cmplx(u.mx(kx), ny, u.mz(kx), i) 

The complete set of such transform functions is

  u.makeSpectral();    // to pure spectral coeffs
  u.makePhysical();    // to pure gridpoint values
  u.makeSpectral_xz(); // to spectral coeffs in x,z
  u.makeSpectral_y();  // to spectral coeffs in y
  u.makePhysical_xz(); // to gridpoint values in x,z
  u.makePhysical_y();  // to gridpoint values in y
  u.makeState(Physical, Spectral); // to x,z Physical and y Spectral
  u.makeState(..., ...); // and the other three combinations of (Physical,Spectral);

The FlowField keeps track of its spectral/physical states in x,z and y performs the desired transform only if it's in the opposite state. You can query the state of a FlowField like this

  fieldstate xzstate = u.xzstate();
  if (xzstate == Physical)
    ....
 
2  fieldstate ystate = u.ystate();
  if (ystate == Spectral)
    ....

The FlowField class has quite a few other member functions and operators. For a complete description, see the header file flowfield.h.

2)
We could eliminate the issue entirely, but at the cost of run-time efficiency
docs/classes/flowfield.1234888216.txt.gz · Last modified: 2009/02/17 08:30 by gibson