burgers 2D works

examples/fourier2d/burgers2d.jl

@@ -18,24 +18,25 @@ end
 using OrdinaryDiffEq, Plots
 using ComponentArrays, CUDA
 
-T = Float32
+T = Float64
 nx = ny = 256
 N = nx * ny
 ν = 1e-3 |> T
 p = nothing
 
-Nmodes = 16
-tol = 1e-6 |> T
+abstol = reltol = 1e-4 |> T
 
 odealg = Tsit5()
 odealg = SSPRK43()
 
 """ spatial discr """
-V = FourierSpace(nx, ny) |> T
+interval = IntervalDomain(0.0, 1.0; periodic = true)
+domain = interval × interval
+V = FourierSpace(nx, ny; domain) |> T
 discr = Collocation()
 x, y = points(V)
 
-odecb = begin
+callback = begin
     function affect!(int)
         println(
                 "[$(int.iter)] \t Time $(round(int.t; digits=8))s"
@@ -46,14 +47,17 @@ odecb = begin
 end
 
 """ IC """
-u0 = begin
-    X = truncationOp(V, (Nmodes//nx, Nmodes//ny))
-    vx0 = X * rand(T, size(x)...)
-    vy0 = X * rand(T, size(x)...)
 
-    ComponentArray(;vx=vx0, vy=vy0)
+function uIC(x, y; μ = 0.9) # 0.9 - 1.1
+    u = @. μ * sin(2π * x) * sin(2π * y)
+    u[x .> 0.5] .= 0
+    u[y .> 0.5] .= 0
+
+    ComponentArray(;vx = u, vy = copy(u))
 end
 
+u0 = uIC(x,y) .|> T
+
 ps = ComponentArray(vel=u0)
 V = make_transform(V, u0.vx; p=ps)
 
@@ -69,43 +73,55 @@ Ax = -diffusionOp(ν, V, discr)
 Ay = -diffusionOp(ν, V, discr)
 
 Cx = advectionOp((zero(x), zero(x)), V, discr;
-                 vel_update_funcs=(
-                                   (v,u,p,t) -> copy!(v, p.vel.vx),
-                                   (v,u,p,t) -> copy!(v, p.vel.vy),
-                                  )
-                )
+    vel_update_funcs = (
+        (v,u,p,t) -> copy!(v, p.vel.vx),
+        (v,u,p,t) -> copy!(v, p.vel.vy),
+
+        # (v,u,p,t) -> fill!(v, true),
+        # (v,u,p,t) -> fill!(v, true),
+    ),
+)
 
 Cy = advectionOp((zero(x), zero(x)), V, discr;
-                 vel_update_funcs=(
-                                   (v,u,p,t) -> copy!(v, p.vel.vx),
-                                   (v,u,p,t) -> copy!(v, p.vel.vy),
-                                  )
-                )
+    vel_update_funcs = (
+        (v,u,p,t) -> copy!(v, p.vel.vx),
+        (v,u,p,t) -> copy!(v, p.vel.vy),
 
-Fx = NullOperator(V)
-Fy = NullOperator(V)
+        # (v,u,p,t) -> fill!(v, true),
+        # (v,u,p,t) -> fill!(v, true),
+    ),
+)
 
-Dtx = cache_operator(Ax-Cx+Fx, x)
-Dty = cache_operator(Ay-Cy+Fy, x)
+Dtx = cache_operator(Ax - Cx, x)
+Dty = cache_operator(Ay - Cy, x)
 
-function ddt(du, u, p, t)
+# function ddt(du, u, p, t)
+#     ps = ComponentArray(vel=u)
+#
+#     Dtx(du.vx, u.vx, ps, t)
+#     Dty(du.vy, u.vy, ps, t)
+#
+#     du
+# end
+
+function ddt(u, p, t)
     ps = ComponentArray(vel=u)
 
-    Dtx(du.vx, u.vx, ps, t)
-    Dty(du.vy, u.vy, ps, t)
+    dvx = Dtx(u.vx, ps, t)
+    dvy = Dty(u.vy, ps, t)
 
-    du
+    ComponentArray(; vx = dvx, vy = dvy)
 end
 
 """ time discr """
-tspan = (0.0, 10.0) .|> T
-tsave = range(tspan...; length=100)
+tspan = (0.0, 1.0) .|> T
+tsave = range(tspan...; length=100) .|> T
 prob  = ODEProblem(ddt, u0, tspan, p)
 
-#@time sol = solve(prob, odealg, saveat=tsave, abstol=1f-2, reltol=1f-2, callback=odecb)
-@time sol = solve(prob, odealg, saveat=tsave, abstol=tol, reltol=tol, callback=odecb)
-#@time sol = solve(prob, odealg, saveat=tsave, abstol=tol, reltol=tol)
-#@time sol = solve(prob, odealg, saveat=tsave, abstol=tol, reltol=tol)
+# plt = reshape(uIC(x, y).vx, nx, ny) |> heatmap
+# display(plt)
+
+@time sol = solve(prob, odealg; saveat = tsave, abstol, reltol, callback)
 
 pred = Array(sol)
 vx = @views pred[:vx, :]