Merge pull request #79 from FourierFlows/FixTracerAdvDiff

Fix TracerAdvDiff module

Author: navidcy

.gitignore

@@ -11,3 +11,4 @@
 
 docs/build/
 docs/site/
+coverage/

docs/src/modules/traceradvdiff.md

@@ -3,13 +3,24 @@
 ### Basic Equations
 
 This module solves the advection diffusion equation for a passive tracer
-concentration $c(x,y,t)$ in two-dimensions:
+concentration $c(x, y, t)$ in two-dimensions:
 
-$$\partial_t c + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} c = \kappa \nabla^2 c\ ,$$
+$$\partial_t c + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} c = \underbrace{\eta \partial_x^2 c + \kappa \partial_y^2 c}_{\textrm{diffusivity}} + \underbrace{\kappa_h (-1)^{n_{h}} \nabla^{2n_{h}}c}_{\textrm{hyper-diffusivity}}\ ,$$
+
+where $\boldsymbol{u} = (u,v)$ is the two-dimensional advecting flow, $\eta$ the $x$-diffusivity and $\kappa$ is the $y$-diffusivity. If $\eta$ is not defined then the code uses isotropic diffusivity, i.e., $\eta \partial_x^2 c + \kappa \partial_y^2 c\mapsto\kappa\nabla^2$. The advecting flow could be either compressible or incompressible. 
 
-where $\mathbf{u} = (u,v)$ is the two-dimensional advecting velocity and $\kappa$
-is the diffusivity.
 
 ### Implementation
 
-Coming soon.
+The equation is time-stepped forward in Fourier space:
+
+$$\partial_t \widehat{c} = - \widehat{\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} c} - \left[ (\eta k_x^2 + \kappa k_y^2) +\kappa_h k^{2\nu_h} \right]\widehat{c}\ .$$
+
+Thus:
+
+```math
+\begin{align*}
+\mathcal{L} &= -\eta k_x^2 - \kappa k_y^2 - \kappa_h k^{2\nu_h}\ , \\
+\mathcal{N}(\widehat{c}) &= - \mathrm{FFT}(u \partial_x c + \upsilon \partial_y c)\ .
+\end{align*}
+```

examples/cudaexamples/cu_randomdecay.jl

@@ -15,7 +15,7 @@ TwoDTurb.set_q!(prob, q0)
 fig = figure(); tic()
 for i = 1:10
   stepforward!(prob, nt)
-  TwoDTurb.updatevars!(prob)  
+  TwoDTurb.updatevars!(prob)
 
   cfl = maximum(prob.vars.U)*prob.grid.dx/prob.ts.dt
   @printf("step: %04d, t: %6.1f, cfl: %.2f, ", prob.step, prob.t, cfl)

examples/traceradvdiff/cellularflow.jl

@@ -0,0 +1,66 @@
+using FourierFlows, PyPlot, JLD2
+
+import FourierFlows.TracerAdvDiff
+
+# Numerical parameters and time-stepping parameters
+nx  = 128         # 2D resolution = nx^2
+stepper = "RK4"   # timestepper
+dt  = 0.02        # timestep
+nsubs  = 200      # number of time-steps for plotting
+nsteps = 4nsubs   # total number of time-steps (must be multiple of nsubs)
+
+# Physical parameters
+Lx  = 2π      # domain size
+kap = 0.002   # diffusivity
+
+gr = TwoDGrid(nx, Lx)
+X, Y = gr.X, gr.Y
+
+# streamfunction (for plotting) and (u,v) flow field
+psiampl = 0.2
+m, n = 1, 1
+psiin = @. psiampl * cos(m*X) * cos(n*Y)
+uvel(x, y) = +psiampl * n * cos(m*x) * sin(n*y)
+vvel(x, y) = -psiampl * m * sin(m*x) * cos(n*y)
+
+prob = TracerAdvDiff.ConstDiffProblem(; steadyflow=true,
+    nx=nx, Lx=Lx, kap=kap, u=uvel, v=vvel, dt=dt, stepper=stepper)
+
+s, v, p, g, eq, ts = prob.state, prob.vars, prob.params, prob.grid, prob.eqn, prob.ts;
+
+# Initial condition c0 = c(x, y, t=0)
+amplc0, sigc0 = 0.1, 0.1
+c0func(x, y) = amplc0*exp(-(x^2+y^2)/(2sigc0^2))
+c0 = c0func.(g.X-1.2, g.Y)
+
+TracerAdvDiff.set_c!(prob, c0)
+
+"Plot the concentration field and the (u, v) streamlines."
+function plotoutput(prob, fig, axs; drawcolorbar=false)
+  s, v, p, g = prob.state, prob.vars, prob.params, prob.grid
+  t = round(prob.state.t, 2)
+
+  TracerAdvDiff.updatevars!(prob)
+
+  cla()
+  pcolormesh(g.X, g.Y, v.c)
+
+  if drawcolorbar; colorbar(); end
+
+  contour(g.X, g.Y, psiin, 15, colors="k", linewidths=0.7)
+  axis("equal")
+  axis("square")
+  title("shading: \$c(x, y, t= $t )\$, contours: \$\\psi(x, y)\$")
+
+  pause(0.001)
+
+  nothing
+end
+
+fig, axs = subplots(ncols=1, nrows=1, figsize=(8, 8))
+plotoutput(prob, fig, axs; drawcolorbar=true)
+
+while prob.step < nsteps
+  stepforward!(prob, nsubs)
+  plotoutput(prob, fig, axs)
+end