Fields
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:
dual nodes,
dual edges,
primal nodes, and
primal edges.
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
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
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
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
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
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
Primal nodes: $\Delta i = 0$, $\Delta j = 0$
Dual nodes: $\Delta i = 1/2$, $\Delta j = 1/2$
Primal edges u: $\Delta i = 1/2$, $\Delta j = 0$
Primal edges v: $\Delta i = 0$, $\Delta j = 1/2$
Dual edges u: $\Delta i = 0$, $\Delta j = 1/2$
Dual edges v: $\Delta i = 1/2$, $\Delta j = 0$
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:
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
or it can be performed volumetrically, corresponding to
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
Whirl.Fields.CircularConvolution
— Type.CircularConvolution{M, N}
A preplanned, circular convolution operator on an M × N matrix.
Fields
Ĝ
: DFT coefficients of the convolution kernelF
: preplanned rFFT operatorF⁻¹
: preplanned irFFT operatorpaddedSpace
: 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
Whirl.Fields.DDF
— Method.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
Whirl.Fields.Edges
— Type.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 typeC
(whereC
is eitherDual
orPrimal
), on a grid of dimensionsdims
. Note thatdims
represent the number of dual cells on the grid.Edges(C,w)
performs the same construction, but uses existing field dataw
ofNodes
type to determine the size of the grid.
Whirl.Fields.InterpolationMatrix
— Type.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
.
Whirl.Fields.Nodes
— Type.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 typeC
(whereC
is eitherDual
orPrimal
), on a grid of dimensionsdims
. Note thatdims
represent the number of dual cells on the grid, even ifC
isPrimal
.Nodes(C,w)
performs the same construction, but uses existing field dataw
ofNodes
type to determine the size of the grid.
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.
Whirl.Fields.Regularize
— Method.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
Whirl.Fields.ScalarData
— Type.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 datad
ScalarData(n::Int)
constructs a wrapper for an array of zeros of lengthn
.ScalarData(x::ScalarData)
constructs a wrapper for an array of zeros of the same length as that wrapped byx
.ScalarData(x::VectorData)
constructs a wrapper for an array of zeros of the same length as that wrapped byx
.
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
Whirl.Fields.VectorData
— Type.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 datau
andv
.VectorData(n::Int)
constructs a wrapper with zeros of lengthn
for both components.VectorData(x::ScalarData)
constructs a wrapper for zero components of the same length as that wrapped byx
.VectorData(x::VectorData)
constructs a wrapper for zero components of the same length as that wrapped byx
.
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
Whirl.Fields.coordinates
— Function.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)
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
Whirl.Fields.curl
— Method.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
Whirl.Fields.divergence!
— Method.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
Whirl.Fields.divergence
— Method.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
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
Whirl.Fields.grad
— Method.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
Whirl.Fields.laplacian!
— Method.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
Whirl.Fields.laplacian
— Method.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
Whirl.Fields.plan_intfact
— Function.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
Whirl.Fields.plan_intfact!
— Function.plan_intfact!(a::Real,dims::Tuple,[fftw_flags=FFTW.ESTIMATE])
Same as plan_intfact
, but operates in-place on data.
Whirl.Fields.plan_laplacian
— Function.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
Whirl.Fields.plan_laplacian!
— Function.plan_laplacian!(dims::Tuple,[with_inverse=false],[fftw_flags=FFTW.ESTIMATE],
[dx=1.0])
Same as plan_laplacian
, but operates in-place on data.
Whirl.Fields.product!
— Method.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
Whirl.Fields.product
— Method.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
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
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
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
Index
Whirl.Fields.CircularConvolution
Whirl.Fields.DDF
Whirl.Fields.Edges
Whirl.Fields.InterpolationMatrix
Whirl.Fields.Nodes
Whirl.Fields.RegularizationMatrix
Whirl.Fields.Regularize
Whirl.Fields.ScalarData
Whirl.Fields.VectorData
Base.shift!
Base.shift!
Base.shift!
Whirl.Fields.coordinates
Whirl.Fields.curl
Whirl.Fields.curl!
Whirl.Fields.divergence
Whirl.Fields.divergence!
Whirl.Fields.grad
Whirl.Fields.grad!
Whirl.Fields.laplacian
Whirl.Fields.laplacian!
Whirl.Fields.plan_intfact
Whirl.Fields.plan_intfact!
Whirl.Fields.plan_laplacian
Whirl.Fields.plan_laplacian!
Whirl.Fields.product
Whirl.Fields.product!