hvlm tighter tolerance + cont print rms0

Author: samayala22

lib/libhbvlm3.cpp

@@ -413,7 +413,7 @@ void HBVLM::run(f64 t_start, f64 omega, const Vars& v, Assembly<f64>& assembly,
         gamma_coeffs_v.reshape(gamma_coeffs_v.shape(0)*gamma_coeffs_v.shape(1)),
         hb_vlm_iter,
         100, // max iterations
-        1e-9, // tolerance
+        1e-12, // tolerance
         10 // history
     );
 

scripts/continuation.py

@@ -355,6 +355,9 @@ def continuation(X0, motion, metadata: Metadata):
             X0 = Xtmp.copy()
             X = X0 * Dscale # unscaled X
             X_mat[:, iteration] = X
+            X_h = X[0:omega_idx:metadata.dims.n_d]
+            A = np.sqrt(X_h[1::2]**2 + X_h[2::2]**2)
+            rms_0 = np.sqrt(X_h[0]**2 + 0.5 * np.sum(A**2, axis=0))
 
             # Stability analysis
             is_stable, floquet_exponents = stability_analysis(J @ np.diag(1.0 / Dscale), motion, X[omega_idx], omega_idx, metadata.dims)
@@ -374,8 +377,8 @@ def continuation(X0, motion, metadata: Metadata):
                     if np.copysign(1.0, metadata.bifurcation_test[i, iteration-1]) != np.copysign(1.0, metadata.bifurcation_test[i, iteration]):
                         print(f"{name} bifurcation detected")
                         metadata.bifurcation[name].append(iteration-1)
-            
-            print(f"{iteration} | ds: {ds * Dscale[-1]:.4f}, param: {X[-1]:.3f}, omega: {X[omega_idx]:.3f}, stable: {is_stable}, nfev: {nfev}, njev: {njev}, timing: {timing:.2f}s")
+
+            print(f"{iteration} | ds: {ds * Dscale[-1]:.4f}, param: {X[-1]:.3f}, omega: {X[omega_idx]:.3f}, rms0: {rms_0:.3e}, stable: {is_stable}, nfev: {nfev}, njev: {njev}, timing: {timing:.2f}s")
 
             if metadata.step_adapt:
                 n_f = nfev + J.shape[1] * njev
@@ -389,10 +392,6 @@ def continuation(X0, motion, metadata: Metadata):
                 print("Continuation reached the end")
                 break
 
-            X_h = X[0:omega_idx:metadata.dims.n_d]
-            A = np.sqrt(X_h[1::2]**2 + X_h[2::2]**2)
-            rms_0 = np.sqrt(X_h[0]**2 + 0.5 * np.sum(A**2, axis=0))
-            # print(f"rms_0: {rms_0:.6e}")
             if rms_0 < 1e-10 and iteration > 1:
                 print("Continuation reached trivial solution, stopping.")
                 break

scripts/dof3_hb.py

@@ -120,7 +120,7 @@ def dKdU(U): return np.zeros((3, 3))
     flutter_speed = 23.9
     # param_start = flutter_speed * 0.3
     # param_end = flutter_speed * 0.6
-    param_start = 6.0
+    param_start = 8.0
     param_end = 20.0
 
     v = dof3.Vars()
@@ -189,7 +189,7 @@ def dKdU(U): return np.zeros((3, 3))
     metadata.name = f"3dof_hbvlm_{torsional_spring_names[torsional_spring]}"
     metadata.param_start = param_start
     metadata.param_end = param_end
-    metadata.max_steps = 1 if INITIAL_ONLY else 10000
+    metadata.max_steps = 1 if INITIAL_ONLY else 5000
     metadata.scaling = True
     metadata.step_adapt = True
     metadata.ds = 0.02
@@ -202,21 +202,18 @@ def dKdU(U): return np.zeros((3, 3))
 
     if helpers.getenv("PLOT"):
         X_mat = metadata.X
-        if X_mat.shape[1] == 1:
-            hb_sol_t, hb_sol0 = hb.to_timedomain(vec_t, dims.n_d, X_mat[:-2, 0], X_mat[-2, 0], dims.n_h)
-            aero_forces = system.aero_forces(sol.y)
-            fig = plot.create_dofs_figure(["Heave", "Pitch", "Control"])
-            dof3.plot_solution(fig, aero_forces, sol, v)
-
-            plot.add_data_and_psd(fig, hb_sol_t, hb_sol0[0, :], "HB-VLM", 1, 1, 3)
-            plot.add_data_and_psd(fig, hb_sol_t, np.degrees(hb_sol0[1, :]), "HB-VLM", 3, 1, 3)
-            plot.add_data_and_psd(fig, hb_sol_t, np.degrees(hb_sol0[2, :]), "HB-VLM", 5, 1, 3)
-
-            plot.add_data_and_psd(fig, hb_sol_t, hb_sol0[3, :], "HB-VLM", 1, 2, 3)
-            plot.add_data_and_psd(fig, hb_sol_t, np.degrees(hb_sol0[4, :]), "HB-VLM", 3, 2, 3)
-            plot.add_data_and_psd(fig, hb_sol_t, np.degrees(hb_sol0[5, :]), "HB-VLM", 5, 2, 3)
-            
-            dof3.format_plot(fig)
-            plot.fig_save(fig, f"build/3dof/hbvlm0", html=True, pdf=False)
-        else:
-            cont.plot_hb_continuation(metadata)
+        hb_sol_t, hb_sol0 = hb.to_timedomain(vec_t, dims.n_d, X_mat[:-2, 0], X_mat[-2, 0], dims.n_h)
+        aero_forces = system.aero_forces(sol.y)
+        fig = plot.create_dofs_figure(["Heave", "Pitch", "Control"])
+        dof3.plot_solution(fig, aero_forces, sol, v)
+
+        plot.add_data_and_psd(fig, hb_sol_t, hb_sol0[0, :], "HB-VLM", 1, 1, 3)
+        plot.add_data_and_psd(fig, hb_sol_t, np.degrees(hb_sol0[1, :]), "HB-VLM", 3, 1, 3)
+        plot.add_data_and_psd(fig, hb_sol_t, np.degrees(hb_sol0[2, :]), "HB-VLM", 5, 1, 3)
+
+        plot.add_data_and_psd(fig, hb_sol_t, hb_sol0[3, :], "HB-VLM", 1, 2, 3)
+        plot.add_data_and_psd(fig, hb_sol_t, np.degrees(hb_sol0[4, :]), "HB-VLM", 3, 2, 3)
+        plot.add_data_and_psd(fig, hb_sol_t, np.degrees(hb_sol0[5, :]), "HB-VLM", 5, 2, 3)
+        
+        dof3.format_plot(fig)
+        plot.fig_save(fig, f"build/3dof/hbvlm0", html=True, pdf=False)