diff --git a/Project.toml b/Project.toml index d0b6bde..c3bb500 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,7 @@ license = "MIT" authors = ["Navid C. Constantinou ", "Gregory L. Wagner "] documentation = "https://fourierflows.github.io/PassiveTracerFlowsDocumentation/dev/" repository = "https://github.com/FourierFlows/PassiveTracerFlows.jl" -version = "0.6.1" +version = "0.7.0" [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" @@ -19,8 +19,8 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [compat] CUDA = "^1.3, ^2, ^3" -DocStringExtensions = "^0.8" -FourierFlows = "^0.6.17, ^0.7, 0.8, 0.9" +DocStringExtensions = "0.8, 0.9" +FourierFlows = "0.9.2" GeophysicalFlows = "0.14" JLD2 = "^0.1, ^0.2, ^0.3, ^0.4" Reexport = "^0.2, ^1" diff --git a/README.md b/README.md index 0ab24f4..35c4b2b 100644 --- a/README.md +++ b/README.md @@ -49,6 +49,6 @@ See `examples/` for example scripts. The code is citable via [zenodo](https://zenodo.org). Please cite as: -> Navid C. Constantinou. and Gregory L. Wagner (2022). FourierFlows/PassiveTracerFlows.jl: PassiveTracerFlows v0.6.1 (Version v0.6.1). Zenodo. [https://doi.org/10.5281/zenodo.2535983](https://doi.org/10.5281/zenodo.2535983) +> Navid C. Constantinou, Josef Bisits, and Gregory L. Wagner (2022). FourierFlows/PassiveTracerFlows.jl: PassiveTracerFlows v0.7.0 (Version v0.7.0). Zenodo. [https://doi.org/10.5281/zenodo.2535983](https://doi.org/10.5281/zenodo.2535983) [FourierFlows.jl]: https://github.com/FourierFlows/FourierFlows.jl diff --git a/docs/make.jl b/docs/make.jl index c757116..7ee0138 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -17,6 +17,7 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples") const OUTPUT_DIR = joinpath(@__DIR__, "src/literated") examples = [ + "onedim_gaussiandiffusion.jl", "cellularflow.jl", "turbulent_advection-diffusion.jl" ] @@ -56,6 +57,7 @@ makedocs( "Home" => "index.md", "Examples" => Any[ "TracerAdvectionDiffusion" => Any[ + "literated/onedim_gaussiandiffusion.md", "literated/cellularflow.md", "literated/turbulent_advection-diffusion.md", ] diff --git a/examples/cellularflow.jl b/examples/cellularflow.jl index 353e32f..889d437 100644 --- a/examples/cellularflow.jl +++ b/examples/cellularflow.jl @@ -45,6 +45,8 @@ nothing # hide # ## Set up cellular flow # We create a two-dimensional grid to construct the cellular flow. Our cellular flow is derived # from a streamfunction ``ψ(x, y) = ψ₀ \cos(x) \cos(y)`` as ``(u, v) = (-∂_y ψ, ∂_x ψ)``. +# The cellular flow is then passed into the `TwoDAdvectingFlow` constructor with `steadyflow = true` +# to indicate that the flow is not time dependent. grid = TwoDGrid(n, L) ψ₀ = 0.2 @@ -54,13 +56,13 @@ mx, my = 1, 1 uvel(x, y) = ψ₀ * my * cos(mx * x) * sin(my * y) vvel(x, y) = -ψ₀ * mx * sin(mx * x) * cos(my * y) +advecting_flow = TwoDAdvectingFlow(; u = uvel, v = vvel, steadyflow = true) nothing # hide # ## Problem setup # We initialize a `Problem` by providing a set of keyword arguments. -prob = TracerAdvectionDiffusion.Problem(dev; nx=n, Lx=L, κ=κ, steadyflow=true, u=uvel, v=vvel, - dt=dt, stepper=stepper) +prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; nx=n, Lx=L, κ, dt, stepper) nothing # hide # and define some shortcuts diff --git a/examples/onedim_gaussiandiffusion.jl b/examples/onedim_gaussiandiffusion.jl new file mode 100644 index 0000000..91a19ae --- /dev/null +++ b/examples/onedim_gaussiandiffusion.jl @@ -0,0 +1,130 @@ +# # Advection-diffusion of tracer in one dimension +# +# This is an example demonstrating the advection-diffusion of a passive tracer in +# one dimension. +# +# ## Install dependencies +# +# First let's make sure we have all the required packages installed + +# ```julia +# using Pkg +# pkg.add(["PassiveTracerFlows", "Printf", "Plots", "JLD2"]) +# ``` +# +# ## Let's begin +# First load packages needed to run this example. +using PassiveTracerFlows, Plots, Printf, JLD2, LinearAlgebra + +# ## Choosing a device: CPU or GPU + +dev = CPU() # Device (CPU/GPU) +nothing # hide + + +# ## Numerical parameters and time-stepping parameters + + n = 128 # 2D resolution = n² +stepper = "RK4" # timestepper + dt = 0.02 # timestep + nsteps = 5000 # total number of time-steps +nothing # hide + + +# ## Physical parameters + +L = 2π # domain size +κ = 0.01 # diffusivity +nothing # hide + +# ## Flow +# We set a constant background flow and pass this to `OneDAdvectingFlow` with `steadyflow = true` to indicate the flow is not time dependent. +u(x) = 0.05 +advecting_flow = OneDAdvectingFlow(; u, steadyflow = true) + +# ## Problem setup +# We initialize a `Problem` by providing a set of keyword arguments. +prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; nx=n, Lx=L, κ, dt, stepper) +nothing # hide + +# and define some shortcuts +sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid +x = grid.x + +# ## Initial condition +# We advect-diffuse a concentration field that has an initial concentration set to Gaussian. +gaussian(x, σ) = exp(-x^2 / 2σ^2) + +amplitude, spread = 1, 0.15 +c₀ = [amplitude * gaussian(x[i], spread) for i in 1:grid.nx] + +TracerAdvectionDiffusion.set_c!(prob, c₀) +nothing #hide + +# ## Saving output +# We create the saved output using the `Output` function from `FourierFlows.jl` then +# save the concentration field using the `get_concentration` function every 50 timesteps. +function get_concentration(prob) + ldiv!(prob.vars.c, prob.grid.rfftplan, deepcopy(prob.sol)) + + return prob.vars.c +end + +output = Output(prob, "advection-diffusion1D.jld2", + (:concentration, get_concentration)) + +# By calling `saveproblem(output)` we save information that we will use for plotting later on. +saveproblem(output) + +# ## Stepping the problem forward +# Now we step the problem forward and save output every 50 timesteps. + +save_frequency = 50 # frequency at which output is saved + +startwalltime = time() +while clock.step <= nsteps + if clock.step % save_frequency == 0 + saveoutput(output) + log = @sprintf("Output saved, step: %04d, t: %.2f, walltime: %.2f min", + clock.step, clock.t, (time()-startwalltime) / 60) + + println(log) + end + + stepforward!(prob) +end + +# ## Visualising the output +# We load the `.jld2` file and create a timeseries of the concentration field +file = jldopen(output.path) + +iterations = parse.(Int, keys(file["snapshots/t"])) + +t = [file["snapshots/t/$i"] for i ∈ iterations] +c = [file["snapshots/concentration/$i"] for i ∈ iterations] +nothing # hide + +# Set up the plotting arguments and look at the initial concentration. +x, Lx = file["grid/x"], file["grid/Lx"] + +plot_args = (xlabel = "x", + ylabel = "c", + framestyle = :box, + xlims = (-Lx/2, Lx/2), + ylims = (0, maximum(c[1])), + legend = :false, + color = :balance) + +p = plot(x, Array(c[1]); + title = "concentration, t = " * @sprintf("%.2f", t[1]), + plot_args...) +nothing # hide + +# Now, we create a movie of the tracer concentration being advected and diffused. + +anim = @animate for i ∈ 1:length(t) + p[1][:title] = "concentration, t = " * @sprintf("%.2f", t[i]) + p[1][1][:y] = Array(c[i]) +end + +mp4(anim, "1D_advection-diffusion.mp4", fps = 12) diff --git a/src/traceradvectiondiffusion.jl b/src/traceradvectiondiffusion.jl index 78d8944..25e4a6f 100644 --- a/src/traceradvectiondiffusion.jl +++ b/src/traceradvectiondiffusion.jl @@ -3,7 +3,9 @@ module TracerAdvectionDiffusion export Problem, set_c!, - updatevars! + updatevars!, + OneDAdvectingFlow, + TwoDAdvectingFlow using CUDA, @@ -17,45 +19,127 @@ using GeophysicalFlows.MultiLayerQG: SingleLayerParams, TwoLayerParams, numberof import LinearAlgebra: mul!, ldiv! import GeophysicalFlows.MultiLayerQG +# -- +# AdvectingFlows +# -- + +"Abstract super type for an advecting flow." +abstract type AbstractAdvectingFlow end + +""" + struct OneDAdvectingFlow <: AbstractAdvectingFlow + +A struct containing the advecting flow for a one dimensional `TracerAdvectionDiffusion.Problem`. +Included are + +$(TYPEDFIELDS) +""" +struct OneDAdvectingFlow <: AbstractAdvectingFlow + "function for the x-component of the advecting flow" + u :: Function + "boolean declaring whether or not the flow is steady (i.e., not time dependent)" + steadyflow :: Bool +end + +noflow(args...) = 0.0 # used as defaults for u, v functions in AdvectingFlow constructors + +""" + OneDAdvectingFlow(; u=noflow, steadyflow=true) + +Return a `OneDAdvectingFlow`. By default, there is no advecting flow `u=noflow` hence +`steadyflow=true`. +""" +OneDAdvectingFlow(; u=noflow, steadyflow=true) = OneDAdvectingFlow(u, steadyflow) + +""" + struct TwoDAdvectingFlow <: AbstractAdvectingFlow + +A struct containing the advecting flow for a two dimensional `TracerAdvectionDiffusion.Problem`. +Included are + +$(TYPEDFIELDS) +""" +struct TwoDAdvectingFlow <: AbstractAdvectingFlow + "function for the x-component of the advecting flow" + u :: Function + "function for the y-component of the advecting flow" + v :: Function + "boolean declaring whether or not the flow is steady (i.e., not time dependent)" + steadyflow :: Bool +end + +""" + TwoDAdvectingFlow(; u=noflow, v=noflow, steadyflow=true) + +Constructor for the `TwoDAdvectingFlow`. The default function for the advecting flow components is `noflow` +hence `steadyflow=true`. +""" +TwoDAdvectingFlow(; u=noflow, v=noflow, steadyflow=true) = TwoDAdvectingFlow(u, v, steadyflow) + # -- # Problems # -- """ - Problem(; parameters...) + Problem(dev, advecting_flow; parameters...) -Construct a constant diffusivity problem with steady or time-varying flow. +Construct a constant diffusivity problem with steady or time-varying `advecting_flow` on device `dev`. +The dimensionality of the problem is inferred from the `advecting_flow`: +* `advecting_flow::OneDAdvectingFlow` for 1D advection-diffusion problem, +* `advecting_flow::TwoDAdvectingFlow` for 2D advection-diffusion problem. """ -noflow(args...) = 0.0 # used as defaults for u, v functions in Problem() +function Problem(dev, advecting_flow::OneDAdvectingFlow; + nx = 128, + Lx = 2π, + κ = 0.1, + dt = 0.01, + stepper = "RK4", + T = Float64 + ) -function Problem(dev; - nx = 128, - Lx = 2π, - ny = nx, - Ly = Lx, - κ = 0.1, - η = κ, - u = noflow, - v = noflow, - dt = 0.01, - stepper = "RK4", - steadyflow = false, - T = Float64 - ) + grid = OneDGrid(dev, nx, Lx; T) - grid = TwoDGrid(dev, nx, Lx, ny, Ly; T=T) - params = steadyflow==true ? ConstDiffSteadyFlowParams(η, κ, u, v, grid) : ConstDiffParams(η, κ, u, v) - vars = Vars(dev, grid) + params = advecting_flow.steadyflow==true ? + ConstDiffSteadyFlowParams(κ, advecting_flow.u, grid::OneDGrid) : + ConstDiffTimeVaryingFlowParams(κ, advecting_flow.u) + + vars = Vars(dev, grid; T) + + equation = Equation(dev, params, grid) + + return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev) +end + +function Problem(dev, advecting_flow::TwoDAdvectingFlow; + nx = 128, + Lx = 2π, + ny = nx, + Ly = Lx, + κ = 0.1, + η = κ, + dt = 0.01, + stepper = "RK4", + T = Float64 + ) + + grid = TwoDGrid(dev, nx, Lx, ny, Ly; T) + + params = advecting_flow.steadyflow==true ? + ConstDiffSteadyFlowParams(η, κ, advecting_flow.u, advecting_flow.v, grid::TwoDGrid) : + ConstDiffTimeVaryingFlowParams(η, κ, advecting_flow.u, advecting_flow.v) + + vars = Vars(dev, grid; T) + equation = Equation(dev, params, grid) return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev) end """ - Problem(dev, MQGprob; parameters...) + Problem(dev, MQGprob::FourierFlows.Problem; parameters...) Construct a constant diffusivity problem on device `dev` using the flow from a -`GeophysicalFlows.MultiLayerQG` problem. +`GeophysicalFlows.MultiLayerQG` problem as the advecting flow. """ function Problem(dev, MQGprob::FourierFlows.Problem; κ = 0.1, @@ -73,7 +157,7 @@ function Problem(dev, MQGprob::FourierFlows.Problem; step_until!(MQGprob, tracer_release_time) end - params = TurbulentFlowParams(η, κ, tracer_release_time, MQGprob) + params = ConstDiffTurbulentFlowParams(η, κ, tracer_release_time, MQGprob) vars = Vars(dev, grid, MQGprob) equation = Equation(dev, params, grid) @@ -86,19 +170,38 @@ end # Params # -- -abstract type AbstractConstDiffParams <: AbstractParams end +abstract type AbstractTimeVaryingFlowParams <: AbstractParams end abstract type AbstractSteadyFlowParams <: AbstractParams end abstract type AbstractTurbulentFlowParams <: AbstractParams end """ - struct ConstDiffParams{T} <: AbstractConstDiffParams + struct ConstDiffTimeVaryingFlowParams1D{T} <: AbstractTimeVaryingFlowParams -A struct containing the parameters for a constant diffusivity problem with time-varying flow. -Included are: +A struct containing the parameters for a constant diffusivity problem with time-varying flow in one +dimension. Included are: $(TYPEDFIELDS) """ -struct ConstDiffParams{T} <: AbstractConstDiffParams +struct ConstDiffTimeVaryingFlowParams1D{T} <: AbstractTimeVaryingFlowParams + "isotropic horizontal diffusivity coefficient" + κ :: T + "isotropic hyperdiffusivity coefficient" + κh :: T + "isotropic hyperdiffusivity order" + nκh :: Int + "function returning the x-component of advecting flow" + u :: Function +end + +""" + struct ConstDiffTimeVaryingFlowParams2D{T} <: AbstractTimeVaryingFlowParams + +A struct containing the parameters for a constant diffusivity problem with time-varying flow in two +dimensions. Included are: + +$(TYPEDFIELDS) +""" +struct ConstDiffTimeVaryingFlowParams2D{T} <: AbstractTimeVaryingFlowParams "isotropic horizontal diffusivity coefficient" η :: T "isotropic vertical diffusivity coefficient" @@ -114,21 +217,42 @@ struct ConstDiffParams{T} <: AbstractConstDiffParams end """ - ConstDiffParams(η, κ, u, v) + ConstDiffTimeVaryingFlowParams(κ, u) + ConstDiffTimeVaryingFlowParams(η, κ, u, v) The constructor for the `params` struct for constant diffusivity problem and time-varying flow. """ -ConstDiffParams(η, κ, u, v) = ConstDiffParams(η, κ, 0η, 0, u, v) +ConstDiffTimeVaryingFlowParams(κ, u) = ConstDiffTimeVaryingFlowParams1D(κ, 0κ, 0, u) +ConstDiffTimeVaryingFlowParams(η, κ, u, v) = ConstDiffTimeVaryingFlowParams2D(η, κ, 0η, 0, u, v) + +""" + struct ConstDiffSteadyFlowParams1D{T} <: AbstractSteadyFlowParams + +A struct containing the parameters for a constant diffusivity problem with steady flow in one dimensions. +Included are: + +$(TYPEDFIELDS) +""" +struct ConstDiffSteadyFlowParams1D{T, A} <: AbstractSteadyFlowParams + "isotropic horizontal diffusivity coefficient" + κ :: T + "isotropic hyperdiffusivity coefficient" + κh :: T + "isotropic hyperdiffusivity order" + nκh :: Int + "x-component of advecting flow" + u :: A +end """ - struct ConstDiffParams{T} <: AbstractConstDiffParams + struct ConstDiffSteadyFlowParams2D{T} <: AbstractSteadyFlowParams -A struct containing the parameters for a constant diffusivity problem with steady flow. +A struct containing the parameters for a constant diffusivity problem with steady flow in two dimensions. Included are: $(TYPEDFIELDS) """ -struct ConstDiffSteadyFlowParams{T, A} <: AbstractSteadyFlowParams +struct ConstDiffSteadyFlowParams2D{T, A} <: AbstractSteadyFlowParams "isotropic horizontal diffusivity coefficient" η :: T "isotropic vertical diffusivity coefficient" @@ -144,30 +268,41 @@ struct ConstDiffSteadyFlowParams{T, A} <: AbstractSteadyFlowParams end """ - ConstDiffSteadyFlowParams(η, κ, κh, nκh, u::Function, v::Function, grid) - ConstDiffSteadyFlowParams(η, κ, u, v, grid) + ConstDiffSteadyFlowParams(κ, κh, nκh, u::Function, grid::OneDGrid) + ConstDiffSteadyFlowParams(κ, u, grid::OneDGrid) + ConstDiffSteadyFlowParams(η, κ, κh, nκh, u::Function, v::Function, grid::TwoDGrid) + ConstDiffSteadyFlowParams(η, κ, u, v, grid::TwoDGrid) -The constructor for the `params` struct for constant diffusivity problem and steady flow. +Return the parameters `params` for a constant diffusivity problem and steady flow. """ -function ConstDiffSteadyFlowParams(η, κ, κh, nκh, u::Function, v::Function, grid) +function ConstDiffSteadyFlowParams(κ, κh, nκh, u::Function, grid::OneDGrid) + x = gridpoints(grid) + ugrid = u.(x) + + return ConstDiffSteadyFlowParams1D(κ, κh, nκh, ugrid) + end + + ConstDiffSteadyFlowParams(κ, u, grid::OneDGrid) = ConstDiffSteadyFlowParams(κ, 0κ, 0, u, grid) + +function ConstDiffSteadyFlowParams(η, κ, κh, nκh, u::Function, v::Function, grid::TwoDGrid) x, y = gridpoints(grid) ugrid = u.(x, y) vgrid = v.(x, y) - return ConstDiffSteadyFlowParams(η, κ, κh, nκh, ugrid, vgrid) + return ConstDiffSteadyFlowParams2D(η, κ, κh, nκh, ugrid, vgrid) end -ConstDiffSteadyFlowParams(η, κ, u, v, grid) = ConstDiffSteadyFlowParams(η, κ, 0η, 0, u, v, grid) +ConstDiffSteadyFlowParams(η, κ, u, v, grid::TwoDGrid) = ConstDiffSteadyFlowParams(η, κ, 0η, 0, u, v, grid) """ - struct TurbulentFlowParams{T} <: AbstractTurbulentFlowParams + struct ConstDiffTurbulentFlowParams{T} <: AbstractTurbulentFlowParams A struct containing the parameters for a constant diffusivity problem with flow obtained from a `GeophysicalFlows.MultiLayerQG` problem. $(TYPEDFIELDS) """ -struct TurbulentFlowParams{T} <: AbstractTurbulentFlowParams +struct ConstDiffTurbulentFlowParams{T} <: AbstractTurbulentFlowParams "isotropic horizontal diffusivity coefficient" η :: T "isotropic vertical diffusivity coefficient" @@ -185,17 +320,17 @@ struct TurbulentFlowParams{T} <: AbstractTurbulentFlowParams end """ - TurbulentFlowParams(η, κ, tracer_release_time, MQGprob) + ConstDiffTurbulentFlowParams(η, κ, tracer_release_time, MQGprob) -The constructor for the `params` for a constant diffusivity problem with flow obtained +Return the parameters `params` for a constant diffusivity problem with flow obtained from a `GeophysicalFlows.MultiLayerQG` problem. """ -function TurbulentFlowParams(η, κ, tracer_release_time, MQGprob) +function ConstDiffTurbulentFlowParams(η, κ, tracer_release_time, MQGprob) nlayers = numberoflayers(MQGprob.params) MultiLayerQG.updatevars!(MQGprob) - return TurbulentFlowParams(η, κ, 0η, 0, nlayers, tracer_release_time, MQGprob) + return ConstDiffTurbulentFlowParams(η, κ, 0η, 0, nlayers, tracer_release_time, MQGprob) end # -- @@ -207,21 +342,35 @@ end Return the equation for constant diffusivity problem with `params` and `grid` on device `dev`. """ -function Equation(dev, params::ConstDiffParams, grid) +function Equation(dev, params::ConstDiffTimeVaryingFlowParams1D, grid) + L = zeros(dev, eltype(grid), (grid.nkr)) + @. L = - params.κ * grid.kr^2 - params.κh * (grid.kr^2)^params.nκh + + return FourierFlows.Equation(L, calcN!, grid) +end + +function Equation(dev, params::ConstDiffTimeVaryingFlowParams2D, grid) L = zeros(dev, eltype(grid), (grid.nkr, grid.nl)) @. L = - params.η * grid.kr^2 - params.κ * grid.l^2 - params.κh * grid.Krsq^params.nκh return FourierFlows.Equation(L, calcN!, grid) end -function Equation(dev, params::ConstDiffSteadyFlowParams, grid) +function Equation(dev, params::ConstDiffSteadyFlowParams1D, grid) + L = zeros(dev, eltype(grid), (grid.nkr)) + @. L = - params.κ * grid.kr^2 - params.κh * (grid.kr^2)^params.nκh + + return FourierFlows.Equation(L, calcN_steadyflow!, grid) +end + +function Equation(dev, params::ConstDiffSteadyFlowParams2D, grid) L = zeros(dev, eltype(grid), (grid.nkr, grid.nl)) @. L = - params.η * grid.kr^2 - params.κ * grid.l^2 - params.κh * grid.Krsq^params.nκh return FourierFlows.Equation(L, calcN_steadyflow!, grid) end -function Equation(dev, params::TurbulentFlowParams, grid) +function Equation(dev, params::ConstDiffTurbulentFlowParams, grid) L = zeros(dev, eltype(grid), (grid.nkr, grid.nl, params.nlayers)) for j in 1:params.nlayers @@ -236,50 +385,74 @@ end # -- """ - struct Vars{Aphys, Atrans} <: AbstractVars + struct Vars1D{Aphys, Atrans} <: AbstractVars -The variables for TracerAdvectionDiffussion problems. +The variables for a 1D `TracerAdvectionDiffussion` problem. $(FIELDS) """ -struct Vars{Aphys, Atrans} <: AbstractVars +struct Vars1D{Aphys, Atrans} <: AbstractVars "tracer concentration" c :: Aphys - "tracer concentration x-derivative, ∂c/∂x" + "tracer concentration ``x``-derivative, ``∂c/∂x``" cx :: Aphys - "tracer concentration y-derivative, ∂c/∂y" + "Fourier transform of tracer concentration" + ch :: Atrans + "Fourier transform of tracer concentration ``x``-derivative, ``∂c/∂x``" + cxh :: Atrans +end + +""" + struct Vars2D{Aphys, Atrans} <: AbstractVars + +The variables for a 2D `TracerAdvectionDiffussion` problem. + +$(FIELDS) +""" +struct Vars2D{Aphys, Atrans} <: AbstractVars + "tracer concentration" + c :: Aphys + "tracer concentration ``x``-derivative, ``∂c/∂x``" + cx :: Aphys + "tracer concentration ``y``-derivative, ``∂c/∂y``" cy :: Aphys "Fourier transform of tracer concentration" ch :: Atrans - "Fourier transform of tracer concentration x-derivative, ∂c/∂x" + "Fourier transform of tracer concentration ``x``-derivative, ``∂c/∂x``" cxh :: Atrans - "Fourier transform of tracer concentration y-derivative, ∂c/∂y" + "Fourier transform of tracer concentration ``y``-derivative, ``∂c/∂y``" cyh :: Atrans end """ - Vars(dev, grid) + Vars(dev, grid; T=Float64) Return the variables `vars` for a constant diffusivity problem on `grid` and device `dev`. """ -function Vars(::Dev, grid::AbstractGrid{T}) where {Dev, T} +function Vars(::Dev, grid::OneDGrid; T=Float64) where Dev + @devzeros Dev T (grid.nx) c cx + @devzeros Dev Complex{T} (grid.nkr) ch cxh + + return Vars1D(c, cx, ch, cxh) +end +function Vars(::Dev, grid::TwoDGrid; T=Float64) where Dev @devzeros Dev T (grid.nx, grid.ny) c cx cy @devzeros Dev Complex{T} (grid.nkr, grid.nl) ch cxh cyh - return Vars(c, cx, cy, ch, cxh, cyh) + return Vars2D(c, cx, cy, ch, cxh, cyh) end function Vars(dev::Dev, grid::AbstractGrid{T}, MQGprob::FourierFlows.Problem) where {Dev, T} nlayers = numberoflayers(MQGprob.params) if nlayers == 1 - return Vars(dev, grid) + return Vars(dev, grid; T=T) else @devzeros Dev T (grid.nx, grid.ny, nlayers) c cx cy @devzeros Dev Complex{T} (grid.nkr, grid.nl, nlayers) ch cxh cyh - return Vars(c, cx, cy, ch, cxh, cyh) + return Vars2D(c, cx, cy, ch, cxh, cyh) end end @@ -294,7 +467,19 @@ end Calculate the advective terms for a tracer equation with constant diffusivity and time-varying flow. """ -function calcN!(N, sol, t, clock, vars, params::AbstractConstDiffParams, grid) +function calcN!(N, sol, t, clock, vars, params::AbstractTimeVaryingFlowParams, grid::OneDGrid) + @. vars.cxh = im * grid.kr * sol + + ldiv!(vars.cx, grid.rfftplan, vars.cxh) # destroys vars.cxh when using fftw + + x = grid.x + @. vars.cx = -params.u(x, clock.t) * vars.cx # copies over vars.cx so vars.cx = N in physical space + mul!(N, grid.rfftplan, vars.cx) + + return nothing +end + +function calcN!(N, sol, t, clock, vars, params::AbstractTimeVaryingFlowParams, grid::TwoDGrid) @. vars.cxh = im * grid.kr * sol @. vars.cyh = im * grid.l * sol @@ -314,7 +499,18 @@ end Calculate the advective terms for a tracer equation with constant diffusivity and time-constant flow. """ -function calcN_steadyflow!(N, sol, t, clock, vars, params::AbstractSteadyFlowParams, grid) +function calcN_steadyflow!(N, sol, t, clock, vars, params::AbstractSteadyFlowParams, grid::OneDGrid) + @. vars.cxh = im * grid.kr * sol + + ldiv!(vars.cx, grid.rfftplan, vars.cxh) # destroys vars.cxh when using fftw + + @. vars.cx = -params.u * vars.cx # copies over vars.cx so vars.cx = N in physical space + mul!(N, grid.rfftplan, vars.cx) + + return nothing +end + +function calcN_steadyflow!(N, sol, t, clock, vars, params::AbstractSteadyFlowParams, grid::TwoDGrid) @. vars.cxh = im * grid.kr * sol @. vars.cyh = im * grid.l * sol @@ -368,7 +564,7 @@ end """ updatevars!(params::AbstractTurbulentFlowParams, vars, grid, sol) -Update the `vars`` on the `grid` with the solution in `sol` for a problem `prob` +Update the `vars` on the `grid` with the solution in `sol` for a problem `prob` that is being advected by a turbulent flow. """ function updatevars!(params::AbstractTurbulentFlowParams, vars, grid, sol) @@ -382,11 +578,11 @@ end updatevars!(prob) = updatevars!(prob.params, prob.vars, prob.grid, prob.sol) """ - set_c!(sol, params::Union{AbstractConstDiffParams, AbstractSteadyFlowParams}, grid, c) + set_c!(sol, params::Union{AbstractTimeVaryingFlowParams, AbstractSteadyFlowParams}, grid, c) Set the solution `sol` as the transform of `c` and update variables `vars`. """ -function set_c!(sol, params::Union{AbstractConstDiffParams, AbstractSteadyFlowParams}, vars, grid, c) +function set_c!(sol, params::Union{AbstractTimeVaryingFlowParams, AbstractSteadyFlowParams}, vars, grid, c) mul!(sol, grid.rfftplan, c) updatevars!(params, vars, grid, sol) diff --git a/test/runtests.jl b/test/runtests.jl index a06d735..b1c7c79 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,10 +25,18 @@ for dev in devices stepper = "RK4" dt, nsteps = 1e-2, 40 + @test test_constvel1D(stepper, dt, nsteps, dev) + dt, tfinal = 0.002, 0.1 + @test test_timedependentvel1D(stepper, dt, tfinal, dev) + dt, nsteps = 1e-2, 40 @test test_constvel(stepper, dt, nsteps, dev) dt, tfinal = 0.002, 0.1 @test test_timedependentvel(stepper, dt, tfinal, dev) dt, tfinal = 0.005, 0.1 + @test test_diffusion1D(stepper, dt, tfinal, dev; steadyflow=true) + dt, tfinal = 0.005, 0.1 + @test test_diffusion1D(stepper, dt, tfinal, dev; steadyflow=false) + dt, tfinal = 0.005, 0.1 @test test_diffusion(stepper, dt, tfinal, dev; steadyflow=true) dt, tfinal = 0.005, 0.1 @test test_diffusion(stepper, dt, tfinal, dev; steadyflow=false) diff --git a/test/test_traceradvectiondiffusion.jl b/test/test_traceradvectiondiffusion.jl index 7f69757..9b0ee84 100644 --- a/test/test_traceradvectiondiffusion.jl +++ b/test/test_traceradvectiondiffusion.jl @@ -1,3 +1,74 @@ +""" + test_constvel1D(; kwargs...) + +Advect a gaussian concentration `c0(x, t)` with a constant velocity flow +`u(x) = uvel` and compare the final state with +`cfinal = c0(x - uvel * tfinal)`. +""" +function test_constvel1D(stepper, dt, nsteps, dev::Device=CPU()) + + nx, Lx = 128, 2π + uvel = 0.05 + u(x) = uvel + advecting_flow = OneDAdvectingFlow(; u) + + prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; nx, Lx, κ=0.0, dt, stepper) + sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid + x = gridpoints(gr) + + σ = 0.1 + c0ampl = 0.1 + c0func(x) = c0ampl * exp(-x^2 / 2σ^2) + + c0 = c0func.(x) + tfinal = nsteps * dt + cfinal = @. c0func(x - uvel * tfinal) + + TracerAdvectionDiffusion.set_c!(prob, c0) + + stepforward!(prob, nsteps) + TracerAdvectionDiffusion.updatevars!(prob) + + return isapprox(cfinal, vs.c, rtol = gr.nx*nsteps*1e-12) +end + +""" + test_timedependenttvel1D(; kwargs...) + +Advect a gaussian concentration `c0(x, t)` with a time-varying velocity flow +`u(x, t) = uvel * sign(-t + tfinal/2)` and compares the final +state with `cfinal = c0(x - uvel * tfinal)`. +""" +function test_timedependentvel1D(stepper, dt, tfinal, dev::Device=CPU(); uvel = 0.05) + + nx, Lx = 128, 2π + nsteps = round(Int, tfinal/dt) + + if !isapprox(tfinal, nsteps*dt, rtol=rtol_traceradvectiondiffusion) + error("tfinal is not multiple of dt") + end + + u(x, t) = uvel * t + uvel * dt/2 + advecting_flow = OneDAdvectingFlow(; u, steadyflow = false) + + prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; nx, Lx, κ=0.0, dt, stepper) + sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid + x = gridpoints(gr) + + σ = 0.2 + c0func(x) = 0.1 * exp(-x^2 / 2σ^2) + c0 = @. c0func(x) + tfinal = nsteps * dt + cfinal = @. c0func(x - 0.5uvel * tfinal^2) + + TracerAdvectionDiffusion.set_c!(prob, c0) + + stepforward!(prob, nsteps) + TracerAdvectionDiffusion.updatevars!(prob) + + return isapprox(cfinal, vs.c, rtol = gr.nx*nsteps*1e-12) +end + """ test_constvel(; kwargs...) @@ -11,15 +82,16 @@ function test_constvel(stepper, dt, nsteps, dev::Device=CPU()) uvel, vvel = 0.2, 0.1 u(x, y) = uvel v(x, y) = vvel + advecting_flow = TwoDAdvectingFlow(; u, v) - prob = TracerAdvectionDiffusion.Problem(dev; nx, Lx, κ=0.0, u, v, dt, stepper, steadyflow=true) + prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; nx, Lx, κ=0.0, dt, stepper) sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid x, y = gridpoints(gr) σ = 0.1 c0ampl = 0.1 - c0func(x, y) = c0ampl * exp(-(x^2 + y^2) / (2σ^2)) + c0func(x, y) = c0ampl * exp(-(x^2 + y^2) / 2σ^2) c0 = c0func.(x, y) tfinal = nsteps * dt @@ -52,13 +124,14 @@ function test_timedependentvel(stepper, dt, tfinal, dev::Device=CPU(); uvel = 0. u(x, y, t) = uvel v(x, y, t) = αv * t + αv * dt/2 + advecting_flow = TwoDAdvectingFlow(; u, v, steadyflow = false) - prob = TracerAdvectionDiffusion.Problem(dev; nx, Lx, κ=0.0, u, v, dt, stepper) + prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; nx, Lx, κ=0.0, dt, stepper) sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid x, y = gridpoints(gr) σ = 0.2 - c0func(x, y) = 0.1 * exp(-(x^2 + y^2) / (2σ^2)) + c0func(x, y) = 0.1 * exp(-(x^2 + y^2) / 2σ^2) c0 = @. c0func(x, y) tfinal = nsteps * dt @@ -72,6 +145,45 @@ function test_timedependentvel(stepper, dt, tfinal, dev::Device=CPU(); uvel = 0. return isapprox(cfinal, vs.c, rtol = gr.nx*gr.ny*nsteps*1e-12) end +""" + test_diffusion1D(; kwargs...) + +Diffuses a gaussian concentration c0(x, t) and compares the final state with +the analytic solution of the heat equation, cfinal +""" +function test_diffusion1D(stepper, dt, tfinal, dev::Device=CPU(); steadyflow = true) + + nx = 128 + Lx = 2π + κ = 0.01 + nsteps = round(Int, tfinal/dt) + + if !isapprox(tfinal, nsteps*dt, rtol=rtol_traceradvectiondiffusion) + error("tfinal is not multiple of dt") + end + + #advecting_flow = steadyflow==true ? u(x) = 0.0 : ut(x, t) = 0.0 + advecting_flow = OneDAdvectingFlow(; steadyflow) + + prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; nx, Lx, κ, dt, stepper) + sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid + x = gridpoints(gr) + + c0ampl, σ₀ = 0.1, 0.1 + σ(t) = sqrt(2κ * t + σ₀) + c0func(x, t) = (c0ampl / σ(t)) * exp(-x^2 / 2σ(t)^2) + + c0 = @. c0func(x, 0) + tfinal = nsteps * dt + cfinal = @. c0func(x, tfinal) + + TracerAdvectionDiffusion.set_c!(prob, c0) + + stepforward!(prob, nsteps) + TracerAdvectionDiffusion.updatevars!(prob) + + return isapprox(cfinal, vs.c, rtol=gr.nx*nsteps*1e-12) +end """ test_diffusion(; kwargs...) @@ -90,18 +202,18 @@ function test_diffusion(stepper, dt, tfinal, dev::Device=CPU(); steadyflow = tru error("tfinal is not multiple of dt") end - prob = TracerAdvectionDiffusion.Problem(dev; steadyflow=steadyflow, nx=nx, - Lx=Lx, κ=κ, dt=dt, stepper=stepper) + advecting_flow = TwoDAdvectingFlow(; steadyflow) + prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; nx, Lx, κ, dt, stepper) sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid x, y = gridpoints(gr) c0ampl, σ = 0.1, 0.1 - c0func(x, y) = c0ampl * exp(-(x^2 + y^2) / (2σ^2)) + c0func(x, y) = c0ampl * exp(-(x^2 + y^2) / 2σ^2) c0 = @. c0func(x, y) tfinal = nsteps * dt σt = sqrt(2κ * tfinal + σ^2) - cfinal = @. c0ampl * (σ^2 / σt^2) * exp(-(x^2 + y^2) / (2σt^2)) + cfinal = @. c0ampl * (σ^2 / σt^2) * exp(-(x^2 + y^2) / 2σt^2) TracerAdvectionDiffusion.set_c!(prob, c0) @@ -113,7 +225,6 @@ end function test_diffusion_multilayerqg(stepper, dt, tfinal, dev::Device=CPU()) # Set up MQGprob to generate zero flow and diffuse concentration field - nx = 128 Lx = 2π @@ -150,11 +261,11 @@ function test_diffusion_multilayerqg(stepper, dt, tfinal, dev::Device=CPU()) x, y = gridpoints(gr) c0ampl, σ = 0.1, 0.1 - c0func(x, y) = c0ampl * exp(-(x^2 + y^2) / (2σ^2)) + c0func(x, y) = c0ampl * exp(-(x^2 + y^2) / 2σ^2) c0 = @. c0func(x, y) tfinal = nsteps * dt σt = sqrt(2κ * tfinal + σ^2) - cfinal = @. c0ampl * (σ^2 / σt^2) * exp(-(x^2 + y^2) / (2σt^2)) + cfinal = @. c0ampl * (σ^2 / σt^2) * exp(-(x^2 + y^2) / 2σt^2) TracerAdvectionDiffusion.set_c!(ADprob, c0) @@ -191,20 +302,22 @@ function test_hyperdiffusion(stepper, dt, tfinal, dev::Device=CPU(); steadyflow gr = TwoDGrid(dev, nx, Lx) x, y = gridpoints(gr) - u, v = zero(x), zero(x) #0*x, 0*x + #u, v = zero(x), zero(x) #0*x, 0*x + u(x, y) = 0.0 + v(x, y) = 0.0 vs = TracerAdvectionDiffusion.Vars(dev, gr) - pr = TracerAdvectionDiffusion.ConstDiffSteadyFlowParams(η, κ, κh, nκh, u, v) + pr = TracerAdvectionDiffusion.ConstDiffSteadyFlowParams(η, κ, κh, nκh, u, v, gr) eq = TracerAdvectionDiffusion.Equation(dev, pr, gr) prob = FourierFlows.Problem(eq, stepper, dt, gr, vs, pr, dev) c0ampl, σ = 0.1, 0.1 - c0func(x, y) = c0ampl * exp(-(x^2 + y^2) / (2σ^2)) + c0func(x, y) = c0ampl * exp(-(x^2 + y^2) / 2σ^2) c0 = @. c0func(x, y) tfinal = nsteps * dt σt = sqrt(2κh * tfinal + σ^2) - cfinal = @. c0ampl * σ^2 / σt^2 * exp(-(x^2 + y^2) / (2σt^2)) + cfinal = @. c0ampl * σ^2 / σt^2 * exp(-(x^2 + y^2) / 2σt^2) TracerAdvectionDiffusion.set_c!(prob, c0)