Postprocessing: write postprocess_statistics and postprocess_bladeloa…

…ding

src/FLOWUnsteady.jl

@@ -53,7 +53,7 @@ for header_name in ["vehicle", "vehicle_vlm",
                     "maneuver", "rotor",
                     "simulation_types", "simulation", "utils",
                     "processing", "processing_force", "monitors",
-                    "noise_wopwop", "noise_bpm"]
+                    "noise_wopwop", "noise_bpm", "postprocessing"]
 
     include("FLOWUnsteady_"*header_name*".jl")
 

src/FLOWUnsteady_postprocessing.jl

@@ -0,0 +1,169 @@
+#=##############################################################################
+# DESCRIPTION
+    Tools for postprocessing simulations.
+
+# ABOUT
+  * Created   : Mar 2023
+  * License   : MIT
+=###############################################################################
+
+
+"""
+    postprocess_statistics(read_path, save_path, nums;
+                            # PROCESSING OPTIONS
+                            idens           = [""],     # Use this to agglomerate multiple simulations
+                            to_exclude      = [],       # Exclude file names containing these words
+                            cyl_axial_dir   = nothing,  # Calculate cylindrical statistics if given an axial axis (vector)
+                            # OUTPUT OPTIONS
+                            prompt=true, debug=false,
+                            verbose=true, v_lvl=0)
+
+Calculate statistics of the simulation VTK outputs over a given range of time
+steps given by `nums`. Use this to calculate, for instance, the mean load
+distribution and fluctuations on wings and rotors.
+"""
+function postprocess_statistics(read_path, save_path, nums;
+                                    # PROCESSING OPTIONS
+                                    idens           = [""],     # Use this to agglomerate multiple simulations
+                                    to_exclude      = [],       # Exclude file names containing these words
+                                    cyl_axial_dir::Union{Nothing, Vector{Float64}} = nothing,  # Calculate cylindrical statistics if given an axial axis (vector)
+                                    # OUTPUT OPTIONS
+                                    prompt=true, debug=false,
+                                    verbose=true, v_lvl=0)
+
+
+    read_paths      = [read_path*iden for iden in idens]  # Path to simulations to process
+
+
+    # -------------- PROCESSING OPTIONS --------------------------------------------
+    # Cartesian statistical operations to perform
+    operations_cart = [
+                        ("mean", gt._agglomerate_mean, 1),
+                        ("cart-variance", gt._agglomerate_variance, 2),
+                        # ("cart-corr1", gt._agglomerate_correlation1, 2),
+                        # ("cart-corr2", gt._agglomerate_correlation2, 2),
+                        # ("cart-corr3", gt._agglomerate_correlation3, 2),
+                      ]
+
+    # Cylindrical statistical operations to perform
+    if cyl_axial_dir != nothing
+        operations_cyl = [
+                          ("cyl-mean", (args...; optargs...) -> gt._agglomerate_cyl_mean(cyl_axial_dir, args...; optargs...), 1),
+                          ("cyl-variance", (args...; optargs...) -> gt._agglomerate_cyl_variance(cyl_axial_dir, args...; optargs...), 2),
+                          # ("cyl-corr1", (args...; optargs...) -> gt._agglomerate_cyl_correlation1(axial_dir, args...; optargs...), 2),
+                          # ("cyl-corr2", (args...; optargs...) -> gt._agglomerate_cyl_correlation2(axial_dir, args...; optargs...), 2),
+                          # ("cyl-corr3", (args...; optargs...) -> gt._agglomerate_cyl_correlation3(axial_dir, args...; optargs...), 2),
+                        ]
+    else
+        operations_cyl = []
+    end
+
+    # Statistical operations
+    operations = vcat(operations_cart, operations_cyl)
+
+    # --------------  PROCESSING SETUP ---------------------------------------------
+    if verbose; println("\t"^(v_lvl)*"Getting ready to process statistics on $(read_paths)"); end;
+
+    # Create save path
+    if !(save_path in read_paths)
+        gt.create_path(save_path, prompt)
+    end
+
+    # Identify vtks to average together
+    to_average_prefixes = unique([f[1:findfirst(".", f)[1]-1] for f in readdir(read_paths[1])
+                                        if contains(f, ".vtk") && prod(.!contains.(f, to_exclude))])
+
+
+    # -------------- PROCESS SIMULATION --------------------------------------------
+    if verbose; println("\n"*"\t"^(v_lvl)*"Iterating over VTK files"); end;
+
+    # Threads.@threads for pref in to_average_prefixes
+    for pref in to_average_prefixes
+
+        if verbose
+            println("\t"^(v_lvl+1)*"Averaging files with prefix $(pref)")
+        end
+
+        # this_operations = contains(pref, "Line-Axial") ? operations_cart : operations
+        this_operations = operations
+
+        # files = vcat([ [joinpath(read_path, f) for f in readdir(read_path) if contains(f, ".vtk") && contains(f, pref)] for read_path in read_paths]...)
+        files = vcat([ [joinpath(read_path, f) for f in readdir(read_path)
+                                    if contains(f, ".vtk") && contains(f, pref) && parse(Int, f[findfirst(".", f)[1]+1:findlast(".", f)[1]-1]) in nums]
+                                        for read_path in read_paths]...)
+
+        t = @elapsed gt.calculate_statistics_vtk(files, save_path;
+                                                    operations=this_operations,
+                                                    # read_path=read_path,
+                                                    read_path="",
+                                                    out_filename=pref*"-statistics",
+                                                    verbose=verbose, v_lvl=v_lvl+2)
+
+        if verbose
+            println("\t"^(v_lvl+2)*"Process statistics:\t$(round(t, digits=1)) s")
+        end
+    end
+
+    println("Statistics saved under $save_path")
+end
+
+
+"""
+    postprocess_bladeloading(read_path;
+                                O           = zeros(3),     # Rotor center
+                                rotor_axis  = [-1, 0, 0],   # Rotor centerline axis
+                                Ftot_axis   = nothing,      # Use a different centerline axis for forces if given
+                                filename    = "singlerotor_Rotor_Blade1_vlm-statistics.vtk", # File name
+                                fieldsuff   = "-mean"       # Suffix of fields "Gamma" and "Ftot", if any
+                                )
+
+Read a blade VTK file `filename` under directory `read_path` and returns
+the circulation and force components of the load distribution along the blade.
+
+Return: `rs, Gamma, Np, Tp, Rp, Zhat, Rhat, That, Ftot`
+"""
+function postprocess_bladeloading(read_path;
+                                    O           = zeros(3),     # Rotor center
+                                    rotor_axis  = [-1, 0, 0],   # Rotor centerline axis
+                                    Ftot_axis   = nothing,      # Use a different centerline axis for forces if given
+                                    filename    = "singlerotor_Rotor_Blade1_vlm-statistics.vtk", # File name
+                                    fieldsuff   = "-mean"       # Suffix of fields "Gamma" and "Ftot", if any
+                                    )
+
+    @assert isapprox(norm(rotor_axis), 1.0) "Non-unitary rotor axis $(rotor_axis)"
+
+    _Ftot_axis = Ftot_axis==nothing ? rotor_axis : Ftot_axis
+
+    points, cells, cell_types, data = gt.read_vtk(filename; path=read_path)
+
+    maxind = Int(size(cells, 1)/2)
+    Xs = [mean([points[:, pi+1] for pi in cell]) for cell in cells[1:maxind]]
+    Rs = [(X-O) - rotor_axis*dot(X-O, rotor_axis) for X in Xs]
+
+    # Blade centerline
+    Rhat = mean(Rs)
+    Rhat /= norm(Rhat)
+
+    rs = [dot(R, Rhat) for R in Rs]
+
+    # Tangent vector
+    That = cross(_Ftot_axis, Rhat)
+
+    # Grabs only the rectangular cells
+    Gamma = data["CELL_DATA"]["Gamma"*fieldsuff][1:maxind]
+    Ftot = data["CELL_DATA"]["Ftot"*fieldsuff][:, 1:maxind]
+
+    nr = length(rs)
+    Np = zeros(nr)              # Normal component
+    Rp = zeros(nr)              # Radial component
+    Tp = zeros(nr)              # Tangential component
+
+    for i in 1:nr
+        F = Ftot[:, i]
+        Np[i] = dot(F, _Ftot_axis)
+        Rp[i] = dot(F, Rhat)
+        Tp[i] = dot(F, That)
+    end
+
+    return rs, Gamma, Np, Tp, Rp, _Ftot_axis, Rhat, That, Ftot
+end