From 642c64622889ce90342a97a7d0c41408e739b98a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 1 Jan 2023 16:02:51 -0500 Subject: [PATCH 1/3] Fixed Python driver bugs --- python/swiftest/swiftest/simulation_class.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index e4fa80049..612aeb4b0 100644 --- 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,8 +2866,8 @@ 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: @@ -2880,6 +2875,7 @@ def save(self, 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) From b63ca31aa8421e0a4ae99201d7479a2b2e337efe Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 1 Jan 2023 17:26:39 -0500 Subject: [PATCH 2/3] Fixed bug that was throwing a segfault when not in debug mode --- src/swiftest/swiftest_util.f90 | 179 ++++++++++++++++----------------- 1 file changed, 88 insertions(+), 91 deletions(-) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 45b91aa73..73b4ef5b6 100644 --- 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 From 86a707f99edff45c25617da9c8e21745ea262aca Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 2 Jan 2023 11:02:15 -0500 Subject: [PATCH 3/3] Fixed bugs prevening SyMBA variables from being written to dump parameter file --- python/swiftest/swiftest/init_cond.py | 4 ++-- src/swiftest/swiftest_io.f90 | 13 ++++++++++--- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index 8dcf08cf2..02ff6f8e1 100644 --- 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/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index fe4ecba18..8ae50ba2d 100644 --- 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,6 +2274,7 @@ 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) @@ -2282,6 +2282,14 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i 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 !!