fix optimisation of poles settings

examples/BUG/modelgeneration.org

@@ -295,13 +295,13 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output
   bdf = BDF()
   bdf.read_bdf(bdf_file)
   node_ids = bdf.node_ids
-  assert len(modes.shape[1]) == len(node_ids), "the modes size does not match the node_ids"
+  assert modes.shape[1] == len(node_ids), "the modes size does not match the node_ids"
   sorted_nodeids = sorted(node_ids)
   asets_ids = bdf.asets[0].node_ids
   asets_ids_sorted = sorted(asets_ids)
   asets_idsfull = np.array([sorted_nodeids.index(ai) for ai in asets_ids_sorted])
   asets_indexes = np.hstack([[6*i + j for j in range(6)] for i in asets_idsfull])
-  modes4simulations = modes[asets_indexes, :]
+  #modes4simulations = modes[asets_indexes, :]
   SAVE = False
   if SAVE:
       np.save(f"./FEM/eigenvecs_{num_modes}.npy", modes4simulations.T)
@@ -327,7 +327,6 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output
       np.save(f"./FEM/eigenvals_{num_modes}.npy", eigs)
 #+end_src
 
-
 *** Plot VTK modes
 Plot the modal shapes in Paraview
 #+begin_src python
@@ -1520,72 +1519,77 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output
   import feniax.unastran.aero as nasaero
   import feniax.unastran.op4handler as op4handler
   import pickle
-  
+
   #op4m = op4.OP4()
   #Qop4 = op4m.read_op4(file_name)
 
   Qhh = op4handler.read_data('./NASTRAN/data_out/Qhh<<parameters_gafs(output="label_gaf")>>.op4',
                              'Q_HH')
   Qhj = op4handler.read_data('./NASTRAN/data_out/Qhj<<parameters_gafs(output="label_gaf")>>.op4',
                              'Q_HJ')
-
+  num_poles = 5
 #+end_src
 
+**** RFA testing for Qhh
 #+begin_src python  :results none  :noweb yes :tangle rogerRFA.py  
   import importlib
   import feniax.aeromodal.roger as roger
   importlib.reload(roger)
 
   optpoles = roger.OptimisePoles(reduced_freqs, Qhh,
-                                 num_poles_=9,
+                                 num_poles_=num_poles,
                                  poles_step_=0.1,
                                  poles_range_=[0.05,1],
-                                 rfa_method_=2)
+                                 rfa_method_=1
+                                 )
+  # optpoles.set_errsettings(#error_name="max",
+  #                          rfa_method=2,
+  #                          norm_order=None)
   optpoles.run(show_info=True)
+
   qhhr1 = optpoles.get_model(label='m1')
-  # optpoles.set_errsettings(error_name="average", rfa_method=2, norm_order=None)
-  # optpoles.run(show_info=True)
-  # qhhr2 = optpoles.get_model(label='m2')
-  # optpoles.set_errsettings(error_name="max", rfa_method=2, norm_order=None)
-  # optpoles.run(show_info=True)
-  # optpoles.save("./AERO", f"A{label_gaf}", f"Poles{label_gaf}")
-  # qhhr3 = optpoles.get_model(label='m3')
+  # optimize with method 2
+  optpoles.set_errsettings(error_name="average", rfa_method=2, norm_order=None)
+  optpoles.run(show_info=True)
+  qhhr2 = optpoles.get_model(label='m2')
+  # optimize for max function
+  optpoles.set_errsettings(error_name="max", rfa_method=2, norm_order=None)
+  optpoles.run(show_info=True)
+  optpoles.save("./AERO", f"A{label_gaf}", f"Poles{label_gaf}")
+  qhhr3 = optpoles.get_model(label='m3')
 
 #+end_src
 
 #+begin_src python  :results none  :noweb yes :tangle rogerRFA.py  
-  poles = jnp.load("./AERO/PolesDd1c7F1Scao-50.npy")
-  rogerhj = roger.ComputeRoger(Qhj, reduced_freqs, poles, 2)
-  jnp.save(f"./AERO/D{label_gaf}.npy", rogerhj.roger_matrices)
-  rogerhjeval = roger.EvaluateRoger.create(rogerhj)
-  Qrogerhj = rogerhjeval.eval_array(reduced_freqs)
-#+end_src
-
-
-#+begin_src python  :results none  :noweb yes :tangle rogerRFA.py  
+  import importlib
+  import feniax.aeromodal.roger as roger
+  importlib.reload(roger)
 
   optpoles = roger.OptimisePoles(reduced_freqs, Qhh,
-                                 num_poles_=4,
-                                 poles_step_=0.05,
-                                 poles_range_=[0.05,1])
-  optpoles.set_errsettings(error_name="max", rfa_method=2, norm_order=None)
+                                 num_poles_=num_poles*2,
+                                 poles_step_=0.1,
+                                 poles_range_=[0.05,1],
+                                 rfa_method_=1
+                                 )
+  optpoles.set_errsettings(error_name="max",
+                           rfa_method=2, norm_order=None)
   optpoles.run(show_info=True)
-  qhhr4 = optpoles.get_model(label='m1')
-  
+  qhhr4 = optpoles.get_model(label='m4')
+
 #+end_src
 
 
 #+begin_src python  :results none  :noweb yes :tangle rogerRFA.py  
   Qroger1 = qhhr1.eval_array(reduced_freqs)
-  # Qroger2 = qhhr2.eval_array(reduced_freqs)
-  # Qroger3 = qhhr3.eval_array(reduced_freqs)
-  # Qroger4 = qhhr4.eval_array(reduced_freqs)
-
+  Qroger2 = qhhr2.eval_array(reduced_freqs)
+  Qroger3 = qhhr3.eval_array(reduced_freqs)
+  Qroger4 = qhhr4.eval_array(reduced_freqs)
 #+end_src
 
+**** Plot Qhh
 #+begin_src python  :results none  :noweb yes :tangle rogerRFA.py
 
-  roger.plot_gafs(0, 1, Qhh, [Qroger, Qroger2, Qroger3, Qroger4])
+  roger.plot_gafs(0, 1, Qhh, [Qroger1, Qroger2, Qroger3, Qroger4])
 
   #roger.plot_gafs(20, 2, Qhh, [Qroger, Qroger2, Qroger3])
 #+end_src
@@ -1598,7 +1602,6 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output
           roger.plot_gafs(li[0], li[1], Qhh, [Qroger1])
 #+end_src
 
-
 #+begin_src python  :results none  :noweb yes :tangle rogerRFA.py
   iterate_vect = list(range(10))
   plot_prod = list(itertools.product(iterate_vect,iterate_vect))
@@ -1607,6 +1610,29 @@ Running Nastran using the tailored functions in run_nastra.sh which moves output
           roger.plot_gafs(li[0], li[1], Qhh, [Qroger, Qroger2, Qroger3, Qroger4])
 #+end_src
 
+**** RFA for Qhw
+
+#+begin_src python  :results none  :noweb yes :tangle rogerRFA.py  
+  poles = jnp.load("./AERO/PolesDd1c7F1Scao-50.npy")
+  rogerhj = roger.ComputeRoger(Qhj, reduced_freqs, poles, 2)
+  jnp.save(f"./AERO/D{label_gaf}.npy", rogerhj.roger_matrices)
+  rogerhjeval = roger.EvaluateRoger.create(rogerhj)
+  Qrogerhj = rogerhjeval.eval_array(reduced_freqs)
+#+end_src
+
+
+#+begin_src python  :results none  :noweb yes :tangle rogerRFA.py  
+
+  optpoles = roger.OptimisePoles(reduced_freqs, Qhh,
+                                 num_poles_=8,
+                                 poles_step_=0.1,
+                                 poles_range_=[0.05,1])
+  optpoles.set_errsettings(#error_name="max",
+                           rfa_method=2, norm_order=None)
+  optpoles.run(show_info=True)
+  qhhr4 = optpoles.get_model(label='m1')
+
+#+end_src
 
 ** Steady
 *** bdf models

feniax/aeromodal/roger.py

@@ -84,9 +84,9 @@ def __init__(self, redfreqs_, sampling_aeromatrices_,
         self.num_poles = num_poles_
         self.poles_step = poles_step_
         self.poles_range = poles_range_
-        self.rfa_method = rfa_method_
+        # self.rfa_method = rfa_method_
         self._build_polesgrid()
-        self.set_errsettings()
+        self.set_errsettings(rfa_method=rfa_method_)
 
     def get_model(self, label='m1'):