Merge branch 'debug' into fraggle_simplification

MintonGroup · Jan 2, 2023 · f36e565 · f36e565
2 parents 613e5d4 + 86a707f
commit f36e565
Show file tree

Hide file tree

Showing 4 changed files with 107 additions and 107 deletions.
diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py
@@ -285,8 +285,8 @@ def vec2xr(param: Dict, **kwargs: Any):
     if param['ROTATION']:
         if "rot" not in kwargs and "Gmass" in kwargs:
             kwargs['rot'] = np.zeros((len(kwargs['Gmass']),3))
-        if "Ip" not in kwargs and "rot" in kwargs:
-            kwargs['Ip'] = np.full_like(kwargs['rot'], 0.4)
+        if "Ip" not in kwargs and "Gmass" in kwargs:
+            kwargs['Ip'] = np.full_like(kwargs['Gmass'], 0.4)
 
     if "time" not in kwargs:
         kwargs["time"] = np.array([0.0])

diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py
@@ -328,8 +328,6 @@ def __init__(self,read_param: bool = False, read_old_output: bool = False, simdi
         else:
             if read_old_output or read_param:
                 raise NotADirectoryError(f"Cannot find directory {self.simdir.resolve()} ")
-            else:
-                self.simdir.mkdir(parents=True, exist_ok=False)
 
         # Set the location of the parameter input file, choosing the default if it isn't specified.
         param_file = kwargs.pop("param_file",Path.cwd() / self.simdir / "param.in")
@@ -1279,9 +1277,8 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N
                 key = valid_var[arg]
                 if key in feature_dict:
                     if arg == "minimum_fragment_gmass":
-                        if self.param['FRAGMENTATION']:
-                           print(f"{arg:<{self._getter_column_width}} {feature_dict[key]} {self.DU_name}^3 / {self.TU_name}^2 ")
-                           print(f"{'minimum_fragment_mass':<{self._getter_column_width}} {feature_dict[key] / self.GU} {self.MU_name}")
+                        print(f"{arg:<{self._getter_column_width}} {feature_dict[key]} {self.DU_name}^3 / {self.TU_name}^2 ")
+                        print(f"{'minimum_fragment_mass':<{self._getter_column_width}} {feature_dict[key] / self.GU} {self.MU_name}")
                     else:
                         print(f"{arg:<{self._getter_column_width}} {feature_dict[key]}")
                 else:
@@ -2174,7 +2171,6 @@ def add_solar_system_body(self,
 
         return
 
-
     def set_ephemeris_date(self,
                           ephemeris_date: str | None = None,
                           verbose: bool | None = None,
@@ -2638,7 +2634,8 @@ def write_param(self,
         if verbose:
             print(f"Writing parameter inputs to file {param_file}")
         param['! VERSION'] = f"{codename} input file"
-
+
+        self.simdir.mkdir(parents=True, exist_ok=True)
         # Check to see if the parameter type matches the output type. If not, we need to convert
         if codename == "Swifter" or codename == "Swiftest":
             if param['IN_TYPE'] == "ASCII":
@@ -2781,7 +2778,6 @@ def _preprocess(ds, param):
 
         return
 
-
     def read_collisions(self):
         col_files = glob(f"{self.simdir}{os.path.sep}collision_*.nc")
         if len(col_files) == 0:
@@ -2801,7 +2797,6 @@ def _preprocess(ds, param):
 
         return
 
-
     def follow(self, codestyle="Swifter"):
         """
         An implementation of the Swift tool_follow algorithm. Under development. Currently only for Swift simulations. 
@@ -2871,15 +2866,16 @@ def save(self,
         if "verbose" in kwargs:
             verbose = kwargs['verbose']
         else:
-            verbose = self%verbose
-
+            verbose = self.verbose
+        
         if codename is None:
             codename = self.codename
         if param_file is None:
             param_file = self.param_file
         if param is None:
             param = self.param
 
+        self.simdir.mkdir(parents=True, exist_ok=True)
         if codename == "Swiftest":
             infile_name = Path(self.simdir) / param['NC_IN']
             io.swiftest_xr2infile(ds=self.data, param=param, in_type=self.param['IN_TYPE'], infile_name=infile_name, framenum=framenum, verbose=verbose)

diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90
@@ -2094,7 +2094,6 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i
             return
          end if
 
-
          if ((param%collision_model /= "MERGE")       .and. &
              (param%collision_model /= "BOUNCE")    .and. &
              (param%collision_model /= "FRAGGLE")) then
@@ -2127,7 +2126,6 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i
             param%inv_c2 = (param%inv_c2)**(-2)
          end if
 
-
          select case(trim(adjustl(param%interaction_loops)))
          case("ADAPTIVE")
             param%ladaptive_interactions = .true.
@@ -2230,6 +2228,7 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i
       character(*),parameter :: Ifmt  = '(I0)'         !! Format label for integer values
       character(*),parameter :: Rfmt  = '(ES25.17)'    !! Format label for real values 
       character(*),parameter :: Lfmt  = '(L1)'         !! Format label for logical values 
+      integer(I4B)                   :: nseeds
 
       associate(param => self)
          call io_param_writer_one("T0", param%t0, unit)
@@ -2275,13 +2274,22 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i
          call io_param_writer_one("INTERACTION_LOOPS", param%interaction_loops, unit)
          call io_param_writer_one("ENCOUNTER_CHECK_PLPL", param%encounter_check_plpl, unit)
          call io_param_writer_one("ENCOUNTER_CHECK_PLTP", param%encounter_check_pltp, unit)
+         call io_param_writer_one("ENCOUNTER_SAVE", param%encounter_save, unit)
 
          if (param%lenergy) then
             call io_param_writer_one("FIRSTENERGY", param%lfirstenergy, unit)
          end if
          call io_param_writer_one("FIRSTKICK",param%lfirstkick, unit)
          call io_param_writer_one("MAXID",param%maxid, unit)
          call io_param_writer_one("MAXID_COLLISION",param%maxid_collision, unit)
+
+         if (param%GMTINY > 0.0_DP) call io_param_writer_one("GMTINY",param%GMTINY, unit)
+         if (param%min_GMfrag > 0.0_DP) call io_param_writer_one("MIN_GMFRAG",param%min_GMfrag, unit)
+         call io_param_writer_one("COLLISION_MODEL",param%collision_model, unit)
+         if (param%collision_model == "FRAGGLE" ) then
+            nseeds = size(param%seed)
+            call io_param_writer_one("SEED", [nseeds, param%seed(:)], unit)
+         end if
 
          iostat = 0
          iomsg = "UDIO not implemented"
@@ -2837,7 +2845,6 @@ module subroutine swiftest_io_toupper(string)
    end subroutine swiftest_io_toupper
 
 
-
    module subroutine swiftest_io_write_discard(self, param)
       !! author: David A. Minton
       !!

diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90
@@ -2552,104 +2552,101 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param)
       type(encounter_storage)  :: encounter_history
       type(collision_storage)  :: collision_history
 
-      select type(param)
-      class is (swiftest_parameters)
-         allocate(swiftest_storage(param%dump_cadence) :: param%system_history)
-         allocate(swiftest_netcdf_parameters :: param%system_history%nc)
-         call param%system_history%reset()
-
-         select case(param%integrator)
-         case (INT_BS)
-            write(*,*) 'Bulirsch-Stoer integrator not yet enabled'
-         case (INT_HELIO)
-            allocate(helio_nbody_system :: nbody_system)
-            select type(nbody_system)
-            class is (helio_nbody_system)
-               allocate(helio_cb :: nbody_system%cb)
-               allocate(helio_pl :: nbody_system%pl)
-               allocate(helio_tp :: nbody_system%tp)
-               allocate(helio_tp :: nbody_system%tp_discards)
-            end select
-            param%collision_model = "MERGE"
-         case (INT_RA15)
-            write(*,*) 'Radau integrator not yet enabled'
-         case (INT_TU4)
-            write(*,*) 'INT_TU4 integrator not yet enabled'
-         case (INT_WHM)
-            allocate(whm_nbody_system :: nbody_system)
-            select type(nbody_system)
-            class is (whm_nbody_system)
-               allocate(whm_cb :: nbody_system%cb)
-               allocate(whm_pl :: nbody_system%pl)
-               allocate(whm_tp :: nbody_system%tp)
-               allocate(whm_tp :: nbody_system%tp_discards)
-            end select
-            param%collision_model = "MERGE"
-         case (INT_RMVS)
-            allocate(rmvs_nbody_system :: nbody_system)
-            select type(nbody_system)
-            class is (rmvs_nbody_system)
-               allocate(rmvs_cb :: nbody_system%cb)
-               allocate(rmvs_pl :: nbody_system%pl)
-               allocate(rmvs_tp :: nbody_system%tp)
-               allocate(rmvs_tp :: nbody_system%tp_discards)
-            end select
-            param%collision_model = "MERGE"
-         case (INT_SYMBA)
-            allocate(symba_nbody_system :: nbody_system)
-            select type(nbody_system)
-            class is (symba_nbody_system)
-               allocate(symba_cb :: nbody_system%cb)
-               allocate(symba_pl :: nbody_system%pl)
-               allocate(symba_tp :: nbody_system%tp)
-
-               allocate(symba_tp :: nbody_system%tp_discards)
-               allocate(symba_pl :: nbody_system%pl_adds)
-               allocate(symba_pl :: nbody_system%pl_discards)
-
-               allocate(symba_list_pltp :: nbody_system%pltp_encounter)
-               allocate(symba_list_plpl :: nbody_system%plpl_encounter)
-               allocate(collision_list_plpl :: nbody_system%plpl_collision)
-
-               if (param%lenc_save_trajectory .or. param%lenc_save_closest) then
-                  allocate(encounter_netcdf_parameters :: encounter_history%nc)
-                  call encounter_history%reset()
-                  select type(nc => encounter_history%nc)
-                  class is (encounter_netcdf_parameters)
-                     nc%file_number = param%iloop / param%dump_cadence
-                  end select
-                  allocate(nbody_system%encounter_history, source=encounter_history)
-               end if
-
-               allocate(collision_netcdf_parameters :: collision_history%nc)
-               call collision_history%reset()
-               select type(nc => collision_history%nc)
-               class is (collision_netcdf_parameters)
+      allocate(swiftest_storage(param%dump_cadence) :: param%system_history)
+      allocate(swiftest_netcdf_parameters :: param%system_history%nc)
+      call param%system_history%reset()
+
+      select case(param%integrator)
+      case (INT_BS)
+         write(*,*) 'Bulirsch-Stoer integrator not yet enabled'
+      case (INT_HELIO)
+         allocate(helio_nbody_system :: nbody_system)
+         select type(nbody_system)
+         class is (helio_nbody_system)
+            allocate(helio_cb :: nbody_system%cb)
+            allocate(helio_pl :: nbody_system%pl)
+            allocate(helio_tp :: nbody_system%tp)
+            allocate(helio_tp :: nbody_system%tp_discards)
+         end select
+         param%collision_model = "MERGE"
+      case (INT_RA15)
+         write(*,*) 'Radau integrator not yet enabled'
+      case (INT_TU4)
+         write(*,*) 'INT_TU4 integrator not yet enabled'
+      case (INT_WHM)
+         allocate(whm_nbody_system :: nbody_system)
+         select type(nbody_system)
+         class is (whm_nbody_system)
+            allocate(whm_cb :: nbody_system%cb)
+            allocate(whm_pl :: nbody_system%pl)
+            allocate(whm_tp :: nbody_system%tp)
+            allocate(whm_tp :: nbody_system%tp_discards)
+         end select
+         param%collision_model = "MERGE"
+      case (INT_RMVS)
+         allocate(rmvs_nbody_system :: nbody_system)
+         select type(nbody_system)
+         class is (rmvs_nbody_system)
+            allocate(rmvs_cb :: nbody_system%cb)
+            allocate(rmvs_pl :: nbody_system%pl)
+            allocate(rmvs_tp :: nbody_system%tp)
+            allocate(rmvs_tp :: nbody_system%tp_discards)
+         end select
+         param%collision_model = "MERGE"
+      case (INT_SYMBA)
+         allocate(symba_nbody_system :: nbody_system)
+         select type(nbody_system)
+         class is (symba_nbody_system)
+            allocate(symba_cb :: nbody_system%cb)
+            allocate(symba_pl :: nbody_system%pl)
+            allocate(symba_tp :: nbody_system%tp)
+
+            allocate(symba_tp :: nbody_system%tp_discards)
+            allocate(symba_pl :: nbody_system%pl_adds)
+            allocate(symba_pl :: nbody_system%pl_discards)
+
+            allocate(symba_list_pltp :: nbody_system%pltp_encounter)
+            allocate(symba_list_plpl :: nbody_system%plpl_encounter)
+            allocate(collision_list_plpl :: nbody_system%plpl_collision)
+
+            if (param%lenc_save_trajectory .or. param%lenc_save_closest) then
+               allocate(encounter_netcdf_parameters :: encounter_history%nc)
+               call encounter_history%reset()
+               select type(nc => encounter_history%nc)
+               class is (encounter_netcdf_parameters)
                   nc%file_number = param%iloop / param%dump_cadence
                end select
-               allocate(nbody_system%collision_history, source=collision_history)
-
+               allocate(nbody_system%encounter_history, source=encounter_history)
+            end if
+
+            allocate(collision_netcdf_parameters :: collision_history%nc)
+            call collision_history%reset()
+            select type(nc => collision_history%nc)
+            class is (collision_netcdf_parameters)
+               nc%file_number = param%iloop / param%dump_cadence
             end select
-         case (INT_RINGMOONS)
-            write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled'
-         case default
-            write(*,*) 'Unkown integrator',param%integrator
-            call util_exit(FAILURE)
-         end select
+            allocate(nbody_system%collision_history, source=collision_history)
 
-         allocate(swiftest_particle_info :: nbody_system%cb%info)
-
-         select case(param%collision_model)
-         case("MERGE")
-            allocate(collision_basic :: nbody_system%collider)
-         case("BOUNCE")
-            allocate(collision_bounce :: nbody_system%collider)
-         case("FRAGGLE")
-            allocate(collision_fraggle :: nbody_system%collider)
          end select
-         call nbody_system%collider%setup(nbody_system)
+      case (INT_RINGMOONS)
+         write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled'
+      case default
+         write(*,*) 'Unkown integrator',param%integrator
+         call util_exit(FAILURE)
+      end select
 
+      allocate(swiftest_particle_info :: nbody_system%cb%info)
+
+      select case(param%collision_model)
+      case("MERGE")
+         allocate(collision_basic :: nbody_system%collider)
+      case("BOUNCE")
+         allocate(collision_bounce :: nbody_system%collider)
+      case("FRAGGLE")
+         allocate(collision_fraggle :: nbody_system%collider)
       end select
+      call nbody_system%collider%setup(nbody_system)
+
 
       return
    end subroutine swiftest_util_setup_construct_system