fix fgs for panel and seminfinite wake

Author: rymanderson

src/FLOWPanel.jl

@@ -38,6 +38,7 @@ import GeometricTools: Meshes
 import ImplicitAD as IAD
 import ImplicitAD: ForwardDiff as FD, ReverseDiff as RD
 import FLOWMath as math
+using WriteVTK
 
 # ------------ GLOBAL VARIABLES AND DATA STRUCTURES ----------------------------
 const module_path = splitdir(@__FILE__)[1]      # Path to this module

src/FLOWPanel_abstractliftingbody.jl

@@ -232,58 +232,33 @@ function _G_phi_wake!(self::AbstractLiftingBody{<:Any,<:Any,TF}, kernel, G, CPs,
     Das, Dbs = self.Das, self.Dbs
     derivatives_switch = FastMultipole.DerivativesSwitch(true,false,false)
 
-    # for chunk in chunks        # Distribute wake panel iteration among all CPU threads
-    # Threads.@threads for chunk in chunks        # Distribute wake panel iteration among all CPU threads
-    for chunk in chunks        # Distribute wake panel iteration among all CPU threads
-
-        # for (ei, (pi, nia, nib, pj, nja, njb)) in enumerate(eachcol(self.shedding))
-        for i_source in chunk                          # Iterate over wake-shedding panels
-
-            pi, nia, nib, pj, nja, njb = view(self.shedding, :, i_source)
-
-            # Fetch nodes of upper wake panel
-            nodes_idx = (self.cells[1, pi], self.cells[2, pi], self.cells[3, pi])
-
-            TE1 = nodes_idx[nia]
-            TE2 = nodes_idx[nib]
-            v1x = self.grid._nodes[1, TE1]
-            v1y = self.grid._nodes[2, TE1]
-            v1z = self.grid._nodes[3, TE1]
-            v2x = self.grid._nodes[1, TE2]
-            v2y = self.grid._nodes[2, TE2]
-            v2z = self.grid._nodes[3, TE2]
-
-            # direction of trailing semi-infinite wake
-            da1, da2, da3 = Das[1, i_source], Das[2, i_source], Das[3, i_source]
-            db1, db2, db3 = Dbs[1, i_source], Dbs[2, i_source], Dbs[3, i_source]
-            @assert isapprox(da1, db1) && isapprox(da2, db2) && isapprox(da3, db3) "Inconsistent wake directions in _G_phi_wake!"
-
-            for i_target in axes(CPs, 2)
-                # get target
-                tx, ty, tz = CPs[1, i_target], CPs[2, i_target], CPs[3, i_target]
-                target = FastMultipole.StaticArrays.SVector{3,TF}(tx, ty, tz)
-
-                # compute influence
-                phi, _ = induced_semiinfinite(target, kernel, v1x, v1y, v1z, v2x, v2y, v2z, da1, da2, da3, 1.0, derivatives_switch; kerneloffset)
-
-                # update G
-                G[i_target, pi] += phi
-            end
+    for isurf in eachindex(self.shedding)
+        # for chunk in chunks        # Distribute wake panel iteration among all CPU threads
+        # Threads.@threads for chunk in chunks        # Distribute wake panel iteration among all CPU threads
+        for chunk in chunks        # Distribute wake panel iteration among all CPU threads
+
+            # for (ei, (pi, nia, nib, pj, nja, njb)) in enumerate(eachcol(self.shedding))
+            for i_source in chunk                          # Iterate over wake-shedding panels
+
+                pi, nia, nib, pj, nja, njb = view(self.shedding[isurf], :, i_source)
 
-            # lower wake panel (if it exists)
-            if pj != -1
-                # Fetch nodes of lower wake panel
-                nodes_idx = (self.cells[1, pj], self.cells[2, pj], self.cells[3, pj])
+                # Fetch nodes of upper wake panel
+                nodes_idx = (self.cells[1, pi], self.cells[2, pi], self.cells[3, pi])
 
-                TE1 = nodes_idx[nja]
-                TE2 = nodes_idx[njb]
+                TE1 = nodes_idx[nia]
+                TE2 = nodes_idx[nib]
                 v1x = self.grid._nodes[1, TE1]
                 v1y = self.grid._nodes[2, TE1]
                 v1z = self.grid._nodes[3, TE1]
                 v2x = self.grid._nodes[1, TE2]
                 v2y = self.grid._nodes[2, TE2]
                 v2z = self.grid._nodes[3, TE2]
 
+                # direction of trailing semi-infinite wake
+                da1, da2, da3 = Das[isurf][1, i_source], Das[isurf][2, i_source], Das[isurf][3, i_source]
+                db1, db2, db3 = Dbs[isurf][1, i_source], Dbs[isurf][2, i_source], Dbs[isurf][3, i_source]
+                @assert isapprox(da1, db1) && isapprox(da2, db2) && isapprox(da3, db3) "Inconsistent wake directions in _G_phi_wake!"
+
                 for i_target in axes(CPs, 2)
                     # get target
                     tx, ty, tz = CPs[1, i_target], CPs[2, i_target], CPs[3, i_target]
@@ -293,9 +268,36 @@ function _G_phi_wake!(self::AbstractLiftingBody{<:Any,<:Any,TF}, kernel, G, CPs,
                     phi, _ = induced_semiinfinite(target, kernel, v1x, v1y, v1z, v2x, v2y, v2z, da1, da2, da3, 1.0, derivatives_switch; kerneloffset)
 
                     # update G
-                    G[i_target, pj] += phi
+                    G[i_target, pi] += phi
                 end
 
+                # lower wake panel (if it exists)
+                if pj != -1
+                    # Fetch nodes of lower wake panel
+                    nodes_idx = (self.cells[1, pj], self.cells[2, pj], self.cells[3, pj])
+
+                    TE1 = nodes_idx[nja]
+                    TE2 = nodes_idx[njb]
+                    v1x = self.grid._nodes[1, TE1]
+                    v1y = self.grid._nodes[2, TE1]
+                    v1z = self.grid._nodes[3, TE1]
+                    v2x = self.grid._nodes[1, TE2]
+                    v2y = self.grid._nodes[2, TE2]
+                    v2z = self.grid._nodes[3, TE2]
+
+                    for i_target in axes(CPs, 2)
+                        # get target
+                        tx, ty, tz = CPs[1, i_target], CPs[2, i_target], CPs[3, i_target]
+                        target = FastMultipole.StaticArrays.SVector{3,TF}(tx, ty, tz)
+
+                        # compute influence
+                        phi, _ = induced_semiinfinite(target, kernel, v1x, v1y, v1z, v2x, v2y, v2z, da1, da2, da3, 1.0, derivatives_switch; kerneloffset)
+
+                        # update G
+                        G[i_target, pj] += phi
+                    end
+
+                end
             end
         end
     end

src/FLOWPanel_liftingbody.jl

@@ -536,18 +536,21 @@ function solve2!(self::RigidWakeBody{Union{ConstantSource, ConstantDoublet}, 2,
     self.strength[:, 2] .= 0.0
 
     # calculate potential due to source strengths
-    _phi!(self, CPs, self.potential, backend; include_wake=false, optargs...) # NOTE: wake is included in the `direct!` function in this line
+    self.potential .= 0.0
+    _phi!(self, CPs, self.potential, backend; include_wake=false, optargs...)
 
     # zero source strengths while we solver for the doublet strengths
     self.strength[:, 1] .= 0.0
 
     # run fgs solver
     FastMultipole.solve!(self, solver.fgs; 
-        max_iterations=solver.max_iterations, 
+        max_iterations=solver.max_iterations,
+        inner_iterations=solver.inner_iterations,
         tolerance=solver.tolerance, 
         rlx=solver.rlx,
         derivatives_switches=FastMultipole.DerivativesSwitch(true, false, false, to_tuple(_unpack_fmm(self))),
-        reverse_pass=false
+        reverse_pass=solver.reverse_pass,
+        verbose=solver.verbose,
     )
 
     # update source strength
@@ -2063,7 +2066,6 @@ function FastMultipole.body_to_multipole!(system::RigidWakeBody{Union{ConstantSo
             # which vertices connect to the wake
             idx1 = Int(buffer[end-7, i_body])
             if idx1 > 0
-                println("wake...")
                 # first vertex and wake
                 Dax, Day, Daz = buffer[end-6, i_body], buffer[end-5, i_body], buffer[end-4, i_body]
                 Da = FastMultipole.SVector{3}(Dax, Day, Daz)
@@ -2090,7 +2092,6 @@ function FastMultipole.body_to_multipole!(system::RigidWakeBody{Union{ConstantSo
 
                 # doublet strength
                 strength_doublet = FastMultipole.SVector{1,TF}(strength[2])
-                @show Da, Db, v1, v2, v1w, v2w, normal
 
                 # update values
                 element = FastMultipole.Panel{FastMultipole.Dipole}
@@ -2253,7 +2254,7 @@ function additional_source_system_to_buffer!(buffer, i_buffer, system::RigidWake
     end
 
     if !system.semiinfinite_wake
-        buffer[4] += update_radius # conservative update to accomadate wake size
+        buffer[4, i_buffer] += update_radius # conservative update to accomadate wake size
     end
 
     return ilast+8
@@ -2585,41 +2586,41 @@ end
 
 
 ##### INTERNAL FUNCTIONS  ######################################################
-function _get_wakestrength_mu(self::RigidWakeBody, i; stri=1)
+function _get_wakestrength_mu(self::RigidWakeBody, i, isurf=1; stri=1)
     # if self.use_wake_strength
     #     return self.wake_strength[i], zero(self.wake_strength[i])
     # else
-        strength1 = self.strength[self.shedding[1, i], stri]
-        strength2 = self.shedding[4, i] != -1 ? self.strength[self.shedding[4, i], stri] : 0.0
+        strength1 = self.strength[self.shedding[isurf][1, i], stri]
+        strength2 = self.shedding[isurf][4, i] != -1 ? self.strength[self.shedding[isurf][4, i], stri] : 0.0
         return strength1, strength2
     # end
 end
-function _get_wakestrength_mu(self::RigidWakeBody{Union{ConstantSource,ConstantDoublet},2,<:Any}, i)
+function _get_wakestrength_mu(self::RigidWakeBody{Union{ConstantSource,ConstantDoublet},2,<:Any}, i, isurf=1)
     # if self.use_wake_strength
     #     return self.wake_strength[i], zero(self.wake_strength[i])
     # else
         stri = 2
-        strength1 = self.strength[self.shedding[1, i], stri]
-        strength2 = self.shedding[4, i] != -1 ? self.strength[self.shedding[4, i], stri] : 0.0
+        strength1 = self.strength[self.shedding[isurf][1, i], stri]
+        strength2 = self.shedding[isurf][4, i] != -1 ? self.strength[self.shedding[isurf][4, i], stri] : 0.0
         return strength1, strength2
     # end
 end
-function _get_wakestrength_Gamma(self::RigidWakeBody, i; stri=1)
+function _get_wakestrength_Gamma(self::RigidWakeBody, i, isurf=1; stri=1)
     # if self.use_wake_strength
     #     return self.wake_strength[i]
     # else
-        strength1 = self.strength[self.shedding[1, i], stri]
-        strength2 = self.shedding[4, i] != -1 ? self.strength[self.shedding[4, i], stri] : 0.0
+        strength1 = self.strength[self.shedding[isurf][1, i], stri]
+        strength2 = self.shedding[isurf][4, i] != -1 ? self.strength[self.shedding[isurf][4, i], stri] : 0.0
         return strength1 - strength2
     # end
 end
-function _get_wakestrength_Gamma(self::RigidWakeBody{Union{ConstantSource,ConstantDoublet},2,<:Any}, i)
+function _get_wakestrength_Gamma(self::RigidWakeBody{Union{ConstantSource,ConstantDoublet},2,<:Any}, i, isurf=1)
     # if self.use_wake_strength
     #     return self.wake_strength[i]
     # else
         stri = 2
-        strength1 = self.strength[self.shedding[1, i], stri]
-        strength2 = self.shedding[4, i] != -1 ? self.strength[self.shedding[4, i], stri] : 0.0
+        strength1 = self.strength[self.shedding[isurf][1, i], stri]
+        strength2 = self.shedding[isurf][4, i] != -1 ? self.strength[self.shedding[isurf][4, i], stri] : 0.0
         return strength1 - strength2
     # end
 end
@@ -2733,3 +2734,72 @@ function _savewake(self::RigidWakeBody, filename::String;
     end
 end
 #### END OF LIFTING BODY  ######################################################
+
+function write_vtk(name::String, body::RigidWakeBody; t=0.0, overwrite::Bool=false, semiinfinite_length=5.0)
+    # set wake length for visualization (if semi-infinite wake is enabled)
+    if body.semiinfinite_wake
+        for (Das, Dbs) in zip(body.Das, body.Dbs)
+            for i in 1:size(Das, 2)
+                # norm_Da = norm(view(Das, :, i))
+                # norm_Db = norm(view(Dbs, :, i))
+                Das[:, i] .*= semiinfinite_length # / norm_Da
+                Dbs[:, i] .*= semiinfinite_length # / norm_Db
+            end
+        end
+    end
+
+    WriteVTK.paraview_collection(name; append=!overwrite) do pvd
+        vtm = WriteVTK.vtk_multiblock(name*"_wake")
+        for i_surf in eachindex(body.shedding)
+            shedding = body.shedding[i_surf]
+            Das = body.Das[i_surf]
+            Dbs = body.Dbs[i_surf]
+            
+            n_wakes = size(shedding, 2)
+            points = zeros(typeof(body.grid._nodes[1]), 3, 4 * n_wakes)
+            cells = Vector{WriteVTK.MeshCell}(undef, n_wakes)
+            strengths = zeros(typeof(body.strength[1]), n_wakes)
+            
+            for i in 1:n_wakes
+                pi = shedding[1, i]
+                nia, nib = shedding[2, i], shedding[3, i]
+                
+                idx1 = body.cells[nia, pi]
+                idx2 = body.cells[nib, pi]
+                
+                p1 = body.grid._nodes[:, idx1]
+                p2 = body.grid._nodes[:, idx2]
+                
+                p3 = p2 + Dbs[:, i]
+                p4 = p1 + Das[:, i]
+                
+                points[:, 1 + 4*(i-1)] = p1
+                points[:, 2 + 4*(i-1)] = p2
+                points[:, 3 + 4*(i-1)] = p3
+                points[:, 4 + 4*(i-1)] = p4
+                
+                cells[i] = WriteVTK.MeshCell(WriteVTK.VTKCellTypes.VTK_QUAD, [1, 2, 3, 4] .+ 4*(i-1))
+
+                mu_upper, mu_lower = _get_wakestrength_mu(body, i)
+                strengths[i] = mu_upper - mu_lower
+            end
+            
+            WriteVTK.vtk_grid(vtm, name * "_wake_$i_surf", points, cells) do vtk
+                vtk["mu"] = strengths
+            end
+        end
+        pvd[t] = vtm
+    end
+
+    # restore original values of Da and Db if they were modified for visualization
+    if body.semiinfinite_wake
+        for (Das, Dbs) in zip(body.Das, body.Dbs)
+            for i in 1:size(Das, 2)
+                # norm_Da = norm(view(Das, :, i))
+                # norm_Db = norm(view(Dbs, :, i))
+                Das[:, i] ./= semiinfinite_length # / norm_Da
+                Dbs[:, i] ./= semiinfinite_length # / norm_Db
+            end
+        end
+    end
+end

src/FLOWPanel_simulate.jl

@@ -121,11 +121,11 @@ function simulate!(system::AbstractBody{TK,NK,TF}, wake::PanelWake, frames#=::Ab
 
         add_field(system, "U", "vector", collect(eachcol(system.velocity)), "cell")
         if has_grad_mu(system)
-            pnl.calcfield_Ugradmu(body; Gammai=get_Gammai(system))
-            add_fields(system, "Ugradmu", "U")
+            calcfield_Ugradmu(system; Gammai=get_Gammai(system))
+            addfields(system, "Ugradmu", "U")
         end
-        pnl.calcfield_Cp(body, norm(uinf))
-        pnl.calcfield_F(body, norm(uinf), rho)
+        calcfield_Cp(system, norm(uinf))
+        calcfield_F(system, norm(uinf), rho)
 
         #------- other solvers -------#
 
@@ -138,12 +138,12 @@ function simulate!(system::AbstractBody{TK,NK,TF}, wake::PanelWake, frames#=::Ab
 
         if !isnothing(path)
             # panel body
-            # write_vtk(joinpath(path, name * "_step_$i_step"), system; vtk_args...) # trailing_edge_list=.!shedding_surfaces, vtk_args...)
+            save(system, joinpath(path, name * "_vehicle_step_$i_step"))
 
             # particle field
 
             # visualize wake
-            write_vtk(name, wake, i_step, t)
+            write_vtk(joinpath(path, name), wake, i_step, t)
 
         end
 

src/FLOWPanel_solver.jl

@@ -283,19 +283,25 @@ struct FGSSolver{TFGS,TF<:Number} <: AbstractMatrixFreeSolver
     leaf_size::Int
     multipole_acceptance::Float64
     max_iterations::Int
+    inner_iterations::Int
     tolerance::Float64
     rlx::Float64
+    reverse_pass::Bool
+    verbose::Bool
 end
 
 function FGSSolver(body::AbstractBody; 
         max_iterations::Int=100,         # Maximum number of iterations
+        inner_iterations::Int=1,
         tolerance::Real=1e-6,            # Convergence tolerance
-        rlx::Real=1.0,                  # Relaxation factor
+        rlx::Real=1.0,                   # Relaxation factor
+        reverse_pass::Bool=true,         # Whether to perform reverse pass for adjoint sensitivities
         expansion_order=7,
         multipole_acceptance=0.4,
         leaf_size=10,
         shrink=false,
         recenter=false,
+        verbose=false,
     )
 
     # ensure CPoffset is negative (we'll solve this in the interior)
@@ -309,7 +315,7 @@ function FGSSolver(body::AbstractBody;
     # restore CPoffset
     body.CPoffset = CPoffset_old
 
-    return FGSSolver{typeof(fgs), TF}(fgs, Int(expansion_order), Int(leaf_size), Float64(multipole_acceptance), max_iterations, Float64(tolerance), Float64(rlx))
+    return FGSSolver{typeof(fgs), TF}(fgs, Int(expansion_order), Int(leaf_size), Float64(multipole_acceptance), max_iterations, Int(inner_iterations), Float64(tolerance), Float64(rlx), Bool(reverse_pass), Bool(verbose))
 end
 
 #--- test solve! ---#
@@ -326,7 +332,7 @@ function solve2!(self::AbstractBody{<:Any,1,<:Any}, Uinfs::Array{<:Real, 2}, sol
     self.velocity .= Uinfs
 
     # solve system
-    FastMultipole.solve!(self, solver.fgs; max_iterations=solver.max_iterations, tolerance=solver.tolerance, rlx=solver.rlx)
+    FastMultipole.solve!(self, solver.fgs; max_iterations=solver.max_iterations, inner_iterations=solver.inner_iterations, tolerance=solver.tolerance, rlx=solver.rlx, verbose=solver.verbose)
 
     # store solution
     set_solution(self, self.strength, self.strength, Tuple{Int,Float64}[], Uinfs)

src/FLOWPanel_wake.jl

@@ -44,9 +44,6 @@ function solve!(body::AbstractBody{TB}, wake::AbstractFreeWake, uinf::AbstractAr
     end
 
     # solve body
-    check_nans(body)
-    @show true in body_probes.gradient
-    @show true in body_probes.scalar_potential
     solve2!(body, body_probes.gradient, body_solver; backend)
 
     # body-on-all influence
@@ -110,7 +107,9 @@ function get_probes(body::AbstractBody{TK,NK,TF}) where {TK,NK,TF}
     normals = _calc_normals(body)
     CPs = _calc_controlpoints(body, normals; off=-1e-10)
     @assert !requires_hessian(body) "`get_probes` must be overloaded to support Hessian output for body type $(typeof(body))"
-    body_probes = FastMultipole.ProbeSystemArray(CPs, body.potential, body.velocity, hessian)
+    potential = zeros(TF, size(CPs, 2))
+    velocity = zeros(TF, 3, size(CPs, 2))
+    body_probes = FastMultipole.ProbeSystemArray(CPs, potential, velocity, hessian)
     return body_probes
 end