|
| 1 | +using Test |
| 2 | +using LinearAlgebra |
| 3 | +import ClimaComms |
| 4 | +ClimaComms.@import_required_backends |
| 5 | +using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 |
| 6 | +using ClimaTimeSteppers |
| 7 | + |
| 8 | +import ClimaCore: |
| 9 | + Fields, |
| 10 | + Domains, |
| 11 | + Topologies, |
| 12 | + Meshes, |
| 13 | + DataLayouts, |
| 14 | + Operators, |
| 15 | + Geometry, |
| 16 | + Spaces |
| 17 | + |
| 18 | + |
| 19 | +# Advection Equation, with constant advective velocity (so advection form = flux form) |
| 20 | +# ∂_t y + w ∂_z y = 0 |
| 21 | +# the solution translates to the right at speed w, |
| 22 | +# so at time t, the solution is y(z - w * t) |
| 23 | + |
| 24 | +# visualization artifacts |
| 25 | +ENV["GKSwstype"] = "nul" |
| 26 | +using ClimaCorePlots, Plots |
| 27 | +Plots.GRBackend() |
| 28 | +dir = "tvd_advection" |
| 29 | +path = joinpath(@__DIR__, "output", dir) |
| 30 | +mkpath(path) |
| 31 | + |
| 32 | + |
| 33 | +function tendency!(yₜ, y, parameters, t) |
| 34 | + (; w, Δt, limiter_method) = parameters |
| 35 | + FT = Spaces.undertype(axes(y.q)) |
| 36 | + bcvel = pulse(-π, t, z₀, zₕ, z₁) |
| 37 | + divf2c = Operators.DivergenceF2C( |
| 38 | + bottom = Operators.SetValue(Geometry.WVector(FT(0))), |
| 39 | + top = Operators.SetValue(Geometry.WVector(FT(0))), |
| 40 | + ) |
| 41 | + upwind1 = Operators.UpwindBiasedProductC2F( |
| 42 | + bottom = Operators.Extrapolate(), |
| 43 | + top = Operators.Extrapolate(), |
| 44 | + ) |
| 45 | + upwind3 = Operators.Upwind3rdOrderBiasedProductC2F( |
| 46 | + bottom = Operators.ThirdOrderOneSided(), |
| 47 | + top = Operators.ThirdOrderOneSided(), |
| 48 | + ) |
| 49 | + FCTZalesak = Operators.FCTZalesak( |
| 50 | + bottom = Operators.FirstOrderOneSided(), |
| 51 | + top = Operators.FirstOrderOneSided(), |
| 52 | + ) |
| 53 | + TVDSlopeLimited = Operators.TVDLimitedFluxC2F( |
| 54 | + bottom = Operators.FirstOrderOneSided(), |
| 55 | + top = Operators.FirstOrderOneSided(), |
| 56 | + method = limiter_method, |
| 57 | + ) |
| 58 | + |
| 59 | + If = Operators.InterpolateC2F() |
| 60 | + |
| 61 | + if limiter_method == "Zalesak" |
| 62 | + @. yₜ.q = |
| 63 | + -divf2c( |
| 64 | + upwind1(w, y.q) + FCTZalesak( |
| 65 | + upwind3(w, y.q) - upwind1(w, y.q), |
| 66 | + y.q / Δt, |
| 67 | + y.q / Δt - divf2c(upwind1(w, y.q)), |
| 68 | + ), |
| 69 | + ) |
| 70 | + else |
| 71 | + Δfluxₕ = @. w * If(y.q) |
| 72 | + Δfluxₗ = @. upwind1(w, y.q) |
| 73 | + @. yₜ.q = |
| 74 | + -divf2c( |
| 75 | + upwind1(w, y.q) + |
| 76 | + TVDSlopeLimited(upwind3(w, y.q) - upwind1(w, y.q), y.q / Δt, w), |
| 77 | + ) |
| 78 | + end |
| 79 | +end |
| 80 | + |
| 81 | +# Define a pulse wave or square wave |
| 82 | + |
| 83 | +FT = Float64 |
| 84 | +t₀ = FT(0.0) |
| 85 | +t₁ = FT(6) |
| 86 | +z₀ = FT(0.0) |
| 87 | +zₕ = FT(2π) |
| 88 | +z₁ = FT(1.0) |
| 89 | +speed = FT(-1.0) |
| 90 | +pulse(z, t, z₀, zₕ, z₁) = abs(z - speed * t) ≤ zₕ ? z₁ : z₀ |
| 91 | + |
| 92 | +n = 2 .^ 8 |
| 93 | +elemlist = 2 .^ [3, 4, 5, 6, 7, 8, 9, 10] |
| 94 | +Δt = FT(0.3) * (20π / n) |
| 95 | +@info "Timestep Δt[s]: $(Δt)" |
| 96 | + |
| 97 | +domain = Domains.IntervalDomain( |
| 98 | + Geometry.ZPoint{FT}(-10π), |
| 99 | + Geometry.ZPoint{FT}(10π); |
| 100 | + boundary_names = (:bottom, :top), |
| 101 | +) |
| 102 | + |
| 103 | +stretch_fns = [Meshes.Uniform()] |
| 104 | +plot_string = ["uniform"] |
| 105 | + |
| 106 | +for (i, stretch_fn) in enumerate(stretch_fns) |
| 107 | + limiter_methods = ( |
| 108 | + Operators.RZeroLimiter(), |
| 109 | + Operators.RMaxLimiter(), |
| 110 | + Operators.KorenLimiter(), |
| 111 | + Operators.SuperbeeLimiter(), |
| 112 | + Operators.MonotonizedCentralLimiter(), |
| 113 | + "Zalesak", |
| 114 | + ) |
| 115 | + for (j, limiter_method) in enumerate(limiter_methods) |
| 116 | + @info (limiter_method, stretch_fn) |
| 117 | + mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n) |
| 118 | + cent_space = Spaces.CenterFiniteDifferenceSpace(mesh) |
| 119 | + face_space = Spaces.FaceFiniteDifferenceSpace(cent_space) |
| 120 | + z = Fields.coordinate_field(cent_space).z |
| 121 | + O = ones(FT, face_space) |
| 122 | + |
| 123 | + # Initial condition |
| 124 | + q_init = pulse.(z, 0.0, z₀, zₕ, z₁) |
| 125 | + y = Fields.FieldVector(q = q_init) |
| 126 | + |
| 127 | + # Unitary, constant advective velocity |
| 128 | + w = Geometry.WVector.(speed .* O) |
| 129 | + |
| 130 | + # Solve the ODE |
| 131 | + parameters = (; w, Δt, limiter_method) |
| 132 | + prob = ODEProblem( |
| 133 | + ClimaODEFunction(; T_exp! = tendency!), |
| 134 | + y, |
| 135 | + (t₀, t₁), |
| 136 | + parameters, |
| 137 | + ) |
| 138 | + sol = solve( |
| 139 | + prob, |
| 140 | + ExplicitAlgorithm(SSP33ShuOsher()), |
| 141 | + dt = Δt, |
| 142 | + saveat = Δt, |
| 143 | + ) |
| 144 | + |
| 145 | + q_final = sol.u[end].q |
| 146 | + q_analytic = pulse.(z, t₁, z₀, zₕ, z₁) |
| 147 | + |
| 148 | + err = norm(q_final .- q_analytic) |
| 149 | + rel_mass_err = norm((sum(q_final) - sum(q_init)) / sum(q_init)) |
| 150 | + |
| 151 | + if j == 1 |
| 152 | + fig = Plots.plot(q_analytic; label = "Exact", color = :red) |
| 153 | + end |
| 154 | + linstyl = [:dash, :dot, :dashdot, :dashdotdot, :dash, :solid] |
| 155 | + clrs = [:orange, :gray, :green, :maroon, :pink, :blue] |
| 156 | + if limiter_method == "Zalesak" |
| 157 | + fig = plot!( |
| 158 | + q_final; |
| 159 | + label = "Zalesak", |
| 160 | + linestyle = linstyl[j], |
| 161 | + color = clrs[j], |
| 162 | + dpi = 400, |
| 163 | + xlim = (-0.5, 1.1), |
| 164 | + ylim = (-15, 10), |
| 165 | + ) |
| 166 | + else |
| 167 | + fig = plot!( |
| 168 | + q_final; |
| 169 | + label = "$(typeof(limiter_method))"[21:end], |
| 170 | + linestyle = linstyl[j], |
| 171 | + color = clrs[j], |
| 172 | + dpi = 400, |
| 173 | + xlim = (-0.5, 1.1), |
| 174 | + ylim = (-15, 10), |
| 175 | + ) |
| 176 | + end |
| 177 | + fig = plot!(legend = :outerbottom, legendcolumns = 2) |
| 178 | + if j == length(limiter_methods) |
| 179 | + Plots.png( |
| 180 | + fig, |
| 181 | + joinpath( |
| 182 | + path, |
| 183 | + "SlopeLimitedFluxSolution_" * |
| 184 | + "$(typeof(limiter_method))"[21:end] * |
| 185 | + ".png", |
| 186 | + ), |
| 187 | + ) |
| 188 | + end |
| 189 | + @test err ≤ 0.25 |
| 190 | + @test rel_mass_err ≤ 10eps() |
| 191 | + end |
| 192 | +end |
0 commit comments