Merge pull request #98 from FourierFlows/FixesDiffusionTestModule+RK4

Fixes the RK4 time-stepper and the Diiffusion test module

Project.toml

@@ -6,8 +6,8 @@ authors = ["Gregory L. Wagner <wagner.greg@gmail.com>", "Navid C. Constantinou <
 description = "Solvers for partial differential equations on periodic domains using Fourier-based pseudospectral methods."
 documentation = "https://fourierflows.github.io/FourierFlows.jl/latest/"
 repository = "https://github.com/FourierFlows/FourierFlows.jl"
-version = "0.3.0"
-versions = ["0.0.1", "0.0.2", "0.1.0", "0.1.1", "0.1.2", "0.1.3", "0.2.0", "0.2.1", "0.2.2", "0.3.0"]
+version = "0.2.2"
+versions = ["0.0.1", "0.0.2", "0.1.0", "0.1.1", "0.1.2", "0.1.3", "0.2.0", "0.2.1", "0.2.2"]
 
 [deps]
 Coverage = "a2441757-f6aa-5fb2-8edb-039e3f45d037"

src/diffusion.jl

@@ -5,7 +5,7 @@ export
   updatevars!,
   set_c!
 
-using 
+using
   FFTW,
   Reexport
 
@@ -50,7 +50,7 @@ function Equation(kappa::T, g) where T<:Number
 end
 
 function Equation(kappa::T, g::AbstractGrid{Tg}) where {T<:AbstractArray,Tg}
-  FourierFlows.Equation(0, calcN!, g; dims=(g.nkr,), T=cxtype(Tg))
+  FourierFlows.Equation(0g.kr, calcN!, g; dims=(g.nkr,), T=cxtype(Tg))
 end
 
 # Construct Vars types
@@ -64,7 +64,7 @@ eval(varsexpression(:Vars, physicalvars, fouriervars))
 
 Returns the vars for constant diffusivity problem on grid g.
 """
-function Vars(g::AbstractGrid{T}) where T 
+function Vars(g::AbstractGrid{T}) where T
   @zeros T g.nx c cx
   @zeros Complex{T} g.nkr ch cxh
   Vars(c, cx, ch, cxh)
@@ -75,13 +75,16 @@ end
 
 Calculate the nonlinear term for the 1D heat equation.
 """
-calcN!(N, sol, t, cl, v, p::Params{T}, g) where T<:Number = nothing
+function calcN!(N, sol, t, cl, v, p::Params{T}, g) where T<:Number
+  @. N = 0
+  nothing
+end
 
 function calcN!(N, sol, t, cl, v, p::Params{T}, g) where T<:AbstractArray
   @. v.cxh = im * g.kr * sol
   ldiv!(v.cx, g.rfftplan, v.cxh)
   @. v.cx *= p.kappa
-  ldiv!(v.cxh, g.rfftplan, v.cx)
+  mul!(v.cxh, g.rfftplan, v.cx)
   @. N = im*g.kr*v.cxh
   nothing
 end

src/timesteppers.jl

@@ -12,7 +12,7 @@ end
 
 Step forward `prob` for `nsteps`.
 """
-function stepforward!(prob::Problem, nsteps::Int) 
+function stepforward!(prob::Problem, nsteps::Int)
   for step = 1:nsteps
     stepforward!(prob)
   end
@@ -46,12 +46,12 @@ isexplicit(stepper) = any(Symbol(stepper) .== fullyexplicitsteppers)
 
 """
     TimeStepper(stepper, eq, dt=nothing)
-    
+
 Generalized timestepper constructor. If `stepper` is explicit, `dt` is not used.
 """
 function TimeStepper(stepper, eq, dt=nothing)
   fullsteppername = Symbol(stepper, :TimeStepper)
-  if isexplicit(stepper) 
+  if isexplicit(stepper)
     return eval(Expr(:call, fullsteppername, eq))
   else
     return eval(Expr(:call, fullsteppername, eq, dt))
@@ -116,7 +116,7 @@ end
 Construct a forward Euler timestepper with spectral filtering.
 """
 struct FilteredForwardEulerTimeStepper{T,Tf} <: AbstractTimeStepper{T}
-  N::T       
+  N::T
   filter::Tf
 end
 
@@ -147,7 +147,7 @@ const third = 1/3
 Construct a 4th-order Runge-Kutta time stepper.
 """
 struct RK4TimeStepper{T} <: AbstractTimeStepper{T}
-  sol₁::T 
+  sol₁::T
   RHS₁::T
   RHS₂::T
   RHS₃::T
@@ -160,7 +160,7 @@ end
 Construct a 4th-order Runge-Kutta time stepper with spectral filtering for the equation `eq`.
 """
 struct FilteredRK4TimeStepper{T,Tf} <: AbstractTimeStepper{T}
-  sol₁::T 
+  sol₁::T
   RHS₁::T
   RHS₂::T
   RHS₃::T
@@ -185,7 +185,7 @@ function addlinearterm!(RHS, L, sol)
 end
 
 function substepsol!(newsol, sol, RHS, dt)
-  @. newsol = sol + 0.5dt*RHS
+  @. newsol = sol + dt*RHS
   nothing
 end
 
@@ -194,15 +194,15 @@ function RK4substeps!(sol, cl, ts, eq, v, p, g, t, dt)
   eq.calcN!(ts.RHS₁, sol, t, cl, v, p, g)
   addlinearterm!(ts.RHS₁, eq.L, sol)
   # Substep 2
-  substepsol!(ts.sol₁, sol, ts.RHS₁, cl.dt)
+  substepsol!(ts.sol₁, sol, ts.RHS₁, 0.5dt)
   eq.calcN!(ts.RHS₂, ts.sol₁, t+0.5dt, cl, v, p, g)
   addlinearterm!(ts.RHS₂, eq.L, ts.sol₁)
   # Substep 3
-  substepsol!(ts.sol₁, sol, ts.RHS₂, cl.dt)
+  substepsol!(ts.sol₁, sol, ts.RHS₂, 0.5dt)
   eq.calcN!(ts.RHS₃, ts.sol₁, t+0.5dt, cl, v, p, g)
   addlinearterm!(ts.RHS₃, eq.L, ts.sol₁)
   # Substep 4
-  substepsol!(ts.sol₁, sol, ts.RHS₃, cl.dt)
+  substepsol!(ts.sol₁, sol, ts.RHS₃, dt)
   eq.calcN!(ts.RHS₄, ts.sol₁, t+dt, cl, v, p, g)
   addlinearterm!(ts.RHS₄, eq.L, ts.sol₁)
   nothing
@@ -213,15 +213,15 @@ function RK4substeps!(sol::AbstractArray{T}, cl, ts, eq, v, p, g) where T<:Abstr
   eq.calcN!(ts.RHS₁, sol, t, cl, v, p, g)
   addlinearterm!.(ts.RHS₁, eq.L, sol)
   # Substep 2
-  substepsol!.(ts.sol₁, sol, ts.RHS₁, cl.dt)
+  substepsol!.(ts.sol₁, sol, ts.RHS₁, 0.5dt)
   eq.calcN!(ts.RHS₂, ts.sol₁, t+0.5dt, cl, v, p, g)
   addlinearterm!.(ts.RHS₂, eq.L, ts.sol₁)
   # Substep 3
-  substepsol!.(ts.sol₁, sol, ts.RHS₂, cl.dt)
+  substepsol!.(ts.sol₁, sol, ts.RHS₂, 0.5dt)
   eq.calcN!(ts.RHS₃, ts.sol₁, t+0.5dt, cl, v, p, g)
   addlinearterm!.(ts.RHS₃, eq.L, ts.sol₁)
   # Substep 4
-  substepsol!.(ts.sol₁, sol, ts.RHS₃, cl.dt)
+  substepsol!.(ts.sol₁, sol, ts.RHS₃, dt)
   eq.calcN!(ts.RHS₄, ts.sol₁, t+dt, cl, v, p, g)
   addlinearterm!.(ts.RHS₄, eq.L, ts.sol₁)
   nothing
@@ -272,15 +272,15 @@ Construct a 4th-order exponential-time-differencing Runge-Kutta time stepper. Th
 """
 struct ETDRK4TimeStepper{T,TL} <: AbstractTimeStepper{T}
   # ETDRK4 coefficents
-  ζ::TL 
+  ζ::TL
   α::TL
   β::TL
   Γ::TL
   expLdt::TL
   expLdt2::TL
   sol₁::T
   sol₂::T
-  N₁::T 
+  N₁::T
   N₂::T
   N₃::T
   N₄::T
@@ -293,15 +293,15 @@ Construct a 4th-order exponential-time-differencing Runge-Kutta time stepper wit
 """
 struct FilteredETDRK4TimeStepper{T,TL,Tf} <: AbstractTimeStepper{T}
   # ETDRK4 coefficents:
-  ζ::TL   
+  ζ::TL
   α::TL
   β::TL
   Γ::TL
   expLdt::TL
   expLdt2::TL
-  sol₁::T 
+  sol₁::T
   sol₂::T
-  N₁::T 
+  N₁::T
   N₂::T
   N₃::T
   N₄::T
@@ -502,14 +502,14 @@ function getexpLs(dt, eq)
   expLdt, expLdt2
 end
 
-function getetdcoeffs(dt, L::AbstractArray{T}; kwargs...) where T<:AbstractArray 
+function getetdcoeffs(dt, L::AbstractArray{T}; kwargs...) where T<:AbstractArray
   ζαβΓ = [ getetdcoeffs(dt, Li) for Li in L ]
   ζ = map(x->x[1], ζαβΓ)
   α = map(x->x[2], ζαβΓ)
   β = map(x->x[3], ζαβΓ)
   Γ = map(x->x[4], ζαβΓ)
   ζ, α, β, Γ
-end 
+end
 
 """
     getetdcoeffs(dt, L; ncirc=32, rcirc=1)
@@ -545,6 +545,6 @@ function getetdcoeffs(dt, L; ncirc=32, rcirc=1)
     β = real.(β)
     Γ = real.(Γ)
   end
- 
+
   ζ, α, β, Γ
 end

test/runtests.jl

@@ -13,7 +13,7 @@ using FourierFlows: parsevalsum2
 using LinearAlgebra: mul!, ldiv!, norm
 
 const rtol_fft = 1e-12
-const rtol_timesteppers = 1e-6
+const rtol_timesteppers = 1e-12
 
 const steppers = [
   "ForwardEuler",
@@ -173,7 +173,8 @@ end
   include("test_timesteppers.jl")
 
   for stepper in steppers
-    @test diffusiontest(stepper)
+    @test constantdiffusiontest(stepper)
+    @test varyingdiffusiontest(stepper)
   end
 
 end

test/test_timesteppers.jl

@@ -1,24 +1,64 @@
-function diffusionproblem(stepper; nx=6, Lx=2π, kappa=1e-2, nsteps=1000)
-  k1 = 2π/Lx
-   τ = 1/(kappa*k1^2) # time-scale for diffusive decay
+function constantdiffusionproblem(stepper; nx=128, Lx=2π, kappa=1e-2, nsteps=1000)
+   τ = 1/kappa  # time-scale for diffusive decay
   dt = 1e-9 * τ # dynamics are resolved.
 
+  kappa = kappa*ones(nx)
   prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper=stepper)
   g = prob.grid
 
-  t1 = dt*nsteps
-  c0 = @. sin(k1*g.x)
-  c1 = @. exp(-kappa*k1^2*t1) * c0 # analytical solution
+  # a gaussian initial condition c(x, t=0)
+  c0ampl, σ = 0.01, 0.2
+  c0func(x) = @. c0ampl*exp(-x^2/(2σ^2))
+  c0 = c0func.(g.x)
+
+  # analytic solution for for 1D heat equation with constant κ
+  tfinal = nsteps*dt
+  σt = sqrt(2*kappa[1]*tfinal + σ^2)
+  cfinal = @. c0ampl*σ/σt * exp(-g.x^2/(2*σt^2))
 
   set_c!(prob, c0)
   tcomp = @elapsed stepforward!(prob, nsteps)
   updatevars!(prob)
 
-  prob, c0, c1, nsteps, tcomp
+  prob, c0, cfinal, nsteps, tcomp
+end
+
+function varyingdiffusionproblem(stepper; nx=128, Lx=2π, kappa=1e-2, nsteps=1000)
+   τ = 1/kappa  # time-scale for diffusive decay
+  dt = 1e-9 * τ # dynamics are resolved.
+
+  kappa = kappa*ones(nx) # this is actually a constant diffusion but defining it as an array will force stepforward to use calcN! function instead of just the linear coefficients L*sol
+
+  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper=stepper)
+  g = prob.grid
+
+  # a gaussian initial condition c(x, t=0)
+  c0ampl, σ = 0.01, 0.2
+  c0func(x) = @. c0ampl*exp(-x^2/(2σ^2))
+  c0 = c0func.(g.x)
+
+  # analytic solution for for 1D heat equation with constant κ
+  tfinal = nsteps*dt
+  σt = sqrt(2*kappa[1]*tfinal + σ^2)
+  cfinal = @. c0ampl*σ/σt * exp(-g.x^2/(2*σt^2))
+
+  set_c!(prob, c0)
+  tcomp = @elapsed stepforward!(prob, nsteps)
+  updatevars!(prob)
+
+  prob, c0, cfinal, nsteps, tcomp
+end
+
+
+function constantdiffusiontest(stepper; kwargs...)
+  prob, c0, c1, nsteps, tcomp = constantdiffusionproblem(stepper; kwargs...)
+  normmsg = "$stepper: relative error ="
+  @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c1-prob.vars.c)/norm(c1), tcomp)
+  isapprox(c1, prob.vars.c, rtol=nsteps*rtol_timesteppers)
 end
 
-function diffusiontest(stepper; kwargs...)
-  prob, c0, c1, nsteps, tcomp = diffusionproblem(stepper; kwargs...)
+function varyingdiffusiontest(stepper; kwargs...)
+  prob, c0, c1, nsteps, tcomp = varyingdiffusionproblem(stepper; kwargs...)
   normmsg = "$stepper: relative error ="
   @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c1-prob.vars.c)/norm(c1), tcomp)
   isapprox(c1, prob.vars.c, rtol=nsteps*rtol_timesteppers)