Towards a High-Order Embedded Boundary Finite Volume Method for the Incompressible Navier-Stokes Equations with Complex Geometries

In this paper


I. Introduction
The finite volume method is inherently conservative in its construction, which is a valuable property when solving for fluid flows with combustion or shock waves present.When operating on structured grids, finite volume methods have been developed to achieve high solution accuracy and code performance.This is demonstrated by the Chombo framework [1], where adaptive mesh refinement is one of the key features for improving solution accuracy with a high degree of parallelism.High-order methods have successfully been created with structured grid solvers [2] due to the convenience of polynomial reconstructions with regular data.Despite these advantages, representation of complex geometries on structured grids is a significant challenge.Technologies such as mapped multi-block methods allow for moderately complex structured grid geometries by combination of several curvilinear grids.However, mapped multi-block grids struggle to represent rough surfaces, and creating quality grids is generally time-consuming.
Alternatively, embedded boundary grids have little restriction in geometries that may be represented, while maintaining the advantages of structured grid solvers.Embedded boundary grids are created by first generating a Cartesian grid independent of boundary geometry.Then, the boundary geometry is overlain on the grid, and cells that intersect the boundaries are cut (Fig. 1).The resulting grid is one that is structured everywhere away from the boundary, while near the boundary the grid is made up of partial, or "cut" cells.The resulting method retains the advantages of solving on structured grids on the interior of the domain, while having relatively little restriction in geometries that may be represented.Additionally, this grid generation process can be done efficiently and quickly.However, the presence of cut cells brings about additional challenges that must be overcome.Elaborate reconstruction schemes are required near the boundaries of embedded boundary methods, since regular or grid-aligned stencils are not accurate or stable in the presence of cut cells.Volumes of cut cells may also be arbitrarily small with respect to the Cartesian cells they originate from, and maintaining stability of these small cells requires special care.
The embedded boundary method has been used with great success in a number of complex fluid dynamics applications [3][4][5][6].Previous time dependent applications have only achieved up to second-order accurate solutions, with first-order or inconsistent results near embedded boundaries [7].A high-order finite volume embedded boundary method for smooth and kinked ( 0 ) domains was demonstrated for Poisson's equation by Devendran et al [8].The goal of this work is to extend the fourth-order embedded boundary method to solve the time dependent incompressible Navier-Stokes equations.Advancing the diffusion terms in time is straightforward using a fourth-order implicit time marching method.Advection can be treated with explicit time marching methods, but doing so requires special consideration for arbitrarily small cut cells.These two schemes can be combined using an additive Runge-Kutta (ARK) scheme coupled with a projection method to solve the incompressible Navier-Stokes equations.

II. The Finite-Volume Method for Embedded Boundaries
Finite volume methods express time dependent PDEs in "flux-divergence" form of based on conservation laws over discrete control volumes.To apply the finite volume method, the system of PDEs is converted to integral form over volumes V  , and the divergence theorem is applied to the flux term ì F yielding The flux term may be broken into separate regions over the volume surface V  as where  indexes a particular face of the cell indexed by , and n  is the corresponding outward normal vector.Averaged quantities are defined as Using cell averaged notation, the finite volume scheme is expressed in terms of averages as This is an exact formulation, as no approximations have been made.In practice, the fluxes and time derivatives are not known exactly, so a numerical solution is created by producing approximations for each.Cut cells can potentially be arbitrarily small in a limit where V  approaches zero.This is numerically troublesome, since the flux term is potentially divided by zero, and the other geometric quantities and stencils must have similar limits to avoid this.Thus numerical stability for small cells requires extra care, and is dealt with by scaling with a cell volume fraction, as detailed in Devendran et al [8].

A. Projection Formulation for the Incompressible Navier-Stokes Equations
The incompressible Navier-Stokes equations with constant density are given by where u is the flow velocity,  the pressure, and  the kinematic viscosity.
A Hodge projection operator P has been used in the finite volume literature [9,10] to enforce a divergence-free velocity field:

B. Finite Volume Projection Formulation
We choose discretizations for each of the spatial operators in Eqs.(3)(4), and write the resulting discrete equations at grid locations  as where A, D, G, and L are fourth-order finite volume approximations of advection, divergence, gradient, and Laplacian terms, respectively.(From this point forward, we will drop the subscript  except for clarity.)The goal is to discretize these operators so that in the regular interior of the domain, with some potential loss of accuracy near boundaries and in cut cells.Because we are using co-located cell-average velocity and pressure, we take the approach of an approximate projection [11].Instead of a strictly zero discrete divergence, we allow u to have a divergence that is at the level of the discretization error.The equivalent discrete projection is which requires the inversion of the Laplacian operator over the whole domain.Using the projection operator, the incompressible flow equations can be approximated by The general procedure for solving this system separates into an intermediate update that integrates the right-hand side of Eq. ( 16), and then projects that update using Eq. ( 14), to maintain an approximately-divergence free stage value for the time integrator.This is essentially a higher-order accurate version of the projection operator described in [7].

Advection-Diffusion
In the projected form of the incompressible Navier-Stokes equations given by Eq. ( 16) there are two components required for a solution: the projection operator, and an advection-diffusion time integration scheme.As building blocks for an incompressible Navier-Stokes solver, the projection and advection-diffusion schemes can be developed and tested separately.Showing stability and correctness of the approximate projection is required for the incompressible Navier-Stokes equations.To show algorithm verification, the scalar form of the advection-diffusion equation is sufficient to prove the order of accuracy of the high-order embedded boundary algorithm developed.Once the components are proven correct, the process of combining these pieces to form an incompressible Navier-Stokes solver should be straightforward.The only remaining challenge is introducing the non-linear vector advection term, A(u).

Projection Formulation for Open Boundaries
We must specify boundary conditions for the projection operators for open domain boundaries to split a given w into three separate components: a scalar potential flow solution, , which only satisfies the boundary conditions, and the remainder which is split into pure gradient, , and divergence-free parts: w = G + v + G, such that In the end, our desired divergence-free velocity field that satisfies the correct boundary conditions is just the components u = v + G.

III. Embedded Boundary Spatial Discretization
In the embedded boundary (EB) approach, cells fall into one of three categories; regular cells, irregular cells, and invalid cells.This distinction between cell types is illustrated in Fig. 1.Regular cells are those that are full Cartesian cells.Irregular cells, or cut cells, are those which are partial cells because they intersect with the boundary geometry.Invalid cells, as the name indicates, are cells that fall outside the domain boundaries and are thus not in the solution domain.To denote the embedded boundary regions, let Ω be the irregular domain and Υ  be any regular cell, then a particular cell can be denoted by The challenge of embedded boundary methods is to approximate the flux terms ⟨F⟩  when regular grid stencils can not be used due to nearby cut cells.To solve this requires knowledge of the flux as a function, so that face averages may be computed appropriately.A general reconstruction is done by creating local polynomials and evaluating their value and derivatives on the face as needed.For a structured grid, reconstructed polynomials lead to grid-aligned regular stencils.When using an embedded boundary method, the cells near the boundary require a consistent approach to generate stencils that depend on the local geometry.On irregular regions, our approach is to use a weighted least-squares polynomial approximation from the cell average values in a local region of neighboring cells.The scheme presented theoretically can produce any order spatial discretization desirable, but for the context of this paper, fourth-order is the focus.The present study continues and extends the work developed by Devendran et al. [8], and the procedure is briefly summarized in this section.

A. Multi-Dimensional Taylor Expansion
A multi-dimensional polynomial can be defined using where q is a multi-index or D dimensional non-negative integer vector, and x is a given point in space.For a sufficiently smooth scalar function , its multi-dimensional Taylor series of order Q can be written as where the convention 0! = 1 is used, and  (q) is the multi-index partial derivative notation and any q  = 0 implies no derivative.In formulating multi-dimensional polynomials to fit cell averaged data, integration over cells is needed.Moments are defined to be the integration of basis polynomials over a region, as for volume moments, while for face moments

B. Flux Reconstruction
To reconstruct a high-order solution from cell averaged data, we use a Taylor expansion about cell centers x .By defining the multi-dimensional polynomial coefficients ( x ), approximated cell averages can be represented as The number of unknown coefficients is determined by  = (D+Q)!D! Q! .Determining the  coefficients requires solving a linear system of equations including at least the same number of cell values ⟨u  ⟩  .If the values ⟨u  ⟩  are more than enough linearly independent values to determine each coefficient, this is an over-determined linear least-squares system.Note that there are degenerate cases where values are almost linearly dependent, in which case the least-squares system will typically distribute the solution between those values.Thus, given a vector  that is a vector of neighboring cell-average velocities ⟨u  ⟩  , using matrix notation we can write the approximation as  ≈   , where the geometric moment matrix  comes from Eq. ( 24) for each of the neighboring cells , and the multi-dimensional coefficients vector  is to be determined.
A diagonal weighting matrix  is added to enable further control of the least-squares fit with the goal to improve stability as in [8]; further discussion is in section III.C.With the weighting matrix, this yields the system and solution for the coefficients as where ( ) † indicates the pseudo-inverse of ( ), which is equivalent to a true inverse when the matrix is square.

Higher-order EB Viscous Flux Stencil
Using this method of reconstruction gives a way to evaluate fluxes for Eq. ( 2).In the convenient case when the flux is linear, such as the diffusion term, a Taylor expansion is applied as This too is a linear equation, which can be expressed in matrix form as in Eq. (25) to derive an expression for the corresponding stencil as in [8]: where  is the same coefficient as before, and and The viscous flux stencil   depends only on the local neighbors selected, their moments, and flux function, so it may be computed once per grid at initialization and cached as a sparse matrix operator.

Higher-order EB Advection Flux Stencil
Unlike the diffusion term, the scalar advection term does not use boundary conditions since the solution is defined only by velocity characteristics.However, we can still evaluate the velocity average on faces as Where   is the row vector of moments of the face  , and  up is an upwind-weight matrix that is evaluated using the face-centroid velocity field for simplicity.The reconstruction vector   may still be cached, since it only depends upon the grid.
For the advection term ⟨∇ • (uu)⟩, we combine upwind polynomial approximations for u to create face-average fluxes using a convolution formula (see [12]): avg =  up  up †  avg (28) avg = vector that evaluates coefficients to get fourth-order ⟨u⟩  (30) 1 = vector that evaluates coefficients to get second-order ⟨  ′ u⟩  , tangential to the face.
In 2D, Eq. ( 27) holds because odd moments of the face around the centroid x are 0, so the convolution formula only involves even moments, which are simply powers of A  .With some algebraic manipulation (see [12]) for  (ℎ 4 ), it reduces to this formula.However, in 3D or higher-order in 2D, there are additional terms that must be included for the general case without this cancellation.

Approximate Projection
The higher-order projection operator must start with cell average velocities, ⟨u⟩  , and correct them to be approximately divergence-free, Eq. ( 16).To accomplish this, we require discretizations and boundary conditions for each operator in Eq. ( 14).
First, the Laplacian operator L is essentially the same as the viscous flux in section III.B.1, but with different boundary conditions based on Eq. ( 18), that are homogeneous Neumann for inflow or solid walls, and homogeneous Dirichlet for outflow.
For the cell-average gradient operator, G, the least-squares polynomial fit is again used to evaluate a stencil just as in Eq. ( 26).In this case, as it is not a finite volume flux, cell moments are used instead of face moments: The gradient operator uses the same boundary conditions as the Laplacian.
Finally, we define a cell-average Divergence operator using the divergence theorem, so the cell average of the divergence of the velocity field can be determined from face average quantities as These face averages are similar to those in the advection flux, however, without any upwind weighting for the reconstruction.The terms on regular faces have single component normal vectors n  =   , and as a result ⟨u • n  ⟩  = ⟨u  ⟩  .Boundary conditions for the divergence are u • n = 0 at any solid walls, including EB boundaries.At inflow boundaries, the velocity is specified, while on outflow no boundary conditions are used, as in [7].

Physical Boundary Conditions
For regions where stencils encompass a Dirichlet boundary condition, the boundary value is defined as u bc = ⟨u⟩  when  corresponds to a face with a prescribed value.A polynomial can be reconstructed to fit the face value using Including this additional equation and value into the stencil system will match the boundary condition in a least-squares sense, with a similar method for Neumann boundaries.These extra boundary condition equations are used when any neighboring cell included in the reconstruction contains a portion of the boundary, and so can accommodate different parts of the boundary (such as corners, etc.).

Interpolation Neighborhoods
An important aspect of the weighted least squares method is neighborhood selection, choosing which neighboring cells to include in the stencil, which becomes important when dealing with embedded boundaries.In this paper, the focus is on fourth-order accurate stencils in two dimensions, which require a minimum of 10 neighbors for flux reconstruction in the most general case, which must span both the  and  directions to estimate higher-order derivatives.To achieve this with the developed weighted least squares approach, stencils of radius 3 cells from the reconstructed face are used (see Fig. 2).These stencils encompass 18 cells in regular regions.For stencils in embedded boundary regions, the same radius of 3 cells is used along with the inclusion of boundary conditions.

C. The Problem of Small Cut Cells
To reconstruct the stencils required for the spatial discretization, the weighted least-squares method Eq. ( 25) is used to simplify generation of stencils near the embedded boundaries.First, a large neighborhood is created to interpolate from, which includes many potential neighboring cells, some of which may be very small cut cells.These cut cells with potential small cell volumes can result in a poorly conditioned system of equations.Rather than exhaustively choosing an interpolation neighborhood which is better conditioned, a weighting scheme is chosen to assign relative importance to each cell's entry in the WLS system.Smaller relative weights mean that cells has less importance and a smaller absolute stencil value.An effective weighting for cell  with target construction face  of where  ,  is the Euclidean distance, has been shown to increase solution stability and spectral properties [8], especially when interpolation neighborhoods are large and there are many possible consistent stencils.

IV. Time Marching Method
When solving the Navier-Stokes equations, time constraints for the advection and diffusion portions of the fluxes can be significantly different.The advection terms are hyperbolic in nature and generally non-linear and less stiff, making them well suited for explicit time marching methods.In contrast, the diffusion fluxes are parabolic, more stiff, and linear, making them well suited for implicit time marching methods using fast linear solvers.In the context of embedded boundary methods, where cut cells can become arbitrarily small, this difference of time step stability constraints between advection and diffusion terms is can be orders of magnitude different.To maintain reasonable time step sizes, a hybrid implicit-explicit (ImEx) Runge-Kutta method is chosen for the embedded boundary algorithm [13].This allows for the advection terms to be updated using an explicit method, and the diffusion terms with an implicit method.As a result, the time evolution is generally not limited by the small time steps required for the diffusion physics.At each stage of the implicit RK time marching, a large sparse matrix system must be solved.To do this efficiently, we use a geometric multigrid solver included in Chombo [14] with PETSC [15,16] as a bottom solver.
To achieve the desired fourth-order accuracy of the presented algorithm, the ARK4(3)6L [2]SA time marching method [17] is used.This method has successfully been demonstrated for stiff, higher-order finite volume methods that use AMR for advection-diffusion problems [18].

V. Algorithm Verification
To show verification of the algorithm developed, the method of manufactured solutions is used.By producing an artificial exact solution, solution error can be computed and used for testing grid convergence.The method of manufactured solutions is explained in detail by Salari et al. [19].First, a manufactured solution form is chosen.Once the manufactured solution is chosen, a source term and boundary conditions are derived by substituting the manufactured solution into the governing equations.
For algorithm verification with more complicated solutions is to use Richardson extrapolation to evaluate the order of convergence, by comparing coarser solutions to the finest solution in a series, and evaluate the relative errors for convergence rates.
To verify our algorithm, convergence tests are used to demonstrate the claimed solution order of accuracy.The solution error is computed, and convergence rates evaluated using each of the  ∞ ,  2 , and  1 norms.The solution error is computed from cell averages as The solution norms are computed as where   is the number of cells in the domain Ω.Convergence can be tested by the expectation on a grid having cell size ℎ that error  ℎ =  (ℎ Q ), and thus Q = log   ℎ  ℎ /log  should be observed when  is the refinement ratio.

VI. Results
To demonstrate the developed algorithm, solutions of projections, scalar advection, and diffusion were evaluated on a set of simple embedded boundary geometries.To begin with, verification is performed for each part of the developed algorithm.First, a potential flow case is examined to verify correctness of the projection operator, verify solution accuracy in space, and determine stability.Next, manufactured solutions for the scalar diffusion equation are tested to verify solution order of accuracy in space and time in the presence of boundary conditions.Following, analytic solutions for the scalar advection equations are used to verify the solution order of accuracy in space and time, and stability.Finally, the scalar advection and scalar diffusion equations are joined together, and again verified for solution order of accuracy.

A. Potential Flow for a Circle in a Channel
The approximate projection operator P(u) is demonstrated to be both fourth-order accurate and stable by solving for bounded potential flow.A geometry is generated by creating an embedded boundary grid with a circle of radius 0.15 placed in the center of a square domain of unit length, shown with a cell size of ℎ = 1/128 in Fig. 3.The left boundary is specified as an inlet condition with unit velocity in the x-direction, the right boundary is specified by an outlet condition, and the remaining boundaries are slip-walls with zero normal velocity.The initial solution is a potential flow in an infinitely long domain evaluated for cell averages.
Stability of the approximate projection operator is demonstrated by showing that repeatedly applications of the projection on a velocity field reduces the divergence towards zero.The velocity divergence Du and pressure gradient G from the projection are computed on a grid with cells size ℎ = 1/128 for 100 iterations and plotted in Fig. 4.These quantities are shown to strictly decrease, showing stability of the high-order approximate projection.
Convergence tests are performed to ensure the spacial discretization achieves the targeted fourth-order accuracy.The projection is applied once for grids of decreasing refinement, and the divergence field Du is used to evaluate the solution error.The finest level is chosen with cell size ℎ = 1/128, and subsequent coarser levels each double ℎ.The  1 ,  2 , and  ∞ errors and convergence rates are calculated and listed in Table 2.The tabulated convergence rates approach the expected fourth-order accuracy, but fail to maintain it consistently.Error in the solution is found to be focused around both the regular and embedded boundaries.The failure to reach consistent fourth-order accuracy is attributed to using an initial condition which does not exactly match the boundary conditions, but further investigation is required.

B. Manufactured Solution for Diffusion inside a Circle
To demonstrate the accuracy of the diffusion equation    = L, a manufactured solution is employed.The algorithm targets fourth-order in space and time using the high-order stencils described in section IV and the ImEx scheme in section III.A circular domain is created of radius 0.3 centered about the point (0.5, 0.5) The manufactured solution is defined on this domain by where  = 0.3 matches the domain radius, and x 0 = (0.5, 0.5) to match the domain center.The embedded boundaries are specified by Dirichlet conditions with values determined by the manufactured solution.Following the method of manufactured solutions, a source term is added to balance the equation.The source term is evaluated explicitly in time, while the Laplacian term is evaluated implicitly.The solution is initialized to the fourth-order cell averaged form of the manufactured solution at time 0.125, and uses a dissipation rate of  = 1.On the finest level, of grid size ℎ = 1/128, a time step of  = 0.1 is taken to advance the solution forward for 128 steps.Subsequent coarser levels double the cell size and time step, while halving the time steps to reach the same end time.Using the chosen exact solution in Eq. (34), the  1 ,  2 , and  ∞ errors and convergence rates are calculated and compiled in Table 2.The expected fourth-order accuracy is demonstrated in all error norms, verifying the algorithm for scalar diffusion.

C. Advection inside a Periodic Channel
To show verification of the scalar advection equation    = −A(), where in this case A = ∇ • (u), comparison to an analytic solution is used.Analytic solutions to scalar advection are defined purely by shifts in the initial profile as a function of time.The algorithm targets fourth-order in space and time using the high-order stencils described in section IV and the explicit portion of the ImEx scheme in section III.The channel geometry is generated in a rectangle of unit height.The left plane is defined by the normal vector (−4, 1) and intersects the point (0.5, 0.5).The right plane uses the same normal vector, but intersects the point (1, 0.5).The top and bottom boundaries are periodic, and the left and right boundaries are slip walls.The initial and exact solution is a Gaussian profile is defined by The profile in Eq. ( 35), and shown in Fig. 5, is initialized along the left boundary about the point x 0 = (0.5, 0.5) with  = 1 and  2 = 1 − 2, with velocity vector u set to unit magnitude aligned with the channel.In the finest case, the grid size is ℎ = 1/256, and the time step corresponds to a CFL number of 1 based on the grid spacing.Each coarser level has the grid size and time step doubled from the finest.The profile is advected for a total of one unit in time, after which the solution should return to the initial profile.Using the exact solution,  1 ,  2 , and  ∞ errors and convergence rates are calculated and compiled in Table 3.The expected fourth-order accuracy is observed in all error norms, even in small cells, verifying the algorithm for scalar advection.Of particular note, this solution achieves a stable high-order explicit advection scheme at CFL of 1 without the use of cell merging or redistribution.This is a significant achievement for embedded boundary methods, where small cell problems are typically problematic for explicit advection schemes.

D. Advection-Diffusion for Solid Body Rotation inside a Circle
Prior results in sections VI.B and VI.C verify scalar advection and diffusion separately.Approaching the incompressible Navier-Stokes equations, these two equations must be coupled and the resulting algorithm verified.The scalar advection-diffusion equation is given by    = −A() + L.Since few analytic solutions are available, the Richardson extrapolation method is used with a "best" solution to compare against.The algorithm targets fourth-order accuracy in space and time using the high-order stencils described in section IV and the ImEx scheme in section III where the advection term is solved explicitly and the diffusion term implicitly.A circular geometry centered about (0.5, 0.5) with radius of 0.3183 is specified, and zero valued Neumann boundary conditions are enforced.The initial profile is in defined by the Gaussian profile in Eq. ( 35), and shown in Fig. 6, is initialized along the left most boundary centered about the point x 0 = (0.1817, 0.5) with  = 1 and  2 = 1e-2.The velocity field for the solid body rotation is specified as u 1 = (x 1 − 0.5) u 2 = −(x 0 − 0.5).The rotation rate  and diffusion rate  are both with unit values.At the finest case, the grid size is ℎ = 1/128, and the time step is set corresponding to a CFL number of 1.Each coarser level has the grid size and time step doubled from the finest.The profile is advanced in time for 3.14159 time units, or half of a rotation by the velocity field.In Fig. 6 the initial and final solution profiles are shown.It is seen that the initial profile diffuses while it travels along the boundary.Using the Richardson interpolated solution at resolution ℎ = 1/256 as a reference, the  1 ,  2 , and  ∞ errors and convergence rates are calculated and compiled in Table 4.The expected fourth-order accuracy is shown for in all error norms, verifying the algorithm for scalar advection-diffusion.Conservation of the scalar with Neumann boundary conditions can be confirmed near round-off.As in the pure advection case, it is particularly noteworthy that this algorithm achieves a stable high-order semi-explicit scheme at a CFL of 1 without the use of cell merging or redistribution.

VII. Conclusions and Future Work
In this work, we demonstrate that our embedded boundary algorithm is fourth-order accurate.Additionally, our algorithm is shown to be stable when encountering small cells, without cell merging, redistribution, or grid remediation to avoid small cells.These results demonstrate the feasibility of high-order embedded boundary methods, and provide confidence that the developed code can be extended to correctly solve engineering problems.
The immediate next step is to extend this algorithm for the incompressible Navier-Stokes equations.The primary challenge of this will be implementing and testing the non-linear advection term, and then coupling that with the existing diffusion and projection methods.Although this work only demonstrates two-dimensional results, the algorithm should be dimension independent and we plan to verify it next in three-dimensions.We also hope to extend this higher-order algorithm to the compressible Navier-Stokes equations, building on second-order methods previously published [3].This introduces a number of additional challenges, such as computing more complicated non-linear upwind fluxes as well as introducing limiters into weighted least-squares stencils.

Fig. 1
Fig. 1 Illustration of notation used for embedded boundary methods.The shaded region lies outside the problem domain of interest.

Fig. 2
Fig. 2 An illustration of stencil neighborhoods for flux construction.The stencils are for the red highlighted faces, with cells numbered according to the neighbor distance from the face.

( a )Fig. 3
Fig. 3 Potential flow around a circle in a channel, where the left boundary is a uniform inlet, the right an outlet, and all other boundaries slip walls.The streamlines are shown in black, and the contours plot the velocity magnitude.
Error norms of the pressure gradient.

Fig. 4 The
Fig. 4 The  1 ,  2 , and  ∞ norms from repeated projections of the channel bounded potential flow over a cylinder of grid size ℎ = 1/128.

( a )
Solution  of the advected Gaussian profile after one periodic cycle in the angled channel, which is visually indistinguishable from the initial condition.(b) A close up view of the angled channel mesh and cell volume fraction .

Fig. 5
Fig. 5 Advection of a Gaussian profile in an angled channel.The channel is periodic in the vertical direction, and imposes slip walls on the left and right boundaries.

( a )Fig. 6
Fig. 6 Advection-diffusion of a Gaussian profile in a counter-clockwise rotating velocity field.The boundary conditions are specified as zero valued Neumann boundaries.