docs:classes:flowfield

The FlowField class represents vector-valued fields of the form

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 and gridpoint values
- 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.

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);

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

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

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.

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

with the action of σ on a velocity field u as

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.

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 and the 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) .... 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.

We could eliminate the issue entirely,
but at the cost of run-time efficiency

docs/classes/flowfield.txt · Last modified: 2010/02/02 07:55 (external edit)