channelflow.org

====== Differences ====== This shows you the differences between two versions of the page.

--- docs:classes:flowfield [2009/02/16 10:20]
gibson created
+++ docs:classes:flowfield [2010/02/02 07:55] (current)
@@ Line 1: / Line 1: @@
 ====== FlowField ======
-The %%FlowField%% class represents scalar, vector, and tensor-valued fields of the form
+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}_{n_y}(y) \; e^{2 \pi i (k_x x/L_x + k_z z/L_z)} {\bf e}_j
+{\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>
+   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}}.

channelflow.org

User Tools

Site Tools

Page Tools