From 1a20d66632caf79d3fd21367d122ff3c9dc49690 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 24 Oct 2022 01:50:18 +1100 Subject: [PATCH 1/5] minor cleanup of timesteppers --- src/timesteppers.jl | 72 ++++++++++++++++----------------------------- 1 file changed, 25 insertions(+), 47 deletions(-) diff --git a/src/timesteppers.jl b/src/timesteppers.jl index c409bd17..56fd8651 100644 --- a/src/timesteppers.jl +++ b/src/timesteppers.jl @@ -73,17 +73,11 @@ end # * Filtered Forward Euler # * RK4 # * Filtered RK4 -# * ETDRK4 # * LSRK54 +# * ETDRK4 # * Filtered ETDRK4 # * AB3 # * Filtered AB3 -# -# Explicit time-steppers are constructed with the signature -# ts = ExplicitTimeStepper(equation::Equation) -# -# Implicit time-steppers are constructed with the signature -# ts = ImplicitTimeStepper(equation::Equation, dt) # -- # Forward Euler @@ -98,7 +92,7 @@ uⁿ⁺¹ = uⁿ + dt * RHS(uⁿ, tⁿ) ``` """ struct ForwardEulerTimeStepper{T} <: AbstractTimeStepper{T} - N::T # Explicit linear and nonlinear terms + N :: T # Explicit linear and nonlinear terms ForwardEulerTimeStepper(N::T) where T = new{T}(0N) end @@ -113,6 +107,7 @@ ForwardEulerTimeStepper(equation::Equation, dev::Device=CPU()) = function stepforward!(sol, clock, ts::ForwardEulerTimeStepper, equation, vars, params, grid) equation.calcN!(ts.N, sol, clock.t, clock, vars, params, grid) @. sol += clock.dt * (equation.L * sol + ts.N) + clock.t += clock.dt clock.step += 1 @@ -143,6 +138,7 @@ end function stepforward!(sol, clock, ts::FilteredForwardEulerTimeStepper, equation, vars, params, grid) equation.calcN!(ts.N, sol, clock.t, clock, vars, params, grid) @. sol = ts.filter * (sol + clock.dt * (ts.N + equation.L * sol)) + clock.t += clock.dt clock.step += 1 @@ -188,7 +184,7 @@ end Construct a 4th-order Runge-Kutta timestepper for `equation` on device `dev`. """ function RK4TimeStepper(equation::Equation, dev::Device=CPU()) - @devzeros typeof(dev) equation.T equation.dims N sol₁ RHS₁ RHS₂ RHS₃ RHS₄ + @devzeros typeof(dev) equation.T equation.dims sol₁ RHS₁ RHS₂ RHS₃ RHS₄ return RK4TimeStepper(sol₁, RHS₁, RHS₂, RHS₃, RHS₄) end @@ -255,13 +251,7 @@ function RK4substeps!(sol, clock, ts, equation, vars, params, grid, t, dt) end function RK4update!(sol, RHS₁, RHS₂, RHS₃, RHS₄, dt) - @. sol += dt*(RHS₁ / 6 + RHS₂ / 3 + RHS₃ / 3 + RHS₄ / 6) - - return nothing -end - -function RK4update!(sol, RHS₁, RHS₂, RHS₃, RHS₄, filter, dt) - @. sol = filter * (sol + dt*(RHS₁ / 6 + RHS₂ / 3 + RHS₃ / 3 + RHS₄ / 6)) + @. sol += dt * (RHS₁ / 6 + RHS₂ / 3 + RHS₃ / 3 + RHS₄ / 6) return nothing end @@ -269,6 +259,7 @@ end function stepforward!(sol, clock, ts::RK4TimeStepper, equation, vars, params, grid) RK4substeps!(sol, clock, ts, equation, vars, params, grid, clock.t, clock.dt) RK4update!(sol, ts.RHS₁, ts.RHS₂, ts.RHS₃, ts.RHS₄, clock.dt) + clock.t += clock.dt clock.step += 1 @@ -277,7 +268,9 @@ end function stepforward!(sol, clock, ts::FilteredRK4TimeStepper, equation, vars, params, grid) RK4substeps!(sol, clock, ts, equation, vars, params, grid, clock.t, clock.dt) - RK4update!(sol, ts.RHS₁, ts.RHS₂, ts.RHS₃, ts.RHS₄, ts.filter, clock.dt) + RK4update!(sol, ts.RHS₁, ts.RHS₂, ts.RHS₃, ts.RHS₄, clock.dt) + @. sol *= ts.filter + clock.t += clock.dt clock.step += 1 @@ -353,26 +346,23 @@ function LSRK54TimeStepper(equation::Equation, dev::Device=CPU()) return LSRK54TimeStepper(S², RHS, Tuple(A), Tuple(B), Tuple(C)) end -function LSRK54!(sol, clock, ts, equation, vars, params, grid, t, dt) - T = equation.T - A, B, C = ts.A, ts.B, ts.C - - # initialize the S² term +function LSRK54update!(sol, clock, ts, equation, vars, params, grid, t, dt) @. ts.S² = 0 for i = 1:5 - equation.calcN!(ts.RHS, sol, t + C[i] * dt , clock, vars, params, grid) + equation.calcN!(ts.RHS, sol, t + ts.C[i] * dt , clock, vars, params, grid) addlinearterm!(ts.RHS, equation.L, sol) - @. ts.S² = A[i] * ts.S² + dt * ts.RHS - @. sol += B[i] * ts.S² + @. ts.S² = ts.A[i] * ts.S² + dt * ts.RHS + @. sol += ts.B[i] * ts.S² end return nothing end function stepforward!(sol, clock, ts::LSRK54TimeStepper, equation, vars, params, grid) - LSRK54!(sol, clock, ts, equation, vars, params, grid, clock.t, clock.dt) + LSRK54update!(sol, clock, ts, equation, vars, params, grid, clock.t, clock.dt) + clock.t += clock.dt clock.step += 1 @@ -466,15 +456,7 @@ end function ETDRK4update!(sol, expLdt, α, β, Γ, N₁, N₂, N₃, N₄) @. sol = (expLdt * sol + α * N₁ + 2β * (N₂ + N₃) - + Γ * N₄ ) - - return nothing -end - -function ETDRK4update!(sol, ts, filter) - @. sol = filter * (ts.expLdt * sol + ts.α * ts.N₁ - + 2 * ts.β * (ts.N₂ + ts.N₃) - + ts.Γ * ts.N₄ ) + + Γ * N₄) return nothing end @@ -515,6 +497,7 @@ end function stepforward!(sol, clock, ts::ETDRK4TimeStepper, equation, vars, params, grid) ETDRK4substeps!(sol, clock, ts, equation, vars, params, grid) ETDRK4update!(sol, ts.expLdt, ts.α, ts.β, ts.Γ, ts.N₁, ts.N₂, ts.N₃, ts.N₄) + clock.t += clock.dt clock.step += 1 @@ -523,7 +506,9 @@ end function stepforward!(sol, clock, ts::FilteredETDRK4TimeStepper, equation, vars, params, grid) ETDRK4substeps!(sol, clock, ts, equation, vars, params, grid) - ETDRK4update!(sol, ts, ts.filter) # update + ETDRK4update!(sol, ts.expLdt, ts.α, ts.β, ts.Γ, ts.N₁, ts.N₂, ts.N₃, ts.N₄) + @. sol *= ts.filter + clock.t += clock.dt clock.step += 1 @@ -606,21 +591,12 @@ function AB3update!(sol, ts, clock) return nothing end -function AB3update!(sol, ts::FilteredAB3TimeStepper, clock) - if clock.step < 3 # forward Euler steps to initialize AB3 - @. sol = ts.filter * (sol + clock.dt * ts.RHS) # Update - else # Otherwise, stepforward with 3rd order Adams Bashforth: - @. sol = ts.filter * (sol + clock.dt * (ab3h1 * ts.RHS - ab3h2 * ts.RHS₋₁ + ab3h3 * ts.RHS₋₂)) - end - - return nothing -end - function stepforward!(sol, clock, ts::AB3TimeStepper, equation, vars, params, grid) equation.calcN!(ts.RHS, sol, clock.t, clock, vars, params, grid) addlinearterm!(ts.RHS, equation.L, sol) - + AB3update!(sol, ts, clock) + clock.t += clock.dt clock.step += 1 @@ -635,6 +611,8 @@ function stepforward!(sol, clock, ts::FilteredAB3TimeStepper, equation, vars, pa addlinearterm!(ts.RHS, equation.L, sol) AB3update!(sol, ts, clock) + @. sol *= ts.filter + clock.t += clock.dt clock.step += 1 From 44b63e6a0df00b33c20a85b34e9bf4e76853e66b Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 24 Oct 2022 01:57:52 +1100 Subject: [PATCH 2/5] code alignment --- src/problem.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/problem.jl b/src/problem.jl index a1a614d9..3fa113ae 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -111,18 +111,18 @@ function Problem(eqn::Equation, stepper, dt, grid::AbstractGrid{T}, end show(io::IO, clock::FourierFlows.Clock) = - print(io, "Clock\n", - " ├─── timestep dt: ", clock.dt, "\n", - " ├────────── step: ", clock.step, "\n", - " └──────── time t: ", clock.t) + print(io, "Clock\n", + " ├─── timestep dt: ", clock.dt, "\n", + " ├────────── step: ", clock.step, "\n", + " └──────── time t: ", clock.t) show(io::IO, eqn::FourierFlows.Equation) = - print(io, "Equation\n", - " ├──────── linear coefficients: L", "\n", - " │ ├───type: ", eltype(eqn.L), "\n", - " │ └───size: ", size(eqn.L), "\n", - " ├───────────── nonlinear term: calcN!()", "\n", - " └─── type of state vector sol: ", eqn.T) + print(io, "Equation\n", + " ├──────── linear coefficients: L", "\n", + " │ ├───type: ", eltype(eqn.L), "\n", + " │ └───size: ", size(eqn.L), "\n", + " ├───────────── nonlinear term: calcN!()", "\n", + " └─── type of state vector sol: ", eqn.T) show(io::IO, problem::FourierFlows.Problem) = print(io, "Problem\n", From de2ff9c7288c5cd4f15c92d1fae1c8adaec55324 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 24 Oct 2022 01:58:10 +1100 Subject: [PATCH 3/5] don't default to CPU on tests --- test/test_timesteppers.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_timesteppers.jl b/test/test_timesteppers.jl index e5c3c6c8..ee69dfb7 100644 --- a/test/test_timesteppers.jl +++ b/test/test_timesteppers.jl @@ -39,7 +39,7 @@ dt = 1e-9 * 1/κ # make sure diffusive decay dynamics are resolved function constantdiffusiontest_stepforward(stepper, dev::Device=CPU(); kwargs...) nsteps = 1000 - prob = construct_diffusion_problem(stepper, κ, dt; nx, Lx=2π, dev=CPU()) + prob = construct_diffusion_problem(stepper, κ, dt; nx, Lx=2π, dev) c_initial, c_final, prob, t_compute = constantdiffusion_stepforward(prob, nsteps; c₀=0.01, σ=0.2, κ=κ) normmsg = "$stepper: relative error =" @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c_final-Array(prob.vars.c))/norm(c_final), t_compute) @@ -50,7 +50,7 @@ end function varyingdiffusiontest_stepforward(stepper, dev::Device=CPU(); kwargs...) nsteps = 1000 - prob = construct_diffusion_problem(stepper, κ*ones(nx), dt; nx=nx, Lx=2π, dev=CPU()) + prob = construct_diffusion_problem(stepper, κ * ones(nx), dt; nx=nx, Lx=2π, dev) c_initial, c_final, prob, t_compute = constantdiffusion_stepforward(prob, nsteps; c₀=0.01, σ=0.2, κ=κ) normmsg = "$stepper: relative error =" @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c_final-Array(prob.vars.c))/norm(c_final), t_compute) @@ -61,7 +61,7 @@ end function constantdiffusiontest_step_until(stepper, dev::Device=CPU(); kwargs...) t_final = 1000*dt + 1e-6/π # make sure t_final is not integer multiple of dt - prob = construct_diffusion_problem(stepper, κ, dt; nx, Lx=2π, dev=CPU()) + prob = construct_diffusion_problem(stepper, κ, dt; nx, Lx=2π, dev) c_initial, c_final, prob, t_compute = constantdiffusion_step_until(prob, t_final; c₀=0.01, σ=0.2, κ=κ) normmsg = "$stepper: relative error =" @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c_final-Array(prob.vars.c))/norm(c_final), t_compute) From deab911a329236c92b63ee7bb1b9a2fb1735fa30 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 24 Oct 2022 02:01:19 +1100 Subject: [PATCH 4/5] minor cleanup --- test/test_timesteppers.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_timesteppers.jl b/test/test_timesteppers.jl index ee69dfb7..bfb23bff 100644 --- a/test/test_timesteppers.jl +++ b/test/test_timesteppers.jl @@ -40,7 +40,7 @@ function constantdiffusiontest_stepforward(stepper, dev::Device=CPU(); kwargs... nsteps = 1000 prob = construct_diffusion_problem(stepper, κ, dt; nx, Lx=2π, dev) - c_initial, c_final, prob, t_compute = constantdiffusion_stepforward(prob, nsteps; c₀=0.01, σ=0.2, κ=κ) + c_initial, c_final, prob, t_compute = constantdiffusion_stepforward(prob, nsteps; c₀=0.01, σ=0.2, κ) normmsg = "$stepper: relative error =" @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c_final-Array(prob.vars.c))/norm(c_final), t_compute) @@ -50,8 +50,8 @@ end function varyingdiffusiontest_stepforward(stepper, dev::Device=CPU(); kwargs...) nsteps = 1000 - prob = construct_diffusion_problem(stepper, κ * ones(nx), dt; nx=nx, Lx=2π, dev) - c_initial, c_final, prob, t_compute = constantdiffusion_stepforward(prob, nsteps; c₀=0.01, σ=0.2, κ=κ) + prob = construct_diffusion_problem(stepper, κ * ones(nx), dt; nx, Lx=2π, dev) + c_initial, c_final, prob, t_compute = constantdiffusion_stepforward(prob, nsteps; c₀=0.01, σ=0.2, κ) normmsg = "$stepper: relative error =" @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c_final-Array(prob.vars.c))/norm(c_final), t_compute) @@ -62,7 +62,7 @@ function constantdiffusiontest_step_until(stepper, dev::Device=CPU(); kwargs...) t_final = 1000*dt + 1e-6/π # make sure t_final is not integer multiple of dt prob = construct_diffusion_problem(stepper, κ, dt; nx, Lx=2π, dev) - c_initial, c_final, prob, t_compute = constantdiffusion_step_until(prob, t_final; c₀=0.01, σ=0.2, κ=κ) + c_initial, c_final, prob, t_compute = constantdiffusion_step_until(prob, t_final; c₀=0.01, σ=0.2, κ) normmsg = "$stepper: relative error =" @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c_final-Array(prob.vars.c))/norm(c_final), t_compute) From 12d4565e9ef601fea671f45339323b77bbaeb584 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 24 Oct 2022 02:40:56 +1100 Subject: [PATCH 5/5] export LSRK54TimeStepper --- src/FourierFlows.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/FourierFlows.jl b/src/FourierFlows.jl index a0f75737..03ef1038 100644 --- a/src/FourierFlows.jl +++ b/src/FourierFlows.jl @@ -51,6 +51,7 @@ export FilteredForwardEulerTimeStepper, RK4TimeStepper, FilteredRK4TimeStepper, + LSRK54TimeStepper, ETDRK4TimeStepper, FilteredETDRK4TimeStepper, AB3TimeStepper,