Implementing SAV bound

Author: thomas-risse

python/base_model.py

@@ -1,9 +1,11 @@
 import numpy as np
 
+
 class Model():
     """General Model class to pass to the SAVSolver.
     """
-    def __init__(self, N = 10):
+
+    def __init__(self, N=10):
         # Number of dofs
         self.N = N
         # Number of inputs
@@ -22,9 +24,9 @@ def __init__(self, N = 10):
         self.Enl = lambda q: (0.25 * np.sum(q**4)) * 1
         self.Fnl = lambda q: (q**3) * 1
         # Both functions are called at the same time in the solver,
-        # in some cases it is then computationally interesting 
+        # in some cases it is then computationally interesting
         # to compute both in the same function.
         self.EandFnl = lambda q: (self.Enl(q), self.Fnl(q))
 
     def setting(self):
-        return {"Name": self.__class__.__name__, "N": self.N} 
+        return {"Name": self.__class__.__name__, "N": self.N}

python/plot_drift_JAES.py

@@ -21,7 +21,7 @@
 Generates figures of section "to drift or not to drift".
 """
 
-#%%
+# %%
 # System description
 
 # Physical parameters
@@ -37,30 +37,32 @@
 T_60_1000 = 3
 
 # Deduce missing physical parameters
-StringParams["T"], StringParams["l0"] = get_T_and_l0_from_f0_beta(f0, beta, StringParams)
-StringParams["eta_0"], StringParams["eta_1"] = get_etas_from_decays(T_60_0, T_60_1000, StringParams)
+StringParams["T"], StringParams["l0"] = get_T_and_l0_from_f0_beta(
+    f0, beta, StringParams)
+StringParams["eta_0"], StringParams["eta_1"] = get_etas_from_decays(
+    T_60_0, T_60_1000, StringParams)
 
 print(StringParams)
 
 
 model = FD_string_model(44100, **StringParams)
 
-modes = ["KC", "GE", "GE4"] # Nonlinear modes
+modes = ["KC", "GE", "GE4"]  # Nonlinear modes
 
-#%% 
+# %%
 # Simulation parameters
 sr = 44100
 duration = 10
 kappa = 0.9
 lambda0s = [0, 1000]
-OF = 2 # Over-sampling factor for reference
+OF = 2  # Over-sampling factor for reference
 
 # Deduce discretization from stability condition
 dt = 1 / sr
-model.recompute_stability(sr, kappa = kappa)
+model.recompute_stability(sr, kappa=kappa)
 
 
-#%%
+# %%
 # Initial conditions and excitation
 
 # External force (applied at the middle of the string)
@@ -69,51 +71,64 @@ def Fext(t):
     width = 2e-3
     period = 500 * width
     out = np.zeros(1)
-    out[0] = Amp * np.sin(np.pi * t / width) * (t%period < width)
+    out[0] = Amp * np.sin(np.pi * t / width) * (t % period < width)
     return out
+
+
 q0 = np.zeros(model.N)
 u0 = np.zeros(model.N)
 
 
-#%% Run simulations and plot results
-fig, axs = plt.subplots(1 + 2*len(lambda0s), 1, figsize=set_size("JAES", height_ratio=0.6), sharex=True)
+# %% Run simulations and plot results
+fig, axs = plt.subplots(1 + 2*len(lambda0s), 1,
+                        figsize=set_size("JAES", height_ratio=0.6), sharex=True)
 linestyles = [":", "--", "-."]
 for i, lambda0 in enumerate(lambda0s):
 
-  for j, mode in enumerate(modes):
-    model.NL_type = mode
-    # Compute SAV solution
-    solver = SAVSolver(model, sr, lambda0)
-
-    storage = STATE_STORAGE_CONFIG
-    storage["Drift"] = True
-    storage["q_idx"] = np.array([model.N//2 + 1])
-    storage["p_idx"] = None
-
-    solver.integrate(q0, u0, Fext, duration, ConstantRmid=True, plotter_config=NO_PLOTTER_CONFIG, storage_config=storage)
-    solver.storage.write(os.path.join(result_folder, f"{mode}/sr{sr}_lambda{lambda0}.h5"))
-
-    f0, _, _ = librosa.pyin(solver.storage.q[:, 0], fmin = 40, fmax = 200, sr = 44100, frame_length = 2048 * 4)
-    write(os.path.join(result_folder, f"{mode}/sr{sr}_lambda{lambda0}.wav"), sr, solver.storage.q[:, 0] / np.max(np.abs(solver.storage.q[:, 0])))
-
-    if mode==modes[0] and i==0:
-        axs[0].plot(solver.storage.t, [Fext(t) for t in solver.storage.t], color = "black", ls = linestyles[j])
-        axs[0].set_ylabel(r"$f_{in}$ [N]")
-    axs[2 * i+1].plot(np.linspace(0, duration, len(f0)), f0, label = mode, ls = linestyles[j])
-    # Here, we divide espilon by the max observed nonlinear energy to get a relative measure
-    print(solver.maxEnl)
-    axs[2 * i+2].semilogy(solver.storage.t, np.abs(solver.storage.epsilon / solver.maxEnl), label = mode)
-
-  axs[2*i+1].set_ylabel(r"$f_0$ [Hz]")
-  axs[2*i+1].set_ylim(80, 110)
-  axs[2*i+2].set_ylabel(r"$\vert\epsilon_{rel}\vert$")
-  axs[2*i+2].set_ylim([1e-4, 1.5e3])
-  axs[2*i+2].set_yticks([1e-4, 1e-1, 1e2])
-  axs[2*i+2].set_yticklabels([1e-4, 1e-1, 1e2])
-  axs[2*i+1].text(0.8, 0.6, fr"$\lambda_0 = {lambda0} s^{-1}$", transform = axs[2*i+1].transAxes, color="red", bbox=dict(facecolor='white', edgecolor='black', boxstyle='round'))
-  axs[2*i+2].text(0.8, 0.6, fr"$\lambda_0 = {lambda0} s^{-1}$", transform = axs[2*i+2].transAxes, color="red", bbox=dict(facecolor='white', edgecolor='black', boxstyle='round'))
-
-axs[1].legend(loc = "lower center", frameon = True, fancybox = True, bbox_to_anchor = (0.5, 2.1), ncol=3)
+    for j, mode in enumerate(modes):
+        model.NL_type = mode
+        # Compute SAV solution
+        solver = SAVSolver(model, sr, lambda0)
+
+        storage = STATE_STORAGE_CONFIG
+        storage["Drift"] = True
+        storage["q_idx"] = np.array([model.N//2 + 1])
+        storage["p_idx"] = None
+
+        solver.integrate(q0, u0, Fext, duration, ConstantRmid=True,
+                         plotter_config=NO_PLOTTER_CONFIG, storage_config=storage, BoundG=True)
+        solver.storage.write(os.path.join(
+            result_folder, f"{mode}/sr{sr}_lambda{lambda0}.h5"))
+
+        f0, _, _ = librosa.pyin(
+            solver.storage.q[:, 0], fmin=40, fmax=200, sr=44100, frame_length=2048 * 4)
+        write(os.path.join(result_folder, f"{mode}/sr{sr}_lambda{lambda0}.wav"),
+              sr, solver.storage.q[:, 0] / np.max(np.abs(solver.storage.q[:, 0])))
+
+        if mode == modes[0] and i == 0:
+            axs[0].plot(solver.storage.t, [Fext(t)
+                        for t in solver.storage.t], color="black", ls=linestyles[j])
+            axs[0].set_ylabel(r"$f_{in}$ [N]")
+        axs[2 * i+1].plot(np.linspace(0, duration, len(f0)),
+                          f0, label=mode, ls=linestyles[j])
+        # Here, we divide espilon by the max observed nonlinear energy to get a relative measure
+        print(solver.maxEnl)
+        axs[2 * i+2].semilogy(solver.storage.t,
+                              np.abs(solver.storage.epsilon / solver.maxEnl), label=mode)
+
+    axs[2*i+1].set_ylabel(r"$f_0$ [Hz]")
+    axs[2*i+1].set_ylim(80, 110)
+    axs[2*i+2].set_ylabel(r"$\vert\epsilon_{rel}\vert$")
+    axs[2*i+2].set_ylim([1e-4, 1.5e3])
+    axs[2*i+2].set_yticks([1e-4, 1e-1, 1e2])
+    axs[2*i+2].set_yticklabels([1e-4, 1e-1, 1e2])
+    axs[2*i+1].text(0.8, 0.6, fr"$\lambda_0 = {lambda0} s^{-1}$", transform=axs[2*i+1].transAxes,
+                    color="red", bbox=dict(facecolor='white', edgecolor='black', boxstyle='round'))
+    axs[2*i+2].text(0.8, 0.6, fr"$\lambda_0 = {lambda0} s^{-1}$", transform=axs[2*i+2].transAxes,
+                    color="red", bbox=dict(facecolor='white', edgecolor='black', boxstyle='round'))
+
+axs[1].legend(loc="lower center", frameon=True, fancybox=True,
+              bbox_to_anchor=(0.5, 2.1), ncol=3)
 axs[4].set_xlim(0, duration)
 axs[4].set_ylim(1e-8, 10)
 axs[4].set_xlabel(r"Time [s]")
@@ -124,4 +139,5 @@ def Fext(t):
 fig.align_ylabels(axs)
 fig.subplots_adjust(hspace=0.1, wspace=0.4)
 # Save figure
-fig.savefig(os.path.join(result_folder, f"test_drift_nl_force.pdf"), bbox_inches='tight')
+fig.savefig(os.path.join(result_folder, f"test_drift_nl_force.pdf"),
+            bbox_inches='tight')