This shows you the differences between two versions of the page.
docs:classes:flowfield [2009/02/17 08:36] gibson |
docs:classes:flowfield [2010/02/02 07:55] |
||
---|---|---|---|
Line 1: | Line 1: | ||
- | ====== FlowField ====== | ||
- | The FlowField class represents vector-valued fields of the form | ||
- | |||
- | <latex> | ||
- | {\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 | ||
- | </latex> | ||
- | |||
- | 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]. ((<latex> | ||
- | \bar{T}_{k_y}(y) = T_{k_y}\left(\frac{2}{b-a}(y - \frac{b+a}{2})\right)</latex>)) 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, +/-, +=, ∇, ∇<sup>2</sup>, norms, inner products, etc. | ||
- | * transforming back and forth between spectral coefficients <latex> \hat{u}_{k_x k_y k_z j}</latex> and gridpoint values <latex>u_j (x_{n_x}, y_{n_y}, z_{n_z})</latex> | ||
- | * 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 {{:librarycode:flowfield.h}}. | ||
- | |||
- | ===== Constructors / Initialization ===== | ||
- | |||
- | FlowFields are initialized with gridsize and cellsize parameters, read from disk, | ||
- | or assigned from computations. Examples: | ||
- | |||
- | <code c++> | ||
- | 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); | ||
- | </code> | ||
- | |||
- | ===== Algebraic and differential operators ===== | ||
- | |||
- | Assume f,g,h etc. are FlowField variables with compatible cell and grid sizes. Examples of | ||
- | possible operations | ||
- | |||
- | <code> | ||
- | <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); | ||
- | </code> | ||
- | |||
- | The latter functions are defined as | ||
- | |||
- | <latex> $ \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*} $ </latex> | ||
- | |||
- | 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 | ||
- | |||
- | <code> | ||
- | // A sequence that results in f = 0.5*(g + h); | ||
- | f = g; | ||
- | f += h; | ||
- | f *= 0.5; | ||
- | </code> | ||
- | | ||
- | 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|FieldSymmetry]] | ||
- | class. Briefly, the symmetries of 3D FlowFields can be parameterized as | ||
- | |||
- | <latex> $ \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*} $ </latex> | ||
- | |||
- | with the action of σ on a velocity field u as | ||
- | |||
- | <latex> | ||
- | \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) | ||
- | </latex> | ||
- | |||
- | A FieldSymmetry can be constructed and applied to a FlowField as follows | ||
- | |||
- | <code c++> | ||
- | 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 | ||
- | </code> | ||
- | |||
- | Or, the symmetric component of a field can be obtained by | ||
- | |||
- | <code c++> | ||
- | FlowField Pu = u; | ||
- | Pu += sigma(u); // Pu now equals u + sigma u | ||
- | Pu *= 0.5; // Pu now equals (u + sigma u)/2 | ||
- | </code> | ||
- | |||
- | For more examples of FlowField and FieldSymmetry usages, see | ||
- | [[:docs:classes:fieldsymmetry|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 {{docs:chflowguide.pdf|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 | ||
- | |||
- | <code> | ||
- | u.makePhysical(); // transform spectral coeffs to gridpt values | ||
- | u.makeSpectral(); // transform gridpt values to spectral coeffs | ||
- | </code> | ||
- | |||
- | 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. | ||
- | <code> | ||
- | 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; | ||
- | </code> | ||
- | |||
- | To print out its spectral coefficients, you need to make it | ||
- | Spectral first | ||
- | <code> | ||
- | 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; | ||
- | } | ||
- | } | ||
- | </code> | ||
- | 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%%. ((We could eliminate the issue entirely, | ||
- | but at the cost of run-time efficiency)). | ||
- | 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 | ||
- | <code> | ||
- | 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; | ||
- | } | ||
- | } | ||
- | </code> | ||
- | |||
- | 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 | ||
- | |||
- | <code c++> | ||
- | // 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) | ||
- | </code> | ||
- | |||
- | The complete set of such transform functions is | ||
- | |||
- | <code c++> | ||
- | 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); | ||
- | </code> | ||
- | 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 | ||
- | |||
- | <code c++> | ||
- | fieldstate xzstate = u.xzstate(); | ||
- | if (xzstate == Physical) | ||
- | .... | ||
- | |||
- | fieldstate ystate = u.ystate(); | ||
- | if (ystate == Spectral) | ||
- | .... | ||
- | </code> | ||
- | |||
- | The FlowField class has quite a few other member functions and operators. | ||
- | For a complete description, see the header file {{:librarycode:flowfield.h}}. |