User Tools

Site Tools


docs:classes:flowfield

Differences

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

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
docs:classes:flowfield [2009/02/16 10:29]
gibson
docs:classes:flowfield [2009/02/17 08:36]
gibson
Line 1: Line 1:
 ====== FlowField ====== ====== FlowField ======
  
-The %%FlowField%% class represents vector-valued fields of the form+The FlowField class represents vector-valued fields of the form
  
 <​latex>​ <​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>​ </​latex>​
  
-and also scalar and tensor fields with appropriate changes in the dimensionality of the coefficients. The barred %%T%% +and also scalar and tensor fields with appropriate changes in the dimensionality of the coefficients. ​ 
-function is a Chebyshev polynomial scaled to fit the domain y ∈ [a,b].+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>​ <​latex>​
-\bar{T}_{k_y}(y) = T_{k_y}\left(y')y' = \frac{2){b-a}\left(y- \frac{b+a}{2}\right)+  ​\sigma [u, v, w](x,y,z) = (s_x us_y v, s_z w)(s_x x + a_x L_x, s_y y, s_z z a_z L_z)
 </​latex>​ </​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}}.
docs/classes/flowfield.txt · Last modified: 2010/02/02 07:55 (external edit)