fix syntax and clean up

Author: rymanderson

examples/duct.jl

@@ -104,21 +104,23 @@ points = hcat(xs, ys)
 # Generate body of revolution
 body = pnl.generate_revolution_liftbody(bodytype, points, NDIVS_theta;
                                         bodyoptargs = (
-                                                        CPoffset=1e-14,
+                                                        CPoffset=1e-10,
                                                         kerneloffset=1e-8,
                                                         kernelcutoff=1e-14,
                                                         characteristiclength=(args...)->d*aspectratio,
-                                                        semiinfinite_wake=true
+                                                        semiinfinite_wake=false
                                             )
                                         )
 
-# body.shedding = body.shedding[:,1:1]
+# body.shedding[1] = body.shedding[1][:,1:1]
 # body.nsheddings = 1
 # body.shedding_full .= -1
 # idx_shed = body.shedding[1, 1]
-# body.shedding_full[:, idx_shed] .= view(body.shedding, 2:3, 1)
-# idx_shed = body.shedding[4, 1]
-# body.shedding_full[:, idx_shed] .= view(body.shedding, 5:6, 1)
+# body.shedding_full[1:2, idx_shed] .= view(body.shedding[1], 2:3, 1)
+# body.shedding_full[3:4, idx_shed] .= 1
+# idx_shed = body.shedding[1][4, 1]
+# body.shedding_full[1:2, idx_shed] .= view(body.shedding[1], 5:6, 1)
+# body.shedding_full[3:4, idx_shed] .= 1
 
 println("Number of panels:\t$(body.ncells)")
 
@@ -141,19 +143,19 @@ AOA = AOAs[i]
 
     # Unitary direction of semi-infinite vortex at points `a` and `b` of each
     # trailing edge panel
-    body.Das[1] .= repeat(Vinf/magVinf, 1, body.nsheddings)
-    body.Dbs[1] .= repeat(Vinf/magVinf, 1, body.nsheddings)
+    body.Das[1] .= repeat(Vinf/magVinf, 1, size(body.Das[1], 2))
+    body.Dbs[1] .= repeat(Vinf/magVinf, 1, size(body.Dbs[1], 2))
 
     # Solve body (panel strengths) giving `Uinfs` as boundary conditions and
     # `Das` and `Dbs` as trailing edge rigid wake direction
     # @time pnl.solve(body, Uinfs, Das, Dbs)
-    leaf_size = 50
+    leaf_size = 100
     expansion_order = 10
     multipole_acceptance = 0.4
     backend = pnl.FastMultipoleBackend(;
                                     expansion_order,
                                     multipole_acceptance,
-                                    leaf_size
+                                    leaf_size=20
                                 )
     # backend = pnl.DirectBackend()
     # global solver = pnl.Backslash(body; least_squares=true)
@@ -169,65 +171,62 @@ AOA = AOAs[i]
     #                 leaf_size=10
     #             )
     # )
-    println("Initializaing solver...")
-    @time solver = pnl.FGSSolver(body;
-        max_iterations=500,         # Maximum number of iterations
-        tolerance=1.0e-6,            # Convergence tolerance
-        rlx=1.0,                  # Relaxation factor
-        expansion_order,
-        multipole_acceptance,
-        leaf_size,
-        shrink=true,
-        recenter=false,
-    )
-    # solver = pnl.BackslashDirichlet(body)
-
-    println("\nSolving...")
-
-    # profile
-    # using Profile, PProf
-
-    pnl.solve2!(body, Uinfs, solver; backend)
-    # @profile pnl.solve2!(body, Uinfs, solver; backend)
-    # Profile.clear()
-    # @profile pnl.solve2!(body, Uinfs, solver; backend)
-    # pprof()
+    # function test_solver(inner_iterations, reverse_pass)
+    #     println("Initializing solver with inner_iterations=$inner_iterations, reverse_pass=$reverse_pass...")
+        println("Initializaing solver...")
+        @time solver = pnl.FGSSolver(body;
+            max_iterations=500,         # Maximum number of iterations
+            tolerance=1.0e-6,            # Convergence tolerance
+            rlx=1.0,                  # Relaxation factor
+            expansion_order,
+            multipole_acceptance,
+            leaf_size,
+            shrink=true,
+            recenter=false,
+            inner_iterations=20,
+            reverse_pass=false,
+            verbose=false
+        )
+        # solver = pnl.BackslashDirichlet(body)
+
+        println("\nSolving...")
+
+        # profile
+        # using Profile, PProf
+
+        @time pnl.solve2!(body, Uinfs, solver; backend)
+        # @profile pnl.solve2!(body, Uinfs, solver; backend)
+        # Profile.clear()
+        # @profile pnl.solve2!(body, Uinfs, solver; backend)
+        # pprof()
+
+    # end
+
+    # for reverse_pass in [true, false]
+    #     for inner_iterations in [1, 2, 4, 8, 16, 32]
+    #         test_solver(inner_iterations, reverse_pass)
+    #     end
+    # end
 
     # ----------------- POST PROCESSING ----------------------------------------
     println("\nPost processing...")
 
     # Calculate surface velocity U on the body
-    @time Us = pnl.calcfield_U(body, body; backend)
+    @time Us = pnl.calcfield_U!(body, body; backend)
+    pnl.apply_freestream!(body, Uinfs[:,1])
 
-    # NOTE: Since the boundary integral equation of the potential flow has a
-    #       discontinuity at the boundary, we need to add the gradient of the
-    #       doublet strength to get an accurate surface velocity
-
-    # Calculate surface velocity U_∇μ due to the gradient of the doublet strength
-    Gammai = kernel == pnl.VortexRing ? 1 : kernel == Union{pnl.ConstantSource, pnl.ConstantDoublet} ? 2 : 0
-    UDeltaGamma = pnl.calcfield_Ugradmu(body; Gammai)
-    # UDeltaGamma = pnl.calcfield_Ugradmu(body; sharpTE=true, force_cellTE=false)
-
-    # Add both velocities together
-    pnl.addfields(body, "Ugradmu", "U")
-
-    # Calculate pressure coefficient (based on U + U_∇μ)
-    @time Cps = pnl.calcfield_Cp(body, magVinf)
+    # Calculate pressure coefficient (based on body.velocity)
+    @time Cps = pnl.calcfield_Cp!(body.Cp, body, magVinf; correct_kuttacondition=true)
 
     # Calculate the force of each panel (based on Cp)
-    @time Fs = pnl.calcfield_F(body, magVinf, rho)
+    @time Fs = pnl.calcfield_F!(body, magVinf, rho)
 
     # check normal flow condition
-    Us_tot = pnl.get_field(body, "U")["field_data"] # total velocity
-    Us_tot = [Us_tot[j][i] for i=1:3, j=1:size(Us_tot,1)]
-    normals = pnl._calc_normals(body)
+    Us_tot = body.velocity
     Udotn = sum(Us_tot .* normals, dims=1)
     resid = maximum(abs.(Udotn))
     println("Max flow tangency residual: $resid")
 
-    normals = pnl.calc_normals(body)
-    CPs_inside = pnl._calc_controlpoints(body, normals; off=-1e-10)
-
     # ----------------- COMPARISON TO EXPERIMENTAL DATA ------------------------
     # Plot surface pressure along slices of the duct
     fig, axs = plot_Cp(body, AOA)
@@ -250,8 +249,7 @@ AOA = AOAs[i]
         name *= kernel == pnl.VortexRing ? "_vortexring" :
                 kernel == Union{pnl.ConstantSource, pnl.ConstantDoublet} ? "_source_doublet" :
                 ""
-        global vtks *= pnl.save(body, name; path=save_path, num=i,
-                                        wake_panel=false, debug=false, out_wake=false)
+        global vtks *= pnl.write_vtk(name, body, i, 0.0; overwrite=true)
     end
 
 # end

examples/duct_postprocessing.jl

@@ -26,10 +26,7 @@ function plot_Cp(body, AOA; plot_experimental=true)
 
         # Data filter to get only +z or only -z data points
         slicefilter(x, i) = (x[3] >= 0) == upperside
-
-        slicepoints, sliceCps = pnl.slicefield(body, "Cp",
-                                                position, direction, row;
-                                                filter=slicefilter)
+        slicepoints, sliceCps = pnl.slice_scalarfield(body, :Cp, 2, 0.0, 0.1; filter=slicefilter)
 
         side = upperside ? "upper" : "lower"
 

examples/duct_unsteady.jl

@@ -112,14 +112,6 @@ body = pnl.generate_revolution_liftbody(bodytype, points, NDIVS_theta;
                                             )
                                         )
 
-# body.shedding = body.shedding[:,1:1]
-# body.nsheddings = 1
-# body.shedding_full .= -1
-# idx_shed = body.shedding[1, 1]
-# body.shedding_full[:, idx_shed] .= view(body.shedding, 2:3, 1)
-# idx_shed = body.shedding[4, 1]
-# body.shedding_full[:, idx_shed] .= view(body.shedding, 5:6, 1)
-
 println("Number of panels:\t$(body.ncells)")
 
 vtks = save_path*"/"                        # String with VTK output files
@@ -137,8 +129,20 @@ AOA = AOAs[i]
     # Freestream at every control point
     Uinf(t) = Vinf
 
+    # prepare other inputs
+    eta = 0.3
+    frames = nothing
+    maneuver = (args...; optargs...) -> nothing
+    l = d * aspectratio
+    dt = magVinf / l / (n_rfl * 500)
+    t_range = range(0.0, step=dt, length=201)
+
+    # update wake directions
+    body.Das[1] .= repeat(Vinf*dt*eta, 1, size(body.Das, 2))
+    body.Dbs[1] .= repeat(Vinf*dt*eta, 1, size(body.Dbs, 2))
+
     # select backend for N-body interactions
-    leaf_size = 150
+    leaf_size = 20
     expansion_order = 10
     multipole_acceptance = 0.4
     backend = pnl.FastMultipoleBackend(;
@@ -163,27 +167,21 @@ AOA = AOAs[i]
     # )
     println("Initializaing solver...")
     @time body_solver = pnl.FGSSolver(body;
-        max_iterations=100,         # Maximum number of iterations
-        tolerance=1.0e-7,            # Convergence tolerance
-        rlx=0.1,                  # Relaxation factor
+        max_iterations=50,         # Maximum number of iterations
+        inner_iterations=10,       # Maximum number of inner iterations
+        reverse_pass=true,        # Whether to do reverse sweeps or not
+        tolerance=1.0e-6,            # Convergence tolerance
+        rlx=1.0,                  # Relaxation factor
         expansion_order,
         multipole_acceptance,
-        leaf_size,
+        leaf_size=150,
         shrink=true,
         recenter=false,
     )
     # solver = pnl.BackslashDirichlet(body)
 
     # initialize wake
-    wake = pnl.PanelWake(body; nwakerows=10)
-
-    # prepare other inputs
-    eta = 0.3
-    frames = nothing
-    maneuver = (args...; optargs...) -> nothing
-    l = d * aspectratio
-    dt = magVinf / l / (n_rfl * 5)
-    t_range = range(0.0, step=dt, length=15)
+    wake = pnl.PanelWake(body; nwakerows=201)
 
     println("\nBegin simulation...")
     @time begin

examples/single_triangle_body.jl

@@ -5,6 +5,8 @@ using FLOWPanel
 const pnl = FLOWPanel
 using FLOWPanel.gt.Meshes
 import FLOWPanel.PyPlot as plt
+colormap = plt.get_cmap("RdBu",15)
+
 
 function make_single_triangle_body(; semiinfinite_wake=false)
     pnl = FLOWPanel
@@ -31,15 +33,15 @@ function make_single_triangle_body(; semiinfinite_wake=false)
     # body = pnl.NonLiftingBody{pnl.ConstantDoublet}(grid)
     # body = pnl.NonLiftingBody{Union{pnl.ConstantSource, pnl.ConstantDoublet}}(grid)
 
-    shedding = zeros(Int, 6, 1)
-    shedding[1,1] = 1
-    shedding[2,1] = 2
-    shedding[3,1] = 3
-    shedding[4,1] = -1
+    shedding = zeros(Int, 6, 0)
+    # shedding[1,1] = 1
+    # shedding[2,1] = 2
+    # shedding[3,1] = 3
+    # shedding[4,1] = -1
     body = pnl.RigidWakeBody{Union{pnl.ConstantSource, pnl.ConstantDoublet}}(grid, [shedding];
         semiinfinite_wake)
-    body.Das[1][2,:] .= -1
-    body.Dbs[1][2,:] .= -1
+    # body.Das[1][2,:] .= -1
+    # body.Dbs[1][2,:] .= -1
 
     return body
 end
@@ -52,13 +54,13 @@ str = 1.0
 # body.strength[1,1] = str
 body.strength[1,2] = str
 
-normals = pnl._calc_normals(body)
+normals = pnl.calc_normals!(body)
 @show normals # outward facing
-cps = pnl._calc_controlpoints(body, normals; off=1e-14)
+cps = pnl.calc_controlpoints!(body, normals; off=1e-14)
 
 # compute velocity at a line of points above the triangle
 npoints = 1000
-dz = range(-1.0, stop=1.0, length=npoints) .* 5.0
+dz = range(-1.0, stop=1.0, length=npoints) .* 1.0
 points = hcat([ cps[:,1] .+ [0.0,0.0,z] for z in dz ]...)  # points above centroid of triangle
 # points = hcat([ cps[:,1] .+ z .* normals[:,1] for z in dz ]...)  # points above centroid of triangle
 # points = hcat([[0.33; 0.33; 0.0] .+ [z; z; z] for z in dz ]...)  # points in y from the centroid
@@ -127,6 +129,21 @@ axs[1].plot(dz, phi_semiinfinite2, ":", label="semi-infinite panel2")
 axs[1].legend()
 plt.tight_layout()
 
+# plot results with flipped axes
+fig3 = plt.figure("verify panel flipped")
+fig3.clf()
+fig3.add_subplot(121, xlabel=L"\phi", ylabel="z")
+fig3.add_subplot(122, xlabel=L"\vec{u} \cdot \hat{n}", ylabel="z")
+axs3 = fig3.get_axes()
+axs3[1].plot(phi_direct, dz, label="direct", color=colormap(14))
+axs3[2].plot(out_direct_z, dz, label="direct", color=colormap(14))
+for ax in axs3
+    ax.spines["top"].set_visible(false)
+    ax.spines["right"].set_visible(false)
+end
+
+plt.tight_layout()
+
 fig2 = plt.figure("error")
 fig2.clf()
 axs2 = fig2.add_subplot(111, xlabel="z", ylabel="error")

src/FLOWPanel.jl

@@ -22,7 +22,7 @@ export  solve, save, Uind!, phi!,
 # ------------ GENERIC MODULES -------------------------------------------------
 import Dierckx
 import LinearAlgebra as LA
-import LinearAlgebra: I
+import LinearAlgebra: I, lu!, ldiv!
 import Krylov
 import LinearOperators
 import SimpleNonlinearSolve
@@ -67,12 +67,15 @@ const Im = Array(1.0I, 3, 3)
 # Shedding matrix for a RigidWakeBody without shedding
 const noshedding = zeros(Int, 6, 0)
 
+const SEMIINFINITE_LENGTH = Array{Float64,0}(undef)
+SEMIINFINITE_LENGTH[] = 10.0
+
 # ------------ HEADERS ---------------------------------------------------------
 for header_name in ["elements", "linearsolver", "fmm",
                     "abstractbody", "solver",
                     "nonliftingbody",
                     "abstractliftingbody", "liftingbody",
-                    "multibody",  "elements_fmm",
+                    "multibody",  "elements_fmm", "frames",
                     "liftingline",
                     "utils", "postprocess",
                     "wake", "simulate",