Fields

Fields

\[\def\ddt#1{\frac{\mathrm{d}#1}{\mathrm{d}t}} \renewcommand{\vec}{\boldsymbol} \newcommand{\uvec}[1]{\vec{\hat{#1}}} \newcommand{\utangent}{\uvec{\tau}} \newcommand{\unormal}{\uvec{n}} \renewcommand{\d}{\,\mathrm{d}}\]

In whirl, field data, such as velocity, vorticity and pressure, are stored on a staggered uniform grid. Such a grid is divided into cells, with edges (which, on a two-dimensional grid, are the same as faces) and nodes (cell centers). Nodes hold scalar-valued data. Edges, on the other hand, hold the components of vector-valued data that are normal to the respective edges; one component lies on the vertical edges, while the other is on the horizontal edges.

Furthermore, there are two different cell types: primal and dual. On the physical grid, these cell types are offset with respect to each other by half a cell spacing in each direction. In other words, the four corners of the primal (resp. dual) cell are the nodes of four dual (resp. primary) cells.

Thus, on a two-dimensional staggered grid, there are four distinct vector spaces, associated with where the data are held on the grid:

In whirl, these are each distinct data types. Furthermore, the relationships between these types are defined by an underlying grid shared by all. By convention, this grid is defined by the number of dual cells NX and NY in each direction; we will often refer to it as the dual grid. For example, Nodes{Dual,NX,NY} is the type for dual node data on this grid; Edges{Primal,NX,NY} is the type for edge data on the primal cells within this same NX by NY dual grid. Note that, even though this latter type is parameterized by NX and NY, these values do not correspond to the number of primal edges in each direction on this dual grid. These values always correspond to the number of dual cells on the grid, for any data type. This makes it clear the grid is shared by all data.

Setting up field data

Let's see an example of creating a blank set of dual node data and filling it with something:

julia> w = Nodes(Dual,(5,4))
Whirl.Fields.Nodes{Whirl.Fields.Dual,5,4} data
Printing in grid orientation (lower left is (1,1)):
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> w .= reshape(1:20,5,4)
Whirl.Fields.Nodes{Whirl.Fields.Dual,5,4} data
Printing in grid orientation (lower left is (1,1)):
 16.0  17.0  18.0  19.0  20.0
 11.0  12.0  13.0  14.0  15.0
  6.0   7.0   8.0   9.0  10.0
  1.0   2.0   3.0   4.0   5.0

Other data types on the same grid can be set up in similar fashion. To ensure that they have a size that is consistent with the dual node data w, we can use this in place of the size:

julia> q = Edges(Primal,w);

julia> q.u[2,3] = 1;

julia> q
Whirl.Fields.Edges{Whirl.Fields.Primal,5,4} data
u (in grid orientation):
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

Field differencing operations

Field operations transform one data type to another. Some of these are differencing operations, analogous to differential counterparts in continuum calculus: curl, divergence, and gradient. For example, a curl operation can act upon dual nodal data (like streamfunction) and return primal edge data (i.e. velocity); a divergence operation acts on edge data (primal or dual) and returns nodal data of the same cell type. Note that these operations are mimetic: they maintain some of the same properties as the continuous counterparts. For example, the divergence of the curl of any dual nodal data is exactly zero. The curl of the gradient of primal nodal data is also zero.

Let's take the curl of the dual nodal data we constructed:

julia> curl(w)
Whirl.Fields.Edges{Whirl.Fields.Primal,5,4} data
u (in grid orientation):
 5.0  5.0  5.0  5.0  5.0
 5.0  5.0  5.0  5.0  5.0
 5.0  5.0  5.0  5.0  5.0
v (in grid orientation):
 -1.0  -1.0  -1.0  -1.0
 -1.0  -1.0  -1.0  -1.0
 -1.0  -1.0  -1.0  -1.0
 -1.0  -1.0  -1.0  -1.0

We could also make this a little more cute by giving the curl operator a symbol and then acting upon the data as though it were a matrix-vector operation:

julia> C = Curl()
Whirl.Fields.Curl()

julia> C*w
Whirl.Fields.Edges{Whirl.Fields.Primal,5,4} data
u (in grid orientation):
 5.0  5.0  5.0  5.0  5.0
 5.0  5.0  5.0  5.0  5.0
 5.0  5.0  5.0  5.0  5.0
v (in grid orientation):
 -1.0  -1.0  -1.0  -1.0
 -1.0  -1.0  -1.0  -1.0
 -1.0  -1.0  -1.0  -1.0
 -1.0  -1.0  -1.0  -1.0

Note that C is not actually a matrix. Rather, it is simply another name for the curl operator, and * is defined in this context to apply curl to whatever is to the right of it. The other operators have similar constructs.

Suppose we wish to apply the curl operation over and over. The curl() function allocates memory for the result whenever it is used; this would become expensive if it is done often. So it makes sense to preallocate space for this result and use the curl!() function, which simply fills in the elements:

julia> q = Edges(Primal,w)
Whirl.Fields.Edges{Whirl.Fields.Primal,5,4} data
u (in grid orientation):
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> curl!(q,w)
Whirl.Fields.Edges{Whirl.Fields.Primal,5,4} data
u (in grid orientation):
 5.0  5.0  5.0  5.0  5.0
 5.0  5.0  5.0  5.0  5.0
 5.0  5.0  5.0  5.0  5.0
v (in grid orientation):
 -1.0  -1.0  -1.0  -1.0
 -1.0  -1.0  -1.0  -1.0
 -1.0  -1.0  -1.0  -1.0
 -1.0  -1.0  -1.0  -1.0

Note that we used a convenience function for setting up primal edge data q of a size that corresponds with w.

Let's check that divergence of the curl is indeed zero:

julia> D = Divergence()
Whirl.Fields.Divergence()

julia> D*(C*w)
Whirl.Fields.Nodes{Whirl.Fields.Primal,5,4} data
Printing in grid orientation (lower left is (1,1)):
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

The Laplacian and its inverse

Whirl also makes heavy use of the discrete Laplacian operator, $L$. This mimics the continuous operator, $\nabla^2$, and acts upon data of any type. Let's apply this to the original data:

julia> laplacian(w)
Whirl.Fields.Nodes{Whirl.Fields.Dual,5,4} data
Printing in grid orientation (lower left is (1,1)):
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

As with the other operators, we can also construct a shorthand of the discrete Laplacian operator,

julia> L = plan_laplacian(size(w))
Discrete Laplacian on a (nx = 5, ny = 4) grid with spacing 1.0

julia> L*w
Whirl.Fields.Nodes{Whirl.Fields.Dual,5,4} data
Printing in grid orientation (lower left is (1,1)):
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

An important part of whirl is the inverse of this operator. That is, we need the ability to solve the discrete Poisson system

\[Ls = w\]

for $s$, for given data $w$. We achieve this in whirl with the lattice Green's function. To outfit the operator with its inverse, we simply set the optional flag:

julia> L = plan_laplacian(size(w),with_inverse=true)
Discrete Laplacian (and inverse) on a (nx = 5, ny = 4) grid with spacing 1.0

Then, the Poisson system is solved with the backslash (\),

julia> s = L\w
Whirl.Fields.Nodes{Whirl.Fields.Dual,5,4} data
Printing in grid orientation (lower left is (1,1)):
 24.082   15.9349  12.2076  13.0464  19.5282
 23.5931  15.0969  11.0184  11.7015  18.2956
 29.2015  21.841   18.0677  18.4457  23.9041
 37.5698  31.9978  28.9659  29.1093  33.016

julia> L*s
Whirl.Fields.Nodes{Whirl.Fields.Dual,5,4} data
Printing in grid orientation (lower left is (1,1)):
 0.0   0.0   0.0   0.0  0.0
 0.0  12.0  13.0  14.0  0.0
 0.0   7.0   8.0   9.0  0.0
 0.0   0.0   0.0   0.0  0.0

It should be observed that the cells on the perimeter have not recovered the original values of w. These are the ghost cells, and the Laplacian operation does not apply to these.

It is also important to note that, although it looks as though we've constructed a matrix L and performed various matrix-vector operations with it, this is not actually the case. In fact, the \ operation associated with L is significantly faster than a matrix inversion. Internally, it carries out a fast convolution between the data in w and the lattice Green's function, via fast Fourier transform. The lattice Green's function (LGF) table is pre-computed and pre-transformed in the original construction of L. (In fact, because this table is not dependent on the size of the grid, it is actually computed once for all time and stored in a file; subsequent applications of it just load it in and use the portion of it necessary for a certain grid.)

The lattice Green's function has the advantage that it is independent of the grid size. Let's solve the Poisson system when $w$ is a unit field, i.e. a field of zeros, except for a single $1$ entry at one node. The solution $s$ represents the influence of this point on all nodes. To see that the LGF does not depend on the grid size, let's use a grid that is long and skinny and plot the solution on it

julia> w = Nodes(Dual,(50,10));

julia> w[20,5] = 1.0
1.0

julia> L = plan_laplacian(w,with_inverse=true)
Discrete Laplacian (and inverse) on a (nx = 50, ny = 10) grid with spacing 1.0

julia> plot(L\w)
Plot{Plots.PyPlotBackend() n=1}

The influence is not affected by the narrow grid dimensions.

The integrating factor

An operator related to the lattice Green's function is the integrating factor. Suppose we have the system of ODEs

\[\ddt u = L u + f(u,t), \quad u(0) = u_0,\]

where $L$ is the discrete Laplacian (on an infinite uniform grid), and $u$ are nodal data (and $f$ is a nodal-valued function acting on this nodal data). The exact solution of this problem is

\[u(t) = E(t)u_0 + \int_0^t E(t-\tau) f(u(\tau),\tau)\,\mathrm{d}\tau,\]

where $E(t)$ is the integrating factor (or matrix exponential) for the system. The easiest way to understand the role of $E(t)$ is to consider its behavior when $f$ is zero and $u_0$ contains a field of zeros except for a single $1$ entry at one cell. Let's set up this initial data:

julia> u0 = Nodes(Dual,(100,100));

julia> u0[40,50] = 1.0
1.0

julia> plot(u0)
Plot{Plots.PyPlotBackend() n=1}

Then, $E(t)u_0$ diffuses this initial unit perturbation in each direction. Here, we apply it with $t = 5$:

julia> E = plan_intfact(5,u0)
Integrating factor with parameter 5 on a (nx = 100, ny = 100) grid

julia> plot(E*u0)
Plot{Plots.PyPlotBackend() n=1}

Note that $E(0) = I$, where $I$ is the identity. Also, the integrating factor has the useful property that $E(t+\tau) = E(t)E(\tau)$. From these properties, it follows that $E^{-1}(t) = E(-t)$. Let's suppose we wish to advance $u$ from time $t = \tau-h$ to time $t = \tau$. Then we can define an auxiliary quantity, $v(t;\tau) = E(\tau-t)u(t)$, and this new quantity satisfies the modified set of ODEs

\[\ddt v = E(\tau-t) f\left[ E(t-\tau) v(t;\tau),t\right],\quad v(\tau-h;\tau) = E(h)u(\tau-h)\]

The result of integrating this set of ODEs to $t = \tau$ is $v(\tau;\tau) = u(\tau)$. In other words, the integrating factor allows us to solve a somewhat reduced set of ODEs.

Other field operations

Other field operations shift the data, by local averaging, from one data type to another. These operations are all called shift!, and they require that the target data be preallocated. For example, to shift dual node data to the dual edges,

julia> w = Nodes(Dual,(5,4));

julia> w .= reshape(1:20,5,4)
Whirl.Fields.Nodes{Whirl.Fields.Dual,5,4} data
Printing in grid orientation (lower left is (1,1)):
 16.0  17.0  18.0  19.0  20.0
 11.0  12.0  13.0  14.0  15.0
  6.0   7.0   8.0   9.0  10.0
  1.0   2.0   3.0   4.0   5.0

julia> Ww = Edges(Dual,w);

julia> shift!(Ww,w)
Whirl.Fields.Edges{Whirl.Fields.Dual,5,4} data
u (in grid orientation):
  0.0   0.0   0.0   0.0
 11.5  12.5  13.5  14.5
  6.5   7.5   8.5   9.5
  0.0   0.0   0.0   0.0
v (in grid orientation):
 0.0  14.5  15.5  16.5  0.0
 0.0   9.5  10.5  11.5  0.0
 0.0   4.5   5.5   6.5  0.0

Note that the edges in the ghost cells are 0; these edges are not assigned any values in the shift operation.

We can then shift this to primal edges:

julia> q = Edges(Primal,w);

julia> shift!(q,Ww)
Whirl.Fields.Edges{Whirl.Fields.Primal,5,4} data
u (in grid orientation):
 0.0  6.0   6.5   7.0  0.0
 0.0  9.5  10.5  11.5  0.0
 0.0  3.5   4.0   4.5  0.0
v (in grid orientation):
 0.0   0.0   0.0  0.0
 6.0  12.5  13.5  7.0
 3.5   7.5   8.5  4.5
 0.0   0.0   0.0  0.0

We can also compute the Hadamard (i.e. element by element) product of any data of the same type, e.g.,

julia> q∘q
Whirl.Fields.Edges{Whirl.Fields.Primal,5,4} data
u (in grid orientation):
 0.0  36.0    42.25   49.0   0.0
 0.0  90.25  110.25  132.25  0.0
 0.0  12.25   16.0    20.25  0.0
v (in grid orientation):
  0.0     0.0     0.0    0.0
 36.0   156.25  182.25  49.0
 12.25   56.25   72.25  20.25
  0.0     0.0     0.0    0.0

The grid in physical space

Thus far, we have not had to consider the relationship between the grid's index space and some physical space. All of the operations thus far have acted on the entries in the discrete fields, based only on their relative indices, and not on their physical coordinates. In this section, we will discuss the relationship between the grid's index space and physical space, and then in the next section we'll discuss how we can transfer data between these spaces.

Generically, we can write the relationship between the physical coordinates $x$ and $y$, and the indices $i$ and $j$ of any grid point as

\[x(i) = (i - \Delta i - i_0)\Delta x, \quad y(j) = (j - \Delta j - j_0)\Delta x\]

The scaling between these spaces is controlled by $\Delta x$, which represents the uniform size of each grid cell; note that grid cells are presumed to be square in whirl. The indices $I_0 = (i_0,j_0)$ represent the location of the origin in the index space for primal nodes. Why primal nodes? Since the underlying grid is composed of dual cells, then primal nodes sit at the corners of the domain, so it is the most convenient for anchoring the grid to a specific point. But, since some field data of the same index are shifted by half a cell in one or both directions, then $\Delta i$ and $\Delta j$ are included for such purposes; these are either $0$ or $1/2$, depending on the field type. For example, for a primal node, $\Delta i = 0$, so that $x(i_0) = 0$; for a dual node, $\Delta i = 1/2$, so that $x(i_0) = -\Delta x/2$.

In particular, for our four different data types and their components

Regularization and interpolation

Based on this relationship between the physical space and the index space, we can now construct a means of transferring data between a point $(x,y)$ in the physical space and the grid points in its immediate vicinity. We say that such a point is immersed in the grid. The process of transferring from the point to the grid is called regularization, since we are effectively smearing this data over some extended neighborhood; the opposite operation, transferring grid field data to an arbitrary point, is interpolation. In whirl, both operations are carried out with the discrete delta function (DDF), which is a discrete analog of the Dirac delta function. The DDF generally has compact support, so that it only interacts with a small number of grid points in the vicinity of a given physical location. Since each of the different field types reside at slightly different locations, the range of indices invoked in this interaction will be different for each field type.

Regularization can actually take different forms. It can be a simple point-wise interpolation, the discrete analog of simply multiplying by the Dirac delta function:

\[f_i \delta(\mathbf{x} - \mathbf{x}_i)\]

to immerse a value $f_i$ based at point $\mathbf{x}_i = (x_i,y_i)$.

Alternatively, regularization can be carried out over a curve $\mathbf{X}(s)$, the analog of

\[\int f(s) \delta(\mathbf{x} - \mathbf{X}(s))\mathrm{d}s\]

or it can be performed volumetrically, corresponding to

\[\int f(\mathbf{y}) \delta(\mathbf{x} - \mathbf{y})\mathrm{d}\mathbf{y}\]

In this case, the function $f$ is distributed over some region of space. In each of these cases, the discrete version is simply a sum over data at a finite number of discrete points, and the type of regularization is specified by providing an optional argument specifying the arclength, area or volume associated with each discrete point. These arguments are used to weight the sum.

Let's see the regularization and interpolation in action. We will set up a ring of 100 points on a circle of radius $1/4$ centered at $(1/2,1/2)$. This curve- type regularization will be weighted by the arclength, $ds$, associated with each of the 100 points. On these points, we will set vector-valued data in which the $x$ component is uniformly equal to 1.0, while the $y$ component is set equal to the vertical position relative to the circle center. We will regularize these vector data to a primal edge field on the grid in which these points are immersed.

julia> n = 100;

julia> θ = linspace(0,2π,n+1);

julia> x = 0.5 + 0.25*cos.(θ[1:n]);

julia> y = 0.5 + 0.25*sin.(θ[1:n]);

julia> ds = 2π/n*0.25;

julia> X = VectorData(x,y);

The variable X now holds the coordinates of the immersed points. Now we will set up the vector-valued data on these points

julia> f = VectorData(X);

julia> fill!(f.u,1.0);

julia> f.v .= X.v.-0.5;

Note that we have ensured that f has the correct dimensions by supplying the coordinate data X. This first step also initializes the data to zeros.

Now, let's set up the grid. The physical domain will be of size $1.0 \times 1.0$, and we will use $100$ dual grid cells in each direction. Allowing a single layer of ghost cells surrounding the domain, we use $102$ cells, and set the cell size to 0.01. Also, we will set the $(x,y)$ origin to coincide with the lower left corner of the domain.

julia> nx = 102; ny = 102;

julia> q = Edges(Primal,(nx,ny));

julia> Lx = 1.0;

julia> dx = Lx/(nx-2)
0.01

Now we set up the regularization operator. To set it up, it needs to know the coordinate data of the set of immersed points, the grid cell size, and the weight to apply to each immersed point. Since this is a regularization of a curve, this weight is the differential arc length ds associated with each point. (This last argument is supplied as a scalar, since it is uniform.)

julia> H = Regularize(X,dx,weights=ds)
Regularization/interpolation operator with non-filtered interpolation
  100 points in grid with cell area 0.0001

We have omitted some optional arguments. For example, it chooses a default DDF kernel (the Roma kernel); this can be changed with the ddftype argument. Also, the lower left corner, where we've set the origin, is the location of the $(1,1)$ primal node; this is the default choice for I0 (the tuple $I_0$ of coordinates in index space discussed in the previous section).

Now we can apply the regularization operator. We supply the target field q as the first argument and the source data f as the second argument.

julia> H(q,f);

julia> plot(q)
Plot{Plots.PyPlotBackend() n=2}

We could also regularize this to a field of dual edges.

julia> p = Edges(Dual,(nx,ny));

julia> H(p,f);

julia> plot(p)
Plot{Plots.PyPlotBackend() n=2}

Scalar-valued data on the immersed points can only be regularized to nodal fields; the syntax is similar, and the regularization operator does not need to be reconstructed:

julia> g = ScalarData(X);

julia> fill!(g,1.0);

julia> w = Nodes(Dual,(nx,ny));

julia> H(w,g);

julia> plot(w)
Plot{Plots.PyPlotBackend() n=1}

For a given regularization operator, $H$, there is a companion interpolation operator, $E$. In whirl, this interpolation is also carried out with the same constructed operator, but with the arguments reversed: the grid field data are the source and the immersed points are the target. Note that interpolation is always a volumetric operation, so the weights assigned during the construction of the operator are not used in interpolation. Let's interpolate our regularized field back onto the immersed points.

julia> f2 = VectorData(X);

julia> H(f2,q);

julia> plot(f2.u,lab="u")
Plot{Plots.PyPlotBackend() n=1}

julia> plot!(f2.v,lab="v")
Plot{Plots.PyPlotBackend() n=2}

Note that interpolation is not the inverse of regularization; we don't recover the original data when we regularize and then interpolate. However, there is generally a way to scale the quantities on the immersed points and on the grid so that $H = E^T$. If we want to force these operations to be transposes of each other, we can supply the issymmetric flag:

julia> H = Regularize(X,dx,issymmetric=true)
Symmetric regularization/interpolation operator with non-filtered interpolation
  100 points in grid with cell area 0.0001

This flag will override any supplied weights.

If we expect to carry out the regularization and interpolation a lot, then it is often sensible to construct matrix versions of these operators. This construction is sometimes a bit slow, but the resulting operators perform their operations much faster than the matrix-free operators described above. To generate these matrix operators, we have to supply the data types of the source and target of the operation. For example, for regularization from scalar field data to dual node data,

julia> g = ScalarData(X);

julia> w = Nodes(Dual,(nx,ny));

julia> Hmat = RegularizationMatrix(H,g,w);

julia> fill!(g,1.0);

julia> w .= Hmat*g;
ERROR: MethodError: no method matching *(::Tuple{Whirl.Fields.RegularizationMatrix{Whirl.Fields.Nodes{Whirl.Fields.Dual,102,102},Whirl.Fields.ScalarData{100}},Whirl.Fields.InterpolationMatrix{Whirl.Fields.Nodes{Whirl.Fields.Dual,102,102},Whirl.Fields.ScalarData{100}}}, ::Whirl.Fields.ScalarData{100})
Closest candidates are:
  *(::Any, ::Any, !Matched::Any, !Matched::Any...) at operators.jl:424
  *(!Matched::PyCall.PyObject, ::Any) at /home/julia/JuliaPro-0.6.2.2/JuliaPro/pkgs-0.6.2.2/v0.6/PyCall/src/pyoperators.jl:11
  *(!Matched::RowVector{T<:Real,CV} where CV<:(ConjArray{T,1,V} where V<:(AbstractArray{T,1} where T) where T), ::AbstractArray{T<:Real,1}) where T<:Real at linalg/rowvector.jl:171
  ...

In general, the interpolation matrix is separately constructed, and the source and target are reversed:

julia> Emat = InterpolationMatrix(H,w,g);

julia> g .= Emat*w;

Alternatively, if the regularization and interpolation are symmetric, then we can get them both when we call for the regularization matrix:

julia> H = Regularize(X,dx,issymmetric=true)
Symmetric regularization/interpolation operator with non-filtered interpolation
  100 points in grid with cell area 0.0001

julia> Hmat, Emat = RegularizationMatrix(H,g,w)
(Regularization matrix acting on type Whirl.Fields.ScalarData{100} and returning type Whirl.Fields.Nodes{Whirl.Fields.Dual,102,102}, Interpolation matrix acting on type Whirl.Fields.Nodes{Whirl.Fields.Dual,102,102} and returning type Whirl.Fields.ScalarData{100})

It might seem a bit funny to store them separately if they are just transposes of each other, but it is essential for the method dispatch that they are given separate types.

Methods

CircularConvolution{M, N}

A preplanned, circular convolution operator on an M × N matrix.

Fields

  • : DFT coefficients of the convolution kernel

  • F: preplanned rFFT operator

  • F⁻¹: preplanned irFFT operator

  • paddedSpace: scratch space to zero-pad the input matrix

  • : scratch space to store the DFT coefficients of the zero-padded input matrix

Constructors:

  • CircularConvolution(G::Matrix{Float64})

Example:

julia> G = repmat(1.0:3,1,4)
3×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0
 3.0  3.0  3.0  3.0

julia> C = CircularConvolution(G)
Circular convolution on a 3 × 4 matrix

julia> C*reshape(1:12, 3, 4)
3×4 Array{Int64,2}:
 164  164  164  164
 130  130  130  130
 148  148  148  148
source
Whirl.Fields.DDFMethod.
DDF([ddftype=Fields.Roma],[dx=1.0])

Construct a discrete delta function operator. This is generally only needed internally by the Regularize operator, so the user doesn't have much need for accessing this directly. The default DDF is the Roma function, which has a support of 3 grid cells. Other choices are the Goza operator, which is a truncated Gaussian with 28 cells support, and the Witchhat, which has 2 cells support. The resulting operator is evaluated with one, two or three coordinate arguments, producing, respectively, 1-d, 2-d, or 3-d smeared delta functions. It can also be called with the usual Julia vectorized dot notation with arrays of arguments. The optional cell spacing argument dx rescales the coordinates by this spacing, and the result is also rescaled by this spacing (raised to the number of dimensions). This spacing argument defaults to 1.0.

julia> ddf = DDF(ddftype=Fields.Roma)
Discrete delta function operator of type Whirl.Fields.Roma, with spacing 1.0

julia> ddf(1)
0.16666666666666666

julia> ddf(-1)
0.16666666666666666

julia> ddf.([-1,0,1])
3-element Array{Float64,1}:
 0.166667
 0.666667
 0.166667
source
Edges{Dual/Primal}

Edges is a wrapper for vector-valued data that lie at the faces of either dual cells or primary cells. Edges type data have fields u and v for the components of the vector field. These are the normal components of the vector field on the vertical and horizontal faces of the corresponding cell.

Constructors

  • Edges(C,dims) creates a vector field of zeros in cells of type C (where C is either Dual or Primal), on a grid of dimensions dims. Note that dims represent the number of dual cells on the grid.

  • Edges(C,w) performs the same construction, but uses existing field data w of Nodes type to determine the size of the grid.

source
InterpolationMatrix(H::Regularize,u::Nodes/Edges,f::ScalarData/VectorData) -> Emat

Construct and store a matrix representation of interpolation associated with H for data of type u to data of type f. The resulting matrix Emat can then be used to apply on grid data of type u to interpolate it to point data of type f, using A_mul_B!(f,Emat,u). It can also be used as just Emat*u.

source
Nodes{Dual/Primal}

Nodes is a wrapper for scalar-valued data that lie at the centers of either dual cells or primary cells. A Nodes type can be accessed by indexing like any other array, and allows the use of [size].

Constructors

  • Nodes(C,dims) creates a field of zeros in cells of type C (where C is either Dual or Primal), on a grid of dimensions dims. Note that dims represent the number of dual cells on the grid, even if C is Primal.

  • Nodes(C,w) performs the same construction, but uses existing field data w of Nodes type to determine the size of the grid.

source
RegularizationMatrix(H::Regularize,f::ScalarData/VectorData,u::Nodes/Edges) -> Hmat

Construct and store a matrix representation of regularization associated with H for data of type f to data of type u. The resulting matrix Hmat can then be used to apply on point data of type f to regularize it to grid data of type u, using A_mul_B!(u,Hmat,f). It can also be used as just Hmat*f.

If H is a symmetric regularization and interpolation operator, then this actually returns a tuple Hmat, Emat, where Emat is the interpolation matrix.

source
Regularize(x,y,dx,[ddftype=Roma],[I0=(1,1)], [weights=1.0], [filter=false],
                   [issymmetric=false])

Constructor to set up an operator for regularizing and interpolating data from/to points immersed in the grid to/from fields on the grid itself. The supplied x and y represent physical coordinates of the immersed points, and dx denotes a uniform physical cell size of the grid. The separate arguments x and y can be replaced by a single argument X of type VectorData holding the coordinates.

The operations of regularization and interpolation are carried out with a discrete delta function (ddf), which defaults to the type Roma. Others are also possible, such as Goza. The optional tuple I0 represents the indices of the primary node that coincides with (x,y) = (0,0). This defaults to (1,1), which leaves one layer of ghost (dual) cells and sets the physical origin in the lower left corner of the grid of interior dual cells.

Another optional parameter, weights, sets the weight of each point in the regularization. This would generally be set with, say, the differential arc length for regularization of data on a curve. It can be a vector (of the same length as x and y) or a scalar if uniform. It defaults to 1.0.

The optional Boolean parameter filter can be set to true if it is desired to apply filtering (see Goza et al, J Comput Phys 2016) to the grid data before interpolating. This is generally only used in the context of preconditioning the solution for forces on the immersed points.

If the optional Boolean parameter issymmetric is set to true, then the regularization and interpolation are constructed to be transposes of each other. Note that this option overrides any supplied weights. The default of this parameter is false.

The resulting operator can be used in either direction, regularization and interpolation, with the first argument representing the target (the entity to regularize/interpolate to), and the second argument the source (the entity to regularize/interpolate from). The regularization does not use the filtering option.

Example

In the example below, we set up a 12 x 12 grid. Using the default value for I0 and setting dx = 0.1, the physical dimensions of the non-ghost part of the grid are 1.0 x 1.0. Three points are set up in the interior, and a vector field is assigned to them, with the x component of each of them set to 1.0. These data are regularized to a field of primal edges on the grid.

julia> x = [0.25,0.75,0.25]; y = [0.75,0.25,0.25];

julia> X = VectorData(x,y);

julia> q = Edges(Primal,(12,12));

julia> dx = 0.1;

julia> H = Regularize(x,y,dx)
Regularization/interpolation operator with non-filtered interpolation
  3 points in grid with cell area 0.01

julia> f = VectorData(X);

julia> fill!(f.u,1.0);

julia> H(q,f)
Whirl.Fields.Edges{Whirl.Fields.Primal,12,12} data
u (in grid orientation):
 0.0  0.0  0.0       0.0     0.0      …  0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     0.0       0.0     0.0      0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0      …  0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     8.33333  33.3333  8.33333  0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     8.33333  33.3333  8.33333  0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0      …  0.0       0.0     0.0      0.0  0.0
v (in grid orientation):
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
source
ScalarData

A wrapper for a one-dimensional array of scalar-valued data. The resulting wrapper can be indexed in the same way as the array itself.

Constructors

  • ScalarData(d) constructs a wrapper for the one-dimensional array of data d

  • ScalarData(n::Int) constructs a wrapper for an array of zeros of length n.

  • ScalarData(x::ScalarData) constructs a wrapper for an array of zeros of the same length as that wrapped by x.

  • ScalarData(x::VectorData) constructs a wrapper for an array of zeros of the same length as that wrapped by x.

Example

julia> f = ScalarData(10);

julia> f[5] = 1.0;

julia> f
10 points of scalar-valued data
 0.0
 0.0
 0.0
 0.0
 1.0
 0.0
 0.0
 0.0
 0.0
 0.0
source
VectorData

A wrapper for a one-dimensional array of two-component vector-valued data. The resulting wrapper can be indexed as though the first component and second component are stacked on top of each other.

Constructors

  • VectorData(u,v) constructs a wrapper for the vector components data u and v.

  • VectorData(n::Int) constructs a wrapper with zeros of length n for both components.

  • VectorData(x::ScalarData) constructs a wrapper for zero components of the same length as that wrapped by x.

  • VectorData(x::VectorData) constructs a wrapper for zero components of the same length as that wrapped by x.

Example

julia> f = VectorData(10);

julia> f.v[1:5] = 1:5;

julia> f
10 points of vector-valued data
 0.0  1.0
 0.0  2.0
 0.0  3.0
 0.0  4.0
 0.0  5.0
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0

julia> f[7] = 1; f[18] = 0.2;

julia> f
10 points of vector-valued data
 0.0  1.0
 0.0  2.0
 0.0  3.0
 0.0  4.0
 0.0  5.0
 0.0  0.0
 1.0  0.0
 0.0  0.2
 0.0  0.0
 0.0  0.0
source
cooordinates(w::Nodes/Edges;[dx=1.0],[I0=(1,1)])

Return a tuple of the ranges of the physical coordinates in each direction for grid data w. If w is of Nodes type, then it returns a tuple of the form xg,yg. If w is of Edges or NodePair type, then it returns a tuple of the form xgu,ygu,xgv,ygv.

The optional keyword argument dx sets the grid spacing; its default is 1.0. The optional keyword I0 accepts a tuple of integers to set the index pair of the primal nodes that coincide with the origin. The default is (1,1).

Example

julia> w = Nodes(Dual,(12,22));

julia> xg, yg = coordinates(w,dx=0.1)
(-0.05:0.1:1.05, -0.05:0.1:2.0500000000000003)
source
Whirl.Fields.curl!Method.
curl!(q::Edges{Primal},w::Nodes{Dual})

Evaluate the discrete curl of w and return it as q.

Example

julia> w = Nodes(Dual,(8,6));

julia> w[3,4] = 1.0;

julia> q = Edges(Primal,w);

julia> curl!(q,w)
Whirl.Fields.Edges{Whirl.Fields.Primal,8,6} data
u (in grid orientation):
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  -1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  -1.0  1.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
source
Whirl.Fields.curlMethod.
curl(w::Nodes{Dual}) --> Edges{Primal}

Evaluate the discrete curl of w. Another way to perform this operation is to construct a Curl object and apply it with *.

Example

julia> C = Curl();

julia> w = Nodes(Dual,(8,6));

julia> w[3,4] = 1.0;

julia> C*w
Whirl.Fields.Edges{Whirl.Fields.Primal,8,6} data
u (in grid orientation):
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  -1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  -1.0  1.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
source
divergence!(w::Nodes,q::Edges)

Evaluate the discrete divergence of edge data q and return it as nodal data w. Note that q can be either primal or dual edge data, but w must be of the same cell type.

Example

julia> q = Edges(Primal,(8,6));

julia> q.u[3,2] = 1.0;

julia> w = Nodes(Primal,(8,6));

julia> divergence!(w,q)
Whirl.Fields.Nodes{Whirl.Fields.Primal,8,6} data
Printing in grid orientation (lower left is (1,1)):
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  1.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
source
divergence(q::Edges) --> Nodes

Evaluate the discrete divergence of edge data q. Can also perform this operation by creating an object of Divergence type and applying it with *.

Example

julia> D = Divergence();

julia> q = Edges(Primal,(8,6));

julia> q.u[3,2] = 1.0;

julia> D*q
Whirl.Fields.Nodes{Whirl.Fields.Primal,8,6} data
Printing in grid orientation (lower left is (1,1)):
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  1.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
source
Whirl.Fields.grad!Method.
grad!(q::Edges{Primal},w::Nodes{Primal})

Evaluate the discrete gradient of primal nodal data w and return it as primal edge data q.

Example

julia> w = Nodes(Primal,(8,6));

julia> w[3,4] = 1.0;

julia> q = Edges(Primal,(8,6));

julia> grad!(q,w)
Whirl.Fields.Edges{Whirl.Fields.Primal,8,6} data
u (in grid orientation):
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
source
Whirl.Fields.gradMethod.
grad(w::Nodes{Primal}) --> Edges{Primal}

Evaluate the discrete gradient of primal nodal data w. Can also perform this operation by creating an object of Grad type and applying it with *.

Example

julia> w = Nodes(Primal,(8,6));

julia> w[3,4] = 1.0;

julia> G = Grad();

julia> G*w
Whirl.Fields.Edges{Whirl.Fields.Primal,8,6} data
u (in grid orientation):
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
source
laplacian!(v,w)

Evaluate the discrete Laplacian of w and return it as v. The data w can be of type dual/primary nodes or edges; v must be of the same type.

Example

julia> w = Nodes(Dual,(8,6));

julia> v = deepcopy(w);

julia> w[4,3] = 1.0;

julia> laplacian!(v,w)
Whirl.Fields.Nodes{Whirl.Fields.Dual,8,6} data
Printing in grid orientation (lower left is (1,1)):
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  -4.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
source
laplacian(w)

Evaluate the discrete Laplacian of w. The data w can be of type dual/primary nodes or edges. The returned result is of the same type as w.

Example

julia> q = Edges(Primal,(8,6));

julia> q.u[2,2] = 1.0;

julia> laplacian(q)
Whirl.Fields.Edges{Whirl.Fields.Primal,8,6} data
u (in grid orientation):
 0.0   0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0   1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  -4.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
source
plan_intfact(a::Real,dims::Tuple,[fftw_flags=FFTW.ESTIMATE])

Constructor to set up an operator for evaluating the integrating factor with real-valued parameter a. This can then be applied with the * operation on data of the appropriate size.

The dims argument can be replaced with w::Nodes to specify the size of the domain.

Example

julia> w = Nodes(Dual,(6,6));

julia> w[4,4] = 1.0;

julia> E = plan_intfact(1.0,(6,6))
Integrating factor with parameter 1.0 on a (nx = 6, ny = 6) grid

julia> E*w
Whirl.Fields.Nodes{Whirl.Fields.Dual,6,6} data
Printing in grid orientation (lower left is (1,1)):
 0.00268447   0.00869352  0.0200715   0.028765    0.0200715   0.00869352
 0.00619787   0.0200715   0.0463409   0.0664124   0.0463409   0.0200715
 0.00888233   0.028765    0.0664124   0.0951774   0.0664124   0.028765
 0.00619787   0.0200715   0.0463409   0.0664124   0.0463409   0.0200715
 0.00268447   0.00869352  0.0200715   0.028765    0.0200715   0.00869352
 0.000828935  0.00268447  0.00619787  0.00888233  0.00619787  0.00268447
source
plan_intfact!(a::Real,dims::Tuple,[fftw_flags=FFTW.ESTIMATE])

Same as plan_intfact, but operates in-place on data.

source
plan_laplacian(dims::Tuple,[with_inverse=false],[fftw_flags=FFTW.ESTIMATE],
                      [dx=1.0])

Constructor to set up an operator for evaluating the discrete Laplacian on dual or primal nodal data of dimension dims. If the optional keyword with_inverse is set to true, then it also sets up the inverse Laplacian (the lattice Green's function, LGF). These can then be applied, respectively, with * and \ operations on data of the appropriate size. The optional parameter dx is used in adjusting the uniform value of the LGF to match the behavior of the continuous analog at large distances; this is set to 1.0 by default.

Instead of the first argument, one can also supply w::Nodes to specify the size of the domain.

Example

julia> w = Nodes(Dual,(5,5));

julia> w[3,3] = 1.0;

julia> L = plan_laplacian(5,5;with_inverse=true)
Discrete Laplacian (and inverse) on a (nx = 5, ny = 5) grid with spacing 1.0

julia> s = L\w
Whirl.Fields.Nodes{Whirl.Fields.Dual,5,5} data
Printing in grid orientation (lower left is (1,1)):
 0.16707    0.129276     0.106037     0.129276    0.16707
 0.129276   0.0609665   -0.00734343   0.0609665   0.129276
 0.106037  -0.00734343  -0.257343    -0.00734343  0.106037
 0.129276   0.0609665   -0.00734343   0.0609665   0.129276
 0.16707    0.129276     0.106037     0.129276    0.16707

julia> L*s ≈ w
true
source
plan_laplacian!(dims::Tuple,[with_inverse=false],[fftw_flags=FFTW.ESTIMATE],
                      [dx=1.0])

Same as plan_laplacian, but operates in-place on data.

source
product!(out::Edges/Nodes,p::Edges/Nodes,q::Edges/Nodes)

Compute the Hadamard (i.e. element by element) product of edge or nodal (primal or dual) data p and q and return the result in out.

Example

julia> q = Edges(Dual,(8,6));

julia> out = p = deepcopy(q);

julia> q.u[3,2] = 0.3;

julia> p.u[3,2] = 0.2;

julia> product!(out,p,q)
Whirl.Fields.Edges{Whirl.Fields.Dual,8,6} data
u (in grid orientation):
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.06  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
source
product(p::Edges/Nodes,q::Edges/Nodes) --> Edges/Nodes

Compute the Hadamard product of edge or nodal (primal or dual) data p and q and return the result. This operation can also be carried out with the operator:

Example

julia> q = Edges(Dual,(8,6));

julia> p = deepcopy(q);

julia> q.u[3,2] = 0.3;

julia> p.u[3,2] = 0.2;

julia> p∘q
Whirl.Fields.Edges{Whirl.Fields.Dual,8,6} data
u (in grid orientation):
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.06  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
source
Base.shift!Method.
shift!((wx::Nodes,wy::Nodes),q::Edges)

Shift (by linear interpolation) the edge data q (of either dual or primal type) to the dual or primal nodes, and return the result in wx and wy. wx holds the shifted q.u data and wy the shifted q.v data.

Example

julia> q = Edges(Primal,(8,6));

julia> q.u[3,2] = 1.0;

julia> wx = Nodes(Dual,(8,6)); wy = deepcopy(wx);

julia> Fields.shift!((wx,wy),q);

julia> wx
Whirl.Fields.Nodes{Whirl.Fields.Dual,8,6} data
Printing in grid orientation (lower left is (1,1)):
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> wy
Whirl.Fields.Nodes{Whirl.Fields.Dual,8,6} data
Printing in grid orientation (lower left is (1,1)):
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
source
Base.shift!Method.
shift!(v::Edges{Dual/Primal},q::Edges{Primal/Dual})

Shift (by linear interpolation) the primal (resp. dual) edge data q to the edges of the dual (resp. primal) cells, and return the result in v.

Example

julia> q = Edges(Primal,(8,6));

julia> q.u[3,2] = 1.0;

julia> v = Edges(Dual,(8,6));

julia> Fields.shift!(v,q)
Whirl.Fields.Edges{Whirl.Fields.Dual,8,6} data
u (in grid orientation):
 0.0  0.0   0.0   0.0  0.0  0.0  0.0
 0.0  0.0   0.0   0.0  0.0  0.0  0.0
 0.0  0.0   0.0   0.0  0.0  0.0  0.0
 0.0  0.25  0.25  0.0  0.0  0.0  0.0
 0.0  0.25  0.25  0.0  0.0  0.0  0.0
 0.0  0.0   0.0   0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
source
Base.shift!Method.
shift!(q::Edges{Dual},w::Nodes{Dual})

Shift (by linear interpolation) the dual nodal data w to the edges of the dual cells, and return the result in q.

Example

julia> w = Nodes(Dual,(8,6));

julia> w[3,4] = 1.0;

julia> q = Edges(Dual,w);

julia> shift!(q,w)
Whirl.Fields.Edges{Whirl.Fields.Dual,8,6} data
u (in grid orientation):
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.5  0.5  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
v (in grid orientation):
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
source

Index