JuliaIBPM · jdeldre · Dec 2, 2024 · Oct 9, 2024 · Nov 18, 2024 · Nov 18, 2024
diff --git a/Project.toml b/Project.toml
@@ -8,6 +8,7 @@ ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -17,6 +18,7 @@ ColorTypes = "<0.10.1, 0.10, 0.11"
 Interpolations = "0.12.10, 0.13, 0.14, 0.15"
 LinearAlgebra = "1.7"
 OrdinaryDiffEq = "6.50 - 6.87"
+Plots = "1.40.9"
 RecipesBase = "1.0"
 RecursiveArrayTools = "2.2 - 2.38.5, 2.38.6, 3"
 Statistics = "1.7"

diff --git a/README.md b/README.md
@@ -11,5 +11,6 @@ This package contains tools that can be used to post-process and analyze the sol
 * proper orthogonal decomposition (POD)
 * dynamic mode decomposition (DMD)
 * finite-time Lyapunov exponent (FTLE)
+* Lagrangian-averaged vorticity deviation (LAVD)
 
 The examples in the documentation, which have companion notebooks in the `examples` folder, are the best way to get started.   
diff --git a/docs/make.jl b/docs/make.jl
@@ -14,6 +14,7 @@ makedocs(
                      "manual/dmdtest.md",
                      "manual/ftle_continuous.md",
                      "manual/ftle.md",
+                     "manual/lavd.md",
                      "manual/functions.md"
                      ]
         #"Internals" => [ "internals/properties.md"]

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -8,6 +8,7 @@ are
 * Proper orthogonal decomposition (POD)
 * Dynamic mode decomposition (DMD)
 * Finite-time Lyapunov exponent (FTLE)
+* Lagrangian-averaged vorticity deviation (LAVD)
 
 ## Installation
 

diff --git a/docs/src/manual/dmdtest.md b/docs/src/manual/dmdtest.md
@@ -121,11 +121,6 @@ vals, idex = findmin(abs2.(true_evals .- transpose(dmdmodes.evals)),dims=2)
 err = sqrt(sum(vals))
 ````
 
-## DMD functions
-```@docs
-dmd
-```
-
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*

diff --git a/docs/src/manual/ftle.md b/docs/src/manual/ftle.md
@@ -87,7 +87,7 @@ X_MAX = 2.0
 Y_MIN = -2.0
 Y_MAX = 2.0
 dx = 0.01
-ftlegrid = PhysicalGrid((X_MIN,X_MAX),(Y_MIN,Y_MAX),dx,optimize=false)
+ftlegrid = PhysicalGrid((X_MIN,X_MAX),(Y_MIN,Y_MAX),dx)
 ftle_cache = SurfaceScalarCache(ftlegrid)
 x0, y0 = x_grid(ftle_cache), y_grid(ftle_cache)
 ````
@@ -148,7 +148,7 @@ xp_max = 0.0
 yp_min = 0.5
 yp_max = 1.5
 dxp = 0.1
-p_grid = PhysicalGrid((xp_min,xp_max),(yp_min,yp_max),dxp,optimize=false)
+p_grid = PhysicalGrid((xp_min,xp_max),(yp_min,yp_max),dxp)
 p_cache = SurfaceScalarCache(p_grid);
 xp0, yp0 = x_grid(p_cache), y_grid(p_cache);
 
@@ -201,17 +201,6 @@ The code here creates a gif
 
     end every 1 fps = 2
 
-## FTLE functions
-```@docs
-make_interp_fields!
-gen_init_conds
-euler_forward
-euler_backward
-adams_bashforth_2_forward
-adams_bashforth_2_backward
-compute_FTLE!
-```
-
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*

diff --git a/docs/src/manual/ftle_continuous.md b/docs/src/manual/ftle_continuous.md
@@ -25,7 +25,7 @@ X_MIN, X_MAX = 0.0, 2.0
 Y_MIN, Y_MAX = 0.0, 1.0
 dx = 0.002
 
-ftlegrid = PhysicalGrid((X_MIN,X_MAX),(Y_MIN,Y_MAX),dx,optimize=false);
+ftlegrid = PhysicalGrid((X_MIN,X_MAX),(Y_MIN,Y_MAX),dx);
 ftle_cache = SurfaceScalarCache(ftlegrid)
 x0, y0 = x_grid(ftle_cache), y_grid(ftle_cache)
 ````
@@ -102,7 +102,7 @@ X_MIN, X_MAX = -6.0, 50.0
 Y_MIN, Y_MAX = -2, 2
 dx = 0.08
 
-ftlegrid = PhysicalGrid((X_MIN,X_MAX),(Y_MIN,Y_MAX),dx,optimize=false);
+ftlegrid = PhysicalGrid((X_MIN,X_MAX),(Y_MIN,Y_MAX),dx);
 ftle_cache = SurfaceScalarCache(ftlegrid)
 x0, y0 = x_grid(ftle_cache), y_grid(ftle_cache)
 ````

diff --git a/docs/src/manual/lavd.md b/docs/src/manual/lavd.md
@@ -0,0 +1,232 @@
+```@meta
+EditURL = "../../../test/literate/lavd.jl"
+```
+
+# Lagrangian-averaged vorticity deviation (LAVD)
+In this example, we will compute the Lagrangian-averaged vorticity deviation (LAVD) field for a flat plate undergoing a 45-degree pitch-up maneuver. The original simluation of the flat plate was used by Wang and Eldredge 2012 (https://doi.org/10.1007/s00162-012-0279-5). The results are compared with Huang and Green 2016 (https://arc.aiaa.org/doi/10.2514/6.2016-2082). The theory behind LAVD can be found in Haller et al. 2016 (https://doi.org/10.1017/jfm.2016.151). A MATLAB package for computing LAVD and extracting coherent vortices can be found here (https://github.com/Hadjighasem/Lagrangian-Averaged-Vorticity-Deviation-LAVD).
+
+```@meta
+CurrentModule = ILMPostProcessing
+```
+
+````@example lavd
+using ILMPostProcessing
+using ViscousFlow
+using Plots
+using Statistics
+````
+
+# Viscous Flow of Pitching Flat Plate
+## Problem Specification
+For faster compututation and testing purposes, the Reynolds number is set to 100 as opposed to 1000 in Huang's paper. The grid Re is also set to 4.0. If better resolution is desired, try grid Re = 3.0. The domain of interest is from x = -0.5 to x = 5.5, but it is set from x = -3.0 to x = 5.5 since the velocity and vorticity fields ahead of the flat plate are required to compute LAVD.
+
+````@example lavd
+my_params = Dict()
+my_params["Re"] = 100
+my_params["grid Re"] = 4.0
+my_params["freestream speed"] = 1.0
+my_params["freestream angle"] = 0.0
+
+xlim = (-3.0,5.5)
+ylim = (-2.0,1.0)
+g = setup_grid(xlim,ylim,my_params)
+Δs = surface_point_spacing(g,my_params)
+````
+
+## Set up Body
+Create a rectangle of length 1.0 and thickness 0.023.
+
+````@example lavd
+Lp = 1.0
+body = Rectangle(Lp/2,0.023/2,Δs)
+bl = BodyList([body])
+````
+
+## Set the Body Motion
+Create a smooth position ramp of 45 degreess of the flat plate's angle of attack about its leading edge.
+
+````@example lavd
+vel = 45pi/180  ## nominal ramp velocity
+Δx = -45pi/180 ## change in position
+t0 = 1.0 ## time of ramp start
+k = SmoothRampDOF(vel,Δx,t0)
+````
+
+Plot the ramp.
+
+````@example lavd
+t = range(0,3,length=301)
+plot(t,dof_position.(k.(t)),xlims=(0,Inf),label="x")
+````
+
+Create the joint.
+
+````@example lavd
+parent_body, child_body = 0, 1
+Xp = MotionTransform([0,0],0) # location of joint in inertial system
+xpiv = [-0.5,0] # place center of motion at LE
+Xc = MotionTransform(xpiv,0)
+joint1 = Joint(RevoluteJoint,parent_body,Xp,child_body,Xc,[k])
+m = RigidBodyMotion([joint1],bl)
+
+x = init_motion_state(bl,m)
+update_body!(bl,x,m)
+plot(bl,xlim=xlim,ylim=ylim)
+````
+
+Animate the motion
+
+````@example lavd
+macro animate_motion(b,m,dt,tmax,xlim,ylim)
+    return esc(quote
+            bc = deepcopy($b)
+            t0, x0 = 0.0, init_motion_state(bc,$m)
+            dxdt = zero(x0)
+            x = copy(x0)
+
+            @gif for t in t0:$dt:t0+$tmax
+                motion_rhs!(dxdt,x,($m,bc),t)
+                global x += dxdt*$dt
+                update_body!(bc,x,$m)
+                plot(bc,xlim=$xlim,ylim=$ylim)
+            end every 5
+        end)
+end
+
+@animate_motion bl m 0.01 4 (-0.5, 5.5) ylim
+````
+
+## Define the Boundary Condition Functions
+
+````@example lavd
+function my_vsplus(t,x,base_cache,phys_params,motions)
+  vsplus = zeros_surface(base_cache)
+  surface_velocity!(vsplus,x,base_cache,motions,t)
+  return vsplus
+end
+
+function my_vsminus(t,x,base_cache,phys_params,motions)
+  vsminus = zeros_surface(base_cache)
+  surface_velocity!(vsminus,x,base_cache,motions,t)
+  return vsminus
+end
+
+bcdict = Dict("exterior" => my_vsplus, "interior" => my_vsminus)
+````
+
+## Construct the system structure
+
+````@example lavd
+sys = viscousflow_system(g,bl,phys_params=my_params,motions=m,bc=bcdict);
+u0 = init_sol(sys)
+````
+
+Check the effective Reynolds number.
+
+````@example lavd
+Umax, imax, tmax, bmax = maxvelocity(u0,sys)
+L = Lp
+Re_eff = my_params["Re"]*Umax*L
+````
+
+Initialize the solver.
+
+````@example lavd
+tspan = (0.0,10.0)
+integrator = init(u0,tspan,sys,alg=LiskaIFHERK(saddlesolver=CG))
+
+step!(integrator)
+@time step!(integrator,9.0)
+````
+
+Plot the solutions.
+
+````@example lavd
+sol = integrator.sol
+plt = plot(layout = (5,2), size = (1000, 1200), legend=:false)
+tsnap = 0.0:1.0:9.0
+for (i, t) in enumerate(tsnap)
+    plot!(plt[i],vorticity(sol,sys,t),sys,layers=false,title="t = $(round(t,digits=2))",clim=(-5,5),levels=range(-5,5,length=30),color = :RdBu, xlim = (-0.5, 5.5))
+    plot!(plt[i],surfaces(sol,sys,t))
+end
+plt
+````
+
+# Compute LAVD
+
+## Generate a Sequence of Velocity and Vorticity Fields
+This step obtains the computed velocity and vorticity fields at a sequence of times, and stores them as a sequence of interpolatable
+fields. This will greatly speed up how we compute the flow properties (i.e. vorticity) along trajectories.
+
+````@example lavd
+t_start = 0.0
+t_end = 8.5
+dt = timestep(u0,sys)
+tr = t_start:dt:t_end
+
+velxy = velocity_xy(sol,sys,tr) # Vector of interpolatable velocities
+velseq = VectorFieldSequence(tr,velxy); # Bundle together with the time array
+vortxy = vorticity_xy(sol,sys,tr)
+vortseq = ScalarFieldSequence(tr,vortxy);
+nothing #hide
+````
+
+## Generate Initial Conditions
+Here, we generate a grid of initial locations from which to integrate trajectories.
+
+````@example lavd
+X_MIN = -0.5
+X_MAX = 5.5
+Y_MIN = -2.0
+Y_MAX = 1.0
+dx = 0.04
+lavdgrid = PhysicalGrid((X_MIN,X_MAX),(Y_MIN,Y_MAX),dx)
+lavd_cache = SurfaceScalarCache(lavdgrid)
+x0, y0 = x_grid(lavd_cache), y_grid(lavd_cache)
+````
+
+## Solve the IVP and Generate LAVD Fields
+
+Compute trajectories
+
+````@example lavd
+T = -2.0
+t0 = 8.5
+t1 = t0 + T
+
+@time traj = compute_trajectory(velseq, (x0, y0), (t0, t1));
+nothing #hide
+````
+
+Compute LAVD
+
+````@example lavd
+LAVD = similar(x0)
+compute_LAVD!(LAVD, traj, X_MIN, Y_MIN, X_MAX, Y_MAX, vortseq)
+plot(LAVD, lavd_cache, colorbar = true, levels = 20)
+plot!(surfaces(sol,sys,t0))
+#savefig("lavd")
+````
+
+Compute IVD
+
+````@example lavd
+IVD = similar(x0)
+compute_IVD!(IVD, traj, X_MIN, Y_MIN, X_MAX, Y_MAX, vortseq)
+plot(IVD, lavd_cache, colorbar = true, levels = 20)
+plot!(surfaces(sol,sys,t0))
+#savefig("ivd")
+````
+
+Plot the vorticity fields
+
+````@example lavd
+plot(vorticity(sol, sys, t0), sys, layers = false,clim = (-5,5),levels = range(-5,5,length=30), colorbar = true, xlim = (-0.5, 5.5))
+plot!(surfaces(sol,sys,t0))
+#savefig("vorticity")
+````
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+
diff --git a/docs/src/manual/pod.md b/docs/src/manual/pod.md
@@ -125,11 +125,6 @@ The energy associated with this mode is
 modes.lambda[end]
 ````
 
-## POD functions
-```@docs
-pod
-```
-
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*

diff --git a/examples/dmdtest.ipynb b/examples/dmdtest.ipynb