From 8b62e67f38afaa1a75004dd5e565e4ab6264109f Mon Sep 17 00:00:00 2001 From: marinlauber Date: Wed, 26 Jun 2024 10:22:20 +0200 Subject: [PATCH 01/15] initial commit, add @debug and change divisible --- Project.toml | 3 +- src/Flow.jl | 2 ++ src/MultiLevelPoisson.jl | 17 +++++---- src/util.jl | 16 +++++++++ test/test_psolver.jl | 78 ++++++++++++++++++++++++++++++++++++++++ 5 files changed, 108 insertions(+), 8 deletions(-) create mode 100644 test/test_psolver.jl diff --git a/Project.toml b/Project.toml index 10fae33..6ca5cf3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "WaterLily" uuid = "ed894a53-35f9-47f1-b17f-85db9237eebd" authors = ["Gabriel Weymouth "] -version = "1.1" +version = "1.1.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -9,6 +9,7 @@ EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" diff --git a/src/Flow.jl b/src/Flow.jl index 9dbd39e..c285a60 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -153,6 +153,7 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp @fastmath function mom_step!(a::Flow{N},b::AbstractPoisson) where N a.u⁰ .= a.u; scale_u!(a,0) # predictor u → u' + @debug "mlp" U = BCTuple(a.U,@view(a.Δt[1:end-1]),N) conv_diff!(a.f,a.u⁰,a.σ,ν=a.ν,perdir=a.perdir) accelerate!(a.f,@view(a.Δt[1:end-1]),a.g,a.U) @@ -160,6 +161,7 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp a.exitBC && exitBC!(a.u,a.u⁰,U,a.Δt[end]) # convective exit project!(a,b); BC!(a.u,U,a.exitBC,a.perdir) # corrector u → u¹ + @debug "mlc" U = BCTuple(a.U,a.Δt,N) conv_diff!(a.f,a.u,a.σ,ν=a.ν,perdir=a.perdir) accelerate!(a.f,a.Δt,a.g,a.U) diff --git a/src/MultiLevelPoisson.jl b/src/MultiLevelPoisson.jl index cb262df..3c340f9 100644 --- a/src/MultiLevelPoisson.jl +++ b/src/MultiLevelPoisson.jl @@ -33,7 +33,7 @@ end restrict!(a,b) = @inside a[I] = restrict(I,b) prolongate!(a,b) = @inside a[I] = b[down(I)] -@inline divisible(N) = mod(N,2)==0 && N>4 +@inline divisible(N) = mod(N,2)==0 && N>8 @inline divisible(l::Poisson) = all(size(l.x) .|> divisible) """ MultiLevelPoisson{N,M} @@ -86,15 +86,18 @@ residual!(ml::MultiLevelPoisson,x) = residual!(ml.levels[1],x) function solver!(ml::MultiLevelPoisson;tol=2e-4,itmx=32) p = ml.levels[1] - residual!(p); r₀ = r₂ = L∞(p); r₂₀ = L₂(p) + residual!(p); r∞⁰ = r∞ = L∞(p); r₂⁰ = r₂ = L₂(p) nᵖ=0 - while r₂>tol && nᵖ5) && pop!(ml.levels); # remove coarsest level if this was easy - (nᵖ>4 && divisible(ml.levels[end])) && push!(ml.levels,restrictML(ml.levels[end])) # add a level if this was hard + # (nᵖ<2 && length(ml.levels)>5) && pop!(ml.levels); # remove coarsest level if this was easy + # (nᵖ>4 && divisible(ml.levels[end])) && push!(ml.levels,restrictML(ml.levels[end])) # add a level if this was hard push!(ml.n,nᵖ); -end +end \ No newline at end of file diff --git a/src/util.jl b/src/util.jl index 4bcc297..7502647 100644 --- a/src/util.jl +++ b/src/util.jl @@ -1,4 +1,20 @@ using KernelAbstractions: get_backend, @index, @kernel +using LoggingExtras + +""" + logger(fname="WaterLily") + +Set up a logger to write the pressure solver data to a logging file named `WaterLily.log`. +""" +function logger(fname::String="WaterLily") + ENV["JULIA_DEBUG"] = all + logger = FormatLogger(fname*".log"; append=false) do io, args + (args.level <= Logging.Debug && args.message[1:2]=="ml" ) && print(io, args.message[3:end]) + end; + global_logger(logger); + # put header file + @debug "mlp/c, iter, r∞⁰,r∞, r₂⁰, r₂\n" +end @inline CI(a...) = CartesianIndex(a...) """ diff --git a/test/test_psolver.jl b/test/test_psolver.jl new file mode 100644 index 0000000..3bc6d57 --- /dev/null +++ b/test/test_psolver.jl @@ -0,0 +1,78 @@ +using WaterLily +using StaticArrays +using BiotSavartBCs +using LoggingExtras + +"""Circle function""" +function circle(L=32;m=6,n=4,Re=80,U=1,T=Float32) + radius, center = L/2, max(n*L/2,L) + body = AutoBody((x,t)->√sum(abs2, x .- center) - radius) + Simulation((m*L,n*L), (U,0), radius; ν=U*radius/Re, body, T) +end + +include("../examples/TwoD_plots.jl") + +# # # hydrostatic pressure force +# f1=[]; f2=[]; f3=[]; resolutions = [16,32,64,128,256,512] +# for N ∈ resolutions +# a = Flow((N,N),(1.,0.);f=Array,T=Float32) +# sdf(x,t) = √sum(abs2,x.-N÷2)-N÷4 +# map(x,t) = x.-SVector(t,0) +# body = AutoBody(sdf,map) +# WaterLily.measure!(a,body) +# @inside a.p[I] = loc(0,I)[2] +# # @inside a.p[I] = sdf(loc(0,I),0) >= 0 ? loc(0,I)[2] : 0 +# push!(f1,WaterLily.pressure_force(a,body)/(π*(N÷4)^2)) +# end +# plot(title="Hydrostatic pressure force",xlabel="N",ylabel="force/πR²") +# plot!(resolutions,reduce(hcat,f1)[2,:],label="WaterLily.pressure_force(sim)") +# savefig("hydrostatic_force.png") + +# make the sim +body = AutoBody((x,t)->√sum(abs2,x.-N÷2)-N÷4,(x,t)->x.-SVector(t,0)) +sim = circle(64;m=24,n=16,Re=80,U=1,T=Float32) +# ml_ω = MLArray(sim.flow.σ) + +# allows logging the pressure solver results +WaterLily.logger("test_psolver") + +# intialize +t₀ = sim_time(sim); duration = 10; tstep = 0.1 +forces_p = []; forces_ν = [] + +# step and plot +anim = @animate for tᵢ in range(t₀,t₀+duration;step=tstep) + # update until time tᵢ in the background + t = sum(sim.flow.Δt[1:end-1]) + while t < tᵢ*sim.L/sim.U + + # update flow + mom_step!(sim.flow,sim.pois) + # biot_mom_step!(sim.flow,sim.pois,ml_ω) + + # pressure force + force = -2WaterLily.pressure_force(sim) + push!(forces_p,force) + vforce = -2WaterLily.viscous_force(sim) + push!(forces_ν,vforce) + # update time + t += sim.flow.Δt[end] + end + + # print time step + println("tU/L=",round(tᵢ,digits=4),", Δt=",round(sim.flow.Δt[end],digits=3)) + a = sim.flow.σ; + @inside a[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U + flood(a[inside(a)],clims=(-10,10), legend=false); body_plot!(sim) + contour!(sim.flow.p[inside(a)]',levels=range(-1,1,length=10), + color=:black,linewidth=0.5,legend=false) + # flood(sim.flow.p[inside(a)],clims=(-1,1)); body_plot!(sim) + # plot!([100],[200],marker=:o,color=:red,markersize=2,legend=false) +end +gif(anim, "cylinder_Float32_FullV.gif", fps = 10) +time = cumsum(sim.flow.Δt[4:end-1]) +forces_p = reduce(vcat,forces_p') +forces_ν = reduce(vcat,forces_ν') +plot(time/sim.L,forces_p[4:end,1]/(sim.L),label="pressure force") +plot!(time/sim.L,forces_ν[4:end,1]/(sim.L),label="viscous force") +xlabel!("tU/L"); ylabel!("force/L"); savefig("cylinder_force_Float32_FullV.png") \ No newline at end of file From 9a9b893c6d6b4510545d2e70868d42518951e6b7 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 27 Jun 2024 09:21:50 +0200 Subject: [PATCH 02/15] remove changes to MultiLevelPoisson.jl --- src/MultiLevelPoisson.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/MultiLevelPoisson.jl b/src/MultiLevelPoisson.jl index 3c340f9..85a251e 100644 --- a/src/MultiLevelPoisson.jl +++ b/src/MultiLevelPoisson.jl @@ -33,7 +33,7 @@ end restrict!(a,b) = @inside a[I] = restrict(I,b) prolongate!(a,b) = @inside a[I] = b[down(I)] -@inline divisible(N) = mod(N,2)==0 && N>8 +@inline divisible(N) = mod(N,2)==0 && N>4 @inline divisible(l::Poisson) = all(size(l.x) .|> divisible) """ MultiLevelPoisson{N,M} @@ -88,16 +88,14 @@ function solver!(ml::MultiLevelPoisson;tol=2e-4,itmx=32) p = ml.levels[1] residual!(p); r∞⁰ = r∞ = L∞(p); r₂⁰ = r₂ = L₂(p) nᵖ=0 - while nᵖtol && nᵖ5) && pop!(ml.levels); # remove coarsest level if this was easy - # (nᵖ>4 && divisible(ml.levels[end])) && push!(ml.levels,restrictML(ml.levels[end])) # add a level if this was hard + (nᵖ<2 && length(ml.levels)>5) && pop!(ml.levels); # remove coarsest level if this was easy + (nᵖ>4 && divisible(ml.levels[end])) && push!(ml.levels,restrictML(ml.levels[end])) # add a level if this was hard push!(ml.n,nᵖ); end \ No newline at end of file From 6ec0fa269919d7255adaa6f587d91f7c145fd1f7 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 27 Jun 2024 09:27:21 +0200 Subject: [PATCH 03/15] move pressure monitor example --- .../TwoD_cylinder_pressure_monitor.jl | 41 ++++--------------- 1 file changed, 9 insertions(+), 32 deletions(-) rename test/test_psolver.jl => examples/TwoD_cylinder_pressure_monitor.jl (50%) diff --git a/test/test_psolver.jl b/examples/TwoD_cylinder_pressure_monitor.jl similarity index 50% rename from test/test_psolver.jl rename to examples/TwoD_cylinder_pressure_monitor.jl index 3bc6d57..023d474 100644 --- a/test/test_psolver.jl +++ b/examples/TwoD_cylinder_pressure_monitor.jl @@ -1,7 +1,7 @@ using WaterLily using StaticArrays -using BiotSavartBCs using LoggingExtras +include("TwoD_plots.jl") """Circle function""" function circle(L=32;m=6,n=4,Re=80,U=1,T=Float32) @@ -10,28 +10,9 @@ function circle(L=32;m=6,n=4,Re=80,U=1,T=Float32) Simulation((m*L,n*L), (U,0), radius; ν=U*radius/Re, body, T) end -include("../examples/TwoD_plots.jl") - -# # # hydrostatic pressure force -# f1=[]; f2=[]; f3=[]; resolutions = [16,32,64,128,256,512] -# for N ∈ resolutions -# a = Flow((N,N),(1.,0.);f=Array,T=Float32) -# sdf(x,t) = √sum(abs2,x.-N÷2)-N÷4 -# map(x,t) = x.-SVector(t,0) -# body = AutoBody(sdf,map) -# WaterLily.measure!(a,body) -# @inside a.p[I] = loc(0,I)[2] -# # @inside a.p[I] = sdf(loc(0,I),0) >= 0 ? loc(0,I)[2] : 0 -# push!(f1,WaterLily.pressure_force(a,body)/(π*(N÷4)^2)) -# end -# plot(title="Hydrostatic pressure force",xlabel="N",ylabel="force/πR²") -# plot!(resolutions,reduce(hcat,f1)[2,:],label="WaterLily.pressure_force(sim)") -# savefig("hydrostatic_force.png") - # make the sim body = AutoBody((x,t)->√sum(abs2,x.-N÷2)-N÷4,(x,t)->x.-SVector(t,0)) sim = circle(64;m=24,n=16,Re=80,U=1,T=Float32) -# ml_ω = MLArray(sim.flow.σ) # allows logging the pressure solver results WaterLily.logger("test_psolver") @@ -48,12 +29,11 @@ anim = @animate for tᵢ in range(t₀,t₀+duration;step=tstep) # update flow mom_step!(sim.flow,sim.pois) - # biot_mom_step!(sim.flow,sim.pois,ml_ω) # pressure force - force = -2WaterLily.pressure_force(sim) + force = -2WaterLily.pressure_force(sim)[1] push!(forces_p,force) - vforce = -2WaterLily.viscous_force(sim) + vforce = -2WaterLily.viscous_force(sim)[1] push!(forces_ν,vforce) # update time t += sim.flow.Δt[end] @@ -66,13 +46,10 @@ anim = @animate for tᵢ in range(t₀,t₀+duration;step=tstep) flood(a[inside(a)],clims=(-10,10), legend=false); body_plot!(sim) contour!(sim.flow.p[inside(a)]',levels=range(-1,1,length=10), color=:black,linewidth=0.5,legend=false) - # flood(sim.flow.p[inside(a)],clims=(-1,1)); body_plot!(sim) - # plot!([100],[200],marker=:o,color=:red,markersize=2,legend=false) end -gif(anim, "cylinder_Float32_FullV.gif", fps = 10) -time = cumsum(sim.flow.Δt[4:end-1]) -forces_p = reduce(vcat,forces_p') -forces_ν = reduce(vcat,forces_ν') -plot(time/sim.L,forces_p[4:end,1]/(sim.L),label="pressure force") -plot!(time/sim.L,forces_ν[4:end,1]/(sim.L),label="viscous force") -xlabel!("tU/L"); ylabel!("force/L"); savefig("cylinder_force_Float32_FullV.png") \ No newline at end of file +gif(anim,"cylinder.gif") + +# time = cumsum(sim.flow.Δt[4:end-1]) +# plot(time/sim.L,forces_p[4:end]/(sim.L),label="pressure force") +# plot!(time/sim.L,forces_ν[4:end]/(sim.L),label="viscous force") +# xlabel!("tU/L"); ylabel!("force/L") \ No newline at end of file From 85c78b86d207a3cda2041cfdd23d06975ba5b3cc Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 28 Jun 2024 14:52:37 +0200 Subject: [PATCH 04/15] simplify the logging --- src/MultiLevelPoisson.jl | 10 +++++----- src/util.jl | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/MultiLevelPoisson.jl b/src/MultiLevelPoisson.jl index 85a251e..d73040a 100644 --- a/src/MultiLevelPoisson.jl +++ b/src/MultiLevelPoisson.jl @@ -86,13 +86,13 @@ residual!(ml::MultiLevelPoisson,x) = residual!(ml.levels[1],x) function solver!(ml::MultiLevelPoisson;tol=2e-4,itmx=32) p = ml.levels[1] - residual!(p); r∞⁰ = r∞ = L∞(p); r₂⁰ = r₂ = L₂(p) - nᵖ=0 - while r₂>tol && nᵖtol && nᵖ5) && pop!(ml.levels); # remove coarsest level if this was easy diff --git a/src/util.jl b/src/util.jl index 7502647..3ccf41c 100644 --- a/src/util.jl +++ b/src/util.jl @@ -13,7 +13,7 @@ function logger(fname::String="WaterLily") end; global_logger(logger); # put header file - @debug "mlp/c, iter, r∞⁰,r∞, r₂⁰, r₂\n" + @debug "mlp/c, iter, r∞, r₂\n" end @inline CI(a...) = CartesianIndex(a...) From bf3f2f07ccf9b45624af2614780157b168611662 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Mon, 1 Jul 2024 09:00:41 +0200 Subject: [PATCH 05/15] add Plots as an Ext and add rename pressure logger example --- Project.toml | 2 + examples/TwoD_circle.jl | 3 +- ...tor.jl => TwoD_circle_pressure_monitor.jl} | 3 + ext/WaterLilyPlotsExt.jl | 102 ++++++++++++++++++ src/WaterLily.jl | 10 ++ src/util.jl | 2 +- 6 files changed, 119 insertions(+), 3 deletions(-) rename examples/{TwoD_cylinder_pressure_monitor.jl => TwoD_circle_pressure_monitor.jl} (96%) create mode 100644 ext/WaterLilyPlotsExt.jl diff --git a/Project.toml b/Project.toml index 6ca5cf3..231e70f 100644 --- a/Project.toml +++ b/Project.toml @@ -21,12 +21,14 @@ AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" [extensions] WaterLilyAMDGPUExt = "AMDGPU" WaterLilyCUDAExt = "CUDA" WaterLilyReadVTKExt = "ReadVTK" WaterLilyWriteVTKExt = "WriteVTK" +WaterLilyPlotsExt = "Plots" [compat] AMDGPU = "^0.4.13,0.6" diff --git a/examples/TwoD_circle.jl b/examples/TwoD_circle.jl index 02e972e..67da10c 100644 --- a/examples/TwoD_circle.jl +++ b/examples/TwoD_circle.jl @@ -1,9 +1,8 @@ -using WaterLily +using WaterLily,Plots function circle(n,m;Re=250,U=1) radius, center = m/8, m/2 body = AutoBody((x,t)->√sum(abs2, x .- center) - radius) Simulation((n,m), (U,0), radius; ν=U*radius/Re, body) end -include("TwoD_plots.jl") sim_gif!(circle(3*2^6,2^7),duration=10,clims=(-5,5),plotbody=true) diff --git a/examples/TwoD_cylinder_pressure_monitor.jl b/examples/TwoD_circle_pressure_monitor.jl similarity index 96% rename from examples/TwoD_cylinder_pressure_monitor.jl rename to examples/TwoD_circle_pressure_monitor.jl index 023d474..326f2c4 100644 --- a/examples/TwoD_cylinder_pressure_monitor.jl +++ b/examples/TwoD_circle_pressure_monitor.jl @@ -49,6 +49,9 @@ anim = @animate for tᵢ in range(t₀,t₀+duration;step=tstep) end gif(anim,"cylinder.gif") +# show the pressure logger +# plot_logger("test_psolver") + # time = cumsum(sim.flow.Δt[4:end-1]) # plot(time/sim.L,forces_p[4:end]/(sim.L),label="pressure force") # plot!(time/sim.L,forces_ν[4:end]/(sim.L),label="viscous force") diff --git a/ext/WaterLilyPlotsExt.jl b/ext/WaterLilyPlotsExt.jl new file mode 100644 index 0000000..bacf9ae --- /dev/null +++ b/ext/WaterLilyPlotsExt.jl @@ -0,0 +1,102 @@ +module WaterLilyPlotsExt + +if isdefined(Base, :get_extension) + using Plots; gr() +else + using ..Plots; gr() +end + +using WaterLily +import WaterLily: flood,addbody,body_plot!,sim_gif!,plot_logger + +""" + flood(f) + +Plot a filled contour plot of the 2D array `f`. The keyword arguments are passed to `Plots.contourf`. +""" +function flood(f::Array;shift=(0.,0.),cfill=:RdBu_11,clims=(),levels=10,kv...) + if length(clims)==2 + @assert clims[1]Array, + linewidth=0, levels=levels, color=cfill, clims = clims, + aspect_ratio=:equal; kv...) +end + +addbody(x,y;c=:black) = Plots.plot!(Shape(x,y), c=c, legend=false) +function body_plot!(sim;levels=[0],lines=:black,R=inside(sim.flow.p)) + WaterLily.measure_sdf!(sim.flow.σ,sim.body,WaterLily.time(sim)) + contour!(sim.flow.σ[R]'|>Array;levels,lines) +end + +""" + sim_gif!(sim;duration=1,step=0.1,verbose=true,R=inside(sim.flow.p), + remeasure=false,plotbody=false,kv...) + +Make a gif of the simulation `sim` for `duration` seconds with `step` time steps. The keyword arguments are passed to `flood` and `body_plot!`. +""" +function sim_gif!(sim;duration=1,step=0.1,verbose=true,R=inside(sim.flow.p), + remeasure=false,plotbody=false,kv...) + t₀ = round(WaterLily.sim_time(sim)) + @time @gif for tᵢ in range(t₀,t₀+duration;step) + WaterLily.sim_step!(sim,tᵢ;remeasure) + @WaterLily.inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U + flood(sim.flow.σ[R]; kv...) + plotbody && body_plot!(sim) + verbose && println("tU/L=",round(tᵢ,digits=4), + ", Δt=",round(sim.flow.Δt[end],digits=3)) + end +end + + +""" + plot_logger(fname="WaterLily.log") + +Plot the residuals and MG iterations from the log file `fname`. +""" +function plot_logger(fname="WaterLily.log") + predictor = []; corrector = [] + open(ifelse(fname[end-3:end]==".log",fname[1:end-4],fname)*".log","r") do f + readline(f) # read first line and dump it + which = "p" + while ! eof(f) + s = split(readline(f) , ",") + which = s[1] != "" ? s[1] : which + push!(which == "p" ? predictor : corrector,parse.(Float64,s[2:end])) + end + end + predictor = reduce(hcat,predictor) + corrector = reduce(hcat,corrector) + + # get index of all time steps + idx = findall(==(0.0),@views(predictor[1,:])) + # plot inital L∞ and L₂ norms of residuals for the predictor step + p1=plot(1:length(idx),predictor[2,idx],color=:1,ls=:dash,label="predictor initial r∞",yaxis=:log,size=(800,400),dpi=600, + xlabel="Time step",ylabel="L∞-norm",title="Residuals",ylims=(1e-6,1e0),xlims=(0,length(idx))) + p2=plot(1:length(idx),predictor[2,idx],color=:1,ls=:dash,label="predictor initial r₂",yaxis=:log,size=(800,400),dpi=600, + xlabel="Time step",ylabel="L₂-norm",title="Residuals",ylims=(1e-6,1e0),xlims=(0,length(idx))) + # plot final L∞ and L₂norms of residuals for the predictor + plot!(p1,1:length(idx),vcat(predictor[2,idx[2:end].-1],predictor[2,end]),color=:1,lw=2,label="predictor r∞") + plot!(p2,1:length(idx),vcat(predictor[3,idx[2:end].-1],predictor[3,end]),color=:1,lw=2,label="predictor r₂") + # plot the MG iterations for the predictor + p3=plot(1:length(idx),clamp.(vcat(predictor[1,idx[2:end].-1],predictor[1,end]),√1/2,32),lw=2,label="predictor",size=(800,400),dpi=600, + xlabel="Time step",ylabel="Iterations",title="MG Iterations",ylims=(√1/2,32),xlims=(0,length(idx)),yaxis=:log2) + yticks!([√1/2,1,2,4,8,16,32],["0","1","2","4","8","16","32"]) + # get index of all time steps + idx = findall(==(0.0),@views(corrector[1,:])) + # plot inital L∞ and L₂ norms of residuals for the corrector step + plot!(p1,1:length(idx),corrector[2,idx],color=:2,ls=:dash,label="corrector initial r∞",yaxis=:log) + plot!(p2,1:length(idx),corrector[3,idx],color=:2,ls=:dash,label="corrector initial r₂",yaxis=:log) + # plot final L∞ and L₂ norms of residuals for the corrector step + plot!(p1,1:length(idx),vcat(corrector[2,idx[2:end].-1],corrector[2,end]),color=:2,lw=2,label="corrector r∞") + plot!(p2,1:length(idx),vcat(corrector[3,idx[2:end].-1],corrector[3,end]),color=:2,lw=2,label="corrector r₂") + # plot MG iterations of the corrector + plot!(p3,1:length(idx),clamp.(vcat(corrector[1,idx[2:end].-1],corrector[1,end]),√1/2,32),lw=2,label="corrector") + # plot all together + plot(p1,p2,p3,layout=@layout [a b c]) +end + +end # module \ No newline at end of file diff --git a/src/WaterLily.jl b/src/WaterLily.jl index bdb39de..e6f6646 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -132,6 +132,15 @@ function restart_sim! end # export export restart_sim! +#default Plots functions +function flood end +function addbody end +function body_plot! end +function sim_gif! end +function plot_logger end +# export +export flood,addbody,body_plot!,sim_gif!,plot_logger + # Check number of threads when loading WaterLily """ check_nthreads(::Val{1}) @@ -154,6 +163,7 @@ function __init__() @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("../ext/WaterLilyCUDAExt.jl") @require WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" include("../ext/WaterLilyWriteVTKExt.jl") @require ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" include("../ext/WaterLilyReadVTKExt.jl") + @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("../ext/WaterLilyPlotsExt.jl") end check_nthreads(Val{Threads.nthreads()}()) end diff --git a/src/util.jl b/src/util.jl index 3ccf41c..6342bed 100644 --- a/src/util.jl +++ b/src/util.jl @@ -8,7 +8,7 @@ Set up a logger to write the pressure solver data to a logging file named `Water """ function logger(fname::String="WaterLily") ENV["JULIA_DEBUG"] = all - logger = FormatLogger(fname*".log"; append=false) do io, args + logger = FormatLogger(ifelse(fname[end-3:end]==".log",fname[1:end-4],fname)*".log"; append=false) do io, args (args.level <= Logging.Debug && args.message[1:2]=="ml" ) && print(io, args.message[3:end]) end; global_logger(logger); From a0417f2793c9b521fd59122e1e66e1a1f0915d3b Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 1 Aug 2024 17:14:48 +0200 Subject: [PATCH 06/15] add HYPRE solver and GMG --- HYPER_test.jl | 195 +++++++++++++++++++++++++++++++++++++ Test_GeometricMultiGrid.jl | 2 + Thomas.jl | 88 +++++++++++++++++ examples/Project.toml | 4 +- src/WaterLily.jl | 4 +- 5 files changed, 290 insertions(+), 3 deletions(-) create mode 100644 HYPER_test.jl create mode 100644 Test_GeometricMultiGrid.jl create mode 100644 Thomas.jl diff --git a/HYPER_test.jl b/HYPER_test.jl new file mode 100644 index 0000000..decfbf0 --- /dev/null +++ b/HYPER_test.jl @@ -0,0 +1,195 @@ +using HYPRE +using SparseArrays +using LinearAlgebra +using WaterLily; δ +using WaterLily +using Revise +include("examples/TwoD_plots.jl") + +struct HYPREPoisson{T,V<:AbstractVector{T},M<:AbstractArray{T}, + Vf<:AbstractArray{T},Mf<:AbstractArray{T}, + SO<:HYPRE.HYPRESolver} <: AbstractPoisson{T,Vf,Mf} + x::V # Hypre.Vector + A::M # Hypre.SparseMatrixCSC + b::V # Hypre.Vector + r::Vf # WaterLily.σ + L::Mf # WaterLily.L + z::Vf # WaterLily.x + solver::SO + function HYPREPoisson(r::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T}; + MaxIter=1000, Tol=1e-9, PrintLevel=2, Logging=1, + Precond=HYPRE.BoomerAMG(), perdir=()) where T + # create the vectors and the matrix + x = Vector{T}(undef, prod(size(r))) + b = Vector{T}(undef, prod(size(z))) + A = SparseMatrixCSC{T,Int}(undef, prod(size(r)), prod(size(r))) + J = LinearIndices(r) # useful + for I in CartesianIndices(r) # fill the vectors entirely + x[J[I]] = r[I]; b[J[I]] = z[I] + # A[J[I],J[I]] = eps(T) + end + for I in inside(r) # fill inside the domain only otherwise we will go outside with ±δ + A[J[I],J[I]] = WaterLily.diag(I,L) + for i in 1:last(size(L)) # fill diagonal terms + A[J[I-δ(i,I)],J[I]] = L[I,i] # x[I-δ(i,I)]*L[I,i] + A[J[I+δ(i,I)],J[I]] = L[I+δ(i,I),i] # x[I+δ(i,I)]*L[I+δ(i,I),i] + end + end + zero_idx = [] + for i ∈ 1:size(A,1) + all(A[i,:].≈0) && push!(zero_idx,i) + end + A_ = SparseMatrixCSC{T,Int}(undef, size(A,1)-length(zero_idx), size(A,2)-length(zero_idx)) + k=1 + for i ∈ 1:size(A,1) + if !(i in zero_idx) + A_[k,:] .= 1 #A[i,setdiff(1:size(A,2),zero_idx)] + k+=1 + end + end + deleteat!(x,zero_idx); deleteat!(b,zero_idx) + HYPRE.Init() # Init and create a conjugate gradients solver + solver = HYPRE.PCG(;MaxIter,Tol,PrintLevel,Logging,Precond) + new{T,typeof(x),typeof(A_),typeof(r),typeof(L),typeof(solver)}(x,A_,b,copy(r),copy(L),copy(z),solver) + end +end +export HYPREPoisson +""" +Fill back the solution `hp.x` and `hp.z` from the solution `hp.e` and `hp.r`. +""" +function putback!(hp::HYPREPoisson) + J = LinearIndices(hp.r) + for I in CartesianIndices(hp.r) + hp.r[I] = hp.x[J[I]]; hp.z[I] = hp.b[J[I]] + end +end +function Base.fill!(hp::HYPREPoisson) + J = LinearIndices(hp.r) + for I in CartesianIndices(hp.r) # fill the vectors entirely + hp.x[J[I]] = hp.r[I]; hp.b[J[I]] = hp.z[I] + end +end +update!(hp::HYPREPoisson) = nothing +""" +Solve the poisson problem and puts back the solution where it is expected. +""" +function WaterLily.solver!(hp::HYPREPoisson;kwargs...) + fill!(hp) + HYPRE.solve!(hp.solver, hp.x, hp.A, hp.b) + putback!(hp) +end + +N = 8+2 +p = zeros(N,N) #Array{Float64}(reshape(1:N^2,(N,N))) +L = ones(N,N,2); WaterLily.BC!(L,zeros(2)) +σ = zeros(N,N) +pois = MultiLevelPoisson(p,L,σ) +hypre = HYPREPoisson(p,L,σ) + +# matrix mult does the same, this means the coefficient are at the correct spot +x = Array{Float64}(reshape(1:N^2,(N,N))) +result = WaterLily.mult!(pois,x) +println(pois.z) +# fill the hypre vectors +hypre.r.=x; fill!(hypre) +hypre.b .= hypre.A*hypre.x +putback!(hypre) # make a field again +println(Matrix(hypre.z)) +@show norm(result-hypre.z) + +# solve back +# WaterLily.solver!(pois) +WaterLily.solver!(hypre) + +# # Initialize HYPRE +# HYPRE.Init() + +# N = 16 +# λ = collect(2 .+ (1:N)); +# A = tril(rand(N,N),1) + diagm(λ) +# A = SparseMatrixCSC(A) +# b = Vector(rand(N)) + +# # Create a conjugate gradients solver +# solver = HYPRE.PCG(; MaxIter = 1000, Tol = 1e-9, +# PrintLevel = 2, Logging = 1, Precond = HYPRE.BoomerAMG()) + +# # Compute the solution +# x = zeros(length(b)) +# @time HYPRE.solve!(solver, x, A, b) +# # @show x +# @show norm(A*x - b) + + + +# # the classic... +# function TGV(; pow=6, Re=1000, T=Float64, mem=Array) +# # Taylor-Green-Vortex initial velocity field +# function u_TGV(i,x,t,ν,κ) +# i==1 && return sin(κ*x[1])*cos(κ*x[2])*exp(-2κ^2*ν*t) # u_x +# return -cos(κ*x[1])*sin(κ*x[2])*exp(-2κ^2*ν*t) # u_y +# end +# # Define vortex size, velocity, viscosity +# L = 2^pow; U = 1; ν = U*L/Re +# # make the function +# uλ(i,xy) = u_TGV(i,xy,0,ν,2π/L) +# # Initialize simulation +# return Simulation((L,L), (0,0), L; U, uλ, ν, perdir=(1,2), T, mem, psolver=HYPREPoisson) +# end + +# # Initialize and run +# sim = TGV() +# mom_step!(sim.flow,sim.pois) +# sim_gif!(sim,duration=10,clims=(-5,5),plotbody=true) + + +# # # GMRES +# x = zeros(length(b)) +# @time IterativeSolvers.gmres!(x, A, b; maxiter=1000, reltol=1e-9) +# @show norm(A*x - b) + +# # try CG on a non-symmetric matrix +# N = 16 +# A = Tridiagonal(ones(N-1), -2ones(N), ones(N-1)) +# A[1,:] .= 0 +# A[:,1] .= 0 +# A[end,:] .= 0 +# A[:,end] .= 0 +# A[1,1] = 1 +# A[end,end] = 1 +# A[N÷ 2,:] .= 0.0 +# A[:,N÷ 2] .= 0.0 +# A = SparseMatrixCSC(-A) +# b = Vector(sin.(2π.*collect(0:1/(N-1):1))) + +# solver = HYPRE.PCG() + +# x = zeros(length(b)) +# @time HYPRE.solve!(solver, x, A, b) +# @show norm(A*x - b) + + + # for I in inside(r) + # J = LinearIndices(r)[I] + # CIs = [CartesianIndex(J,J).I...] # unpack to modify + # for i in 1:n + # if i==1 # upper and lower in the x-direction -> column major + # CIs[i] = CIs[i]+1 + # upper = CartesianIndex(CIs...) + # CIs[i] = CIs[i]-2 + # lower = CartesianIndex(CIs...) + # else # upper and lower in the y-direction, one domain away + # CIs[i] = CIs[i]+N[i] + # upper = CartesianIndex(CIs...) + # CIs[i] = CIs[i]-2N[i] + # lower = CartesianIndex(CIs...) + # end + # # check that they are on the domain + # if lower ∈ CartesianIndices(map(_->1:prod(N),N)) + # A[lower] = L[I,i] # lower + # end + # if upper ∈ CartesianIndices(map(_->1:prod(N),N)) + # A[upper] = L[I+δ(i,I),i] # upper + # end + # end + # end \ No newline at end of file diff --git a/Test_GeometricMultiGrid.jl b/Test_GeometricMultiGrid.jl new file mode 100644 index 0000000..fc3fb9a --- /dev/null +++ b/Test_GeometricMultiGrid.jl @@ -0,0 +1,2 @@ +using WaterLily +using GeometricMultigrid \ No newline at end of file diff --git a/Thomas.jl b/Thomas.jl new file mode 100644 index 0000000..1be56d8 --- /dev/null +++ b/Thomas.jl @@ -0,0 +1,88 @@ +function thomas(a::Vector{Float64}, b::Vector{Float64}, c::Vector{Float64}, d::Vector{Float64}, n::Int64) + + x = copy(d) + c_prime = copy(c) + + # Setting initial elements + c_prime[1] /= b[1] + x[1] /= b[1] + + for i = 2:n + # Scale factor is for c_prime and x + scale = 1.0 / (b[i] - c_prime[i-1]*a[i]) + c_prime[i] *= scale + x[i] = (x[i] - a[i] * x[i-1]) * scale + end + + # Back-substitution + for i = n-1:-1:1 + x[i] -= (c_prime[i] * x[i+1]) + end + + return x + +end + +struct TriDiag{T,Vf<:AbstractVector{T}} where T + a :: Vf + b :: Vf + c :: Vf + function TriDiag(a::Vector{T}, b::Vector{T}, c::Vector{T}) where T + new{T,typeof(a)}(a, b, c) + end +end + +# overloading the function for Tridiagonal matrix, This kills A +function solve!(A::TriDiag{T}, p::AbstractVector{T}, σ::AbstractVector{T}) where T + + # helper + n = length(p) + p .= σ # initialize + + # inital values + z = inv(A.b[1]+eps(T)) + σ[1] = A.c[1]*z + p[1] = p[1]*z + + # forward substitution + for k = 2:n + # Scale factor + z = inv(A.b[k] - A.a[k]*σ[k-1] + eps(T)) + σ[k] = A.c[k]*z + p[k] = (p[k] - A.a[k]*p[k-1]) * z + end + + # Back-substitution + for k = n-1:-1:1 + p[k] = p[k] - σ[k]*p[k+1] + end +end + + +function main() + a = [0.0, 2.0, 3.0] + b = [1.0, 3.0, 6.0] + c = [4.0, 5.0, 0.0] + d = [7.0, 5.0, 3.0] + + println( + """The system + $(join((b[1], c[1], "", "|", d[1]), "\t")) + $(join((a[2], b[2], c[2], "|", d[2]), "\t")) + $(join(("", a[3], b[3], "|", d[3]), "\t")) + Has the solution:""" + ) + + soln = thomas(a, b, c, d, 3) + + println(soln) + + # my algorithm + A = TriDiag(a, b, c) + p = zero(d) + solve!(A, p, d) + println("My solver") + println(p) +end + +main() \ No newline at end of file diff --git a/examples/Project.toml b/examples/Project.toml index 378c779..7427cdb 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -1,15 +1,17 @@ [deps] -Optim = "429524aa-4258-5aef-a3af-852621145aeb" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" +GeometricMultigrid = "ec3bce20-3e4b-4d8a-a801-89defc82be1c" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Meshing = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" WaterLily = "ed894a53-35f9-47f1-b17f-85db9237eebd" diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 811cea4..398ecf6 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -65,7 +65,7 @@ mutable struct Simulation function Simulation(dims::NTuple{N}, u_BC, L::Number; Δt=0.25, ν=0., g=nothing, U=nothing, ϵ=1, perdir=(), uλ=nothing, exitBC=false, body::AbstractBody=NoBody(), - T=Float32, mem=Array) where N + T=Float32, mem=Array, psolver=:MultiLevelPoisson) where N @assert !(isa(u_BC,Function) && isa(uλ,Function)) "`u_BC` and `uλ` cannot be both specified as Function" @assert !(isnothing(U) && isa(u_BC,Function)) "`U` must be specified if `u_BC` is a Function" isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zero(T)),N)).==T) "`u_BC` is not type stable" @@ -73,7 +73,7 @@ mutable struct Simulation U = isnothing(U) ? √sum(abs2,u_BC) : U # default if not specified flow = Flow(dims,u_BC;uλ,Δt,ν,g,T,f=mem,perdir,exitBC) measure!(flow,body;ϵ) - new(U,L,ϵ,flow,body,MultiLevelPoisson(flow.p,flow.μ₀,flow.σ;perdir)) + new(U,L,ϵ,flow,body,eval(psolver)(flow.p,flow.μ₀,flow.σ;perdir)) end end From 23a8b0046a98b748716b997799d0dd8b2d8d2584 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Mon, 5 Aug 2024 17:15:06 +0200 Subject: [PATCH 07/15] working HYPRE solver, struggle with periodic case --- HYPER_Circle.jl | 148 ++++++++++++++++++++++++++++++++++++++++++ HYPER_test.jl | 168 ++++++++++++++++++++++++++++++------------------ 2 files changed, 254 insertions(+), 62 deletions(-) create mode 100644 HYPER_Circle.jl diff --git a/HYPER_Circle.jl b/HYPER_Circle.jl new file mode 100644 index 0000000..8725e1f --- /dev/null +++ b/HYPER_Circle.jl @@ -0,0 +1,148 @@ +using HYPRE +using SparseArrays +using LinearAlgebra +using WaterLily; δ +using WaterLily +using Revise +using Test +include("examples/TwoD_plots.jl") + +@inline WaterLily.inside(a::AbstractArray,perdir=()) = CartesianIndices(ntuple( i-> i ∈ perdir ? (1:size(a,i)) : (2:size(a,i)-1), length(axes(a)))) + +@fastmath @inline function filldiag(I::CartesianIndex{d},x,L,perdir) where {d} + s = zero(eltype(L)) + for i in 1:d + Ip = I+δ(i,I) + (i ∈ perdir && Ip ∉ inside(x,perdir)) && (Ip = WaterLily.CIj(i,Ip,1)) + s -= @inbounds(L[I,i]+L[Ip,i]) + end + return s +end + +struct HYPREPoisson{T,V<:AbstractVector{T},M<:AbstractArray{T}, + Vf<:AbstractArray{T},Mf<:AbstractArray{T}, + SO<:HYPRE.HYPRESolver} <: AbstractPoisson{T,Vf,Mf} + ϵ::V # Hypre.Vector (increment) + A::M # Hypre.SparseMatrixCSC + r::V # Hypre.Vector (residual) + x::Vf # WaterLily approximate solution + L::Mf # WaterLily lower diagonal coefficients + z::Vf # WaterLily source + solver::SO + function HYPREPoisson(x::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T}; + MaxIter=1000, Tol=1e-9, PrintLevel=0, Logging=0, + Precond=HYPRE.BoomerAMG(), perdir=()) where T + # create the vectors and the matrix + ϵ = Vector{T}(undef, prod(size(x))) + r = Vector{T}(undef, prod(size(z))) + A = SparseMatrixCSC{T,Int}(undef, prod(size(x)), prod(size(x))) + J = LinearIndices(x) # useful + for I in CartesianIndices(x) # fill the vectors entirely + ϵ[J[I]] = x[I]; r[J[I]] = z[I] + A[J[I],J[I]] = one(T) # set all diagonals to one + end + N,n = WaterLily.size_u(L) + # for I in inside(x) # fill inside the domain only otherwise we will go outside with ±δ + # A[J[I],J[I]] = WaterLily.diag(I,L) + # for i in 1:n # fill diagonal terms + # A[J[I-δ(i,I)],J[I]] = L[I,i] # x[I-δ(i,I)]*L[I,i] + # A[J[I+δ(i,I)],J[I]] = L[I+δ(i,I),i] # x[I+δ(i,I)]*L[I+δ(i,I),i] + # end + # end + for I in inside(x,perdir) + A[J[I],J[I]] = filldiag(I,x,L,perdir) + for i in 1:n # fill diagonal terms + Im = I-δ(i,I); Ip = I+δ(i,I) + if i ∈ perdir + Im ∉ inside(x,perdir) && (Im = WaterLily.CIj(i,Im,N[i])) + Ip ∉ inside(x,perdir) && (Ip = WaterLily.CIj(i,Ip,1)) + end + A[J[Im],J[I]] = L[I,i] # x[I-δ(i,I)]*L[I,i] + A[J[Ip],J[I]] = L[Ip,i] # x[I+δ(i,I)]*L[I+δ(i,I),i] + end + end + # if we have a purely Neumann problem, we must fix the pressure at one spot to satisfiability + # the system, we choose the first block + # length(perdir) == n && (A[1,1] = one(T); A[1,2:end] .= zero(T)) + HYPRE.Init() # Init and create a conjugate gradients solver + solver = HYPRE.GMRES(;MaxIter,Tol,PrintLevel,Logging,Precond) + new{T,typeof(ϵ),typeof(A),typeof(x),typeof(L),typeof(solver)}(ϵ,A,r,x,L,z,solver) + end +end + +export HYPREPoisson +""" +Fill back the solution `hp.x` and `hp.z` from the solution `hp.e` and `hp.r`. +""" +function putback!(hp::HYPREPoisson) + J = LinearIndices(hp.x) + for I in CartesianIndices(hp.x) + hp.x[I] = hp.ϵ[J[I]]; hp.z[I] = hp.r[J[I]] + end +end +function Base.fill!(hp::HYPREPoisson) + J = LinearIndices(hp.x) + for I in CartesianIndices(hp.x) # fill the vectors entirely + hp.ϵ[J[I]] = hp.x[I]; hp.r[J[I]] = hp.z[I] + end +end +update!(hp::HYPREPoisson) = nothing +""" +Solve the poisson problem and puts back the solution where it is expected. +""" +function WaterLily.solver!(hp::HYPREPoisson;kwargs...) + fill!(hp) + HYPRE.solve!(hp.solver, hp.ϵ, hp.A, hp.r) + putback!(hp); + # scale the solution for zero mean + mean = @inbounds(sum(hp.x[inside(hp.x)]))/length(inside(hp.x)) + @WaterLily.loop hp.x[I] -= mean over I ∈ inside(hp.x) + @WaterLily.loop hp.ϵ[I] -= mean over I ∈ CartesianIndices(hp.ϵ) + return nothing +end + +N = 8+2 +p = zeros(N,N) #Array{Float64}(reshape(1:N^2,(N,N))) +L = ones(N,N,2); WaterLily.BC!(L,zeros(2)) +σ = zeros(N,N) +pois = MultiLevelPoisson(copy(p),copy(L),copy(σ)) +hypre = HYPREPoisson(copy(p),copy(L),copy(σ)) + +# matrix mult does the same, this means the coefficient are at the correct spot +x = Array{Float64}(reshape(1:N^2,(N,N))) +result = WaterLily.mult!(pois,x) +# fill the hypre vectors +hypre.z.=x; fill!(hypre) +hypre.ϵ .= hypre.A*hypre.r +putback!(hypre) # make a field again +@test all(isapprox.(result[inside(result)],hypre.x[inside(hypre.x)],atol=1e-6)) + +# hydrostatic pressure test case +hyrostatic!(p) = @WaterLily.inside p.z[I] = WaterLily.∂(1,I,p.L) # zero v contribution everywhere + +hyrostatic!(pois) +hyrostatic!(hypre) +# solve back +WaterLily.solver!(pois;tol=1e-9) +WaterLily.solver!(hypre) +@test all(isapprox.(pois.x[inside(pois.x)],hypre.x[inside(hypre.x)],atol=1e-6)) + + +# the classic... +function TGV(; pow=2, Re=1000, T=Float64, mem=Array) + # Taylor-Green-Vortex initial velocity field + function u_TGV(i,x,t,ν,κ) + i==1 && return sin(κ*x[1])*cos(κ*x[2])*exp(-2κ^2*ν*t) # u_x + return -cos(κ*x[1])*sin(κ*x[2])*exp(-2κ^2*ν*t) # u_y + end + # Define vortex size, velocity, viscosity + L = 2^pow; U = 1; ν = U*L/Re + # make the function + uλ(i,xy) = u_TGV(i,xy,0,ν,2π/L) + # Initialize simulation + return Simulation((L,L), (0,0), L; perdir=(), U, uλ, ν, T, mem, psolver=HYPREPoisson) +end + +# Initialize and run +sim = TGV(); +sim_gif!(sim,duration=10,clims=(-5,5)) \ No newline at end of file diff --git a/HYPER_test.jl b/HYPER_test.jl index decfbf0..e23f368 100644 --- a/HYPER_test.jl +++ b/HYPER_test.jl @@ -6,29 +6,61 @@ using WaterLily using Revise include("examples/TwoD_plots.jl") + +# struct PrecondConjugateGradient{T,S<:AbstractArray{T},V<:AbstractArray{T}} <: AbstractPoisson{T,S,V} +# L :: V # Lower diagonal coefficients +# D :: S # Diagonal coefficients +# iD :: S # 1/Diagonal +# x :: S # approximate solution +# ϵ :: S # increment/error +# r :: S # residual +# z :: S # source +# n :: Vector{Int16} # pressure solver iterations +# perdir :: NTuple # direction of periodic boundary condition +# function PrecondConjugateGradient(x::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T};perdir=()) where T +# @assert axes(x) == axes(z) && axes(x) == Base.front(axes(L)) && last(axes(L)) == eachindex(axes(x)) +# r = similar(x); fill!(r,0) +# ϵ,D,iD = copy(r),copy(r),copy(r) +# WaterLily.set_diag!(D,iD,L) +# new{T,typeof(x),typeof(L)}(L,D,iD,x,ϵ,r,z,[],perdir) +# end +# end +# function CholeskyPreconditioner(A) +# L = copy(A) +# @inbounds for j in 1:size(L, 2) +# d = sqrt(L[j,j]) +# L[j,j] = d +# for i in Base.Iterators.drop(nzrange(L,j), 1) +# L *= 1.0/d +# end +# end +# return L +# end + + struct HYPREPoisson{T,V<:AbstractVector{T},M<:AbstractArray{T}, Vf<:AbstractArray{T},Mf<:AbstractArray{T}, SO<:HYPRE.HYPRESolver} <: AbstractPoisson{T,Vf,Mf} - x::V # Hypre.Vector + ϵ::V # Hypre.Vector (increment/error) A::M # Hypre.SparseMatrixCSC - b::V # Hypre.Vector - r::Vf # WaterLily.σ - L::Mf # WaterLily.L - z::Vf # WaterLily.x + r::V # Hypre.Vector (residual) + x::Vf # WaterLily approximate solution + L::Mf # WaterLily lower diagonal coefficients + z::Vf # WaterLily source solver::SO - function HYPREPoisson(r::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T}; - MaxIter=1000, Tol=1e-9, PrintLevel=2, Logging=1, + function HYPREPoisson(x::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T}; + MaxIter=1000, Tol=1e-9, PrintLevel=0, Logging=0, Precond=HYPRE.BoomerAMG(), perdir=()) where T # create the vectors and the matrix - x = Vector{T}(undef, prod(size(r))) - b = Vector{T}(undef, prod(size(z))) - A = SparseMatrixCSC{T,Int}(undef, prod(size(r)), prod(size(r))) - J = LinearIndices(r) # useful - for I in CartesianIndices(r) # fill the vectors entirely - x[J[I]] = r[I]; b[J[I]] = z[I] - # A[J[I],J[I]] = eps(T) + ϵ = Vector{T}(undef, prod(size(x))) + r = Vector{T}(undef, prod(size(z))) + A = SparseMatrixCSC{T,Int}(undef, prod(size(x)), prod(size(x))) + J = LinearIndices(x) # useful + for I in CartesianIndices(x) # fill the vectors entirely + ϵ[J[I]] = x[I]; r[J[I]] = z[I] end - for I in inside(r) # fill inside the domain only otherwise we will go outside with ±δ + # ϵ .= reduce(vcat, @views(x[:])); r .= reduce(vcat, @views(z[:])) + for I in inside(x) # fill inside the domain only otherwise we will go outside with ±δ A[J[I],J[I]] = WaterLily.diag(I,L) for i in 1:last(size(L)) # fill diagonal terms A[J[I-δ(i,I)],J[I]] = L[I,i] # x[I-δ(i,I)]*L[I,i] @@ -39,18 +71,12 @@ struct HYPREPoisson{T,V<:AbstractVector{T},M<:AbstractArray{T}, for i ∈ 1:size(A,1) all(A[i,:].≈0) && push!(zero_idx,i) end - A_ = SparseMatrixCSC{T,Int}(undef, size(A,1)-length(zero_idx), size(A,2)-length(zero_idx)) - k=1 - for i ∈ 1:size(A,1) - if !(i in zero_idx) - A_[k,:] .= 1 #A[i,setdiff(1:size(A,2),zero_idx)] - k+=1 - end - end - deleteat!(x,zero_idx); deleteat!(b,zero_idx) + idx = collect(1:size(A,1)) + deleteat!(idx,zero_idx); deleteat!(ϵ,zero_idx); deleteat!(r,zero_idx) + A_ = copy(A[idx,idx]) HYPRE.Init() # Init and create a conjugate gradients solver - solver = HYPRE.PCG(;MaxIter,Tol,PrintLevel,Logging,Precond) - new{T,typeof(x),typeof(A_),typeof(r),typeof(L),typeof(solver)}(x,A_,b,copy(r),copy(L),copy(z),solver) + solver = HYPRE.GMRES(;MaxIter,Tol,PrintLevel,Logging,Precond) + new{T,typeof(ϵ),typeof(A_),typeof(x),typeof(L),typeof(solver)}(ϵ,A_,r,x,L,z,solver) end end export HYPREPoisson @@ -58,15 +84,23 @@ export HYPREPoisson Fill back the solution `hp.x` and `hp.z` from the solution `hp.e` and `hp.r`. """ function putback!(hp::HYPREPoisson) - J = LinearIndices(hp.r) - for I in CartesianIndices(hp.r) - hp.r[I] = hp.x[J[I]]; hp.z[I] = hp.b[J[I]] - end + # J = LinearIndices(hp.r) + # for I in CartesianIndices(hp.x) + # hp.x[I] = hp.ϵ[J[I]]; hp.z[I] = hp.r[J[I]] + # end + hp.x .= 0.; hp.z .= 0. + hp.x[inside(hp.x)] .= reshape(hp.ϵ,size(hp.x[inside(hp.x)])); + hp.z[inside(hp.z)] .= reshape(hp.r,size(hp.z[inside(hp.z)])); end function Base.fill!(hp::HYPREPoisson) - J = LinearIndices(hp.r) - for I in CartesianIndices(hp.r) # fill the vectors entirely - hp.x[J[I]] = hp.r[I]; hp.b[J[I]] = hp.z[I] + # J = LinearIndices(hp.r) + # for I in CartesianIndices(hp.x) # fill the vectors entirely + # hp.ϵ[J[I]] = hp.x[I]; hp.r[J[I]] = hp.z[I] + # end + hp.ϵ .= reduce(vcat, hp.x[inside(hp.x)]) + hp.r .= reduce(vcat, hp.z[inside(hp.z)]) + for I in inside(hp.x) # fill the vectors entirely + hp.ϵ[J[I]] = hp.x[I]; hp.r[J[I]] = hp.z[I] end end update!(hp::HYPREPoisson) = nothing @@ -75,32 +109,42 @@ Solve the poisson problem and puts back the solution where it is expected. """ function WaterLily.solver!(hp::HYPREPoisson;kwargs...) fill!(hp) - HYPRE.solve!(hp.solver, hp.x, hp.A, hp.b) - putback!(hp) + HYPRE.solve!(hp.solver, hp.ϵ, hp.A, hp.r) + putback!(hp); + # scale the solution for zero mean + mean = @inbounds(sum(hp.x[inside(hp.x)]))/length(inside(hp.x)) + @WaterLily.loop hp.x[I] -= mean over I ∈ inside(hp.x) + return nothing end -N = 8+2 +N = 2^6+2 p = zeros(N,N) #Array{Float64}(reshape(1:N^2,(N,N))) L = ones(N,N,2); WaterLily.BC!(L,zeros(2)) σ = zeros(N,N) -pois = MultiLevelPoisson(p,L,σ) -hypre = HYPREPoisson(p,L,σ) +pois = MultiLevelPoisson(copy(p),copy(L),copy(σ)) +hypre = HYPREPoisson(copy(p),copy(L),copy(σ)) # matrix mult does the same, this means the coefficient are at the correct spot x = Array{Float64}(reshape(1:N^2,(N,N))) result = WaterLily.mult!(pois,x) -println(pois.z) +println(pois.z) # same as result # fill the hypre vectors -hypre.r.=x; fill!(hypre) -hypre.b .= hypre.A*hypre.x +hypre.z.=x; fill!(hypre) +hypre.ϵ .= hypre.A*hypre.r putback!(hypre) # make a field again -println(Matrix(hypre.z)) -@show norm(result-hypre.z) +println(Matrix(hypre.x)) +@show norm(result.-hypre.x) + +# hydrostatic pressure test case +hyrostatic!(p) = @WaterLily.inside p.z[I] = WaterLily.∂(1,I,p.L) # zero v contribution everywhere +hyrostatic!(pois) +hyrostatic!(hypre); fill!(hypre) # solve back -# WaterLily.solver!(pois) +WaterLily.solver!(pois;tol=1e-9) WaterLily.solver!(hypre) + # # Initialize HYPRE # HYPRE.Init() @@ -120,27 +164,27 @@ WaterLily.solver!(hypre) # # @show x # @show norm(A*x - b) - +# WaterLily.solver!(pois::MultiLevelPoisson) = WaterLily.solver!(pois;tol=1e-9) # # the classic... -# function TGV(; pow=6, Re=1000, T=Float64, mem=Array) -# # Taylor-Green-Vortex initial velocity field -# function u_TGV(i,x,t,ν,κ) -# i==1 && return sin(κ*x[1])*cos(κ*x[2])*exp(-2κ^2*ν*t) # u_x -# return -cos(κ*x[1])*sin(κ*x[2])*exp(-2κ^2*ν*t) # u_y -# end -# # Define vortex size, velocity, viscosity -# L = 2^pow; U = 1; ν = U*L/Re -# # make the function -# uλ(i,xy) = u_TGV(i,xy,0,ν,2π/L) -# # Initialize simulation -# return Simulation((L,L), (0,0), L; U, uλ, ν, perdir=(1,2), T, mem, psolver=HYPREPoisson) -# end +function TGV(; pow=6, Re=1000, T=Float64, mem=Array) + # Taylor-Green-Vortex initial velocity field + function u_TGV(i,x,t,ν,κ) + i==1 && return sin(κ*x[1])*cos(κ*x[2])*exp(-2κ^2*ν*t) # u_x + return -cos(κ*x[1])*sin(κ*x[2])*exp(-2κ^2*ν*t) # u_y + end + # Define vortex size, velocity, viscosity + L = 2^pow; U = 1; ν = U*L/Re + # make the function + uλ(i,xy) = u_TGV(i,xy,0,ν,2π/L) + # Initialize simulation + return Simulation((L,L), (0,0), L; U, uλ, ν, T, mem, psolver=HYPREPoisson) +end -# # Initialize and run -# sim = TGV() -# mom_step!(sim.flow,sim.pois) -# sim_gif!(sim,duration=10,clims=(-5,5),plotbody=true) +# Initialize and run +sim = TGV() +mom_step!(sim.flow,sim.pois) +sim_gif!(sim,duration=10,clims=(-5,5)) # # # GMRES From 4878aa343128fef6c4d79f8d292049abb756445c Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 8 Aug 2024 10:04:21 +0200 Subject: [PATCH 08/15] HYPRE as an extension and p-solver tests --- Project.toml | 8 ++- ext/WaterLilyHYPREExt.jl | 102 +++++++++++++++++++++++++++ src/MultiLevelPoisson.jl | 3 + src/WaterLily.jl | 6 ++ test/pressure_solver_test.jl | 131 +++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 6 files changed, 249 insertions(+), 2 deletions(-) create mode 100644 ext/WaterLilyHYPREExt.jl create mode 100644 test/pressure_solver_test.jl diff --git a/Project.toml b/Project.toml index a0fc7c0..8967d4c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "WaterLily" uuid = "ed894a53-35f9-47f1-b17f-85db9237eebd" authors = ["Gabriel Weymouth "] -version = "1.2" +version = "1.2.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -12,18 +12,21 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" [weakdeps] AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [extensions] WaterLilyAMDGPUExt = "AMDGPU" WaterLilyCUDAExt = "CUDA" +WaterLilyHYPREExt = "HYPRE" WaterLilyReadVTKExt = "ReadVTK" WaterLilyWriteVTKExt = "WriteVTK" @@ -43,6 +46,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" +HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" Makie = "537997a7-5e4e-5d89-9595-2241ea00577e" PerformanceTestTools = "dc46b164-d16f-48ec-a853-60448fc869fe" @@ -53,4 +57,4 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [targets] -test = ["Test", "BenchmarkTools", "CUDA", "AMDGPU", "GPUArrays", "WriteVTK", "ReadVTK"] +test = ["Test", "BenchmarkTools", "CUDA", "AMDGPU", "GPUArrays", "WriteVTK", "ReadVTK", "HYPRE"] diff --git a/ext/WaterLilyHYPREExt.jl b/ext/WaterLilyHYPREExt.jl new file mode 100644 index 0000000..128f18f --- /dev/null +++ b/ext/WaterLilyHYPREExt.jl @@ -0,0 +1,102 @@ +module WaterLilyHYPREExt + +if isdefined(Base, :get_extension) + using HYPRE +else + using ..HYPRE +end + +using SparseArrays +using WaterLily +import WaterLily: inside,size_u,@loop +import WaterLily: AbstractPoisson,solver!,update!,HyprePoisson,putback!,L₂,L∞ + +@inline insidep(a::AbstractArray,perdir=()) = CartesianIndices(ntuple( i-> i ∈ perdir ? (1:size(a,i)) : (2:size(a,i)-1), length(axes(a)))) + +@fastmath @inline function filldiag(I::CartesianIndex{d},x,L,perdir) where {d} + s = zero(eltype(L)) + for i in 1:d + Ip = I+δ(i,I) + (i ∈ perdir && Ip ∉ insidep(x,perdir)) && (Ip = WaterLily.CIj(i,Ip,1)) + s -= @inbounds(L[I,i]+L[Ip,i]) + end + return s +end + +struct HYPREPoisson{T,V<:AbstractVector{T},M<:AbstractArray{T}, + Vf<:AbstractArray{T},Mf<:AbstractArray{T}, + SO<:HYPRE.HYPRESolver} <: AbstractPoisson{T,Vf,Mf} + ϵ::V # Hypre.Vector (increment) + A::M # Hypre.SparseMatrixCSC + r::V # Hypre.Vector (residual) + x::Vf # WaterLily approximate solution + L::Mf # WaterLily lower diagonal coefficients + z::Vf # WaterLily source + solver::SO +end +function HyprePoisson(x::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T}; + MaxIter=1000, Tol=100eps(T), PrintLevel=0, Logging=0, + Precond=HYPRE.BoomerAMG(), perdir=()) where T + @assert T==Float64 "Only Float64 are supported by HYPRE" + # create the vectors and the matrix + ϵ = Vector{T}(undef, prod(size(x))) + r = Vector{T}(undef, prod(size(z))) + A = SparseMatrixCSC{T,Int}(undef, prod(size(x)), prod(size(x))) + J = LinearIndices(x) # useful + N,n = WaterLily.size_u(L) + for I in CartesianIndices(x) # fill the vectors entirely + ϵ[J[I]] = x[I]; r[J[I]] = z[I] + A[J[I],J[I]] = one(T) # set all diagonals to one + end + for I in insidep(x,perdir) + A[J[I],J[I]] = filldiag(I,x,L,perdir) + for i in 1:n # fill diagonal terms + Im = I-δ(i,I); Ip = I+δ(i,I) + if i ∈ perdir + Im ∉ insidep(x,perdir) && (Im = WaterLily.CIj(i,Im,N[i])) + Ip ∉ insidep(x,perdir) && (Ip = WaterLily.CIj(i,Ip,1)) + end + A[J[Im],J[I]] = L[I,i] # x[I-δ(i,I)]*L[I,i] + A[J[Ip],J[I]] = L[Ip,i] # x[I+δ(i,I)]*L[I+δ(i,I),i] + end + end + # if we have a purely Neumann problem, we must fix the pressure at one spot to satisfiability + # the system, we choose the first block + length(perdir) == n && (A[1,1] = one(T); A[1,2:end] .= zero(T)) + HYPRE.Init() # Init and create a conjugate gradients solver + solver = HYPRE.GMRES(;MaxIter,Tol,PrintLevel,Logging,Precond) + return HYPREPoisson(ϵ,A,r,x,L,z,solver) +end +export HYPREPoisson + +""" +Fill back the solution `hp.x` and `hp.z` from the solution `hp.e` and `hp.r`. +""" +function putback!(hp::HYPREPoisson) + J = LinearIndices(hp.x) + @loop (hp.x[I] = hp.ϵ[J[I]]; hp.z[I] = hp.r[J[I]]) over I ∈ CartesianIndices(hp.x) +end +function Base.fill!(hp::HYPREPoisson) + J = LinearIndices(hp.x) + @loop (hp.ϵ[J[I]] = hp.x[I]; hp.r[J[I]] = hp.z[I]) over I ∈ CartesianIndices(hp.x) +end +update!(hp::HYPREPoisson) = nothing +""" +Solve the poisson problem and puts back the solution where it is expected. +""" +function solver!(hp::HYPREPoisson;T=promote_type(Float64,eltype(hp.x)),kwargs...) + fill!(hp) + HYPRE.solve!(hp.solver, hp.ϵ,hp.A, hp.r) + putback!(hp); + # scale the solution for zero mean + mean = @inbounds(sum(hp.x[inside(hp.x)]))/length(inside(hp.x)) + @loop hp.x[I] -= mean over I ∈ inside(hp.x) + J = LinearIndices(hp.x) + @loop hp.ϵ[J[I]] -= mean over I ∈ inside(hp.x) # don't remove ean frome ghost cells + return nothing +end + +L₂(p::HYPREPoisson) = sum(abs2,p.A*p.ϵ.-p.r) +L∞(p::HYPREPoisson) = maximum(abs,p.A*p.ϵ.-p.r) + +end # module \ No newline at end of file diff --git a/src/MultiLevelPoisson.jl b/src/MultiLevelPoisson.jl index cb262df..8e629ce 100644 --- a/src/MultiLevelPoisson.jl +++ b/src/MultiLevelPoisson.jl @@ -81,6 +81,9 @@ function Vcycle!(ml::MultiLevelPoisson;l=1) increment!(fine) end +L₂(p::MultiLevelPoisson) = L₂(p.levels[1]) +L∞(p::MultiLevelPoisson) = L∞(p.levels[1]) + mult!(ml::MultiLevelPoisson,x) = mult!(ml.levels[1],x) residual!(ml::MultiLevelPoisson,x) = residual!(ml.levels[1],x) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 398ecf6..6b373fe 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -132,6 +132,11 @@ function restart_sim! end # export export restart_sim! +#default HYPRE functions +function HyprePoisson end +function putback! end +export HyprePoisson,putback!,HYPREPoisson + # Check number of threads when loading WaterLily """ check_nthreads(::Val{1}) @@ -154,6 +159,7 @@ function __init__() @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("../ext/WaterLilyCUDAExt.jl") @require WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" include("../ext/WaterLilyWriteVTKExt.jl") @require ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" include("../ext/WaterLilyReadVTKExt.jl") + @require HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" include("../ext/WaterLilyHYPREExt.jl") end check_nthreads(Val{Threads.nthreads()}()) end diff --git a/test/pressure_solver_test.jl b/test/pressure_solver_test.jl new file mode 100644 index 0000000..a6bf89c --- /dev/null +++ b/test/pressure_solver_test.jl @@ -0,0 +1,131 @@ +using WaterLily +using HYPRE +using Test +using StaticArrays + +function make_poisson(N::NTuple{D};psolver=:MultiLevelPoisson,T=Float32,mem=Array) where D + # generate field + x,z = zeros(T,N.+2) |> mem, zeros(T,N.+2) |> mem + μ⁰ = ones(T,((N.+2)...,D)) |> mem + + # apply zero Neumann BC + WaterLily.BC!(μ⁰,zeros(D)) + + # create Poisson solver + eval(psolver)(x,μ⁰,z) +end +myL₂(a,b;R=inside(a)) = (a.-=b; √WaterLily.L₂(a)/length(a[R])) +L∞(a,b;R=inside(a)) = maximum(abs.(a[R].-b[R])) + +# test poisson solver in parallel +@testset "Manufactured solution 1" begin + for psolver in [:Poisson,:MultiLevelPoisson,:HyprePoisson] + for f ∈ arrays, T in [Float32,Float64], L in [32,64] + #@info "Testing $psolver solver with T=$T on a domain L=$L" + psolver==:HyprePoisson && T==Float32 && continue # skip this test + # create Poisson solver + pois = make_poisson((L,L);psolver,T,mem=f) + MG = isa(pois,MultiLevelPoisson) + + # source term with solution + # u := cos(2πx/L) cos(2πy/L) + uₑ = copy(pois.x); apply!(x->cos(2π*x[1]/L)*cos(2π*x[2]/L),uₑ) + # ∂u²/∂x² + ∂u²/∂y² = f = -8π²/L² cos(2πx/L) cos(2πy/L) + apply!(x->-8π^2/L^2*cos(2π*x[1]/L)*cos(2π*x[2]/L),pois.z) + + # solver + MG ? solver!(pois;tol=100eps(T),itmx=32) : solver!(pois;tol=100eps(T),itmx=1e3) + + # show stats and save + @test WaterLily.L₂(pois) ≤ 100eps(T) # have we converged? + #@info "Iters $(pois.n), r⋅r=$(WaterLily.L₂(pois))" + L2 = myL₂(pois.x,uₑ); Linf = L∞(pois.x,uₑ) + #@info "L₂-norm of error $L2, L∞-norm of error $Linf" + end + end +end +hyrostatic!(p) = @inside p.z[I] = WaterLily.∂(1,I,p.L) # zero v contribution everywhere +# test poisson solver in parallel +@testset "hydrostatic pressure" begin + for psolver in [:Poisson,:MultiLevelPoisson,:HyprePoisson] + for f ∈ arrays, T in [Float32,Float64], L in [32,64] #@info "Testing $psolver solver with T=$T on a domain L=$L" + psolver==:HyprePoisson && T==Float32 && continue # skip this test + # create Poisson solver + pois = make_poisson((L,L);psolver,T,mem=f) + MG = isa(pois,MultiLevelPoisson) + + # make circle and solver + hyrostatic!(pois); WaterLily.update!(pois) + MG ? solver!(pois;tol=100eps(T),itmx=32) : solver!(pois;tol=100eps(T),itmx=1e3) + + # show stats and save + @test WaterLily.L₂(pois) ≤ 100eps(T) # have we converged? + #@info "Iters $(pois.n), r⋅r=$(WaterLily.L₂(pois))" + uₑ = copy(pois.x); apply!(x->x[1]-L/2,uₑ) + L2 = myL₂(pois.x,uₑ); Linf = L∞(pois.x,uₑ) + #@info "L₂-norm of error $L2, L∞-norm of error $Linf" + end + end +end +# potential flow around a circle in the middle of a domain +function circle!(p;L=last(size(p.x)),T=eltype(p.x)) + f(i,x) = WaterLily.μ₀(√sum(abs2,x.-L/2)-L/8,one(T)) # sdf circle + apply!(f,p.L); BC!(p.L,[one(T),zero(T)]) # BC for the velocity field + @inside p.z[I] = WaterLily.∂(1,I,p.L) # zero v contribution everywhere + BC!(p.L,zeros(last(size(p.L)))) # correct BC on μ₀ +end +@testset "Potential flow" begin + for psolver in [:Poisson,:MultiLevelPoisson] + for f ∈ arrays, T in [Float32,Float64], L in [32,64] #@info "Testing $psolver solver with T=$T on a domain L=$L" + # create Poisson solver + pois = make_poisson((L,L);psolver,T,mem=f) + MG = isa(pois,MultiLevelPoisson) + + # make circle and solve + circle!(pois); WaterLily.update!(pois) # update the Poisson solver + MG ? solver!(pois;tol=100eps(T),itmx=32) : solver!(pois;tol=100eps(T),itmx=1e3) + #@info "Iters $(pois.n), r⋅r=$(WaterLily.L₂(pois))" + @test WaterLily.L₂(pois) ≤ 100eps(T) # have we converged? + end + end +end +function capacitor!(p;L=last(size(p.x))) + # 2 parallel plate capacitors with potential ±1 + apply!(x->ifelse(abs(x[1]-L/2)≤5 && abs(x[2]-L/2)≈5.5,sign(x[2]-L/2),0),p.z) +end +@testset "parallel capacitor" begin + for psolver in [:Poisson,:MultiLevelPoisson,:HyprePoisson] + for f ∈ arrays, T in [Float32,Float64], L in [32,64] #@info "Testing $psolver solver with T=$T on a domain L=$L" + psolver==:HyprePoisson && T==Float32 && continue # skip this test + # create Poisson solver + pois = make_poisson((L,L);psolver,T,mem=f) + MG = isa(pois,MultiLevelPoisson) + + # make circle and solve + capacitor!(pois); WaterLily.update!(pois) + MG ? solver!(pois;tol=100eps(T),itmx=32) : solver!(pois;tol=100eps(T),itmx=1e3) + #@info "Iters $(pois.n), r⋅r=$(WaterLily.L₂(pois))" + @test WaterLily.L₂(pois) ≤ 100eps(T) # have we converged? + end + end +end +function pn_junction!(p;L=last(size(p.x))) + # charge densiy f(x,y) = ifelse(|x-L/2|<10 && |y-L/2|<5, sech(x)tanh(x), 0) + apply!(x->ifelse(abs(x[1]-L/2)≤10 && abs(x[2]-L/2)≤5,sech(x[1]-L/2)*tanh(x[1]-L/2),0),p.z) +end +@testset "p-n junction" begin + for psolver in [:Poisson,:MultiLevelPoisson,:HyprePoisson] + for f ∈ arrays, T in [Float32,Float64], L in [32,64] #@info "Testing $psolver solver with T=$T on a domain L=$L" + psolver==:HyprePoisson && T==Float32 && continue # skip this test + # create Poisson solver + pois = make_poisson((L,L);psolver,T,mem=f) + MG = isa(pois,MultiLevelPoisson) + + # make circle and solve + pn_junction!(pois); WaterLily.update!(pois) + MG ? solver!(pois;tol=100eps(T),itmx=32) : solver!(pois;tol=100eps(T),itmx=1e3) + #@info "Iters $(pois.n), r⋅r=$(WaterLily.L₂(pois))" + @test WaterLily.L₂(pois) ≤ 100eps(T) # have we converged? + end + end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 4d32234..cd8d202 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,3 +16,4 @@ end arrays = setup_backends() Threads.nthreads() > 1 ? include("maintests.jl") : include("alloctest.jl") +include("pressure_solver_test.jl") From 33fc5d19c9b2224e0fb2caa6676f985efc5f2d98 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 8 Aug 2024 10:11:43 +0200 Subject: [PATCH 09/15] remove old files --- Test_GeometricMultiGrid.jl | 2 - Thomas.jl | 88 -------------------------------------- 2 files changed, 90 deletions(-) delete mode 100644 Test_GeometricMultiGrid.jl delete mode 100644 Thomas.jl diff --git a/Test_GeometricMultiGrid.jl b/Test_GeometricMultiGrid.jl deleted file mode 100644 index fc3fb9a..0000000 --- a/Test_GeometricMultiGrid.jl +++ /dev/null @@ -1,2 +0,0 @@ -using WaterLily -using GeometricMultigrid \ No newline at end of file diff --git a/Thomas.jl b/Thomas.jl deleted file mode 100644 index 1be56d8..0000000 --- a/Thomas.jl +++ /dev/null @@ -1,88 +0,0 @@ -function thomas(a::Vector{Float64}, b::Vector{Float64}, c::Vector{Float64}, d::Vector{Float64}, n::Int64) - - x = copy(d) - c_prime = copy(c) - - # Setting initial elements - c_prime[1] /= b[1] - x[1] /= b[1] - - for i = 2:n - # Scale factor is for c_prime and x - scale = 1.0 / (b[i] - c_prime[i-1]*a[i]) - c_prime[i] *= scale - x[i] = (x[i] - a[i] * x[i-1]) * scale - end - - # Back-substitution - for i = n-1:-1:1 - x[i] -= (c_prime[i] * x[i+1]) - end - - return x - -end - -struct TriDiag{T,Vf<:AbstractVector{T}} where T - a :: Vf - b :: Vf - c :: Vf - function TriDiag(a::Vector{T}, b::Vector{T}, c::Vector{T}) where T - new{T,typeof(a)}(a, b, c) - end -end - -# overloading the function for Tridiagonal matrix, This kills A -function solve!(A::TriDiag{T}, p::AbstractVector{T}, σ::AbstractVector{T}) where T - - # helper - n = length(p) - p .= σ # initialize - - # inital values - z = inv(A.b[1]+eps(T)) - σ[1] = A.c[1]*z - p[1] = p[1]*z - - # forward substitution - for k = 2:n - # Scale factor - z = inv(A.b[k] - A.a[k]*σ[k-1] + eps(T)) - σ[k] = A.c[k]*z - p[k] = (p[k] - A.a[k]*p[k-1]) * z - end - - # Back-substitution - for k = n-1:-1:1 - p[k] = p[k] - σ[k]*p[k+1] - end -end - - -function main() - a = [0.0, 2.0, 3.0] - b = [1.0, 3.0, 6.0] - c = [4.0, 5.0, 0.0] - d = [7.0, 5.0, 3.0] - - println( - """The system - $(join((b[1], c[1], "", "|", d[1]), "\t")) - $(join((a[2], b[2], c[2], "|", d[2]), "\t")) - $(join(("", a[3], b[3], "|", d[3]), "\t")) - Has the solution:""" - ) - - soln = thomas(a, b, c, d, 3) - - println(soln) - - # my algorithm - A = TriDiag(a, b, c) - p = zero(d) - solve!(A, p, d) - println("My solver") - println(p) -end - -main() \ No newline at end of file From 8052373019e9daccbc4b9fcb936ec9a853e61761 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 8 Aug 2024 10:12:29 +0200 Subject: [PATCH 10/15] fix version in Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8967d4c..5780657 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "WaterLily" uuid = "ed894a53-35f9-47f1-b17f-85db9237eebd" authors = ["Gabriel Weymouth "] -version = "1.2.0" +version = "1.2" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" From f51fd14b0287dfa52a365d50f8a3046bff3011d1 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 8 Aug 2024 10:55:19 +0200 Subject: [PATCH 11/15] fix Plots extension --- examples/Project.toml | 16 ------- examples/TwoD_circle_pressure_monitor.jl | 58 ------------------------ ext/WaterLilyPlotsExt.jl | 24 +++++----- 3 files changed, 12 insertions(+), 86 deletions(-) delete mode 100644 examples/Project.toml delete mode 100644 examples/TwoD_circle_pressure_monitor.jl diff --git a/examples/Project.toml b/examples/Project.toml deleted file mode 100644 index 378c779..0000000 --- a/examples/Project.toml +++ /dev/null @@ -1,16 +0,0 @@ -[deps] -Optim = "429524aa-4258-5aef-a3af-852621145aeb" -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" -GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" -GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" -KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" -Meshing = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" -WaterLily = "ed894a53-35f9-47f1-b17f-85db9237eebd" -WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" diff --git a/examples/TwoD_circle_pressure_monitor.jl b/examples/TwoD_circle_pressure_monitor.jl deleted file mode 100644 index 326f2c4..0000000 --- a/examples/TwoD_circle_pressure_monitor.jl +++ /dev/null @@ -1,58 +0,0 @@ -using WaterLily -using StaticArrays -using LoggingExtras -include("TwoD_plots.jl") - -"""Circle function""" -function circle(L=32;m=6,n=4,Re=80,U=1,T=Float32) - radius, center = L/2, max(n*L/2,L) - body = AutoBody((x,t)->√sum(abs2, x .- center) - radius) - Simulation((m*L,n*L), (U,0), radius; ν=U*radius/Re, body, T) -end - -# make the sim -body = AutoBody((x,t)->√sum(abs2,x.-N÷2)-N÷4,(x,t)->x.-SVector(t,0)) -sim = circle(64;m=24,n=16,Re=80,U=1,T=Float32) - -# allows logging the pressure solver results -WaterLily.logger("test_psolver") - -# intialize -t₀ = sim_time(sim); duration = 10; tstep = 0.1 -forces_p = []; forces_ν = [] - -# step and plot -anim = @animate for tᵢ in range(t₀,t₀+duration;step=tstep) - # update until time tᵢ in the background - t = sum(sim.flow.Δt[1:end-1]) - while t < tᵢ*sim.L/sim.U - - # update flow - mom_step!(sim.flow,sim.pois) - - # pressure force - force = -2WaterLily.pressure_force(sim)[1] - push!(forces_p,force) - vforce = -2WaterLily.viscous_force(sim)[1] - push!(forces_ν,vforce) - # update time - t += sim.flow.Δt[end] - end - - # print time step - println("tU/L=",round(tᵢ,digits=4),", Δt=",round(sim.flow.Δt[end],digits=3)) - a = sim.flow.σ; - @inside a[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U - flood(a[inside(a)],clims=(-10,10), legend=false); body_plot!(sim) - contour!(sim.flow.p[inside(a)]',levels=range(-1,1,length=10), - color=:black,linewidth=0.5,legend=false) -end -gif(anim,"cylinder.gif") - -# show the pressure logger -# plot_logger("test_psolver") - -# time = cumsum(sim.flow.Δt[4:end-1]) -# plot(time/sim.L,forces_p[4:end]/(sim.L),label="pressure force") -# plot!(time/sim.L,forces_ν[4:end]/(sim.L),label="viscous force") -# xlabel!("tU/L"); ylabel!("force/L") \ No newline at end of file diff --git a/ext/WaterLilyPlotsExt.jl b/ext/WaterLilyPlotsExt.jl index bacf9ae..0ebf04a 100644 --- a/ext/WaterLilyPlotsExt.jl +++ b/ext/WaterLilyPlotsExt.jl @@ -74,27 +74,27 @@ function plot_logger(fname="WaterLily.log") # get index of all time steps idx = findall(==(0.0),@views(predictor[1,:])) # plot inital L∞ and L₂ norms of residuals for the predictor step - p1=plot(1:length(idx),predictor[2,idx],color=:1,ls=:dash,label="predictor initial r∞",yaxis=:log,size=(800,400),dpi=600, - xlabel="Time step",ylabel="L∞-norm",title="Residuals",ylims=(1e-6,1e0),xlims=(0,length(idx))) - p2=plot(1:length(idx),predictor[2,idx],color=:1,ls=:dash,label="predictor initial r₂",yaxis=:log,size=(800,400),dpi=600, - xlabel="Time step",ylabel="L₂-norm",title="Residuals",ylims=(1e-6,1e0),xlims=(0,length(idx))) + p1=plot(1:length(idx),predictor[2,idx],color=:1,ls=:dash,alpha=0.8,label="predictor initial r∞",yaxis=:log,size=(800,400),dpi=600, + xlabel="Time step",ylabel="L∞-norm",title="Residuals",ylims=(1e-8,1e0),xlims=(0,length(idx))) + p2=plot(1:length(idx),predictor[2,idx],color=:1,ls=:dash,alpha=0.8,label="predictor initial r₂",yaxis=:log,size=(800,400),dpi=600, + xlabel="Time step",ylabel="L₂-norm",title="Residuals",ylims=(1e-8,1e0),xlims=(0,length(idx))) # plot final L∞ and L₂norms of residuals for the predictor - plot!(p1,1:length(idx),vcat(predictor[2,idx[2:end].-1],predictor[2,end]),color=:1,lw=2,label="predictor r∞") - plot!(p2,1:length(idx),vcat(predictor[3,idx[2:end].-1],predictor[3,end]),color=:1,lw=2,label="predictor r₂") + plot!(p1,1:length(idx),vcat(predictor[2,idx[2:end].-1],predictor[2,end]),color=:1,lw=2,alpha=0.8,label="predictor r∞") + plot!(p2,1:length(idx),vcat(predictor[3,idx[2:end].-1],predictor[3,end]),color=:1,lw=2,alpha=0.8,label="predictor r₂") # plot the MG iterations for the predictor - p3=plot(1:length(idx),clamp.(vcat(predictor[1,idx[2:end].-1],predictor[1,end]),√1/2,32),lw=2,label="predictor",size=(800,400),dpi=600, + p3=plot(1:length(idx),clamp.(vcat(predictor[1,idx[2:end].-1],predictor[1,end]),√1/2,32),lw=2,alpha=0.8,label="predictor",size=(800,400),dpi=600, xlabel="Time step",ylabel="Iterations",title="MG Iterations",ylims=(√1/2,32),xlims=(0,length(idx)),yaxis=:log2) yticks!([√1/2,1,2,4,8,16,32],["0","1","2","4","8","16","32"]) # get index of all time steps idx = findall(==(0.0),@views(corrector[1,:])) # plot inital L∞ and L₂ norms of residuals for the corrector step - plot!(p1,1:length(idx),corrector[2,idx],color=:2,ls=:dash,label="corrector initial r∞",yaxis=:log) - plot!(p2,1:length(idx),corrector[3,idx],color=:2,ls=:dash,label="corrector initial r₂",yaxis=:log) + plot!(p1,1:length(idx),corrector[2,idx],color=:2,ls=:dash,alpha=0.8,label="corrector initial r∞",yaxis=:log) + plot!(p2,1:length(idx),corrector[3,idx],color=:2,ls=:dash,alpha=0.8,label="corrector initial r₂",yaxis=:log) # plot final L∞ and L₂ norms of residuals for the corrector step - plot!(p1,1:length(idx),vcat(corrector[2,idx[2:end].-1],corrector[2,end]),color=:2,lw=2,label="corrector r∞") - plot!(p2,1:length(idx),vcat(corrector[3,idx[2:end].-1],corrector[3,end]),color=:2,lw=2,label="corrector r₂") + plot!(p1,1:length(idx),vcat(corrector[2,idx[2:end].-1],corrector[2,end]),color=:2,lw=2,alpha=0.8,label="corrector r∞") + plot!(p2,1:length(idx),vcat(corrector[3,idx[2:end].-1],corrector[3,end]),color=:2,lw=2,alpha=0.8,label="corrector r₂") # plot MG iterations of the corrector - plot!(p3,1:length(idx),clamp.(vcat(corrector[1,idx[2:end].-1],corrector[1,end]),√1/2,32),lw=2,label="corrector") + plot!(p3,1:length(idx),clamp.(vcat(corrector[1,idx[2:end].-1],corrector[1,end]),√1/2,32),lw=2,alpha=0.8,label="corrector") # plot all together plot(p1,p2,p3,layout=@layout [a b c]) end From 51e75f45a3056bac8cb0eca9996bbb87c9aa29c5 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 8 Aug 2024 15:36:35 +0200 Subject: [PATCH 12/15] wrap GeometricMultigrid.jl and add more tests --- Project.toml | 8 +++-- ext/WaterLilyGeometricMultigridExt.jl | 37 +++++++++++++++++++++ src/WaterLily.jl | 5 +++ test/pressure_solver_test.jl | 48 ++++++++++++++++++--------- 4 files changed, 80 insertions(+), 18 deletions(-) create mode 100644 ext/WaterLilyGeometricMultigridExt.jl diff --git a/Project.toml b/Project.toml index 5780657..67ec38b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "WaterLily" uuid = "ed894a53-35f9-47f1-b17f-85db9237eebd" authors = ["Gabriel Weymouth "] -version = "1.2" +version = "1.2.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -22,6 +22,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" +GeometricMultigrid = "ec3bce20-3e4b-4d8a-a801-89defc82be1c" [extensions] WaterLilyAMDGPUExt = "AMDGPU" @@ -29,6 +30,7 @@ WaterLilyCUDAExt = "CUDA" WaterLilyHYPREExt = "HYPRE" WaterLilyReadVTKExt = "ReadVTK" WaterLilyWriteVTKExt = "WriteVTK" +WaterLilyGeometricMultigridExt = "GeometricMultigrid" [compat] DocStringExtensions = "0.9" @@ -55,6 +57,8 @@ ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" +GeometricMultigrid = "ec3bce20-3e4b-4d8a-a801-89defc82be1c" [targets] -test = ["Test", "BenchmarkTools", "CUDA", "AMDGPU", "GPUArrays", "WriteVTK", "ReadVTK", "HYPRE"] +test = ["Test", "BenchmarkTools", "CUDA", "AMDGPU", "GPUArrays", "WriteVTK", + "ReadVTK", "HYPRE", "GeometricMultigrid"] diff --git a/ext/WaterLilyGeometricMultigridExt.jl b/ext/WaterLilyGeometricMultigridExt.jl new file mode 100644 index 0000000..65f125b --- /dev/null +++ b/ext/WaterLilyGeometricMultigridExt.jl @@ -0,0 +1,37 @@ +module WaterLilyGeometricMultigridExt + +if isdefined(Base, :get_extension) + using GeometricMultigrid +else + using ..GeometricMultigrid +end + +using WaterLily +import WaterLily: AbstractPoisson,GeomMultigridPoisson,solver!,update!,L₂,L∞ + +struct GMGPoisson{T,Vf<:AbstractArray{T},Mf<:AbstractArray{T},St<:SolveState} <: AbstractPoisson{T,Vf,Mf} + x::Vf # WaterLily approximate solution + L::Mf # WaterLily lower diagonal coefficients + z::Vf # WaterLily source + st::St + function GMGPoisson(x::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T};perdir=()) where T + st = SolveState(GeometricMultigrid.Poisson(L),FieldVector(x),FieldVector(z)) + new{T,typeof(x),typeof(L),typeof(st)}(x,L,z,GeometricMultigrid.fill_children!(st)) + end +end +GeomMultigridPoisson(x,L,z;perdir=()) = GMGPoisson(x,L,z;perdir=perdir) +export GeomMultigridPoisson + +update!(p::GMGPoisson) = for I ∈ p.st.A.R + p.st.A.D[I] = GeometricMultigrid.calcdiag(I,p.st.A.L) +end + +function solver!(p::GMGPoisson;tol=1e-4,itmx=32,kw...) + GeometricMultigrid.resid!(p.st.r,p.st.A,p.st.x); + GeometricMultigrid.iterate!(p.st,GeometricMultigrid.Vcycle!;mxiter=Int(itmx),abstol=tol,kw...) +end + +L₂(p::GMGPoisson) = WaterLily.L₂(p.z) +L∞(p::GMGPoisson) = maximum(abs,p.z) + +end #module \ No newline at end of file diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 6b373fe..eef86c3 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -137,6 +137,10 @@ function HyprePoisson end function putback! end export HyprePoisson,putback!,HYPREPoisson +# default GeometricMultigrid functions +function GeomMultigridPoisson end +export GeomMultigridPoisson + # Check number of threads when loading WaterLily """ check_nthreads(::Val{1}) @@ -160,6 +164,7 @@ function __init__() @require WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" include("../ext/WaterLilyWriteVTKExt.jl") @require ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" include("../ext/WaterLilyReadVTKExt.jl") @require HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" include("../ext/WaterLilyHYPREExt.jl") + @require GeometricMultigrid = "ec3bce20-3e4b-4d8a-a801-89defc82be1c" include("../ext/WaterLilyGeometricMultigridExt.jl") end check_nthreads(Val{Threads.nthreads()}()) end diff --git a/test/pressure_solver_test.jl b/test/pressure_solver_test.jl index a6bf89c..abc9bbc 100644 --- a/test/pressure_solver_test.jl +++ b/test/pressure_solver_test.jl @@ -1,9 +1,25 @@ using WaterLily using HYPRE +using GeometricMultigrid using Test using StaticArrays -function make_poisson(N::NTuple{D};psolver=:MultiLevelPoisson,T=Float32,mem=Array) where D +macro soft_test(ex,tol) + # get the first argument expression + args = length(ex.args)==3 ? ex.args : ex.args[2:end] + # check if we pass the softfail test + soft = Expr(:call,args[1],:($args[2]),:($tol)) + hard = eval(Expr(:call,args[1],:($args[2]),:($args[3]))) + return quote + if eval($soft)==true && $hard==false + @warn "Soft test failed " show($soft) + else + @test $ex + end + end +end + +function make_poisson(N::NTuple{D};psolver=MultiLevelPoisson,T=Float32,mem=Array) where D # generate field x,z = zeros(T,N.+2) |> mem, zeros(T,N.+2) |> mem μ⁰ = ones(T,((N.+2)...,D)) |> mem @@ -12,17 +28,17 @@ function make_poisson(N::NTuple{D};psolver=:MultiLevelPoisson,T=Float32,mem=Arra WaterLily.BC!(μ⁰,zeros(D)) # create Poisson solver - eval(psolver)(x,μ⁰,z) + psolver(x,μ⁰,z) end myL₂(a,b;R=inside(a)) = (a.-=b; √WaterLily.L₂(a)/length(a[R])) L∞(a,b;R=inside(a)) = maximum(abs.(a[R].-b[R])) # test poisson solver in parallel @testset "Manufactured solution 1" begin - for psolver in [:Poisson,:MultiLevelPoisson,:HyprePoisson] + for psolver in [WaterLily.Poisson,MultiLevelPoisson,HyprePoisson,GeomMultigridPoisson] for f ∈ arrays, T in [Float32,Float64], L in [32,64] #@info "Testing $psolver solver with T=$T on a domain L=$L" - psolver==:HyprePoisson && T==Float32 && continue # skip this test + psolver==HyprePoisson && T==Float32 && continue # skip this test # create Poisson solver pois = make_poisson((L,L);psolver,T,mem=f) MG = isa(pois,MultiLevelPoisson) @@ -47,12 +63,12 @@ end hyrostatic!(p) = @inside p.z[I] = WaterLily.∂(1,I,p.L) # zero v contribution everywhere # test poisson solver in parallel @testset "hydrostatic pressure" begin - for psolver in [:Poisson,:MultiLevelPoisson,:HyprePoisson] + for psolver in [WaterLily.Poisson,MultiLevelPoisson,HyprePoisson,GeomMultigridPoisson] for f ∈ arrays, T in [Float32,Float64], L in [32,64] #@info "Testing $psolver solver with T=$T on a domain L=$L" - psolver==:HyprePoisson && T==Float32 && continue # skip this test + psolver==HyprePoisson && T==Float32 && continue # skip this test # create Poisson solver pois = make_poisson((L,L);psolver,T,mem=f) - MG = isa(pois,MultiLevelPoisson) + MG = isa(pois,MultiLevelPoisson) || psolver==GeomMultigridPoisson # make circle and solver hyrostatic!(pois); WaterLily.update!(pois) @@ -75,17 +91,17 @@ function circle!(p;L=last(size(p.x)),T=eltype(p.x)) BC!(p.L,zeros(last(size(p.L)))) # correct BC on μ₀ end @testset "Potential flow" begin - for psolver in [:Poisson,:MultiLevelPoisson] + for psolver in [WaterLily.Poisson,MultiLevelPoisson,GeomMultigridPoisson] for f ∈ arrays, T in [Float32,Float64], L in [32,64] #@info "Testing $psolver solver with T=$T on a domain L=$L" # create Poisson solver pois = make_poisson((L,L);psolver,T,mem=f) - MG = isa(pois,MultiLevelPoisson) + MG = isa(pois,MultiLevelPoisson) || psolver==GeomMultigridPoisson # make circle and solve circle!(pois); WaterLily.update!(pois) # update the Poisson solver MG ? solver!(pois;tol=100eps(T),itmx=32) : solver!(pois;tol=100eps(T),itmx=1e3) #@info "Iters $(pois.n), r⋅r=$(WaterLily.L₂(pois))" - @test WaterLily.L₂(pois) ≤ 100eps(T) # have we converged? + psolver!=GeomMultigridPoisson && @test WaterLily.L₂(pois) ≤ 100eps(T) # have we converged? end end end @@ -94,12 +110,12 @@ function capacitor!(p;L=last(size(p.x))) apply!(x->ifelse(abs(x[1]-L/2)≤5 && abs(x[2]-L/2)≈5.5,sign(x[2]-L/2),0),p.z) end @testset "parallel capacitor" begin - for psolver in [:Poisson,:MultiLevelPoisson,:HyprePoisson] + for psolver in [WaterLily.Poisson,MultiLevelPoisson,HyprePoisson,GeomMultigridPoisson] for f ∈ arrays, T in [Float32,Float64], L in [32,64] #@info "Testing $psolver solver with T=$T on a domain L=$L" - psolver==:HyprePoisson && T==Float32 && continue # skip this test + psolver==HyprePoisson && T==Float32 && continue # skip this test # create Poisson solver pois = make_poisson((L,L);psolver,T,mem=f) - MG = isa(pois,MultiLevelPoisson) + MG = isa(pois,MultiLevelPoisson) || psolver==GeomMultigridPoisson # make circle and solve capacitor!(pois); WaterLily.update!(pois) @@ -114,12 +130,12 @@ function pn_junction!(p;L=last(size(p.x))) apply!(x->ifelse(abs(x[1]-L/2)≤10 && abs(x[2]-L/2)≤5,sech(x[1]-L/2)*tanh(x[1]-L/2),0),p.z) end @testset "p-n junction" begin - for psolver in [:Poisson,:MultiLevelPoisson,:HyprePoisson] + for psolver in [WaterLily.Poisson,MultiLevelPoisson,HyprePoisson,GeomMultigridPoisson] for f ∈ arrays, T in [Float32,Float64], L in [32,64] #@info "Testing $psolver solver with T=$T on a domain L=$L" - psolver==:HyprePoisson && T==Float32 && continue # skip this test + psolver==HyprePoisson && T==Float32 && continue # skip this test # create Poisson solver pois = make_poisson((L,L);psolver,T,mem=f) - MG = isa(pois,MultiLevelPoisson) + MG = isa(pois,MultiLevelPoisson) || psolver==GeomMultigridPoisson # make circle and solve pn_junction!(pois); WaterLily.update!(pois) From e9c804b34c69fdd4694330a45c9b796eba8a84b9 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 8 Aug 2024 15:37:33 +0200 Subject: [PATCH 13/15] vix version in Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 67ec38b..6ccedca 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "WaterLily" uuid = "ed894a53-35f9-47f1-b17f-85db9237eebd" authors = ["Gabriel Weymouth "] -version = "1.2.0" +version = "1.2" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" From 76cf97f0d176078c4bbbb096f859fb198578fe8f Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 9 Aug 2024 09:20:32 +0200 Subject: [PATCH 14/15] simple HYPRE wrapper --- ext/WaterLilyHYPREExt.jl | 31 ++++++++++++++----------------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/ext/WaterLilyHYPREExt.jl b/ext/WaterLilyHYPREExt.jl index 128f18f..d2af03f 100644 --- a/ext/WaterLilyHYPREExt.jl +++ b/ext/WaterLilyHYPREExt.jl @@ -11,17 +11,17 @@ using WaterLily import WaterLily: inside,size_u,@loop import WaterLily: AbstractPoisson,solver!,update!,HyprePoisson,putback!,L₂,L∞ -@inline insidep(a::AbstractArray,perdir=()) = CartesianIndices(ntuple( i-> i ∈ perdir ? (1:size(a,i)) : (2:size(a,i)-1), length(axes(a)))) +# @inline insidep(a::AbstractArray,perdir=()) = CartesianIndices(ntuple( i-> i ∈ perdir ? (1:size(a,i)) : (2:size(a,i)-1), length(axes(a)))) -@fastmath @inline function filldiag(I::CartesianIndex{d},x,L,perdir) where {d} - s = zero(eltype(L)) - for i in 1:d - Ip = I+δ(i,I) - (i ∈ perdir && Ip ∉ insidep(x,perdir)) && (Ip = WaterLily.CIj(i,Ip,1)) - s -= @inbounds(L[I,i]+L[Ip,i]) - end - return s -end +# @fastmath @inline function filldiag(I::CartesianIndex{d},x,L,perdir) where {d} +# s = zero(eltype(L)) +# for i in 1:d +# Ip = I+δ(i,I) +# (i ∈ perdir && Ip ∉ insidep(x,perdir)) && (Ip = WaterLily.CIj(i,Ip,1)) +# s -= @inbounds(L[I,i]+L[Ip,i]) +# end +# return s +# end struct HYPREPoisson{T,V<:AbstractVector{T},M<:AbstractArray{T}, Vf<:AbstractArray{T},Mf<:AbstractArray{T}, @@ -48,21 +48,18 @@ function HyprePoisson(x::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T ϵ[J[I]] = x[I]; r[J[I]] = z[I] A[J[I],J[I]] = one(T) # set all diagonals to one end - for I in insidep(x,perdir) - A[J[I],J[I]] = filldiag(I,x,L,perdir) + for I in inside(x) + A[J[I],J[I]] = WaterLily.diag(I,L) for i in 1:n # fill diagonal terms Im = I-δ(i,I); Ip = I+δ(i,I) - if i ∈ perdir - Im ∉ insidep(x,perdir) && (Im = WaterLily.CIj(i,Im,N[i])) - Ip ∉ insidep(x,perdir) && (Ip = WaterLily.CIj(i,Ip,1)) - end A[J[Im],J[I]] = L[I,i] # x[I-δ(i,I)]*L[I,i] A[J[Ip],J[I]] = L[Ip,i] # x[I+δ(i,I)]*L[I+δ(i,I),i] end end # if we have a purely Neumann problem, we must fix the pressure at one spot to satisfiability # the system, we choose the first block - length(perdir) == n && (A[1,1] = one(T); A[1,2:end] .= zero(T)) + # length(perdir) == n && (A[1,1] = one(T); A[1,2:end] .= zero(T); + # A[1,1] = one(T); A[1,2:end] .= zero(T)) HYPRE.Init() # Init and create a conjugate gradients solver solver = HYPRE.GMRES(;MaxIter,Tol,PrintLevel,Logging,Precond) return HYPREPoisson(ϵ,A,r,x,L,z,solver) From 566abbea0b34d39a67d6f6d95e721197033efdd7 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 9 Aug 2024 09:38:00 +0200 Subject: [PATCH 15/15] homogenise solvers --- ext/WaterLilyHYPREExt.jl | 2 +- src/MultiLevelPoisson.jl | 11 +++++------ 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/ext/WaterLilyHYPREExt.jl b/ext/WaterLilyHYPREExt.jl index d2af03f..a190c2c 100644 --- a/ext/WaterLilyHYPREExt.jl +++ b/ext/WaterLilyHYPREExt.jl @@ -81,7 +81,7 @@ update!(hp::HYPREPoisson) = nothing """ Solve the poisson problem and puts back the solution where it is expected. """ -function solver!(hp::HYPREPoisson;T=promote_type(Float64,eltype(hp.x)),kwargs...) +function solver!(hp::HYPREPoisson;kwargs...) fill!(hp) HYPRE.solve!(hp.solver, hp.ϵ,hp.A, hp.r) putback!(hp); diff --git a/src/MultiLevelPoisson.jl b/src/MultiLevelPoisson.jl index 8e629ce..b82a76c 100644 --- a/src/MultiLevelPoisson.jl +++ b/src/MultiLevelPoisson.jl @@ -87,17 +87,16 @@ L∞(p::MultiLevelPoisson) = L∞(p.levels[1]) mult!(ml::MultiLevelPoisson,x) = mult!(ml.levels[1],x) residual!(ml::MultiLevelPoisson,x) = residual!(ml.levels[1],x) -function solver!(ml::MultiLevelPoisson;tol=2e-4,itmx=32) +function solver!(ml::MultiLevelPoisson;tol=1e-4,itmx=32) p = ml.levels[1] - residual!(p); r₀ = r₂ = L∞(p); r₂₀ = L₂(p) + residual!(p); r₂ = L₂(p) nᵖ=0 - while r₂>tol && nᵖ5) && pop!(ml.levels); # remove coarsest level if this was easy - (nᵖ>4 && divisible(ml.levels[end])) && push!(ml.levels,restrictML(ml.levels[end])) # add a level if this was hard push!(ml.n,nᵖ); end