diff --git a/CMakeLists.txt b/CMakeLists.txt index 96ce5efa9..5d4c8637f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,18 +29,21 @@ ENDIF(NOT CMAKE_Fortran_COMPILER_SUPPORTS_F90) # Set some options the user may choose # Uncomment the below if you want the user to choose a parallelization library -OPTION(USE_MPI "Use the MPI library for parallelization" OFF) +OPTION(USE_MPI "Use the MPI library for parallelization" ON) OPTION(USE_OPENMP "Use OpenMP for parallelization" ON) -# This INCLUDE statement executes code that sets the compile flags for DEBUG, -# RELEASE, and TESTING. You should review this file and make sure the flags -# are to your liking. -INCLUDE(${CMAKE_MODULE_PATH}/SetFortranFlags.cmake) + # Locate and set parallelization libraries. There are some CMake peculiarities # taken care of here, such as the fact that the FindOpenMP routine doesn't know # about Fortran. INCLUDE(${CMAKE_MODULE_PATH}/SetParallelizationLibrary.cmake) INCLUDE(${CMAKE_MODULE_PATH}/SetUpNetCDF.cmake) +INCLUDE(${CMAKE_MODULE_PATH}/SetMKL.cmake) + +# This INCLUDE statement executes code that sets the compile flags for DEBUG, +# RELEASE, PROFILING, and TESTING. +INCLUDE(${CMAKE_MODULE_PATH}/SetFortranFlags.cmake) + # There is an error in CMAKE with this flag for pgf90. Unset it GET_FILENAME_COMPONENT(FCNAME ${CMAKE_Fortran_COMPILER} NAME) @@ -59,9 +62,10 @@ SET(SWIFTEST_DRIVER swiftest_driver) SET(SRC ${CMAKE_SOURCE_DIR}/src) SET(LIB ${CMAKE_SOURCE_DIR}/lib) SET(BIN ${CMAKE_SOURCE_DIR}/bin) +SET(MOD ${CMAKE_SOURCE_DIR}/include) # Have the .mod files placed in the lib folder -SET(CMAKE_Fortran_MODULE_DIRECTORY ${LIB}) +SET(CMAKE_Fortran_MODULE_DIRECTORY ${MOD}) # The source for the SWIFTEST binary and have it placed in the bin folder ADD_SUBDIRECTORY(${SRC} ${BIN}) diff --git a/cmake/Modules/FindMKL.cmake b/cmake/Modules/FindMKL.cmake new file mode 100644 index 000000000..9e48932c3 --- /dev/null +++ b/cmake/Modules/FindMKL.cmake @@ -0,0 +1,17 @@ +# Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +# This file is part of Swiftest. +# Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +# Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# You should have received a copy of the GNU General Public License along with Swiftest. +# If not, see: https://www.gnu.org/licenses. + +# - Finds the Intel MKL libraries +find_path(MKL_INCLUDE_DIR NAMES mkl.h HINTS ENV MKLROOT PATH_SUFFIXES include) +find_library(MKL_LIBRARY NAMES libmkl_core.a HINTS ENV MKLROOT PATH_SUFFIXES lib lib/intel64 ) + +set(MKL_FOUND TRUE) +set(MKL_INCLUDE_DIRS ${MKL_INCLUDE_DIR}) +set(MKL_LIBRARIES ${MKL_LIBRARY}) +mark_as_advanced(MKL_LIBRARY MKL_INCLUDE_DIR) \ No newline at end of file diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake index 7850fbdb8..d869e89b6 100644 --- a/cmake/Modules/SetFortranFlags.cmake +++ b/cmake/Modules/SetFortranFlags.cmake @@ -222,9 +222,9 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_TESTING "${CMAKE_Fortran_FLAGS_TESTING}" # Unroll loops SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "-funroll-loops" # GNU - "-unroll" # Intel + Fortran "-unroll" # Intel "/unroll" # Intel Windows + "-funroll-loops" # GNU "-Munroll" # Portland Group ) @@ -288,6 +288,12 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" Fortran "-fma" # Intel ) +# Generate fused multiply-add instructions +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" + Fortran "-qmkl=cluster" # Intel + Fortran "-qmkl" # Intel + Fortran "-mkl" # Old Intel + ) ##################### ### MATH FLAGS ### diff --git a/cmake/Modules/SetMKL.cmake b/cmake/Modules/SetMKL.cmake new file mode 100644 index 000000000..e58c9f51a --- /dev/null +++ b/cmake/Modules/SetMKL.cmake @@ -0,0 +1,14 @@ +# Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +# This file is part of Swiftest. +# Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +# Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# You should have received a copy of the GNU General Public License along with Swiftest. +# If not, see: https://www.gnu.org/licenses. + +# Find MKL if not already found +IF(NOT MKL_FOUND) + ENABLE_LANGUAGE(C) # Some libraries need a C compiler to find + FIND_PACKAGE(MKL REQUIRED) +ENDIF(NOT MKL_FOUND) diff --git a/cmake/Modules/SetParallelizationLibrary.cmake b/cmake/Modules/SetParallelizationLibrary.cmake index 03ab970e6..224806406 100644 --- a/cmake/Modules/SetParallelizationLibrary.cmake +++ b/cmake/Modules/SetParallelizationLibrary.cmake @@ -12,9 +12,7 @@ # When one is turned on, the other is turned off # If both are off, we explicitly disable them just in case -IF (USE_OPENMP AND USE_MPI) - MESSAGE (FATAL_ERROR "Cannot use both OpenMP and MPI") -ELSEIF (USE_OPENMP) +IF (USE_OPENMP) # Find OpenMP IF (NOT OpenMP_Fortran_FLAGS) FIND_PACKAGE (OpenMP_Fortran) @@ -23,20 +21,16 @@ ELSEIF (USE_OPENMP) ENDIF (NOT OpenMP_Fortran_FLAGS) ENDIF (NOT OpenMP_Fortran_FLAGS) # Turn of MPI - UNSET (MPI_FOUND CACHE) - UNSET (MPI_COMPILER CACHE) - UNSET (MPI_LIBRARY CACHE) -ELSEIF (USE_MPI) +ENDIF (USE_OPENMP) + +IF (USE_MPI) # Find MPI IF (NOT MPI_Fortran_FOUND) FIND_PACKAGE (MPI REQUIRED) ENDIF (NOT MPI_Fortran_FOUND) - # Turn off OpenMP - SET (OMP_NUM_PROCS 0 CACHE - STRING "Number of processors OpenMP may use" FORCE) - UNSET (OpenMP_C_FLAGS CACHE) - UNSET (GOMP_Fortran_LINK_FLAGS CACHE) -ELSE () +ENDIF (USE_MPI) + +IF (NOT USE_OPENMP AND NOT USE_MPI) # Turn off both OpenMP and MPI SET (OMP_NUM_PROCS 0 CACHE STRING "Number of processors OpenMP may use" FORCE) @@ -45,4 +39,4 @@ ELSE () UNSET (MPI_FOUND CACHE) UNSET (MPI_COMPILER CACHE) UNSET (MPI_LIBRARY CACHE) -ENDIF (USE_OPENMP AND USE_MPI) +ENDIF (NOT USE_OPENMP AND NOT USE_MPI) diff --git a/examples/.gitignore b/examples/.gitignore index 643407776..4aceee79f 100644 --- a/examples/.gitignore +++ b/examples/.gitignore @@ -1 +1,4 @@ */param.* +*/simdata/* +*/*/simdata/* + diff --git a/examples/Basic_Simulation/initial_conditions.ipynb b/examples/Basic_Simulation/initial_conditions.ipynb index cef068fa0..2bcd9cfe7 100644 --- a/examples/Basic_Simulation/initial_conditions.ipynb +++ b/examples/Basic_Simulation/initial_conditions.ipynb @@ -9,7 +9,8 @@ "source": [ "import swiftest\n", "import numpy as np\n", - "from numpy.random import default_rng" + "from numpy.random import default_rng\n", + "%env OMP_NUM_THREADS=4" ] }, { @@ -20,8 +21,7 @@ "outputs": [], "source": [ "# Initialize the simulation object as a variable\n", - "sim = swiftest.Simulation(tstart=0.0, tstop=10.0, dt=0.005, tstep_out=1.0, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8)\n", - "sim.get_parameter()" + "sim = swiftest.Simulation(tstart=0.0, tstop=1.0e6, dt=0.005, tstep_out=1.0e3, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8)" ] }, { @@ -119,6 +119,17 @@ "sim.add_body(name_tp, a_tp, e_tp, inc_tp, capom_tp, omega_tp, capm_tp)" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "135bd3f8-0a56-4d62-b728-f150debc1a76", + "metadata": {}, + "outputs": [], + "source": [ + "# Display the run configuration parameters\n", + "sim.get_parameter()" + ] + }, { "cell_type": "code", "execution_count": null, diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index 78900ce1f..4accb67e7 100644 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -23,8 +23,7 @@ from numpy.random import default_rng # Initialize the simulation object as a variable -sim = swiftest.Simulation(tstart=0.0, tstop=10.0, dt=0.005, tstep_out=1.0, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8) -sim.get_parameter() +sim = swiftest.Simulation(tstart=0.0, tstop=1.0e3, dt=0.005, tstep_out=1.0e0, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8) # Add the modern planets and the Sun using the JPL Horizons Database sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Pluto"]) @@ -64,6 +63,8 @@ capm_tp = default_rng().uniform(0.0, 360.0, ntp) sim.add_body(name_tp, a_tp, e_tp, inc_tp, capom_tp, omega_tp, capm_tp) +# Display the run configuration parameters +sim.get_parameter() # Run the simulation sim.run() diff --git a/examples/Basic_Simulation/run_from_file.ipynb b/examples/Basic_Simulation/run_from_file.ipynb new file mode 100644 index 000000000..d2f4fdeb1 --- /dev/null +++ b/examples/Basic_Simulation/run_from_file.ipynb @@ -0,0 +1,119 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "aa57e957-f141-4373-9a62-a91845203aa3", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "env: OMP_NUM_THREADS=8\n" + ] + } + ], + "source": [ + "import swiftest\n", + "import xarray as xr\n", + "import numpy as np\n", + "import os\n", + "%env OMP_NUM_THREADS=8" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "a9c020aa-0a87-4fa9-9b11-62213edb370c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/simdata/param.in\n", + "\n", + "Creating Dataset from NetCDF file\n", + "Successfully converted 1 output frames.\n" + ] + } + ], + "source": [ + "sim = swiftest.Simulation(read_param=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "209523b6-d7a8-46f0-8687-54f199015c2d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Writing parameter inputs to file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/simdata/param.in\n", + "Running a Swiftest symba run from tstart=0.0 y to tstop=1000.0 y\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "ca7e9529651b46209bc86174955f7a01", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Time: 0.0 / 1000.0 y : 0%| , npl: 14 ntp: 10 nplm: 13" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/daminton/git_debug/swiftest/python/swiftest/swiftest/simulation_class.py:465: UserWarning: Error executing main swiftest_driver program\n", + " self._run_swiftest_driver()\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Creating Dataset from NetCDF file\n", + "Successfully converted 19 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .data\n" + ] + } + ], + "source": [ + "sim.run()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (My debug_env Kernel)", + "language": "python", + "name": "debug_env" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/Basic_Simulation/run_from_file.py b/examples/Basic_Simulation/run_from_file.py index 9c477410a..146286814 100644 --- a/examples/Basic_Simulation/run_from_file.py +++ b/examples/Basic_Simulation/run_from_file.py @@ -1,3 +1,3 @@ import swiftest -sim = swiftest.Simulation() -sim.run(tstop=20.0) \ No newline at end of file +sim = swiftest.Simulation(read_param=True) +sim.run() \ No newline at end of file diff --git a/examples/Basic_Simulation/run_simulation.ipynb b/examples/Basic_Simulation/run_simulation.ipynb index ea78c2690..b35fc9ca6 100644 --- a/examples/Basic_Simulation/run_simulation.ipynb +++ b/examples/Basic_Simulation/run_simulation.ipynb @@ -21,7 +21,7 @@ "metadata": {}, "outputs": [], "source": [ - "sim = swiftest.Simulation()" + "sim = swiftest.Simulation(read_param=True)" ] }, { @@ -31,7 +31,7 @@ "metadata": {}, "outputs": [], "source": [ - "sim.run(tstop=10.0)" + "sim.run()" ] }, { diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 3c23a510a..ab666fb00 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -23,28 +23,67 @@ Returns ------- fragmentation.mp4 : mp4 movie file - Movide of a fragmentation event. + Movie of a fragmentation event. """ import swiftest import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation -from matplotlib.animation import FuncAnimation +from pathlib import Path -# Change this to be the parameter input file correlated with the run that you -# wish to test. Swiftest will pull in the corresponding out.nc file automatically. -param_file = "param.hitandrun.in" +print("Select a fragmentation movie to generate.") +print("1. Head-on disruption") +print("2. Off-axis supercatastrophic") +print("3. Hit and run") +print("4. All of the above") +user_selection = int(input("? ")) -# Change this to an appropriate title and filename to appear on the movie. -movie_title = "Hit and Run" -movie_filename = "hitandrun.mp4" +available_movie_styles = ["disruption_headon", "supercatastrophic_off_axis", "hitandrun"] +movie_title_list = ["Head-on Disrutption", "Off-axis Supercatastrophic", "Hit and Run"] +movie_titles = dict(zip(available_movie_styles, movie_title_list)) -# Pull in the Swiftest output data from the parameter file and store it as a Xarray dataset. -ds = swiftest.Simulation(param_file=param_file).ds +# These initial conditions were generated by trial and error +pos_vectors = {"disruption_headon" : [np.array([1.0, -1.807993e-05, 0.0]), + np.array([1.0, 1.807993e-05 ,0.0])], + "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05, 0.0]), + np.array([1.0, 4.2e-05, 0.0])], + "hitandrun" : [np.array([1.0, -4.2e-05, 0.0]), + np.array([1.0, 4.2e-05, 0.0])] + } -# Calculate the number of frames in the dataset. -nframes = int(ds['time'].size) +vel_vectors = {"disruption_headon" : [np.array([-2.562596e-04, 6.280005, 0.0]), + np.array([-2.562596e-04, -6.280005, 0.0])], + "supercatastrophic_off_axis": [np.array([0.0, 6.28, 0.0]), + np.array([1.0, -6.28, 0.0])], + "hitandrun" : [np.array([0.0, 6.28, 0.0]), + np.array([-1.5, -6.28, 0.0])] + } + +rot_vectors = {"disruption_headon" : [np.array([0.0, 0.0, 0.0]), + np.array([0.0, 0.0, 0.0])], + "supercatastrophic_off_axis": [np.array([0.0, 0.0, -6.0e4]), + np.array([0.0, 0.0, 1.0e5])], + "hitandrun" : [np.array([0.0, 0.0, 6.0e4]), + np.array([0.0, 0.0, 1.0e5])] + } + +body_Gmass = {"disruption_headon" : [1e-7, 1e-10], + "supercatastrophic_off_axis": [1e-7, 1e-8], + "hitandrun" : [1e-7, 7e-10] + } + +density = 3000 * swiftest.AU2M**3 / swiftest.MSun +GU = swiftest.GMSun * swiftest.YR2S**2 / swiftest.AU2M**3 +body_radius = body_Gmass.copy() +for k,v in body_Gmass.items(): + body_radius[k] = [((Gmass/GU)/(4./3.*np.pi*density))**(1./3.) for Gmass in v] + +if user_selection > 0 and user_selection < 4: + movie_styles = [available_movie_styles[user_selection-1]] +else: + print("Generating all movie styles") + movie_styles = available_movie_styles.copy() # Define a function to calculate the center of mass of the system. def center(xhx, xhy, xhz, Gmass): @@ -53,24 +92,19 @@ def center(xhx, xhy, xhz, Gmass): z_com = np.sum(Gmass * xhz) / np.sum(Gmass) return x_com, y_com, z_com -# Calculate the distance along the y-axis between the colliding bodies at the start of the simulation. -# This will be used to scale the axis limits on the movie. -scale_frame = abs(ds['xhy'].isel(time=0).isel(id=1).values) + abs(ds['xhy'].isel(time=0).isel(id=2).values) +def animate(i,ds,movie_title): -# Set up the figure and the animation. -fig, ax = plt.subplots(figsize=(4,4)) -def animate(i): # Calculate the position and mass of all bodies in the system at time i and store as a numpy array. - xhx = ds['xhx'].isel(time=i).dropna(dim='id').values - xhy = ds['xhy'].isel(time=i).dropna(dim='id').values - xhz = ds['xhx'].isel(time=i).dropna(dim='id').values - Gmass = ds['Gmass'].isel(time=i).dropna(dim='id').values[1:] # Drop the Sun from the numpy array. + xhx = ds['xhx'].isel(time=i).dropna(dim='name').values + xhy = ds['xhy'].isel(time=i).dropna(dim='name').values + xhz = ds['xhx'].isel(time=i).dropna(dim='name').values + Gmass = ds['Gmass'].isel(time=i).dropna(dim='name').values[1:] # Drop the Sun from the numpy array. - # Calculate the center of mass of the system at time i. While the center of mass relative to the + # Calculate the center of mass of the system at time i. While the center of mass relative to the # colliding bodies does not change, the center of mass of the collision will move as the bodies # orbit the system center of mass. x_com, y_com, z_com = center(xhx, xhy, xhz, Gmass) - + # Create the figure and plot the bodies as points. fig.clear() ax = fig.add_subplot(111) @@ -82,11 +116,36 @@ def animate(i): ax.grid(False) ax.set_xticks([]) ax.set_yticks([]) - - ax.scatter(xhx, xhy, s = (5000000000 * Gmass)) - + + ax.scatter(xhx, xhy, s=(5000000000 * Gmass)) + plt.tight_layout() -# Generate the movie. -ani = animation.FuncAnimation(fig, animate, interval=1, frames=nframes, repeat=False) -ani.save(movie_filename, fps=60, dpi=300, extra_args=['-vcodec', 'libx264']) \ No newline at end of file +for style in movie_styles: + param_file = Path(style) / "param.in" + + movie_filename = f"{style}.mp4" + + # Pull in the Swiftest output data from the parameter file and store it as a Xarray dataset. + sim = swiftest.Simulation(param_file=param_file, rotation=True, init_cond_format = "XV", compute_conservation_values=True) + sim.add_solar_system_body("Sun") + sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], xh=pos_vectors[style], vh=vel_vectors[style], rot=rot_vectors[style]) + + # Set fragmentation parameters + minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body + gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades + sim.set_parameter(fragmentation = True, gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass) + sim.run(dt=1e-8, tstop=1.e-5) + + # Calculate the number of frames in the dataset. + nframes = int(sim.data['time'].size) + + # Calculate the distance along the y-axis between the colliding bodies at the start of the simulation. + # This will be used to scale the axis limits on the movie. + scale_frame = abs(sim.data['xhy'].isel(time=0, name=1).values) + abs(sim.data['xhy'].isel(time=0, name=2).values) + + # Set up the figure and the animation. + fig, ax = plt.subplots(figsize=(4,4)) + # Generate the movie. + ani = animation.FuncAnimation(fig, animate, fargs=(sim.data, movie_titles[style]), interval=1, frames=nframes, repeat=False) + ani.save(movie_filename, fps=60, dpi=300, extra_args=['-vcodec', 'libx264']) \ No newline at end of file diff --git a/examples/Fragmentation/cb.in b/examples/Fragmentation/cb.in deleted file mode 100644 index 8766a22ae..000000000 --- a/examples/Fragmentation/cb.in +++ /dev/null @@ -1,16 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -0 -39.47841760435743 ! G*Mass -0.005 ! Radius -0.0 ! J2 -0.0 ! J4 -0.4 0.4 0.4 !Ip -0.0 0.0 0.0 !rot !11.2093063 -38.75937204 82.25088158 ! rot (radian / year) \ No newline at end of file diff --git a/examples/Fragmentation/disruption_headon.in b/examples/Fragmentation/disruption_headon.in deleted file mode 100644 index 1f4b208f1..000000000 --- a/examples/Fragmentation/disruption_headon.in +++ /dev/null @@ -1,22 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -2 -1 1e-07 0.0009 -7e-06 -1.0 -1.807993e-05 0.0 --2.562596e-04 6.280005 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 0.0 !rot -2 7e-10 0.0004 -3.25e-06 -1.0 1.807993e-05 0.0 --2.562596e-04 -6.280005 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 0.0 !rot diff --git a/examples/Fragmentation/hitandrun.in b/examples/Fragmentation/hitandrun.in deleted file mode 100644 index 285fc63a2..000000000 --- a/examples/Fragmentation/hitandrun.in +++ /dev/null @@ -1,22 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -2 -1 1e-07 0.0009 -7e-06 -1.0 -4.20E-05 0.0 -0.00 6.28 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 6.0e4 !rot -2 7e-10 0.0004 -3.25e-06 -1.0 4.20E-05 0.0 --1.50 -6.28 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 1.0e5 !rot diff --git a/examples/Fragmentation/param.disruption_headon.in b/examples/Fragmentation/param.disruption_headon.in deleted file mode 100644 index 0fd657831..000000000 --- a/examples/Fragmentation/param.disruption_headon.in +++ /dev/null @@ -1,51 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -T0 0.0e0 -TSTOP 0.00001 -DT 0.00000001 -ISTEP_OUT 1 -ISTEP_DUMP 1 -OUT_FORM XVEL -OUT_TYPE NETCDF_DOUBLE -OUT_STAT REPLACE -IN_FORM XV -IN_TYPE ASCII -NC_IN -1.0 -PL_IN disruption_headon.in -TP_IN tp.in -CB_IN cb.in -BIN_OUT disruption_headon.nc -CHK_QMIN -1.0 -CHK_RMIN 0.005 -CHK_RMAX 1e6 -CHK_EJECT -1.0 -CHK_QMIN_COORD -1.0 -CHK_QMIN_RANGE -1.0 -1.0 -MU2KG 1.98908e30 -TU2S 3.1556925e7 -DU2M 1.49598e11 -EXTRA_FORCE no -PARTICLE_OUT -1.0 -BIG_DISCARD no -CHK_CLOSE yes -GR NO -INTERACTION_LOOPS TRIANGULAR -ENCOUNTER_CHECK TRIANGULAR -RHILL_PRESENT yes -FRAGMENTATION yes -ROTATION yes -ENERGY yes -ENERGY_OUT -1.0 -ENC_OUT -1.0 -GMTINY 1.0e-11 -MIN_GMFRAG 1.0e-11 -TIDES NO -YORP NO -YARKOVSKY NO diff --git a/examples/Fragmentation/param.hitandrun.in b/examples/Fragmentation/param.hitandrun.in deleted file mode 100644 index 1bd02166c..000000000 --- a/examples/Fragmentation/param.hitandrun.in +++ /dev/null @@ -1,51 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -T0 0.0e0 -TSTOP 0.00001 -DT 0.00000001 -ISTEP_OUT 1 -ISTEP_DUMP 1 -OUT_FORM XVEL -OUT_TYPE NETCDF_DOUBLE -OUT_STAT REPLACE -IN_FORM XV -IN_TYPE ASCII -NC_IN -1.0 -PL_IN hitandrun.in -TP_IN tp.in -CB_IN cb.in -BIN_OUT hitandrun.nc -CHK_QMIN -1.0 -CHK_RMIN 0.005 -CHK_RMAX 1e6 -CHK_EJECT -1.0 -CHK_QMIN_COORD -1.0 -CHK_QMIN_RANGE -1.0 -1.0 -MU2KG 1.98908e30 -TU2S 3.1556925e7 -DU2M 1.49598e11 -EXTRA_FORCE no -PARTICLE_OUT -1.0 -BIG_DISCARD no -CHK_CLOSE yes -GR NO -INTERACTION_LOOPS TRIANGULAR -ENCOUNTER_CHECK TRIANGULAR -RHILL_PRESENT yes -FRAGMENTATION yes -ROTATION yes -ENERGY yes -ENERGY_OUT -1.0 -ENC_OUT -1.0 -GMTINY 1.0e-11 -MIN_GMFRAG 1.0e-11 -TIDES NO -YORP NO -YARKOVSKY NO diff --git a/examples/Fragmentation/param.supercatastrophic_off_axis.in b/examples/Fragmentation/param.supercatastrophic_off_axis.in deleted file mode 100644 index 08b5dd71d..000000000 --- a/examples/Fragmentation/param.supercatastrophic_off_axis.in +++ /dev/null @@ -1,51 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -T0 0.0e0 -TSTOP 0.00001 -DT 0.00000001 -ISTEP_OUT 1 -ISTEP_DUMP 1 -OUT_FORM XVEL -OUT_TYPE NETCDF_DOUBLE -OUT_STAT REPLACE -IN_FORM XV -IN_TYPE ASCII -NC_IN -1.0 -PL_IN supercatastrophic_off_axis.in -TP_IN tp.in -CB_IN cb.in -BIN_OUT supercatastrophic_off_axis.nc -CHK_QMIN -1.0 -CHK_RMIN 0.005 -CHK_RMAX 1e6 -CHK_EJECT -1.0 -CHK_QMIN_COORD -1.0 -CHK_QMIN_RANGE -1.0 -1.0 -MU2KG 1.98908e30 -TU2S 3.1556925e7 -DU2M 1.49598e11 -EXTRA_FORCE no -PARTICLE_OUT -1.0 -BIG_DISCARD no -CHK_CLOSE yes -GR NO -INTERACTION_LOOPS TRIANGULAR -ENCOUNTER_CHECK TRIANGULAR -RHILL_PRESENT yes -FRAGMENTATION yes -ROTATION yes -ENERGY yes -ENERGY_OUT -1.0 -ENC_OUT -1.0 -GMTINY 1.0e-11 -MIN_GMFRAG 1.0e-11 -TIDES NO -YORP NO -YARKOVSKY NO diff --git a/examples/Fragmentation/supercatastrophic_off_axis.in b/examples/Fragmentation/supercatastrophic_off_axis.in deleted file mode 100644 index 03315636d..000000000 --- a/examples/Fragmentation/supercatastrophic_off_axis.in +++ /dev/null @@ -1,22 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -2 -1 1e-07 0.0009 -7e-06 -1.0 -4.2e-05 0.0 -0.00 6.28 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 -6.0e4 !rot -2 1e-08 0.0004 -3.25e-06 -1.0 4.2e-05 0.0 -1.00 -6.28 0.0 -0.4 0.4 0.4 !Ip -0.0 0.0 1.0e5 !rot diff --git a/examples/Fragmentation/tp.in b/examples/Fragmentation/tp.in deleted file mode 100644 index 3c6f40630..000000000 --- a/examples/Fragmentation/tp.in +++ /dev/null @@ -1,10 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -0 diff --git a/examples/helio_gr_test/grsim/param.gr.in b/examples/helio_gr_test/grsim/param.gr.in new file mode 100644 index 000000000..0616db203 --- /dev/null +++ b/examples/helio_gr_test/grsim/param.gr.in @@ -0,0 +1,35 @@ +! VERSION Swiftest input file +T0 0.0 +TSTART 0.0 +TSTOP 1000.0 +DT 0.005 +ISTEP_OUT 2000 +ISTEP_DUMP 2000 +NC_IN init_cond.nc +IN_TYPE NETCDF_DOUBLE +IN_FORM EL +BIN_OUT bin.gr.nc +OUT_FORM XVEL +OUT_TYPE NETCDF_DOUBLE +OUT_STAT REPLACE +CHK_QMIN 0.004650467260962157 +CHK_RMIN 0.004650467260962157 +CHK_RMAX 10000.0 +CHK_EJECT 10000.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 10000.0 +MU2KG 1.988409870698051e+30 +TU2S 31557600.0 +DU2M 149597870700.0 +FRAGMENTATION NO +RESTART NO +CHK_CLOSE YES +GR YES +ROTATION NO +ENERGY NO +EXTRA_FORCE NO +BIG_DISCARD NO +RHILL_PRESENT NO +INTERACTION_LOOPS TRIANGULAR +ENCOUNTER_CHECK TRIANGULAR +TIDES NO diff --git a/examples/helio_gr_test/init_cond.py b/examples/helio_gr_test/init_cond.py deleted file mode 100755 index 9e89dc358..000000000 --- a/examples/helio_gr_test/init_cond.py +++ /dev/null @@ -1,48 +0,0 @@ -#!/usr/bin/env python3 -import swiftest - -sim = swiftest.Simulation() -sim.param['PL_IN'] = "pl.swiftest.in" -sim.param['TP_IN'] = "tp.swiftest.in" -sim.param['CB_IN'] = "cb.swiftest.in" -sim.param['BIN_OUT'] = "bin.swiftest.nc" - -sim.param['MU2KG'] = swiftest.MSun -sim.param['TU2S'] = swiftest.YR2S -sim.param['DU2M'] = swiftest.AU2M -sim.param['T0'] = 0.0 -sim.param['DT'] = 0.25 * swiftest.JD2S / swiftest.YR2S -sim.param['TSTOP'] = 1000.0 -sim.param['ISTEP_OUT'] = 1461 -sim.param['ISTEP_DUMP'] = 1461 -sim.param['CHK_QMIN_COORD'] = "HELIO" -sim.param['CHK_QMIN'] = swiftest.RSun / swiftest.AU2M -sim.param['CHK_QMIN_RANGE'] = f"{swiftest.RSun / swiftest.AU2M} 1000.0" -sim.param['CHK_RMIN'] = swiftest.RSun / swiftest.AU2M -sim.param['CHK_RMAX'] = 1000.0 -sim.param['CHK_EJECT'] = 1000.0 -sim.param['OUT_STAT'] = "UNKNOWN" -sim.param['IN_FORM'] = "EL" -sim.param['OUT_FORM'] = "XVEL" -sim.param['OUT_TYPE'] = "NETCDF_DOUBLE" -sim.param['RHILL_PRESENT'] = "YES" -sim.param['GR'] = 'YES' - -bodyid = { - "Sun": 0, - "Mercury": 1, - "Venus": 2, - "Earth": 3, - "Mars": 4, - "Jupiter": 5, - "Saturn": 6, - "Uranus": 7, - "Neptune": 8, -} - -for name, id in bodyid.items(): - sim.add(name, idval=id, date="2027-04-30") - -sim.save("param.swiftest.in") - - diff --git a/examples/helio_gr_test/nogrsim/param.nogr.in b/examples/helio_gr_test/nogrsim/param.nogr.in new file mode 100644 index 000000000..9e2ab0b22 --- /dev/null +++ b/examples/helio_gr_test/nogrsim/param.nogr.in @@ -0,0 +1,35 @@ +! VERSION Swiftest input file +T0 0.0 +TSTART 0.0 +TSTOP 1000.0 +DT 0.005 +ISTEP_OUT 2000 +ISTEP_DUMP 2000 +NC_IN init_cond.nc +IN_TYPE NETCDF_DOUBLE +IN_FORM EL +BIN_OUT bin.nogr.nc +OUT_FORM XVEL +OUT_TYPE NETCDF_DOUBLE +OUT_STAT REPLACE +CHK_QMIN 0.004650467260962157 +CHK_RMIN 0.004650467260962157 +CHK_RMAX 10000.0 +CHK_EJECT 10000.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 10000.0 +MU2KG 1.988409870698051e+30 +TU2S 31557600.0 +DU2M 149597870700.0 +FRAGMENTATION NO +RESTART NO +CHK_CLOSE YES +GR NO +ROTATION NO +ENERGY NO +EXTRA_FORCE NO +BIG_DISCARD NO +RHILL_PRESENT NO +INTERACTION_LOOPS TRIANGULAR +ENCOUNTER_CHECK TRIANGULAR +TIDES NO diff --git a/examples/helio_gr_test/swiftest_relativity.ipynb b/examples/helio_gr_test/swiftest_relativity.ipynb index 6946ef658..4b6d13106 100644 --- a/examples/helio_gr_test/swiftest_relativity.ipynb +++ b/examples/helio_gr_test/swiftest_relativity.ipynb @@ -19,7 +19,7 @@ "metadata": {}, "outputs": [], "source": [ - "sim_gr = swiftest.Simulation(param_file=\"param.gr.in\", output_file_name=\"bin.gr.nc\")\n", + "sim_gr = swiftest.Simulation(param_file=\"grsim/param.gr.in\", output_file_name=\"bin.gr.nc\")\n", "sim_gr.add_solar_system_body([\"Sun\",\"Mercury\",\"Venus\",\"Earth\",\"Mars\",\"Jupiter\",\"Saturn\",\"Uranus\",\"Neptune\"])" ] }, @@ -29,7 +29,7 @@ "metadata": {}, "outputs": [], "source": [ - "sim_nogr = swiftest.Simulation(param_file=\"param.nogr.in\", output_file_name=\"bin.nogr.nc\")\n", + "sim_nogr = swiftest.Simulation(param_file=\"nogrsim/param.nogr.in\", output_file_name=\"bin.nogr.nc\")\n", "sim_nogr.add_solar_system_body([\"Sun\",\"Mercury\",\"Venus\",\"Earth\",\"Mars\",\"Jupiter\",\"Saturn\",\"Uranus\",\"Neptune\"])" ] }, @@ -39,7 +39,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "tstep_out = 10.0\n", "sim_gr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"helio\",general_relativity=True)" ] @@ -50,7 +49,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "sim_nogr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"helio\",general_relativity=False)" ] }, @@ -138,13 +136,6 @@ "print(f'Obs - Swiftest No GR : {np.mean(dvarpi_obs - dvarpi_nogr)}')" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "code", "execution_count": null, diff --git a/python/swiftest/swiftest/constants.py b/python/swiftest/swiftest/constants.py index fe3253c86..2d3f89f7c 100644 --- a/python/swiftest/swiftest/constants.py +++ b/python/swiftest/swiftest/constants.py @@ -13,7 +13,7 @@ import astropy.constants as const # Constants in SI units -GC = const.G.value +GC = const.G.value[()] AU2M = const.au.value GMSun = const.GM_sun.value MSun = const.M_sun.value diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index ea28bcb8c..1d43fb599 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -329,7 +329,6 @@ def vec2xr(param: Dict, infodims = ['id', 'vec'] # The central body is always given id 0 - if GMpl is not None: icb = (~np.isnan(GMpl)) & (idvals == 0) ipl = (~np.isnan(GMpl)) & (idvals != 0) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 4fd797b76..5d1f72ad6 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -109,6 +109,8 @@ def str2bool(input_str): {True, False} """ + if type(input_str) is bool: + return input_str valid_true = ["YES", "Y", "T", "TRUE", ".TRUE."] valid_false = ["NO", "N", "F", "FALSE", ".FALSE."] if input_str.upper() in valid_true: @@ -175,12 +177,13 @@ def read_swiftest_param(param_file_name, param, verbose=True): param[uc] = param[uc].upper() for i in int_param: - if i in param and type(i) != int: - param[i] = int(param[i]) + if i in param and type(param[i]) != int: + param[i] = int(float(param[i])) for f in float_param: - if f in param and type(f) is str: + if f in param and type(param[f]) is str: param[f] = real2float(param[f]) + for b in bool_param: if b in param: param[b] = str2bool(param[b]) @@ -427,7 +430,7 @@ def write_labeled_param(param, param_file_name): 'TU2S', 'DU2M', 'GMTINY', - 'FRAGMENTATION' + 'FRAGMENTATION', 'MIN_GMFRAG', 'RESTART'] ptmp = param.copy() @@ -1067,7 +1070,7 @@ def select_active_from_frame(ds, param, framenum=-1): return frame -def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,framenum=-1): +def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,framenum=-1,verbose=True): """ Writes a set of Swiftest input files from a single frame of a Swiftest xarray dataset @@ -1100,7 +1103,8 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram if infile_name is None: infile_name = param['NC_IN'] frame = unclean_string_values(frame) - print(f"Writing initial conditions to file {infile_name}") + if verbose: + print(f"Writing initial conditions to file {infile_name}") frame.to_netcdf(path=infile_name) return frame diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index c58d68db8..7669c774f 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -15,7 +15,9 @@ from swiftest import tool from swiftest import constants from swiftest import __file__ as _pyfile +import json import os +from pathlib import Path import datetime import xarray as xr import numpy as np @@ -24,10 +26,12 @@ import subprocess import shlex import warnings +from tqdm.auto import tqdm from typing import ( Literal, Dict, List, + Tuple, Any ) @@ -37,7 +41,7 @@ class Simulation: This is a class that defines the basic Swift/Swifter/Swiftest simulation object """ - def __init__(self,read_param: bool = True, **kwargs: Any): + def __init__(self,read_param: bool = False, **kwargs: Any): """ Parameters @@ -55,8 +59,8 @@ def __init__(self,read_param: bool = True, **kwargs: Any): 1. Arguments to Simulation() 2. The parameter input file given by `param_file` under the following conditions: - `read_param` is set to True (default behavior). - - The file given by `param_file` exists. The default file is `param.in` located in the current working - directory, which can be changed by passing `param_file` as an argument. + - The file given by `param_file` exists. The default file is `param.in` located in the `simdata` directory + inside the current working directory, which can be changed by passing `param_file` as an argument. - The argument has an equivalent parameter or set of parameters in the parameter input file. 3. Default values (see below) @@ -274,26 +278,42 @@ def __init__(self,read_param: bool = True, **kwargs: Any): Parameter input file equivalent: None """ - # Width of the column in the printed name of the parameter in parameter getters - self._getter_column_width = '32' + # Configuration parameters will be stored in a json file alongside the Python source scripts. + self._config_file = Path(_pyfile).parent / "swiftest_configuration.json" + config_exists = self._config_file.exists() + if config_exists: + try: + with open(self._config_file, 'r') as f: + self._swiftest_configuration = json.load(f) + except: + config_exists = False + if not config_exists: + self._swiftest_configuration = {"shell" : str(Path(os.environ['SHELL']).name), + "shell_full" : str(Path(os.environ['SHELL'])), + "getter_column_width" : '32'} + self._swiftest_configuration['startup_script'] = str(Path.home() / f".{str(self._swiftest_configuration['shell'])}rc") + config_json = json.dumps(self._swiftest_configuration, indent=4) + with open(self._config_file, 'w') as f: + f.write(config_json) + + self._getter_column_width = self._swiftest_configuration['getter_column_width'] + self._shell = Path(self._swiftest_configuration['shell']) + self._shell_full = Path(self._swiftest_configuration['shell_full']) self.param = {} self.data = xr.Dataset() + # Set the location of the parameter input file, choosing the default if it isn't specified. + param_file = kwargs.pop("param_file",Path.cwd() / "simdata" / "param.in") + self.verbose = kwargs.pop("verbose",True) + # Parameters are set in reverse priority order. First the defaults, then values from a pre-existing input file, # then using the arguments passed via **kwargs. #-------------------------- # Lowest Priority: Defaults #-------------------------- - # Quietly set all parameters to their defaults. - self.verbose = kwargs.pop("verbose",True) - self.set_parameter(verbose=False) - - # Set the location of the parameter input file - param_file = kwargs.pop("param_file",self.param_file) - read_param = kwargs.pop("read_param",False) self.set_parameter(verbose=False,param_file=param_file) #----------------------------------------------------------------- @@ -308,6 +328,8 @@ def __init__(self,read_param: bool = True, **kwargs: Any): # We will add the parameter file to the kwarg list. This will keep the set_parameter method from # overriding everything with defaults when there are no arguments passed to Simulation() kwargs['param_file'] = self.param_file + param_file_found = True + else: param_file_found = False @@ -318,7 +340,7 @@ def __init__(self,read_param: bool = True, **kwargs: Any): # Let the user know that there was a problem reading an old parameter file and we're going to create a new one if read_param and not param_file_found: - warnings.warn(f"{self.param_file} not found. Creating a new file using default values for parameters not passed to Simulation().") + warnings.warn(f"{self.param_file} not found. Creating a new file using default values for parameters not passed to Simulation().",stacklevel=2) self.write_param() # Read in an old simulation file if requested @@ -328,9 +350,90 @@ def __init__(self,read_param: bool = True, **kwargs: Any): if os.path.exists(binpath): self.bin2xr() else: - warnings.warn(f"BIN_OUT file {binpath} not found.") + warnings.warn(f"BIN_OUT file {binpath} not found.",stacklevel=2) return + def _run_swiftest_driver(self): + """ + Internal callable function that executes the swiftest_driver run + """ + + # Get current environment variables + + env = os.environ.copy() + driver_script = os.path.join(self.binary_path, "swiftest_driver.sh") + with open(driver_script, 'w') as f: + f.write(f"#{self._shell_full} -l\n") + f.write(f"source ~/.{self._shell}rc\n") + f.write(f"cd {self.sim_dir}\n") + f.write(f"{str(self.driver_executable)} {self.integrator} {str(self.param_file)} compact\n") + + cmd = f"{env['SHELL']} -l {driver_script}" + + def _type_scrub(output_data): + int_vars = ["ILOOP","NPL","NTP","NPLM"] + for k,v in output_data.items(): + if k in int_vars: + output_data[k] = int(v) + else: + output_data[k] = float(v) + return output_data + + process_output = False + noutput = int((self.param['TSTOP'] - self.param['T0']) / self.param['DT']) + iloop = int((self.param['TSTART'] - self.param['T0']) / self.param['DT']) + twidth = int(np.ceil(np.log10(self.param['TSTOP']/self.param['DT']))) + pre_message = f"Time: {self.param['TSTART']:.{twidth}e} / {self.param['TSTOP']:.{twidth}e} {self.TU_name} " + post_message = f"npl: {self.data['npl'].values[0]} ntp: {self.data['ntp'].values[0]}" + if "nplm" in self.data: + post_message += f" nplm: {self.data['nplm'].values[0]}" + if self.param['ENERGY']: + post_message += f" dL/L0: {0.0:.5e} dE/|E0|: {0.0:+.5e}" + post_message += f" Wall time / step: {0.0:.5e} s" + pbar = tqdm(total=noutput, desc=pre_message, postfix=post_message, bar_format='{l_bar}{bar}{postfix}') + try: + with subprocess.Popen(shlex.split(cmd), + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + env=env, + universal_newlines=True) as p: + + for line in p.stdout: + if "SWIFTEST STOP" in line: + process_output = False + + if process_output: + kvstream=line.replace('\n','').strip().split(';') # Removes the newline character, + output_data = _type_scrub({kv.split()[0]: kv.split()[1] for kv in kvstream[:-1]}) + pre_message = f"Time: {output_data['T']:.{twidth}e} / {self.param['TSTOP']:.{twidth}e} {self.TU_name}" + post_message = f" npl: {output_data['NPL']} ntp: {output_data['NTP']}" + if "NPLM" in output_data: + post_message += f" nplm: {output_data['NPLM']}" + if "LTOTERR" in output_data: + post_message += f" dL/L0: {output_data['LTOTERR']:.5e}" + if "ETOTERR" in output_data: + post_message += f" dE/|E0|: {output_data['ETOTERR']:+.5e}" + post_message += f" Wall time / step: {output_data['WTPS']:.5e} s" + interval = output_data['ILOOP'] - iloop + if interval > 0: + pbar.update(interval) + pbar.set_description_str(pre_message) + pbar.set_postfix_str(post_message) + iloop = output_data['ILOOP'] + + if "SWIFTEST START" in line: + process_output = True + + res = p.communicate() + if p.returncode != 0: + for line in res[1]: + print(line, end='') + warnings.warn("Failure in swiftest_driver", stacklevel=2) + except: + warnings.warn(f"Error executing main swiftest_driver program", stacklevel=2) + + pbar.close() + return def run(self,**kwargs): """ @@ -354,43 +457,18 @@ def run(self,**kwargs): self.write_param() if self.codename != "Swiftest": - warnings.warn(f"Running an integration is not yet supported for {self.codename}") + warnings.warn(f"Running an integration is not yet supported for {self.codename}",stacklevel=2) return if self.driver_executable is None: - warnings.warn("Path to swiftest_driver has not been set!") - warnings.warn(f"Make sure swiftest_driver is compiled and the executable is in {self.binary_path}") + msg = "Path to swiftest_driver has not been set!" + msg += f"\nMake sure swiftest_driver is compiled and the executable is in {str(self.binary_path)}" + warnings.warn(msg,stacklevel=2) return print(f"Running a {self.codename} {self.integrator} run from tstart={self.param['TSTART']} {self.TU_name} to tstop={self.param['TSTOP']} {self.TU_name}") - # Get current environment variables - env = os.environ.copy() - driver_script = os.path.join(self.binary_path,"swiftest_driver.sh") - shell = os.path.basename(env['SHELL']) - with open(driver_script,'w') as f: - f.write(f"#{env['SHELL']} -l {os.linesep}") - f.write(f"source ~/.{shell}rc {os.linesep}") - f.write(f"cd {self.sim_dir} {os.linesep}") - f.write(f"{self.driver_executable} {self.integrator} {self.param_file}") - - try: - cmd = f"{env['SHELL']} -l {driver_script}" - with subprocess.Popen(shlex.split(cmd), - stdout=subprocess.PIPE, - stderr=subprocess.PIPE, - env=env, - universal_newlines=True) as p: - for line in p.stdout: - print(line, end='') - res = p.communicate() - if p.returncode != 0: - for line in res[1]: - print(line, end='') - raise Exception ("Failure in swiftest_driver") - except: - warnings.warn(f"Error executing main swiftest_driver program") - return + self._run_swiftest_driver() # Read in new data self.bin2xr() @@ -512,7 +590,7 @@ def set_simulation_time(self, if tstop is not None: if tstop <= tstart: - warnings.warn("tstop must be greater than tstart.") + warnings.warn("tstop must be greater than tstart.",stacklevel=2) return {} if tstop is not None: @@ -525,8 +603,9 @@ def set_simulation_time(self, if dt is not None and tstop is not None: if dt > (tstop - tstart): - warnings.warn("dt must be smaller than tstop-tstart") - warnings.warn(f"Setting dt = {tstop - tstart} instead of {dt}") + msg = "dt must be smaller than tstop-tstart" + msg +=f"\nSetting dt = {tstop - tstart} instead of {dt}" + warnings.warn(msg,stacklevel=2) dt = tstop - tstart if dt is not None: @@ -535,7 +614,7 @@ def set_simulation_time(self, if istep_out is None and tstep_out is None: istep_out = self.param.pop("ISTEP_OUT", None) elif istep_out is not None and tstep_out is not None: - warnings.warn("istep_out and tstep_out cannot both be set") + warnings.warn("istep_out and tstep_out cannot both be set",stacklevel=2) return {} else: update_list.append("istep_out") @@ -643,17 +722,17 @@ def set_parameter(self, verbose: bool = True, **kwargs): default_arguments = { "codename" : "Swiftest", "integrator": "symba", - "param_file": "param.in", "t0": 0.0, "tstart": 0.0, "tstop": None, "dt": None, - "istep_out": None, + "istep_out": 1, "tstep_out": None, - "istep_dump": None, + "istep_dump": 1, "init_cond_file_type": "NETCDF_DOUBLE", "init_cond_file_name": None, "init_cond_format": "EL", + "read_old_output_file": False, "output_file_type": "NETCDF_DOUBLE", "output_file_name": None, "output_format": "XVEL", @@ -669,13 +748,13 @@ def set_parameter(self, verbose: bool = True, **kwargs): "rmin": constants.RSun / constants.AU2M, "rmax": 10000.0, "qmin_coord": "HELIO", - "gmtiny": None, + "gmtiny": 0.0, "mtiny": None, "close_encounter_check": True, "general_relativity": True, "fragmentation": False, "minimum_fragment_mass": None, - "minimum_fragment_gmass": None, + "minimum_fragment_gmass": 0.0, "rotation": False, "compute_conservation_values": False, "extra_force": False, @@ -686,21 +765,32 @@ def set_parameter(self, verbose: bool = True, **kwargs): "ephemeris_date": "MBCL", "restart": False, } + param_file = kwargs.pop("param_file",None) + + # Extract the simulation directory and create it if it doesn't exist + if param_file is not None: + self.param_file = Path.cwd() / param_file + self.sim_dir = self.param_file.parent + if self.sim_dir.exists(): + if not self.sim_dir.is_dir(): + msg = f"Cannot create the {self.sim_dir} directory: File exists." + msg += "\nDelete the file or change the location of param_file" + warnings.warn(msg,stacklevel=2) + else: + self.sim_dir.mkdir(parents=True, exist_ok=False) # If no arguments (other than, possibly, verbose) are requested, use defaults if len(kwargs) == 0: kwargs = default_arguments + unrecognized = [k for k,v in kwargs.items() if k not in default_arguments] + if len(unrecognized) > 0: + for k in unrecognized: + warnings.warn(f'Unrecognized argument "{k}"',stacklevel=2) + # Add the verbose flag to the kwargs for passing down to the individual setters kwargs["verbose"] = verbose - param_file = kwargs.pop("param_file",None) - - if param_file is not None: - self.param_file = os.path.realpath(param_file) - self.sim_dir = os.path.dirname(self.param_file) - - # Setters returning parameter dictionary values param_dict = {} param_dict.update(self.set_unit_system(**kwargs)) @@ -783,7 +873,7 @@ def set_integrator(self, if codename is not None: valid_codename = ["Swiftest", "Swifter", "Swift"] if codename.title() not in valid_codename: - warnings.warn(f"{codename} is not a valid codename. Valid options are ",",".join(valid_codename)) + warnings.warn(f"{codename} is not a valid codename. Valid options are ",",".join(valid_codename),stacklevel=2) try: self.codename except: @@ -794,10 +884,10 @@ def set_integrator(self, self.param['! VERSION'] = f"{self.codename} input file" update_list.append("codename") if self.codename == "Swiftest": - self.binary_path = os.path.realpath(os.path.join(os.path.dirname(os.path.realpath(_pyfile)),os.pardir,os.pardir,os.pardir,"bin")) - self.driver_executable = os.path.join(self.binary_path,"swiftest_driver") - if not os.path.exists(self.driver_executable): - warnings.warn(f"Cannot find the Swiftest driver in {self.binary_path}") + self.binary_path = Path(_pyfile).parent.parent.parent.parent / "bin" + self.driver_executable = self.binary_path / "swiftest_driver" + if not self.driver_executable.exists(): + warnings.warn(f"Cannot find the Swiftest driver in {str(self.binary_path)}",stacklevel=2) self.driver_executable = None else: self.binary_path = "NOT IMPLEMENTED FOR THIS CODE" @@ -807,7 +897,7 @@ def set_integrator(self, if integrator is not None: valid_integrator = ["symba","rmvs","whm","helio"] if integrator.lower() not in valid_integrator: - warnings.warn(f"{integrator} is not a valid integrator. Valid options are ",",".join(valid_integrator)) + warnings.warn(f"{integrator} is not a valid integrator. Valid options are ",",".join(valid_integrator),stacklevel=2) try: self.integrator except: @@ -818,9 +908,9 @@ def set_integrator(self, if mtiny is not None or gmtiny is not None: if self.integrator != "symba": - warnings.warn("mtiny and gmtiny are only used by SyMBA.") + warnings.warn("mtiny and gmtiny are only used by SyMBA.",stacklevel=2) if mtiny is not None and gmtiny is not None: - warnings.warn("Only set mtiny or gmtiny, not both!") + warnings.warn("Only set mtiny or gmtiny, not both.",stacklevel=2) elif gmtiny is not None: self.param['GMTINY'] = gmtiny update_list.append("gmtiny") @@ -859,19 +949,19 @@ def get_integrator(self,arg_list: str | List[str] | None = None, verbose: bool | valid_instance_vars = {"codename": self.codename, "integrator": self.integrator, - "param_file": self.param_file, - "driver_executable": self.driver_executable} + "param_file": str(self.param_file), + "driver_executable": str(self.driver_executable)} try: self.integrator except: - warnings.warn(f"integrator is not set") + warnings.warn(f"integrator is not set",stacklevel=2) return {} try: self.codename except: - warnings.warn(f"codename is not set") + warnings.warn(f"codename is not set",stacklevel=2) return {} if verbose is None: @@ -1008,19 +1098,19 @@ def set_feature(self, if fragmentation is not None: if self.codename != "Swiftest" and self.integrator != "symba" and fragmentation: - warnings.warn("Fragmentation is only available on Swiftest SyMBA.") + warnings.warn("Fragmentation is only available on Swiftest SyMBA.",stacklevel=2) self.param['FRAGMENTATION'] = False else: self.param['FRAGMENTATION'] = fragmentation update_list.append("fragmentation") if fragmentation: if "MIN_GMFRAG" not in self.param and minimum_fragment_mass is None and minimum_fragment_gmass is None: - warnings.warn("Minimum fragment mass is not set. Set it using minimum_fragment_gmass or minimum_fragment_mass") + warnings.warn("Minimum fragment mass is not set. Set it using minimum_fragment_gmass or minimum_fragment_mass",stacklevel=2) else: update_list.append("minimum_fragment_gmass") if minimum_fragment_gmass is not None and minimum_fragment_mass is not None: - warnings.warn("Only set either minimum_fragment_mass or minimum_fragment_gmass, but not both!") + warnings.warn("Only set either minimum_fragment_mass or minimum_fragment_gmass, but not both!",stacklevel=2) if minimum_fragment_gmass is not None: self.param["MIN_GMFRAG"] = minimum_fragment_gmass @@ -1065,8 +1155,9 @@ def set_feature(self, if interaction_loops is not None: valid_vals = ["TRIANGULAR", "FLAT", "ADAPTIVE"] if interaction_loops not in valid_vals: - warnings.warn(f"{interaction_loops} is not a valid option for interaction loops.") - warnings.warn(f"Must be one of {valid_vals}") + msg = f"{interaction_loops} is not a valid option for interaction loops." + msg += f"\nMust be one of {valid_vals}" + warnings.warn(msg,stacklevel=2) if "INTERACTION_LOOPS" not in self.param: self.param["INTERACTION_LOOPS"] = valid_vals[0] else: @@ -1076,8 +1167,9 @@ def set_feature(self, if encounter_check_loops is not None: valid_vals = ["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] if encounter_check_loops not in valid_vals: - warnings.warn(f"{encounter_check_loops} is not a valid option for interaction loops.") - warnings.warn(f"Must be one of {valid_vals}") + msg = f"{encounter_check_loops} is not a valid option for interaction loops." + msg += f"\nMust be one of {valid_vals}" + warnings.warn(msg,stacklevel=2) if "ENCOUNTER_CHECK" not in self.param: self.param["ENCOUNTER_CHECK"] = valid_vals[0] else: @@ -1208,13 +1300,14 @@ def set_init_cond_files(self, return {} def ascii_file_input_error_msg(codename): - warnings.warn(f"in set_init_cond_files: init_cond_file_name must be a dictionary of the form: ") - warnings.warn('{') + msg = f"in set_init_cond_files: init_cond_file_name must be a dictionary of the form: " + msg += "\n {" if codename == "Swiftest": - warnings.warn('"CB" : *path to central body initial conditions file*,') - warnings.warn('"PL" : *path to massive body initial conditions file*,') - warnings.warn('"TP" : *path to test particle initial conditions file*') - warnings.warn('}') + msg += '\n"CB" : *path to central body initial conditions file*,' + msg += '\n"PL" : *path to massive body initial conditions file*,' + msg += '\n"TP" : *path to test particle initial conditions file*' + msg += '\n}' + warnings.warn(msg,stacklevel=2) return {} if init_cond_format is None: @@ -1234,21 +1327,21 @@ def ascii_file_input_error_msg(codename): else: init_cond_keys = ["PL", "TP"] if init_cond_file_type != "ASCII": - warnings.warn(f"{init_cond_file_type} is not supported by {self.codename}. Using ASCII instead") + warnings.warn(f"{init_cond_file_type} is not supported by {self.codename}. Using ASCII instead",stacklevel=2) init_cond_file_type = "ASCII" if init_cond_format != "XV": - warnings.warn(f"{init_cond_format} is not supported by {self.codename}. Using XV instead") + warnings.warn(f"{init_cond_format} is not supported by {self.codename}. Using XV instead",stacklevel=2) init_cond_format = "XV" valid_formats = {"EL", "XV"} if init_cond_format not in valid_formats: - warnings.warn(f"{init_cond_format} is not a valid input format") + warnings.warn(f"{init_cond_format} is not a valid input format",stacklevel=2) else: self.param['IN_FORM'] = init_cond_format valid_types = {"NETCDF_DOUBLE", "NETCDF_FLOAT", "ASCII"} if init_cond_file_type not in valid_types: - warnings.warn(f"{init_cond_file_type} is not a valid input type") + warnings.warn(f"{init_cond_file_type} is not a valid input type",stackevel=2) else: self.param['IN_TYPE'] = init_cond_file_type @@ -1275,7 +1368,7 @@ def ascii_file_input_error_msg(codename): elif type(init_cond_file_name) is dict: # Oops, accidentally passed a dictionary instead of the expected single string or path-like for NetCDF # input type. - warnings.warn(f"Only a single input file is used for NetCDF files") + warnings.warn(f"Only a single input file is used for NetCDF files",stacklevel=2) else: self.param["NC_IN"] = init_cond_file_name @@ -1416,7 +1509,7 @@ def set_output_files(self, if output_file_type is None: output_file_type = "NETCDF_DOUBLE" elif output_file_type not in ["NETCDF_DOUBLE", "NETCDF_FLOAT"]: - warnings.warn(f"{output_file_type} is not compatible with Swiftest. Setting to NETCDF_DOUBLE") + warnings.warn(f"{output_file_type} is not compatible with Swiftest. Setting to NETCDF_DOUBLE",stacklevel=2) output_file_type = "NETCDF_DOUBLE" elif self.codename == "Swifter": if output_file_type is None: @@ -1424,7 +1517,7 @@ def set_output_files(self, if output_file_type is None: output_file_type = "REAL8" elif output_file_type not in ["REAL4", "REAL8", "XDR4", "XDR8"]: - warnings.warn(f"{output_file_type} is not compatible with Swifter. Setting to REAL8") + warnings.warn(f"{output_file_type} is not compatible with Swifter. Setting to REAL8",stacklevel=2) output_file_type = "REAL8" elif self.codename == "Swift": if output_file_type is None: @@ -1432,7 +1525,7 @@ def set_output_files(self, if output_file_type is None: output_file_type = "REAL4" if output_file_type not in ["REAL4"]: - warnings.warn(f"{output_file_type} is not compatible with Swift. Setting to REAL4") + warnings.warn(f"{output_file_type} is not compatible with Swift. Setting to REAL4",stacklevel=2) output_file_type = "REAL4" self.param['OUT_TYPE'] = output_file_type @@ -1445,7 +1538,7 @@ def set_output_files(self, self.param['BIN_OUT'] = output_file_name if output_format != "XV" and self.codename != "Swiftest": - warnings.warn(f"{output_format} is not compatible with {self.codename}. Setting to XV") + warnings.warn(f"{output_format} is not compatible with {self.codename}. Setting to XV",stacklevel=2) output_format = "XV" self.param["OUT_FORM"] = output_format @@ -1620,7 +1713,7 @@ def set_unit_system(self, self.param['MU2KG'] = 1000.0 self.MU_name = "g" else: - warnings.warn(f"{MU} not a recognized unit system. Using MSun as a default.") + warnings.warn(f"{MU} not a recognized unit system. Using MSun as a default.",stacklevel=2) self.param['MU2KG'] = constants.MSun self.MU_name = "MSun" @@ -1643,7 +1736,7 @@ def set_unit_system(self, self.param['DU2M'] = 100.0 self.DU_name = "cm" else: - warnings.warn(f"{DU} not a recognized unit system. Using AU as a default.") + warnings.warn(f"{DU} not a recognized unit system. Using AU as a default.",stacklevel=2) self.param['DU2M'] = constants.AU2M self.DU_name = "AU" @@ -1663,7 +1756,7 @@ def set_unit_system(self, self.param['TU2S'] = 1.0 self.TU_name = "s" else: - warnings.warn(f"{TU} not a recognized unit system. Using YR as a default.") + warnings.warn(f"{TU} not a recognized unit system. Using YR as a default.",stacklevel=2) self.param['TU2S'] = constants.YR2S self.TU_name = "y" @@ -1846,7 +1939,7 @@ def set_distance_range(self, if qmin_coord is not None: valid_qmin_coord = ["HELIO","BARY"] if qmin_coord.upper() not in valid_qmin_coord: - warnings.warn(f"qmin_coord = {qmin_coord} is not a valid option. Must be one of",','.join(valid_qmin_coord)) + warnings.warn(f"qmin_coord = {qmin_coord} is not a valid option. Must be one of",','.join(valid_qmin_coord),stacklevel=2) self.param['CHK_QMIN_COORD'] = valid_qmin_coord[0] else: self.param['CHK_QMIN_COORD'] = qmin_coord.upper() @@ -1969,7 +2062,7 @@ def add_solar_system_body(self, if type(ephemeris_id) is int: ephemeris_id = [ephemeris_id] if len(ephemeris_id) != len(name): - warnings.warn(f"The length of ephemeris_id ({len(ephemeris_id)}) does not match the length of name ({len(name)})") + warnings.warn(f"The length of ephemeris_id ({len(ephemeris_id)}) does not match the length of name ({len(name)})",stacklevel=2) return None else: ephemeris_id = [None] * len(name) @@ -1982,42 +2075,47 @@ def add_solar_system_body(self, try: datetime.datetime.fromisoformat(date) except: - warnings.warn(f"{date} is not a valid date format. Must be 'YYYY-MM-DD'. Setting to {self.ephemeris_date}") + warnings.warn(f"{date} is not a valid date format. Must be 'YYYY-MM-DD'. Setting to {self.ephemeris_date}",stacklevel=2) date = self.ephemeris_date if source.upper() != "HORIZONS": - warnings.warn("Currently only the JPL Horizons ephemeris service is supported") + warnings.warn("Currently only the JPL Horizons ephemeris service is supported",stacklevel=2) body_list = [] for i,n in enumerate(name): body_list.append(init_cond.solar_system_horizons(n, self.param, date, idval=ephemeris_id[i])) #Convert the list receieved from the solar_system_horizons output and turn it into arguments to vec2xr - name,v1,v2,v3,v4,v5,v6,ephemeris_id,GMpl,Rpl,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4 = tuple(np.squeeze(np.hsplit(np.array(body_list),19))) + if len(body_list) == 1: + name,v1,v2,v3,v4,v5,v6,ephemeris_id,Gmass,radius,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4 = tuple(np.hsplit(np.array(body_list[0]),19)) + else: + name,v1,v2,v3,v4,v5,v6,ephemeris_id,Gmass,radius,rhill,Ip1,Ip2,Ip3,rotx,roty,rotz,J2,J4 = tuple(np.squeeze(np.hsplit(np.array(body_list),19))) + ephemeris_id = ephemeris_id.astype(int) v1 = v1.astype(np.float64) v2 = v2.astype(np.float64) v3 = v3.astype(np.float64) v4 = v4.astype(np.float64) v5 = v5.astype(np.float64) v6 = v6.astype(np.float64) - ephemeris_id = ephemeris_id.astype(int) - GMpl = GMpl.astype(np.float64) - Rpl = Rpl.astype(np.float64) rhill = rhill.astype(np.float64) + J2 = J2.astype(np.float64) + J4 = J4.astype(np.float64) + + Gmass = Gmass.astype(np.float64) + radius = radius.astype(np.float64) Ip1 = Ip1.astype(np.float64) Ip2 = Ip2.astype(np.float64) Ip3 = Ip3.astype(np.float64) rotx = rotx.astype(np.float64) roty = roty.astype(np.float64) rotz = rotz.astype(np.float64) - J2 = J2.astype(np.float64) - J4 = J4.astype(np.float64) - if all(np.isnan(GMpl)): - GMpl = None - if all(np.isnan(Rpl)): - Rpl = None + + if all(np.isnan(Gmass)): + Gmass = None + if all(np.isnan(radius)): + radius = None if all(np.isnan(rhill)): rhill = None if all(np.isnan(Ip1)): @@ -2040,13 +2138,14 @@ def add_solar_system_body(self, t = self.param['TSTART'] dsnew = init_cond.vec2xr(self.param,name,v1,v2,v3,v4,v5,v6,ephemeris_id, - GMpl=GMpl, Rpl=Rpl, rhill=rhill, + GMpl=Gmass, Rpl=radius, rhill=rhill, Ip1=Ip1, Ip2=Ip2, Ip3=Ip3, rotx=rotx, roty=roty, rotz=rotz, J2=J2, J4=J4, t=t) dsnew = self._combine_and_fix_dsnew(dsnew) - self.save() + if dsnew['npl'] > 0 or dsnew['ntp'] > 0: + self.save(verbose=False) return dsnew @@ -2091,8 +2190,9 @@ def set_ephemeris_date(self, datetime.datetime.fromisoformat(ephemeris_date) except: valid_date_args = ['"MBCL"', '"TODAY"', '"YYYY-MM-DD"'] - warnings.warn(f"{ephemeris_date} is not a valid format. Valid options include:", ', '.join(valid_date_args)) - warnings.warn("Using MBCL for date.") + msg = f"{ephemeris_date} is not a valid format. Valid options include:", ', '.join(valid_date_args) + msg += "\nUsing MBCL for date." + warnings.warn(msg,stacklevel=2) ephemeris_date = minton_bcl self.ephemeris_date = ephemeris_date @@ -2127,7 +2227,7 @@ def get_ephemeris_date(self, arg_list: str | List[str] | None = None, verbose: b try: self.ephemeris_date except: - warnings.warn(f"ephemeris_date is not set") + warnings.warn(f"ephemeris_date is not set",stacklevel=2) return valid_arg = {"ephemeris_date": self.ephemeris_date} @@ -2173,23 +2273,28 @@ def _get_instance_var(self, arg_list: str | List[str], valid_arg: Dict, verbose: return tuple(arg_vals) def add_body(self, - name: str | List[str] | npt.NDArray[np.str_], - v1: float | List[float] | npt.NDArray[np.float_], - v2: float | List[float] | npt.NDArray[np.float_], - v3: float | List[float] | npt.NDArray[np.float_], - v4: float | List[float] | npt.NDArray[np.float_], - v5: float | List[float] | npt.NDArray[np.float_], - v6: float | List[float] | npt.NDArray[np.float_], + name: str | List[str] | npt.NDArray[np.str_] | None=None, idvals: int | list[int] | npt.NDArray[np.int_] | None=None, - GMpl: float | List[float] | npt.NDArray[np.float_] | None=None, - Rpl: float | List[float] | npt.NDArray[np.float_] | None=None, + v1: float | List[float] | npt.NDArray[np.float_] | None = None, + v2: float | List[float] | npt.NDArray[np.float_] | None = None, + v3: float | List[float] | npt.NDArray[np.float_] | None = None, + v4: float | List[float] | npt.NDArray[np.float_] | None = None, + v5: float | List[float] | npt.NDArray[np.float_] | None = None, + v6: float | List[float] | npt.NDArray[np.float_] | None = None, + xh: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None = None, + vh: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None = None, + mass: float | List[float] | npt.NDArray[np.float_] | None=None, + Gmass: float | List[float] | npt.NDArray[np.float_] | None=None, + radius: float | List[float] | npt.NDArray[np.float_] | None=None, rhill: float | List[float] | npt.NDArray[np.float_] | None=None, Ip1: float | List[float] | npt.NDArray[np.float_] | None=None, Ip2: float | List[float] | npt.NDArray[np.float_] | None=None, Ip3: float | List[float] | npt.NDArray[np.float_] | None=None, + Ip: List[float] | npt.NDArray[np.float_] | None=None, rotx: float | List[float] | npt.NDArray[np.float_] | None=None, roty: float | List[float] | npt.NDArray[np.float_] | None=None, rotz: float | List[float] | npt.NDArray[np.float_] | None=None, + rot: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None=None, J2: float | List[float] | npt.NDArray[np.float_] | None=None, J4: float | List[float] | npt.NDArray[np.float_] | None=None): """ @@ -2200,32 +2305,45 @@ def add_body(self, Parameters ---------- - name : str or array-like of str - Name or names of - v1 : float or array-like of float + name : str or array-like of str, optional + Name or names of Bodies. If none passed, name will be "Body" + idvals : int or array-like of int, optional + Unique id values. If not passed, an id will be assigned in ascending order starting from the pre-existing + Dataset ids. + v1 : float or array-like of float, optional xhx for param['IN_FORM'] == "XV"; a for param['IN_FORM'] == "EL" - v2 : float or array-like of float + v2 : float or array-like of float, optional xhy for param['IN_FORM'] == "XV"; e for param['IN_FORM'] == "EL" - v3 : float or array-like of float + v3 : float or array-like of float, optional xhz for param['IN_FORM'] == "XV"; inc for param['IN_FORM'] == "EL" - v4 : float or array-like of float + v4 : float or array-like of float, optional vhx for param['IN_FORM'] == "XV"; capom for param['IN_FORM'] == "EL" - v5 : float or array-like of float + v5 : float or array-like of float, optional vhy for param['IN_FORM'] == "XV"; omega for param['IN_FORM'] == "EL" - v6 : float or array-like of float + v6 : float or array-like of float, optional vhz for param['IN_FORM'] == "XV"; capm for param['IN_FORM'] == "EL" - idvals : int or array-like of int, optional - Unique id values. If not passed, this will be computed based on the pre-existing Dataset ids. + xh : (n,3) array-like of float, optional + Position vector array. This can be used instead of passing v1, v2, and v3 sepearately for "XV" input format + vh : (n,3) array-like of float, optional + Velocity vector array. This can be used instead of passing v4, v5, and v6 sepearately for "XV" input format + mass : float or array-like of float, optional + mass values if these are massive bodies (only one of mass or Gmass can be passed) Gmass : float or array-like of float, optional - G*mass values if these are massive bodies + G*mass values if these are massive bodies (only one of mass or Gmass can be passed) radius : float or array-like of float, optional Radius values if these are massive bodies - rhill : float, optional + rhill : float or array-like of float, optional Hill's radius values if these are massive bodies - Ip1,y,z : float, optional - Principal axes moments of inertia these are massive bodies with rotation enabled - rotx,y,z: float, optional + Ip<1,2,3> : float or array-like of float, optional + Principal axes moments of inertia if these are massive bodies with rotation enabled + rot: float or array-like of float, optional Rotation rate vector components if these are massive bodies with rotation enabled + rot: (3) or (n,3) array-like of float, optional + Rotation rate vectors if these are massive bodies with rotation enabled. This can be used instead of passing + rotx, roty, and rotz separately + Ip: (3) or (n,3) array-like of flaot, optional + Principal axes moments of inertia vectors if these are massive bodies with rotation enabled. This can be used + instead of passing Ip1, Ip2, and Ip3 separately Returns ------- @@ -2242,40 +2360,76 @@ def input_to_array(val,t,n=None): t = np.int64 elif t == "s": t = np.str + if val is None: - return None + return None, n + elif isinstance(val, np.ndarray): + pass elif np.isscalar(val): val = np.array([val],dtype=t) - elif type(val) is list: - val = np.array(val,dtype=t) + else: + try: + val = np.array(val,dtype=t) + except: + raise ValueError(f"{val} cannot be converted to a numpy array") if n is None: return val, len(val) else: if n != len(val): - raise ValueError(f"Error! Mismatched array lengths in add_body. Got {len(val)} when expecting {n}") - return val - - - name,nbodies = input_to_array(name,"s") - v1 = input_to_array(v1,"f",nbodies) - v2 = input_to_array(v2,"f",nbodies) - v3 = input_to_array(v3,"f",nbodies) - v4 = input_to_array(v4,"f",nbodies) - v5 = input_to_array(v5,"f",nbodies) - v6 = input_to_array(v6,"f",nbodies) - idvals = input_to_array(idvals,"i",nbodies) - GMpl = input_to_array(GMpl,"f",nbodies) - rhill = input_to_array(rhill,"f",nbodies) - Rpl = input_to_array(Rpl,"f",nbodies) - Ip1 = input_to_array(Ip1,"f",nbodies) - Ip2 = input_to_array(Ip2,"f",nbodies) - Ip3 = input_to_array(Ip3,"f",nbodies) - rotx = input_to_array(rotx,"f",nbodies) - roty = input_to_array(roty,"f",nbodies) - rotz = input_to_array(rotz,"f",nbodies) - J2 = input_to_array(J2,"f",nbodies) - J4 = input_to_array(J4,"f",nbodies) + raise ValueError(f"Mismatched array lengths in add_body. Got {len(val)} when expecting {n}") + return val, n + + def input_to_array_3d(val,n=None): + if val is None: + return None, n + elif isinstance(val, np.ndarray): + pass + else: + try: + val = np.array(val,dtype=np.float64) + except: + raise ValueError(f"{val} cannot be converted to a numpy array") + if n is None: + if val.dim > 2 or val.dim == 0: + raise ValueError(f"Argument must be an (n,3) array. This one is {val.shape}") + else: + if val.shape[-1] != 3: + raise ValueError(f"Argument must be a 3-dimensional vector. This one has {val.shape[0]}!") + if val.dim == 1: + n = 1 + else: + n = val.shape[0] + elif val.shape != (n,3): + raise ValueError(f"Argument is an incorrect shape. Expected {(n,3)}. Got {val.shape} instead") + return val, n + + nbodies = None + name,nbodies = input_to_array(name,"s",nbodies) + v1,nbodies = input_to_array(v1,"f",nbodies) + v2,nbodies = input_to_array(v2,"f",nbodies) + v3,nbodies = input_to_array(v3,"f",nbodies) + v4,nbodies = input_to_array(v4,"f",nbodies) + v5,nbodies = input_to_array(v5,"f",nbodies) + v6,nbodies = input_to_array(v6,"f",nbodies) + idvals,nbodies = input_to_array(idvals,"i",nbodies) + mass,nbodies = input_to_array(mass,"f",nbodies) + Gmass,nbodies = input_to_array(Gmass,"f",nbodies) + rhill,nbodies = input_to_array(rhill,"f",nbodies) + radius,nbodies = input_to_array(radius,"f",nbodies) + Ip1,nbodies = input_to_array(Ip1,"f",nbodies) + Ip2,nbodies = input_to_array(Ip2,"f",nbodies) + Ip3,nbodies = input_to_array(Ip3,"f",nbodies) + rotx,nbodies = input_to_array(rotx,"f",nbodies) + roty,nbodies = input_to_array(roty,"f",nbodies) + rotz,nbodies = input_to_array(rotz,"f",nbodies) + J2,nbodies = input_to_array(J2,"f",nbodies) + J4,nbodies = input_to_array(J4,"f",nbodies) + + xh,nbodies = input_to_array_3d(xh,nbodies) + vh,nbodies = input_to_array_3d(vh,nbodies) + rot,nbodies = input_to_array_3d(rot,nbodies) + Ip,nbodies = input_to_array_3d(Ip,nbodies) if len(self.data) == 0: maxid = -1 @@ -2285,6 +2439,9 @@ def input_to_array(val,t,n=None): if idvals is None: idvals = np.arange(start=maxid+1,stop=maxid+1+nbodies,dtype=int) + if name is None: + name=np.char.mod(f"Body%d",idvals) + if len(self.data) > 0: dup_id = np.in1d(idvals, self.data.id) if any(dup_id): @@ -2292,14 +2449,52 @@ def input_to_array(val,t,n=None): t = self.param['TSTART'] + if xh is not None: + if v1 is not None or v2 is not None or v3 is not None: + raise ValueError("Cannot use xh and v1,v2,v3 inputs simultaneously!") + else: + v1 = xh.T[0] + v2 = xh.T[1] + v3 = xh.T[2] + + if vh is not None: + if v4 is not None or v5 is not None or v6 is not None: + raise ValueError("Cannot use vh and v4,v5,v6 inputs simultaneously!") + else: + v4 = vh.T[0] + v5 = vh.T[1] + v6 = vh.T[2] + + if rot is not None: + if rotx is not None or roty is not None or rotz is not None: + raise ValueError("Cannot use rot and rotx,roty,rotz inputs simultaneously!") + else: + rotx = rot.T[0] + roty = rot.T[1] + rotz = rot.T[2] + + if Ip is not None: + if Ip1 is not None or Ip2 is not None or Ip3 is not None: + raise ValueError("Cannot use Ip and Ip1,Ip2,Ip3 inputs simultaneously!") + else: + Ip1 = Ip.T[0] + Ip2 = Ip.T[1] + Ip3 = Ip.T[2] + + if mass is not None: + if Gmass is not None: + raise ValueError("Cannot use mass and Gmass inputs simultaneously!") + else: + Gmass = self.param['GU'] * mass + dsnew = init_cond.vec2xr(self.param, name, v1, v2, v3, v4, v5, v6, idvals, - GMpl=GMpl, Rpl=Rpl, rhill=rhill, + GMpl=Gmass, Rpl=radius, rhill=rhill, Ip1=Ip1, Ip2=Ip2, Ip3=Ip3, rotx=rotx, roty=roty, rotz=rotz, J2=J2, J4=J4,t=t) dsnew = self._combine_and_fix_dsnew(dsnew) - self.save() + self.save(verbose=False) return dsnew @@ -2322,8 +2517,9 @@ def _combine_and_fix_dsnew(self,dsnew): dsnew = dsnew.swap_dims({"id" : "name"}) dsnew = dsnew.reset_coords("id") else: - warnings.warn("Non-unique names detected for bodies. The Dataset will be dimensioned by integer id instead of name.") - warnings.warn("Consider using unique names instead.") + msg = "Non-unique names detected for bodies. The Dataset will be dimensioned by integer id instead of name." + msg +="\nConsider using unique names instead." + print(msg) self.data = xr.combine_by_coords([self.data, dsnew]) @@ -2342,9 +2538,13 @@ def get_nvals(ds): if "Gmass" in ds: ds['ntp'] = ds[count_dim].where(np.isnan(ds['Gmass'])).count(dim=count_dim) ds['npl'] = ds[count_dim].where(~(np.isnan(ds['Gmass']))).count(dim=count_dim) - 1 + if self.integrator == "symba" and "GMTINY" in self.param and self.param['GMTINY'] is not None: + ds['nplm'] = ds[count_dim].where(ds['Gmass'] > self.param['GMTINY']).count(dim=count_dim) - 1 else: ds['ntp'] = ds[count_dim].count(dim=count_dim) ds['npl'] = xr.full_like(ds['ntp'],0) + if self.integrator == "symba" and "GMTINY" in self.param and self.param['GMTINY'] is not None: + ds['nplm'] = xr.full_like(ds['ntp'],0) return ds dsnew = get_nvals(dsnew) @@ -2377,7 +2577,7 @@ def read_param(self, if param_file is None: param_file = self.param_file - if coename is None: + if codename is None: codename = self.codename if verbose is None: @@ -2387,13 +2587,21 @@ def read_param(self, return False if codename == "Swiftest": - self.param = io.read_swiftest_param(param_file, param, verbose=verbose) + self.param = io.read_swiftest_param(param_file, self.param, verbose=verbose) + if "NETCDF" in self.param['IN_TYPE']: + init_cond_file = self.sim_dir / self.param['NC_IN'] + if os.path.exists(init_cond_file): + param_tmp = self.param.copy() + param_tmp['BIN_OUT'] = init_cond_file + self.data = io.swiftest2xr(param_tmp, verbose=self.verbose) + else: + warnings.warn(f"Initial conditions file file {init_cond_file} not found.", stacklevel=2) elif codename == "Swifter": self.param = io.read_swifter_param(param_file, verbose=verbose) elif codename == "Swift": self.param = io.read_swift_param(param_file, verbose=verbose) else: - warnings.warn(f'{codename} is not a recognized code name. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn(f'{codename} is not a recognized code name. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return False return True @@ -2429,7 +2637,14 @@ def write_param(self, param_file = self.param_file if param is None: param = self.param - print(f"Writing parameter inputs to file {param_file}") + + if "verbose" in kwargs: + verbose = kwargs['verbose'] + else: + verbose = self.verbose + + if verbose: + print(f"Writing parameter inputs to file {param_file}") param['! VERSION'] = f"{codename} input file" # Check to see if the parameter type matches the output type. If not, we need to convert @@ -2444,7 +2659,7 @@ def write_param(self, elif codename == "Swift": io.write_swift_param(param, param_file) else: - warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", tpname="tp.swiftest.in", @@ -2474,10 +2689,10 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t """ oldparam = self.param if self.codename == newcodename: - warnings.warn(f"This parameter configuration is already in {newcodename} format") + warnings.warn(f"This parameter configuration is already in {newcodename} format",stacklevel=2) return oldparam if newcodename != "Swift" and newcodename != "Swifter" and newcodename != "Swiftest": - warnings.warn(f'{newcodename} is an invalid code type. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn(f'{newcodename} is an invalid code type. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return oldparam goodconversion = True if self.codename == "Swifter": @@ -2498,7 +2713,7 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t if goodconversion: self.write_param(param_file) else: - warnings.warn(f"Conversion from {self.codename} to {newcodename} is not supported.") + warnings.warn(f"Conversion from {self.codename} to {newcodename} is not supported.",stacklevel=2) return oldparam def bin2xr(self): @@ -2525,9 +2740,9 @@ def bin2xr(self): self.data = io.swifter2xr(param_tmp, verbose=self.verbose) if self.verbose: print('Swifter simulation data stored as xarray DataSet .data') elif self.codename == "Swift": - warnings.warn("Reading Swift simulation data is not implemented yet") + warnings.warn("Reading Swift simulation data is not implemented yet",stacklevel=2) else: - warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return def follow(self, codestyle="Swifter"): @@ -2558,7 +2773,7 @@ def follow(self, codestyle="Swifter"): i_list = [i for i in line.split(" ") if i.strip()] nskp = int(i_list[0]) except IOError: - warnings.warn('No follow.in file found') + warnings.warn('No follow.in file found',stacklevel=2) ifol = None nskp = None fol = tool.follow_swift(self.data, ifol=ifol, nskp=nskp) @@ -2595,6 +2810,12 @@ def save(self, ------- None """ + + if "verbose" in kwargs: + verbose = kwargs['verbose'] + else: + verbose = self%verbose + if codename is None: codename = self.codename if param_file is None: @@ -2603,17 +2824,18 @@ def save(self, param = self.param if codename == "Swiftest": - io.swiftest_xr2infile(ds=self.data, param=param, in_type=self.param['IN_TYPE'], framenum=framenum) - self.write_param(param_file=param_file) + infile_name = Path(self.sim_dir) / 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) + self.write_param(param_file=param_file,**kwargs) elif codename == "Swifter": if codename == "Swiftest": swifter_param = io.swiftest2swifter_param(param) else: swifter_param = param io.swifter_xr2infile(self.data, swifter_param, framenum) - self.write_param(param_file, param=swifter_param) + self.write_param(param_file, param=swifter_param,**kwargs) else: - warnings.warn(f'Saving to {codename} not supported') + warnings.warn(f'Saving to {codename} not supported',stacklevel=2) return @@ -2658,7 +2880,7 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil elif self.param['OUT_TYPE'] == 'NETCDF_FLOAT': new_param['IN_TYPE'] = 'NETCDF_FLOAT' else: - warnings.warn(f"{self.param['OUT_TYPE']} is an invalid OUT_TYPE file") + warnings.warn(f"{self.param['OUT_TYPE']} is an invalid OUT_TYPE file",stacklevel=2) return if self.param['BIN_OUT'] != new_param['BIN_OUT'] and restart: diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d467332f8..fd0bf10c3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -45,6 +45,7 @@ SET(FAST_MATH_FILES ${SRC}/helio/helio_step.f90 ${SRC}/helio/helio_util.f90 ${SRC}/io/io.f90 + ${SRC}/io/io_progress_bar.f90 ${SRC}/netcdf/netcdf.f90 ${SRC}/obl/obl.f90 ${SRC}/operators/operator_cross.f90 diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index d59e2a9b7..3ec23ef99 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -57,8 +57,12 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa end if f_spin = F_SPIN_FIRST - lk_plpl = allocated(pl%k_plpl) - if (lk_plpl) deallocate(pl%k_plpl) + if (param%lflatten_interactions) then + lk_plpl = allocated(pl%k_plpl) + if (lk_plpl) deallocate(pl%k_plpl) + else + lk_plpl = .false. + end if call frag%set_natural_scale(colliders) @@ -250,19 +254,21 @@ subroutine fraggle_generate_spins(frag, f_spin, lfailure) frag%rot(:,:) = 0.0_DP frag%ke_spin = 0.0_DP - do i = 1, nfrag - ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets - rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) & - * L_remainder(:) / norm2(L_remainder(:)) - rot_L(:) = f_spin * L_remainder(:) / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i)) - if (norm2(rot_ke) < norm2(rot_L)) then - frag%rot(:,i) = rot_ke(:) - else - frag%rot(:, i) = rot_L(:) - end if - frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 & - * dot_product(frag%rot(:, i), frag%rot(:, i)) - end do + if (norm2(L_remainder(:)) > FRAGGLE_LTOL) then + do i = 1, nfrag + ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets + rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) & + * L_remainder(:) / norm2(L_remainder(:)) + rot_L(:) = f_spin * L_remainder(:) / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i)) + if (norm2(rot_ke) < norm2(rot_L)) then + frag%rot(:,i) = rot_ke(:) + else + frag%rot(:, i) = rot_L(:) + end if + frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 & + * dot_product(frag%rot(:, i), frag%rot(:, i)) + end do + end if frag%ke_spin = 0.5_DP * frag%ke_spin lfailure = ((frag%ke_budget - frag%ke_spin - frag%ke_orbit) < 0.0_DP) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 8d9594974..e03e30eb5 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -87,7 +87,6 @@ module subroutine fraggle_util_construct_temporary_system(frag, system, param, t !! Author: David A. Minton !! !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments - !! and optionally including fragments. implicit none ! Arguments class(fraggle_fragments), intent(in) :: frag !! Fraggle fragment system object @@ -147,7 +146,6 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa ! Internals - logical, dimension(:), allocatable :: lexclude class(swiftest_nbody_system), allocatable, save :: tmpsys class(swiftest_parameters), allocatable, save :: tmpparam integer(I4B) :: npl_before, npl_after @@ -162,24 +160,24 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para npl_before = pl%nbody npl_after = npl_before + nfrag - ! Build the exluded body logical mask - allocate(lexclude(npl_after)) if (lbefore) then - lexclude(1:npl_before) = .false. - lexclude(npl_before+1:npl_after) = .true. call fraggle_util_construct_temporary_system(frag, system, param, tmpsys, tmpparam) + ! Build the exluded body logical mask for the *before* case: Only the original bodies are used to compute energy and momentum + tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = ACTIVE + tmpsys%pl%status(npl_before+1:npl_after) = INACTIVE else - lexclude(1:npl_after) = .false. - lexclude(colliders%idx(1:colliders%ncoll)) = .true. if (.not.allocated(tmpsys)) then write(*,*) "Error in fraggle_util_get_energy_momentum. " // & " This must be called with lbefore=.true. at least once before calling it with lbefore=.false." call util_exit(FAILURE) end if + ! Build the exluded body logical mask for the *after* case: Only the new bodies are used to compute energy and momentum call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam) + tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = INACTIVE + tmpsys%pl%status(npl_before+1:npl_after) = ACTIVE end if - call tmpsys%pl%flatten(param) + if (param%lflatten_interactions) call tmpsys%pl%flatten(param) call tmpsys%get_energy_and_momentum(param) diff --git a/src/io/io.f90 b/src/io/io.f90 index 86cf55728..2d7f59de2 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -9,8 +9,107 @@ submodule (swiftest_classes) s_io use swiftest + contains + module subroutine io_compact_output(self, param, timer) + !! author: David Minton + !! + !! Generates the terminal output displayed when display_style is set to COMPACT. This is used by the Python driver to + !! make nice-looking progress reports. + implicit none + + interface fmt + !! author: David Minton + !! + !! Formats a pair of variables and corresponding values for the compact display output. Generic interface for different variable types to format. + procedure :: fmt_I4B, fmt_I8B, fmt_DP + end interface + + ! Arguments + class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Input colleciton of user-defined parameters + class(*), intent(in) :: timer !! Object used for computing elapsed wall time (must be unlimited polymorphic because the walltimer module requires swiftest_classes) + ! Internals + character(len=:), allocatable :: formatted_output + + select type(timer) + class is (walltimer) + formatted_output = fmt("ILOOP",param%iloop) // fmt("T",param%t) // fmt("NPL",self%pl%nbody) // fmt("NTP",self%tp%nbody) + select type(pl => self%pl) + class is (symba_pl) + formatted_output = formatted_output // fmt("NPLM",pl%nplm) + end select + if (param%lenergy) then + formatted_output = formatted_output // fmt("LTOTERR",self%Ltot_error) // fmt("ETOTERR",self%Etot_error) // fmt("MTOTERR",self%Mtot_error) & + // fmt("KEOERR",self%ke_orbit_error) // fmt("PEERR",self%pe_error) // fmt("EORBERR",self%Eorbit_error) & + // fmt("EUNTRERR",self%Euntracked_error) // fmt("LESCERR",self%Lescape_error) // fmt("MESCERR",self%Mescape_error) + if (param%lclose) formatted_output = formatted_output // fmt("ECOLLERR",self%Ecoll_error) + if (param%lrotation) formatted_output = formatted_output // fmt("KESPINERR",self%ke_spin_error) // fmt("LSPINERR",self%Lspin_error) + end if + + if (.not. timer%main_is_started) then ! This is the start of a new run + formatted_output = formatted_output // fmt("WT",0.0_DP) // fmt("IWT",0.0_DP) // fmt("WTPS",0.0_DP) + else + formatted_output = formatted_output // fmt("WT",timer%wall_main) // fmt("IWT",timer%wall_step) // fmt("WTPS",timer%wall_per_substep) + end if + write(*,*) formatted_output + end select + return + + contains + + function fmt_I4B(varname,val) result(pair_string) + implicit none + ! Arguments + character(*), intent(in) :: varname !! The variable name of the pair + integer(I4B), intent(in) :: val !! A 4-byte integer value + ! Result + character(len=:), allocatable :: pair_string + ! Internals + character(len=24) :: str_value + + write(str_value,*) val + pair_string = trim(adjustl(varname)) // " " // trim(adjustl(str_value)) // ";" + + return + end function fmt_I4B + + function fmt_I8B(varname, val) result(pair_string) + implicit none + ! Arguments + character(*), intent(in) :: varname !! The variable name of the pair + integer(I8B), intent(in) :: val !! An 8-byte integer value + ! Result + character(len=:), allocatable :: pair_string + ! Internals + character(len=24) :: str_value + + write(str_value,*) val + pair_string = trim(adjustl(varname)) // " " // trim(adjustl(str_value)) // ";" + + return + end function fmt_I8B + + function fmt_DP(varname, val) result(pair_string) + implicit none + ! Arguments + character(*), intent(in) :: varname !! The variable name of the pair + real(DP), intent(in) :: val !! A double precision floating point value + ! Result + character(len=:), allocatable :: pair_string + ! Internals + character(len=24) :: str_value + + write(str_value,'(ES24.16)') val + pair_string = trim(adjustl(varname)) // " " // trim(adjustl(str_value)) // ";" + + return + end function fmt_DP + + end subroutine io_compact_output + + module subroutine io_conservation_report(self, param, lterminal) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! @@ -23,28 +122,17 @@ module subroutine io_conservation_report(self, param, lterminal) ! Internals real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now - real(DP) :: Eorbit_error, Etotal_error, Ecoll_error + real(DP) :: Eorbit_error, Etot_error, Ecoll_error real(DP) :: GMtot_now real(DP) :: Lerror, Merror character(len=STRMAX) :: errmsg - character(len=*), parameter :: EGYFMT = '(ES23.16,10(",",ES23.16,:))' ! Format code for all simulation output - character(len=*), parameter :: EGYHEADER = '("t,Eorbit,Ecollisions,Lx,Ly,Lz,Mtot")' integer(I4B), parameter :: EGYIU = 72 - character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5 & + character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5 & "; DEcollisions/|E0| = ", ES12.5, & "; D(Eorbit+Ecollisions)/|E0| = ", ES12.5, & "; DM/M0 = ", ES12.5)' - associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody) - if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then - if (param%lfirstenergy .and. (param%out_stat /= "OLD")) then - open(unit=EGYIU, file=param%energy_out, form="formatted", status="replace", action="write", err=667, iomsg=errmsg) - write(EGYIU,EGYHEADER, err=667, iomsg=errmsg) - else - open(unit=EGYIU, file=param%energy_out, form="formatted", status="old", action="write", & - position="append", err=667, iomsg=errmsg) - end if - end if + associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit) call pl%vb2vh(cb) call pl%xh2xb(cb) @@ -60,6 +148,9 @@ module subroutine io_conservation_report(self, param, lterminal) GMtot_now = system%GMtot + system%GMescape if (param%lfirstenergy) then + system%ke_orbit_orig = ke_orbit_now + system%ke_spin_orig = ke_spin_now + system%pe_orig = pe_now system%Eorbit_orig = Eorbit_now system%GMtot_orig = GMtot_now system%Lorbit_orig(:) = Lorbit_now(:) @@ -68,34 +159,31 @@ module subroutine io_conservation_report(self, param, lterminal) param%lfirstenergy = .false. end if - if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then - write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, system%Ecollisions, Ltot_now, GMtot_now - close(EGYIU, err = 667, iomsg = errmsg) - end if - if (.not.param%lfirstenergy) then - Lerror = norm2(Ltot_now(:) - system%Ltot_orig(:)) / norm2(system%Ltot_orig(:)) - Eorbit_error = (Eorbit_now - system%Eorbit_orig) / abs(system%Eorbit_orig) - Ecoll_error = system%Ecollisions / abs(system%Eorbit_orig) - Etotal_error = (Eorbit_now - system%Ecollisions - system%Eorbit_orig - system%Euntracked) / abs(system%Eorbit_orig) - Merror = (GMtot_now - system%GMtot_orig) / system%GMtot_orig - if (lterminal) write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror - if (abs(Merror) > 100 * epsilon(Merror)) then + system%ke_orbit_error = (ke_orbit_now - system%ke_orbit_orig) / abs(system%Eorbit_orig) + system%ke_spin_error = (ke_spin_now - system%ke_spin_orig) / abs(system%Eorbit_orig) + system%pe_error = (pe_now - system%pe_orig) / abs(system%Eorbit_orig) + system%Eorbit_error = (Eorbit_now - system%Eorbit_orig) / abs(system%Eorbit_orig) + system%Ecoll_error = system%Ecollisions / abs(system%Eorbit_orig) + system%Euntracked_error = system%Euntracked / abs(system%Eorbit_orig) + system%Etot_error = (Eorbit_now - system%Ecollisions - system%Eorbit_orig - system%Euntracked) / abs(system%Eorbit_orig) + + system%Lorbit_error = norm2(Lorbit_now(:) - system%Lorbit_orig(:)) / norm2(system%Ltot_orig(:)) + system%Lspin_error = norm2(Lspin_now(:) - system%Lspin_orig(:)) / norm2(system%Ltot_orig(:)) + system%Lescape_error = norm2(system%Lescape(:)) / norm2(system%Ltot_orig(:)) + system%Ltot_error = norm2(Ltot_now(:) - system%Ltot_orig(:)) / norm2(system%Ltot_orig(:)) + system%Mescape_error = system%GMescape / system%GMtot_orig + system%Mtot_error = (GMtot_now - system%GMtot_orig) / system%GMtot_orig + if (lterminal) write(display_unit, EGYTERMFMT) system%Ltot_error, system%Ecoll_error, system%Etot_error,system%Mtot_error + if (abs(system%Mtot_error) > 100 * epsilon(system%Mtot_error)) then write(*,*) "Severe error! Mass not conserved! Halting!" - if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then - write(*,*) "Merror = ", Merror - write(*,*) "GMtot_now : ",GMtot_now - write(*,*) "GMtot_orig: ",system%GMtot_orig - write(*,*) "Difference: ",GMtot_now - system%GMtot_orig - else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - ! Save the frame of data to the bin file in the slot just after the present one for diagnostics - param%ioutput = param%ioutput + 1_I8B - call pl%xv2el(cb) - call self%write_hdr(param%nciu, param) - call cb%write_frame(param%nciu, param) - call pl%write_frame(param%nciu, param) - call param%nciu%close() - end if + ! Save the frame of data to the bin file in the slot just after the present one for diagnostics + param%ioutput = param%ioutput + 1_I8B + call pl%xv2el(cb) + call self%write_hdr(param%nciu, param) + call cb%write_frame(param%nciu, param) + call pl%write_frame(param%nciu, param) + call param%nciu%close() call util_exit(FAILURE) end if end if @@ -140,97 +228,6 @@ module subroutine io_dump_param(self, param_file_name) end subroutine io_dump_param - module subroutine io_dump_particle_info(self, iu) - !! author: David A. Minton - !! - !! Reads in particle information object information from an open file unformatted file - implicit none - ! Arguments - class(swiftest_particle_info), intent(in) :: self !! Particle metadata information object - integer(I4B), intent(in) :: iu !! Open file unit number - ! Internals - character(STRMAX) :: errmsg - - write(iu, err = 667, iomsg = errmsg) self%name - write(iu, err = 667, iomsg = errmsg) self%particle_type - write(iu, err = 667, iomsg = errmsg) self%origin_type - write(iu, err = 667, iomsg = errmsg) self%origin_time - write(iu, err = 667, iomsg = errmsg) self%collision_id - write(iu, err = 667, iomsg = errmsg) self%origin_xh(:) - write(iu, err = 667, iomsg = errmsg) self%origin_vh(:) - - return - - 667 continue - write(*,*) "Error writing particle metadata information from file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine io_dump_particle_info - - - module subroutine io_dump_particle_info_base(self, param, idx) - !! author: David A. Minton - !! - !! Dumps the particle information data to a file. - !! Pass a list of array indices for test particles (tpidx) and/or massive bodies (plidx) to append - implicit none - ! Arguments - class(swiftest_base), intent(inout) :: self !! Swiftest base object (can be cb, pl, or tp) - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I4B), dimension(:), optional, intent(in) :: idx !! Array of test particle indices to append to the particle file - - ! Internals - logical, save :: lfirst = .true. - integer(I4B) :: i - character(STRMAX) :: errmsg - - if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then - if (lfirst) then - select case(param%out_stat) - case('APPEND') - open(unit=LUN, file=param%particle_out, status='OLD', position='APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) - case('NEW', 'UNKNOWN', 'REPLACE') - open(unit=LUN, file=param%particle_out, status=param%out_stat, form='UNFORMATTED', err=667, iomsg=errmsg) - case default - write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) - call util_exit(FAILURE) - end select - - lfirst = .false. - else - open(unit=LUN, file=param%particle_out, status='OLD', position= 'APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) - end if - - select type(self) - class is (swiftest_cb) - write(LUN, err = 667, iomsg = errmsg) self%id - call self%info%dump(LUN) - class is (swiftest_body) - if (present(idx)) then - do i = 1, size(idx) - write(LUN, err = 667, iomsg = errmsg) self%id(idx(i)) - call self%info(idx(i))%dump(LUN) - end do - else - do i = 1, self%nbody - write(LUN, err = 667, iomsg = errmsg) self%id(i) - call self%info(i)%dump(LUN) - end do - end if - end select - - close(unit = LUN, err = 667, iomsg = errmsg) - else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - call self%write_particle_info(param%nciu, param) - end if - - return - - 667 continue - write(*,*) "Error writing particle information file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine io_dump_particle_info_base - - module subroutine io_dump_base(self, param) !! author: David A. Minton !! @@ -294,49 +291,25 @@ module subroutine io_dump_system(self, param) param_file_name = trim(adjustl(DUMP_PARAM_FILE(idx))) dump_param%in_form = XV dump_param%out_stat = 'APPEND' - if ((param%out_type == REAL8_TYPE) .or. (param%out_type == REAL4_TYPE)) then - dump_param%in_type = REAL8_TYPE - dump_param%incbfile = trim(adjustl(DUMP_CB_FILE(idx))) - dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx))) - dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx))) - - dump_param%Eorbit_orig = self%Eorbit_orig - dump_param%GMtot_orig = self%GMtot_orig - dump_param%Ltot_orig(:) = self%Ltot_orig(:) - dump_param%Lorbit_orig(:) = self%Lorbit_orig(:) - dump_param%Lspin_orig(:) = self%Lspin_orig(:) - dump_param%GMescape = self%GMescape - dump_param%Ecollisions = self%Ecollisions - dump_param%Euntracked = self%Euntracked - dump_param%Lescape(:) = self%Lescape - - else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - dump_param%in_type = NETCDF_DOUBLE_TYPE - dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx))) - dump_param%nciu%id_chunk = self%pl%nbody + self%tp%nbody - dump_param%nciu%time_chunk = 1 - end if + dump_param%in_type = NETCDF_DOUBLE_TYPE + dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx))) + dump_param%nciu%id_chunk = self%pl%nbody + self%tp%nbody + dump_param%nciu%time_chunk = 1 dump_param%T0 = param%t call dump_param%dump(param_file_name) dump_param%out_form = XV - if ((param%out_type == REAL8_TYPE) .or. (param%out_type == REAL4_TYPE)) then - call self%cb%dump(dump_param) - call self%pl%dump(dump_param) - call self%tp%dump(dump_param) - else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - dump_param%outfile = trim(adjustl(DUMP_NC_FILE(idx))) - dump_param%ioutput = 0 - call dump_param%nciu%initialize(dump_param) - call self%write_hdr(dump_param%nciu, dump_param) - call self%cb%write_frame(dump_param%nciu, dump_param) - call self%pl%write_frame(dump_param%nciu, dump_param) - call self%tp%write_frame(dump_param%nciu, dump_param) - call dump_param%nciu%close() - ! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) - call param%nciu%flush(param) - end if + dump_param%outfile = trim(adjustl(DUMP_NC_FILE(idx))) + dump_param%ioutput = 0 + call dump_param%nciu%initialize(dump_param) + call self%write_hdr(dump_param%nciu, dump_param) + call self%cb%write_frame(dump_param%nciu, dump_param) + call self%pl%write_frame(dump_param%nciu, dump_param) + call self%tp%write_frame(dump_param%nciu, dump_param) + call dump_param%nciu%close() + ! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) + call param%nciu%flush(param) idx = idx + 1 if (idx > NDUMPFILES) idx = 1 @@ -345,67 +318,78 @@ module subroutine io_dump_system(self, param) end subroutine io_dump_system - module function io_get_args(integrator, param_file_name) result(ierr) + module subroutine io_get_args(integrator, param_file_name, display_style) !! author: David A. Minton !! !! Reads in the name of the parameter file from command line arguments. implicit none ! Arguments - integer(I4B) :: integrator !! Symbolic code of the requested integrator - character(len=:), allocatable :: param_file_name !! Name of the input parameters file - ! Result - integer(I4B) :: ierr !! I/O error code + character(len=:), intent(inout), allocatable :: integrator !! Symbolic code of the requested integrator + character(len=:), intent(inout), allocatable :: param_file_name !! Name of the input parameters file + character(len=:), intent(inout), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" ! Internals - character(len=STRMAX) :: arg1, arg2 - integer :: narg,ierr_arg1, ierr_arg2 + character(len=STRMAX), dimension(:), allocatable :: arg + integer(I4B), dimension(:), allocatable :: ierr + integer :: i,narg character(len=*),parameter :: linefmt = '(A)' - ierr = -1 ! Default is to fail - narg = command_argument_count() ! - if (narg == 2) then - call get_command_argument(1, arg1, status = ierr_arg1) - call get_command_argument(2, arg2, status = ierr_arg2) - if ((ierr_arg1 == 0) .and. (ierr_arg2 == 0)) then - ierr = 0 - call io_toupper(arg1) - select case(arg1) - case('BS') - integrator = BS - case('HELIO') - integrator = HELIO - case('RA15') - integrator = RA15 - case('TU4') - integrator = TU4 - case('WHM') - integrator = WHM - case('RMVS') - integrator = RMVS - case('SYMBA') - integrator = SYMBA - case('RINGMOONS') - integrator = RINGMOONS - case default - integrator = UNKNOWN_INTEGRATOR - write(*,*) trim(adjustl(arg1)) // ' is not a valid integrator.' - ierr = -1 - end select - param_file_name = trim(adjustl(arg2)) - end if - else - call get_command_argument(1, arg1, status = ierr_arg1) - if (ierr_arg1 == 0) then - if (arg1 == '-v' .or. arg1 == '--version') then - call util_version() - else if (arg1 == '-h' .or. arg1 == '--help') then - call util_exit(HELP) - end if + narg = command_argument_count() + if (narg > 0) then + allocate(arg(narg),ierr(narg)) + do i = 1,narg + call get_command_argument(i, arg(i), status = ierr(i)) + end do + if (any(ierr /= 0)) call util_exit(USAGE) + else + call util_exit(USAGE) + end if + + if (narg == 1) then + if (arg(1) == '-v' .or. arg(1) == '--version') then + call util_version() + else if (arg(1) == '-h' .or. arg(1) == '--help') then + call util_exit(HELP) + else + call util_exit(USAGE) end if + else if (narg >= 2) then + call io_toupper(arg(1)) + select case(arg(1)) + case('BS') + integrator = BS + case('HELIO') + integrator = HELIO + case('RA15') + integrator = RA15 + case('TU4') + integrator = TU4 + case('WHM') + integrator = WHM + case('RMVS') + integrator = RMVS + case('SYMBA') + integrator = SYMBA + case('RINGMOONS') + integrator = RINGMOONS + case default + integrator = UNKNOWN_INTEGRATOR + write(*,*) trim(adjustl(arg(1))) // ' is not a valid integrator.' + call util_exit(USAGE) + end select + param_file_name = trim(adjustl(arg(2))) + end if + + if (narg == 2) then + display_style = "STANDARD" + else if (narg == 3) then + call io_toupper(arg(3)) + display_style = trim(adjustl(arg(3))) + else + call util_exit(USAGE) end if - if (ierr /= 0) call util_exit(USAGE) return - end function io_get_args + end subroutine io_get_args module function io_get_old_t_final_system(self, param) result(old_t_final) @@ -554,14 +538,13 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) !! !! Adapted from David E. Kaufmann's Swifter routine io_init_param.f90 !! Adapted from Martin Duncan's Swift routine io_init_param.f - use, intrinsic :: iso_fortran_env implicit none ! Arguments class(swiftest_parameters), intent(inout) :: self !! Collection of parameters integer, intent(in) :: unit !! File unit number character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. !! If you do not include a char-literal-constant, the iotype argument contains only DT. - integer, intent(in) :: v_list(:) !! The first element passes the integrator code to the reader + character(len=*), intent(in) :: v_list(:) !! The first element passes the integrator code to the reader integer, intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 ! Internals @@ -651,12 +634,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) ifirst = ilast + 1 param_value = io_get_token(line, ifirst, ilast, iostat) read(param_value, *, err = 667, iomsg = iomsg) param%qmin_ahi - case ("ENC_OUT") - param%enc_out = param_value - case ("DISCARD_OUT") - param%discard_out = param_value - case ("ENERGY_OUT") - param%energy_out = param_value case ("EXTRA_FORCE") call io_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') param%lextra_force = .true. @@ -745,8 +722,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *, err = 667, iomsg = iomsg) param%maxid case ("MAXID_COLLISION") read(param_value, *, err = 667, iomsg = iomsg) param%maxid_collision - case ("PARTICLE_OUT") - param%particle_out = param_value case ("RESTART") if (param_value == "NO" .or. param_value == 'F') then param%lrestart = .false. @@ -797,8 +772,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) end if param%lrestart = (param%out_stat == "APPEND") if (param%outfile /= "") then - if ((param%out_type /= REAL4_TYPE) .and. (param%out_type /= REAL8_TYPE) .and. & - (param%out_type /= NETCDF_FLOAT_TYPE) .and. (param%out_type /= NETCDF_DOUBLE_TYPE)) then + if ((param%out_type /= NETCDF_FLOAT_TYPE) .and. (param%out_type /= NETCDF_DOUBLE_TYPE)) then write(iomsg,*) 'Invalid out_type: ',trim(adjustl(param%out_type)) iostat = -1 return @@ -939,7 +913,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = 0 ! Print the contents of the parameter file to standard output - ! call param%writer(unit = OUTPUT_UNIT, iotype = "none", v_list = [0], iostat = iostat, iomsg = iomsg) + call param%writer(unit = param%display_unit, iotype = "none", v_list = [0], iostat = iostat, iomsg = iomsg) end associate @@ -992,12 +966,6 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) call io_param_writer_one("OUT_FORM", param%out_form, unit) call io_param_writer_one("OUT_STAT", "APPEND", unit) end if - if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then - call io_param_writer_one("PARTICLE_OUT", param%particle_out, unit) - end if - if (param%enc_out /= "") then - call io_param_writer_one("ENC_OUT", param%enc_out, unit) - end if call io_param_writer_one("CHK_RMIN", param%rmin, unit) call io_param_writer_one("CHK_RMAX", param%rmax, unit) call io_param_writer_one("CHK_EJECT", param%rmaxu, unit) @@ -1011,17 +979,8 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) call io_param_writer_one("DU2M", param%DU2M, unit) call io_param_writer_one("RHILL_PRESENT", param%lrhill_present, unit) call io_param_writer_one("EXTRA_FORCE", param%lextra_force, unit) - if (param%discard_out /= "") then - call io_param_writer_one("DISCARD_OUT", param%discard_out, unit) - end if - if (param%discard_out /= "") then - call io_param_writer_one("BIG_DISCARD", param%lbig_discard, unit) - end if call io_param_writer_one("CHK_CLOSE", param%lclose, unit) call io_param_writer_one("ENERGY", param%lenergy, unit) - if (param%lenergy .and. (param%energy_out /= "")) then - call io_param_writer_one("ENERGY_OUT", param%energy_out, unit) - end if call io_param_writer_one("GR", param%lgr, unit) call io_param_writer_one("ROTATION", param%lrotation, unit) call io_param_writer_one("TIDES", param%ltides, unit) @@ -1031,17 +990,6 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) if (param%lenergy) then call io_param_writer_one("FIRSTENERGY", param%lfirstenergy, unit) - if ((param%out_type == REAL8_TYPE) .or. (param%out_type == REAL4_TYPE)) then - call io_param_writer_one("EORBIT_ORIG", param%Eorbit_orig, unit) - call io_param_writer_one("GMTOT_ORIG", param%GMtot_orig, unit) - call io_param_writer_one("LTOT_ORIG", param%Ltot_orig(:), unit) - call io_param_writer_one("LORBIT_ORIG", param%Lorbit_orig(:), unit) - call io_param_writer_one("LSPIN_ORIG", param%Lspin_orig(:), unit) - call io_param_writer_one("LESCAPE", param%Lescape(:), unit) - call io_param_writer_one("GMESCAPE",param%GMescape, unit) - call io_param_writer_one("ECOLLISIONS",param%Ecollisions, unit) - call io_param_writer_one("EUNTRACKED",param%Euntracked, unit) - end if end if call io_param_writer_one("FIRSTKICK",param%lfirstkick, unit) call io_param_writer_one("MAXID",param%maxid, unit) @@ -1799,8 +1747,7 @@ module subroutine io_read_in_param(self, param_file_name) character(STRMAX) :: errmsg !! Error message in UDIO procedure ! Read in name of parameter file - write(*, *) 'Parameter input file is ', trim(adjustl(param_file_name)) - write(*, *) ' ' + write(self%display_unit, *) 'Parameter input file is ', trim(adjustl(param_file_name)) self%param_file_name = param_file_name !! todo: Currently this procedure does not work in user-defined derived-type input mode @@ -1811,7 +1758,7 @@ module subroutine io_read_in_param(self, param_file_name) if (ierr == 0) return 667 continue - write(*,*) "Error reading parameter file: " // trim(adjustl(errmsg)) + write(self%display_unit,*) "Error reading parameter file: " // trim(adjustl(errmsg)) call util_exit(FAILURE) end subroutine io_read_in_param @@ -1843,73 +1790,39 @@ module subroutine io_read_in_particle_info(self, iu) end subroutine io_read_in_particle_info - module subroutine io_read_particle_info_system(self, param) + module subroutine io_set_display_param(self, display_style) !! author: David A. Minton !! - !! Reads an old particle information file for a restartd run + !! Sets the display style parameters. If display is "STANDARD" then output goes to stdout. If display is "COMPACT" + !! then it is redirected to a log file and a progress-bar is used for stdout implicit none ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: self !! Current run configuration parameters + character(*), intent(in) :: display_style !! Style of the output display ! Internals - integer(I4B) :: id, idx - logical :: lmatch character(STRMAX) :: errmsg - type(swiftest_particle_info), allocatable :: tmpinfo - - if (.not.((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE))) return ! This subroutine is only necessary for classic binary input files - - open(unit = LUN, file = param%particle_out, status = 'OLD', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - - allocate(tmpinfo, mold=self%cb%info) - select type(cb => self%cb) - class is (swiftest_cb) - select type(pl => self%pl) - class is (swiftest_pl) - select type(tp => self%tp) - class is (swiftest_tp) - associate(npl => pl%nbody, ntp => tp%nbody) - do - lmatch = .false. - read(LUN, err = 667, iomsg = errmsg, end = 333) id - - if (id == cb%id) then - call cb%info%read_in(LUN) - lmatch = .true. - else - if (npl > 0) then - idx = findloc(pl%id(1:npl), id, dim=1) - if (idx /= 0) then - call pl%info(idx)%read_in(LUN) - lmatch = .true. - end if - end if - if (.not.lmatch .and. ntp > 0) then - idx = findloc(tp%id(1:ntp), id, dim=1) - if (idx /= 0) then - call tp%info(idx)%read_in(LUN) - lmatch = .true. - end if - end if - end if - if (.not.lmatch) then - call tmpinfo%read_in(LUN) - end if - end do - end associate - close(unit = LUN, err = 667, iomsg = errmsg) - end select - end select + select case(display_style) + case ('STANDARD') + self%display_unit = OUTPUT_UNIT !! stdout from iso_fortran_env + self%log_output = .false. + case ('COMPACT', 'PROGRESS') + open(unit=SWIFTEST_LOG_OUT, file=SWIFTEST_LOG_FILE, status='replace', err = 667, iomsg = errmsg) + self%display_unit = SWIFTEST_LOG_OUT + self%log_output = .true. + case default + write(*,*) display_style, " is an unknown display style" + call util_exit(USAGE) end select - 333 continue + self%display_style = display_style + return 667 continue - write(*,*) "Error reading particle information file: " // trim(adjustl(errmsg)) + write(*,*) "Error opening swiftest log file: " // trim(adjustl(errmsg)) call util_exit(FAILURE) - end subroutine io_read_particle_info_system + end subroutine io_set_display_param module subroutine io_toupper(string) @@ -1979,9 +1892,9 @@ module subroutine io_write_discard(self, param) case('APPEND') open(unit=LUN, file=param%discard_out, status='OLD', position='APPEND', form='FORMATTED', err=667, iomsg=errmsg) case('NEW', 'REPLACE', 'UNKNOWN') - open(unit=LUN, file=param%discard_out, status=param%out_stat, form='FORMATTED', err=667, iomsg=errmsg) + open(unit=LUN, file=param%discard_out, status=out_stat, form='FORMATTED', err=667, iomsg=errmsg) case default - write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) + write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(out_stat)) call util_exit(FAILURE) end select lfirst = .false. @@ -2154,74 +2067,36 @@ module subroutine io_write_frame_system(self, param) allocate(tp, source = self%tp) iu = BINUNIT - if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then - if (lfirst) then - select case(param%out_stat) - case('APPEND') - open(unit=iu, file=param%outfile, status='OLD', position='APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) - case('NEW', 'REPLACE', 'UNKNOWN') - open(unit=iu, file=param%outfile, status=param%out_stat, form='UNFORMATTED', err=667, iomsg=errmsg) - case default - write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) - call util_exit(FAILURE) - end select - - lfirst = .false. - else - open(unit=iu, file=param%outfile, status='OLD', position= 'APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) - end if - else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - - param%nciu%id_chunk = pl%nbody + tp%nbody - param%nciu%time_chunk = max(param%istep_dump / param%istep_out, 1) - if (lfirst) then - inquire(file=param%outfile, exist=fileExists) - - select case(param%out_stat) - case('APPEND') - if (.not.fileExists) then - errmsg = param%outfile // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN" - goto 667 - end if - case('NEW') - if (fileExists) then - errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN" - goto 667 - end if - call param%nciu%initialize(param) - case('REPLACE', 'UNKNOWN') - call param%nciu%initialize(param) - end select + param%nciu%id_chunk = pl%nbody + tp%nbody + param%nciu%time_chunk = max(param%istep_dump / param%istep_out, 1) + if (lfirst) then + inquire(file=param%outfile, exist=fileExists) + + select case(param%out_stat) + case('APPEND') + if (.not.fileExists) then + errmsg = param%outfile // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN" + goto 667 + end if + case('NEW') + if (fileExists) then + errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN" + goto 667 + end if + call param%nciu%initialize(param) + case('REPLACE', 'UNKNOWN') + call param%nciu%initialize(param) + end select - lfirst = .false. - end if + lfirst = .false. end if ! Write out each data type frame - if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then - ! For these data types, do these conversion here before writing the output. - if (param%lgr) then - call pl%pv2v(param) - call tp%pv2v(param) - end if - - if ((param%out_form == EL) .or. (param%out_form == XVEL)) then ! Do an orbital element conversion prior to writing out the frame, as we have access to the central body here - call pl%xv2el(cb) - call tp%xv2el(cb) - end if - - call self%write_hdr(iu, param) - call cb%write_frame(iu, param) - call pl%write_frame(iu, param) - call tp%write_frame(iu, param) - close(iu, err = 667, iomsg = errmsg) - else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - ! For NetCDF output, because we want to store the pseudovelocity separately from the true velocity, we need to do the orbital element conversion internally - call self%write_hdr(param%nciu, param) - call cb%write_frame(param%nciu, param) - call pl%write_frame(param%nciu, param) - call tp%write_frame(param%nciu, param) - end if + ! For NetCDF output, because we want to store the pseudovelocity separately from the true velocity, we need to do the orbital element conversion internally + call self%write_hdr(param%nciu, param) + call cb%write_frame(param%nciu, param) + call pl%write_frame(param%nciu, param) + call tp%write_frame(param%nciu, param) return diff --git a/src/io/io_progress_bar.f90 b/src/io/io_progress_bar.f90 new file mode 100644 index 000000000..9a49ff935 --- /dev/null +++ b/src/io/io_progress_bar.f90 @@ -0,0 +1,98 @@ +module io_progress_bar + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Definition of classes and methods used to determine close encounters + use swiftest_globals + use swiftest_classes + implicit none + public + + character(len=1),parameter, private :: barchar = "#" !! The progress bar character + + type :: progress_bar + !! author: David A. Minton + !! + !! Implements a class for a simple progress bar that can print on the screen. + integer(I4B) :: PBARSIZE = 80 !! Number of characters acros for a whole progress bar + integer(I8B) :: loop_length !! The total number of loops that the progrees bar is executing + character(len=:), allocatable :: barstr !! The string that prints out as the progress bar + integer(I4B) :: bar_pos !! The current position of the progress bar + character(len=32) :: fmt !! The format string that is used to define the progress bar itself + character(len=64) :: message !! The current message displayed at the end of the progress bar + contains + procedure :: reset => io_pbar_reset !! Resets the progress bar to the beginning + procedure :: update => io_pbar_update !! Updates the progress bar with new values + end type progress_bar + +contains + + subroutine io_pbar_reset(self, loop_length) + !! author: David A. Minton + !! + !! Resets the progress bar to the beginning + implicit none + ! Arguments + class(progress_bar),intent(inout) :: self !! The progress bar object + integer(I8B), intent(in) :: loop_length !! The length of the loop that the progress bar is attached to + ! Internals + character(len=2) :: numchar + integer(I4B) :: k + + if (.not.allocated(self%barstr)) then + allocate(character(self%PBARSIZE) :: self%barstr) + end if + do k = 1, self%PBARSIZE + self%barstr(k:k) = " " + end do + write(numchar,'(I2)') self%PBARSIZE + self%fmt = '(A1,"[",A' // numchar // ',"] ",A,$)' + self%loop_length = loop_length + self%bar_pos = 0 + self%message = "" + + write(*,fmt=self%fmt) char(13),self%barstr,trim(adjustl(self%message)) + + return + end subroutine io_pbar_reset + + + subroutine io_pbar_update(self,i,message) + !! author: David A. Minton + !! + !! Updates the progress bar with new values + implicit none + ! Arguments + class(progress_bar), intent(inout) :: self !! Progres bar object + integer(I8B), intent(in) :: i !! The current loop index of the progress loop + character(len=*), intent(in), optional :: message !! An optional message to display to the right of the progress bar + ! Internals + real(DP) :: frac + integer(I4B) :: bar_pos !! The current integer position of the progress bar + logical :: update = .false. + + ! Compute the current position + frac = real(i,kind=DP) / real(self%loop_length,kind=DP) + bar_pos = min(int(ceiling(frac * self%PBARSIZE),kind=I4B),self%PBARSIZE) + + if (bar_pos /= self%bar_pos) then + ! Fill in the bar character up to the current position + self%barstr(bar_pos:bar_pos) = barchar + update = .true. + self%bar_pos = bar_pos + end if + + if (present(message)) then + if (message /= self%message) then + update = .true. + self%message = message + end if + end if + + if (update) write(*,fmt=self%fmt) char(13),self%barstr,trim(adjustl(self%message)) + + + return + end subroutine io_pbar_update + + +end module io_progress_bar diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 467403269..c25566fe0 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -20,10 +20,10 @@ program swiftest_driver class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated class(swiftest_parameters), allocatable :: param !! Run configuration parameters - integer(I4B) :: integrator !! Integrator type code (see swiftest_globals for symbolic names) + character(len=:), allocatable :: integrator !! Integrator type code (see swiftest_globals for symbolic names) character(len=:),allocatable :: param_file_name !! Name of the file containing user-defined parameters + character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" integer(I4B) :: ierr !! I/O error code - integer(I8B) :: iloop !! Loop counter integer(I8B) :: idump !! Dump cadence counter integer(I8B) :: iout !! Output cadence counter integer(I8B) :: ioutput_t0 !! The output frame counter at time 0 @@ -31,16 +31,19 @@ program swiftest_driver real(DP) :: old_t_final = 0.0_DP !! Output time at which writing should start, in order to prevent duplicate lines being written for restarts type(walltimer) :: integration_timer !! Object used for computing elapsed wall time real(DP) :: tfrac - character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & - '"; Number of active pl, tp = ", I5, ", ", I5)' - character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & - '"; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)' - - ierr = io_get_args(integrator, param_file_name) - if (ierr /= 0) then - write(*,*) 'Error reading in arguments from the command line' - call util_exit(FAILURE) - end if + type(progress_bar) :: pbar !! Object used to print out a progress bar + character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & + '"; Number of active pl, tp = ", I6, ", ", I6)' + character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & + '"; Number of active plm, pl, tp = ", I6, ", ", I6, ", ", I6)' + character(*), parameter :: pbarfmt = '("Time = ", ES12.5," of ",ES12.5)' + character(len=64) :: pbarmessage + + character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)' + + + call io_get_args(integrator, param_file_name, display_style) + !> Read in the user-defined parameters file and the initial conditions of the system select case(integrator) case(symba) @@ -49,17 +52,29 @@ program swiftest_driver allocate(swiftest_parameters :: param) end select param%integrator = integrator + call param%set_display(display_style) + + !> Define the maximum number of threads + nthreads = 1 ! In the *serial* case + !$ nthreads = omp_get_max_threads() ! In the *parallel* case + !$ write(param%display_unit,'(a)') ' OpenMP parameters:' + !$ write(param%display_unit,'(a)') ' ------------------' + !$ write(param%display_unit,'(a,i3,/)') ' Number of threads = ', nthreads + !$ if (param%log_output) write(*,'(a,i3)') ' OpenMP: Number of threads = ',nthreads call setup_construct_system(nbody_system, param) call param%read_in(param_file_name) - associate(t => param%t, & - t0 => param%t0, & - dt => param%dt, & - tstop => param%tstop, & - istep_out => param%istep_out, & - istep_dump => param%istep_dump, & - ioutput => param%ioutput) + associate(t => param%t, & + t0 => param%t0, & + dt => param%dt, & + tstop => param%tstop, & + iloop => param%iloop, & + istep_out => param%istep_out, & + istep_dump => param%istep_dump, & + ioutput => param%ioutput, & + display_style => param%display_style, & + display_unit => param%display_unit) call nbody_system%initialize(param) t = t0 @@ -80,14 +95,16 @@ program swiftest_driver if (istep_out > 0) call nbody_system%write_frame(param) end if - !> Define the maximum number of threads - nthreads = 1 ! In the *serial* case - !$ nthreads = omp_get_max_threads() ! In the *parallel* case - !$ write(*,'(a)') ' OpenMP parameters:' - !$ write(*,'(a)') ' ------------------' - !$ write(*,'(a,i3,/)') ' Number of threads = ', nthreads - write(*, *) " *************** Main Loop *************** " + write(display_unit, *) " *************** Main Loop *************** " if (param%lrestart .and. param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) + if (display_style == "PROGRESS") then + call pbar%reset(nloops) + write(pbarmessage,fmt=pbarfmt) t0, tstop + call pbar%update(1,message=pbarmessage) + else if (display_style == "COMPACT") then + write(*,*) "SWIFTEST START " // trim(adjustl(param%integrator)) + call nbody_system%compact_output(param,integration_timer) + end if do iloop = 1, nloops !> Step the system forward in time call integration_timer%start() @@ -98,24 +115,33 @@ program swiftest_driver !> Evaluate any discards or collisional outcomes call nbody_system%discard(param) + if (display_style == "PROGRESS") call pbar%update(iloop) !> If the loop counter is at the output cadence value, append the data file with a single frame if (istep_out > 0) then iout = iout - 1 if (iout == 0) then ioutput = ioutput_t0 + iloop / istep_out - if (t > old_t_final) call nbody_system%write_frame(param) + call nbody_system%write_frame(param) tfrac = (param%t - param%t0) / (param%tstop - param%t0) select type(pl => nbody_system%pl) class is (symba_pl) - write(*, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody + write(display_unit, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody class default - write(*, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody + write(display_unit, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody end select if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) - call integration_timer%report(message="Integration steps:", nsubsteps=istep_out) + call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out) + + if (display_style == "PROGRESS") then + write(pbarmessage,fmt=pbarfmt) t, tstop + call pbar%update(1,message=pbarmessage) + else if (display_style == "COMPACT") then + call nbody_system%compact_output(param,integration_timer) + end if + call integration_timer%reset() iout = istep_out @@ -131,6 +157,7 @@ program swiftest_driver end if end if end do + if (display_style == "COMPACT") write(*,*) "SWIFTEST STOP" // trim(adjustl(param%integrator)) end associate call nbody_system%dealloc() diff --git a/src/modules/swiftest.f90 b/src/modules/swiftest.f90 index 6f84c66b9..edc41f134 100644 --- a/src/modules/swiftest.f90 +++ b/src/modules/swiftest.f90 @@ -23,6 +23,7 @@ module swiftest use lambda_function use walltime_classes use encounter_classes + use io_progress_bar !use advisor_annotate !$ use omp_lib implicit none diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 6c0ac2be3..f04346624 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -32,14 +32,15 @@ module swiftest_classes !> User defined parameters that are read in from the parameters input file. !> Each paramter is initialized to a default values. type :: swiftest_parameters - integer(I4B) :: integrator = UNKNOWN_INTEGRATOR !! Symbolic name of the nbody integrator used + character(STRMAX) :: integrator = UNKNOWN_INTEGRATOR !! Symbolic name of the nbody integrator used character(STRMAX) :: param_file_name = "param.in" !! The default name of the parameter input file integer(I4B) :: maxid = -1 !! The current maximum particle id number - integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number + integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number real(DP) :: t0 = -1.0_DP !! Integration start time real(DP) :: t = -1.0_DP !! Integration current time real(DP) :: tstop = -1.0_DP !! Integration stop time real(DP) :: dt = -1.0_DP !! Time step + integer(I8B) :: iloop = 0_I8B !! Main loop counter integer(I8B) :: ioutput = 0_I8B !! Output counter character(STRMAX) :: incbfile = CB_INFILE !! Name of input file for the central body character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies @@ -52,7 +53,6 @@ module swiftest_classes character(STRMAX) :: out_type = NETCDF_DOUBLE_TYPE !! Binary format of output file character(STRMAX) :: out_form = XVEL !! Data to write to output file character(STRMAX) :: out_stat = 'NEW' !! Open status for output binary file - character(STRMAX) :: particle_out = PARTICLE_OUTFILE !! Name of output particle information file integer(I4B) :: istep_dump = -1 !! Number of time steps between dumps real(DP) :: rmin = -1.0_DP !! Minimum heliocentric radius for test particle real(DP) :: rmax = -1.0_DP !! Maximum heliocentric radius for test particle @@ -72,6 +72,7 @@ module swiftest_classes character(NAMELEN) :: interaction_loops = "ADAPTIVE" !! Method used to compute interaction loops. Options are "TRIANGULAR", "FLAT", or "ADAPTIVE" character(NAMELEN) :: encounter_check_plpl = "ADAPTIVE" !! Method used to compute pl-pl encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" character(NAMELEN) :: encounter_check_pltp = "ADAPTIVE" !! Method used to compute pl-tp encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" + ! The following are used internally, and are not set by the user, but instead are determined by the input value of INTERACTION_LOOPS logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interaction loops logical :: ladaptive_interactions = .false. !! Adaptive interaction loop is turned on (choose between TRIANGULAR and FLAT based on periodic timing tests) @@ -104,6 +105,10 @@ module swiftest_classes logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step logical :: lrestart = .false. !! Indicates whether or not this is a restarted run + character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT"}). Default is "STANDARD" + integer(I4B) :: display_unit !! File unit number for display (either to stdout or to a log file) + logical :: log_output = .false. !! Logs the output to file instead of displaying it on the terminal + ! Future features not implemented or in development logical :: lgr = .false. !! Turn on GR logical :: lyarkovsky = .false. !! Turn on Yarkovsky effect @@ -111,10 +116,11 @@ module swiftest_classes type(netcdf_parameters) :: nciu !! Object containing NetCDF parameters contains - procedure :: reader => io_param_reader - procedure :: writer => io_param_writer - procedure :: dump => io_dump_param - procedure :: read_in => io_read_in_param + procedure :: reader => io_param_reader + procedure :: writer => io_param_writer + procedure :: dump => io_dump_param + procedure :: read_in => io_read_in_param + procedure :: set_display => io_set_display_param end type swiftest_parameters @@ -137,7 +143,6 @@ module swiftest_classes real(DP), dimension(NDIM) :: discard_vh !! The heliocentric velocity vector at the time of the particle's discard integer(I4B) :: discard_body_id !! The id of the other body involved in the discard (0 if no other body involved) contains - procedure :: dump => io_dump_particle_info !! Dumps contents of particle information to file procedure :: read_in => io_read_in_particle_info !! Read in a particle information object from an open file procedure :: copy => util_copy_particle_info !! Copies one set of information object components into another, component-by-component procedure :: set_value => util_set_particle_info !! Sets one or more values of the particle information metadata object @@ -151,12 +156,11 @@ module swiftest_classes contains !! The minimal methods that all systems must have procedure :: dump => io_dump_base !! Dump contents to file - procedure :: dump_particle_info => io_dump_particle_info_base !! Dump contents of particle information metadata to file procedure :: read_in => io_read_in_base !! Read in body initial conditions from a file procedure :: write_frame_netcdf => netcdf_write_frame_base !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format - procedure :: write_particle_info_netcdf => netcdf_write_particle_info_base !! Writes out the particle information metadata to NetCDF file + procedure :: write_particle_info_netcdf => netcdf_write_particle_info_base !! Dump contents of particle information metadata to file generic :: write_frame => write_frame_netcdf !! Set up generic procedure that will switch between NetCDF or Fortran binary depending on arguments - generic :: write_particle_info => write_particle_info_netcdf + generic :: write_particle_info => write_particle_info_netcdf !! Set up generic procedure that will switch between NetCDF or Fortran binary depending on arguments end type swiftest_base !******************************************************************************************************************************** @@ -364,6 +368,9 @@ module swiftest_classes real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! System orbital angular momentum vector real(DP), dimension(NDIM) :: Lspin = 0.0_DP !! System spin angular momentum vector real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector + real(DP) :: ke_orbit_orig = 0.0_DP!! Initial orbital kinetic energy + real(DP) :: ke_spin_orig = 0.0_DP!! Initial spin kinetic energy + real(DP) :: pe_orig = 0.0_DP !! Initial potential energy real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector @@ -373,6 +380,22 @@ module swiftest_classes real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies + + ! Energy, momentum, and mass errors (used in error reporting) + real(DP) :: ke_orbit_error = 0.0_DP + real(DP) :: ke_spin_error = 0.0_DP + real(DP) :: pe_error = 0.0_DP + real(DP) :: Eorbit_error = 0.0_DP + real(DP) :: Ecoll_error = 0.0_DP + real(DP) :: Euntracked_error = 0.0_DP + real(DP) :: Etot_error = 0.0_DP + real(DP) :: Lorbit_error = 0.0_DP + real(DP) :: Lspin_error = 0.0_DP + real(DP) :: Lescape_error = 0.0_DP + real(DP) :: Ltot_error = 0.0_DP + real(DP) :: Mtot_error = 0.0_DP + real(DP) :: Mescape_error = 0.0_DP + logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated !! separately from massive bodies. Massive body variables are saved at half steps, and passed to !! the test particles @@ -382,6 +405,7 @@ module swiftest_classes ! Concrete classes that are common to the basic integrator (only test particles considered for discard) procedure :: discard => discard_system !! Perform a discard step on the system + procedure :: compact_output => io_compact_output !! Prints out out terminal output when display_style is set to COMPACT procedure :: conservation_report => io_conservation_report !! Compute energy and momentum and print out the change with time procedure :: dump => io_dump_system !! Dump the state of the system to a file procedure :: get_old_t_final_bin => io_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the binary output. @@ -394,7 +418,6 @@ module swiftest_classes procedure :: read_hdr_netcdf => netcdf_read_hdr_system !! Read a header for an output frame in NetCDF format procedure :: write_hdr_netcdf => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format procedure :: read_in => io_read_in_system !! Reads the initial conditions for an nbody system - procedure :: read_particle_info_bin => io_read_particle_info_system !! Read in particle metadata from file procedure :: read_particle_info_netcdf => netcdf_read_particle_info_system !! Read in particle metadata from file procedure :: write_discard => io_write_discard !! Write out information about discarded test particles procedure :: obl_pot => obl_pot_system !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body @@ -411,9 +434,10 @@ module swiftest_classes generic :: read_hdr => read_hdr_netcdf !! Generic method call for reading headers generic :: read_frame => read_frame_bin, read_frame_netcdf !! Generic method call for reading a frame of output data generic :: write_frame => write_frame_bin, write_frame_netcdf !! Generic method call for writing a frame of output data - generic :: read_particle_info => read_particle_info_bin, read_particle_info_netcdf !! Genereric method call for reading in the particle information metadata + generic :: read_particle_info => read_particle_info_netcdf !! Genereric method call for reading in the particle information metadata end type swiftest_nbody_system + abstract interface subroutine abstract_accel(self, system, param, t, lbeg) @@ -582,6 +606,13 @@ pure module subroutine gr_vh2pv_body(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine gr_vh2pv_body + module subroutine io_compact_output(self, param, timer) + implicit none + class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Input colleciton of user-defined parameters + class(*), intent(in) :: timer !! Object used for computing elapsed wall time + end subroutine io_compact_output + module subroutine io_conservation_report(self, param, lterminal) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object @@ -595,19 +626,6 @@ module subroutine io_dump_param(self, param_file_name) character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) end subroutine io_dump_param - module subroutine io_dump_particle_info_base(self, param, idx) - implicit none - class(swiftest_base), intent(inout) :: self !! Swiftest base object (can be cb, pl, or tp) - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I4B), dimension(:), optional, intent(in) :: idx !! Array of test particle indices to append to the particle file - end subroutine io_dump_particle_info_base - - module subroutine io_dump_particle_info(self, iu) - implicit none - class(swiftest_particle_info), intent(in) :: self !! Swiftest particle info metadata object - integer(I4B), intent(in) :: iu !! Open unformatted file unit number - end subroutine io_dump_particle_info - module subroutine io_dump_base(self, param) implicit none class(swiftest_base), intent(inout) :: self !! Swiftest base object @@ -620,12 +638,12 @@ module subroutine io_dump_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine io_dump_system - module function io_get_args(integrator, param_file_name) result(ierr) + module subroutine io_get_args(integrator, param_file_name, display_style) implicit none - integer(I4B) :: integrator !! Symbolic code of the requested integrator - character(len=:), allocatable :: param_file_name !! Name of the input parameters file - integer(I4B) :: ierr !! I/O error code - end function io_get_args + character(len=:), allocatable, intent(inout) :: integrator !! Symbolic code of the requested integrator + character(len=:), allocatable, intent(inout) :: param_file_name !! Name of the input parameters file + character(len=:), allocatable, intent(inout) :: display_style !! Style of the output display {"STANDARD", "COMPACT"}). Default is "STANDARD" + end subroutine io_get_args module function io_get_old_t_final_system(self, param) result(old_t_final) implicit none @@ -662,7 +680,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) integer(I4B), intent(in) :: unit !! File unit number character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. !! If you do not include a char-literal-constant, the iotype argument contains only DT. - integer(I4B), intent(in) :: v_list(:) !! The first element passes the integrator code to the reader + character(len=*), intent(in) :: v_list(:) !! The first element passes the integrator code to the reader integer(I4B), intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine io_param_reader @@ -787,11 +805,16 @@ module function io_read_frame_system(self, iu, param) result(ierr) integer(I4B) :: ierr !! Error code: returns 0 if the read is successful end function io_read_frame_system - module subroutine io_read_particle_info_system(self, param) + module subroutine io_set_display_param(self, display_style) implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine io_read_particle_info_system + class(swiftest_parameters), intent(inout) :: self !! Current run configuration parameters + character(*), intent(in) :: display_style !! Style of the output display + end subroutine io_set_display_param + + module subroutine io_toupper(string) + implicit none + character(*), intent(inout) :: string !! String to make upper case + end subroutine io_toupper module subroutine io_write_discard(self, param) implicit none @@ -799,11 +822,6 @@ module subroutine io_write_discard(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine io_write_discard - module subroutine io_toupper(string) - implicit none - character(*), intent(inout) :: string !! String to make upper case - end subroutine io_toupper - module subroutine io_write_frame_body(self, iu, param) implicit none class(swiftest_body), intent(in) :: self !! Swiftest body object diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index f6644a302..49fa0b834 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -44,15 +44,15 @@ module swiftest_globals real(SP), parameter :: VERSION_NUMBER = 1.0_SP !! swiftest version !> Symbolic name for integrator types - integer(I4B), parameter :: UNKNOWN_INTEGRATOR = 1 - integer(I4B), parameter :: BS = 2 - integer(I4B), parameter :: HELIO = 3 - integer(I4B), parameter :: RA15 = 4 - integer(I4B), parameter :: TU4 = 5 - integer(I4B), parameter :: WHM = 6 - integer(I4B), parameter :: RMVS = 7 - integer(I4B), parameter :: SYMBA = 8 - integer(I4B), parameter :: RINGMOONS = 9 + character(*), parameter :: UNKNOWN_INTEGRATOR = "UKNOWN INTEGRATOR" + character(*), parameter :: BS = "Bulirsch-Stoer" + character(*), parameter :: HELIO = "Democratic Heliocentric" + character(*), parameter :: RA15 = "Radau 15th order" + character(*), parameter :: TU4 = "T+U 4th order" + character(*), parameter :: WHM = "Wisdom-Holman Method" + character(*), parameter :: RMVS = "Regularized Mixed Variable Symplectic" + character(*), parameter :: SYMBA = "SyMBA" + character(*), parameter :: RINGMOONS = "SyMBA-RINGMOONS" integer(I4B), parameter :: STRMAX = 512 !! Maximum size of character strings integer(I4B), parameter :: NAMELEN = 32 !! Maximum size of name strings @@ -81,11 +81,6 @@ module swiftest_globals integer(I4B), parameter :: USAGE = -2 !! Symbolic name for function return/flag code for printing the usage message integer(I4B), parameter :: HELP = -3 !! Symbolic name for function return/flag code for printing the usage message - character(*), parameter :: SUCCESS_MSG = '(/, "Normal termination of Swiftest (version ", f3.1, ")")' - character(*), parameter :: FAIL_MSG = '(/, "Terminating Swiftest (version ", f3.1, ") due to error!!")' - character(*), parameter :: USAGE_MSG = '("Usage: swiftest [bs|helio|ra15|rmvs|symba|tu4|whm] ")' - character(*), parameter :: HELP_MSG = USAGE_MSG - integer(I4B), parameter :: ELLIPSE = -1 !! Symbolic names for orbit types - ellipse integer(I4B), parameter :: PARABOLA = 0 !! Symbolic names for orbit types - parabola integer(I4B), parameter :: HYPERBOLA = 1 !! Symbolic names for orbit types - hyperbola @@ -123,11 +118,13 @@ module swiftest_globals !> Standard file names integer(I4B), parameter :: NDUMPFILES = 2 - character(*), dimension(2), parameter :: DUMP_CB_FILE = ['dump_cb1.bin', 'dump_cb2.bin' ] - character(*), dimension(2), parameter :: DUMP_PL_FILE = ['dump_pl1.bin', 'dump_pl2.bin' ] - character(*), dimension(2), parameter :: DUMP_TP_FILE = ['dump_tp1.bin', 'dump_tp2.bin' ] - character(*), dimension(2), parameter :: DUMP_NC_FILE = ['dump_bin1.nc', 'dump_bin2.nc' ] + character(*), dimension(2), parameter :: DUMP_CB_FILE = ['dump_cb1.bin', 'dump_cb2.bin' ] + character(*), dimension(2), parameter :: DUMP_PL_FILE = ['dump_pl1.bin', 'dump_pl2.bin' ] + character(*), dimension(2), parameter :: DUMP_TP_FILE = ['dump_tp1.bin', 'dump_tp2.bin' ] + character(*), dimension(2), parameter :: DUMP_NC_FILE = ['dump_bin1.nc', 'dump_bin2.nc' ] character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.in', 'dump_param2.in'] + character(*), parameter :: SWIFTEST_LOG_FILE = "swiftest.log" !! Name of file to use to log output when using "COMPACT" display style + integer(I4B), parameter :: SWIFTEST_LOG_OUT = 33 !! File unit for log file when using "COMPACT" display style !> Default file names that can be changed by the user in the parameters file character(*), parameter :: CB_INFILE = 'cb.in' @@ -136,7 +133,6 @@ module swiftest_globals character(*), parameter :: NC_INFILE = 'in.nc' character(*), parameter :: BIN_OUTFILE = 'bin.nc' integer(I4B), parameter :: BINUNIT = 20 !! File unit number for the binary output file - character(*), parameter :: PARTICLE_OUTFILE = 'particle.dat' integer(I4B), parameter :: PARTICLEUNIT = 44 !! File unit number for the binary particle info output file integer(I4B), parameter :: LUN = 42 !! File unit number for files that are opened and closed within a single subroutine call, and therefore should not collide diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index cacad4b4d..e016a36b9 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -375,7 +375,7 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms integer, intent(in) :: unit !! File unit number character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. !! If you do not include a char-literal-constant, the iotype argument contains only DT. - integer, intent(in) :: v_list(:) !! The first element passes the integrator code to the reader + character(len=*), intent(in) :: v_list(:) !! The first element passes the integrator code to the reader integer, intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine symba_io_param_reader diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 0dc3f2892..536272b44 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -29,6 +29,8 @@ module walltime_classes integer(I8B) :: count_stop_step !! Value of the clock ticker at the end of a timed step integer(I8B) :: count_pause !! Value of the clock ticker at the end of a timed step real(DP) :: wall_step !! Value of the step elapsed time + real(DP) :: wall_main !! Value of the main clock elapsed time + real(DP) :: wall_per_substep !! Value of time per substep logical :: main_is_started = .false. !! Logical flag indicating whether or not the main timer has been reset or not logical :: is_paused = .false. !! Logical flag indicating whether or not the timer is paused @@ -60,10 +62,11 @@ module walltime_classes end type interaction_timer interface - module subroutine walltime_report(self, message, nsubsteps) + module subroutine walltime_report(self, message, unit, nsubsteps) implicit none class(walltimer), intent(inout) :: self !! Walltimer object character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output + integer(I4B), intent(in) :: unit !! Output file unit for report text to be directed integer(I4B), optional, intent(in) :: nsubsteps !! Number of substeps used to compute the time per step end subroutine walltime_report diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 859e8c6ba..3d0943d95 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -157,9 +157,7 @@ module subroutine setup_initialize_system(self, param) pl%lfirst = param%lfirstkick tp%lfirst = param%lfirstkick - if (param%lrestart) then - call system%read_particle_info(param) - else + if (.not.param%lrestart) then call system%init_particle_info(param) end if end associate diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 7f2792309..d9a48b52a 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -24,7 +24,7 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms integer, intent(in) :: unit !! File unit number character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. !! If you do not include a char-literal-constant, the iotype argument contains only DT. - integer, intent(in) :: v_list(:) !! The first element passes the integrator code to the reader + character(len=*), intent(in) :: v_list(:) !! The first element passes the integrator code to the reader integer, intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 ! internals diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 8d110451f..0087c24e4 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -676,7 +676,7 @@ module subroutine symba_util_rearray_pl(self, system, param) end where end select - call pl%dump_particle_info(param, idx=pack([(i, i=1, npl)], ldump_mask)) + call pl%write_particle_info(param%nciu, param) deallocate(ldump_mask) ! Reindex the new list of bodies diff --git a/src/util/util_exit.f90 b/src/util/util_exit.f90 index 61dacdf99..a7b77c197 100644 --- a/src/util/util_exit.f90 +++ b/src/util/util_exit.f90 @@ -23,6 +23,10 @@ module subroutine util_exit(code) integer(I4B), intent(in) :: code ! Internals character(*), parameter :: BAR = '("------------------------------------------------")' + character(*), parameter :: SUCCESS_MSG = '(/, "Normal termination of Swiftest (version ", f3.1, ")")' + character(*), parameter :: FAIL_MSG = '(/, "Terminating Swiftest (version ", f3.1, ") due to error!!")' + character(*), parameter :: USAGE_MSG = '("Usage: swiftest [bs|helio|ra15|rmvs|symba|tu4|whm] [standard|compact|progress|NONE]")' + character(*), parameter :: HELP_MSG = USAGE_MSG select case(code) case(SUCCESS) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 621ea80a6..ed7119d8b 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -23,7 +23,6 @@ module subroutine util_get_energy_momentum_system(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i - integer(I8B) :: nplpl real(DP) :: kecb, kespincb real(DP), dimension(self%pl%nbody) :: kepl, kespinpl real(DP), dimension(self%pl%nbody) :: Lplorbitx, Lplorbity, Lplorbitz @@ -32,7 +31,6 @@ module subroutine util_get_energy_momentum_system(self, param) real(DP) :: hx, hy, hz associate(system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) - nplpl = pl%nplpl system%Lorbit(:) = 0.0_DP system%Lspin(:) = 0.0_DP system%Ltot(:) = 0.0_DP @@ -89,7 +87,7 @@ module subroutine util_get_energy_momentum_system(self, param) end if if (param%lflatten_interactions) then - call util_get_energy_potential_flat(npl, nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + call util_get_energy_potential_flat(npl, pl%nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) else call util_get_energy_potential_triangular(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) end if diff --git a/src/walltime/walltime.f90 b/src/walltime/walltime.f90 index 6c53e2276..104b63d95 100644 --- a/src/walltime/walltime.f90 +++ b/src/walltime/walltime.f90 @@ -38,7 +38,7 @@ module subroutine walltime_stop(self) end subroutine walltime_stop - module subroutine walltime_report(self, message, nsubsteps) + module subroutine walltime_report(self, message, unit, nsubsteps) !! author: David A. Minton !! !! Prints the elapsed time information to the terminal @@ -46,6 +46,7 @@ module subroutine walltime_report(self, message, nsubsteps) ! Arguments class(walltimer), intent(inout) :: self !! Walltimer object character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output + integer(I4B), intent(in) :: unit !! Output file unit for report text to be directed integer(I4B), optional, intent(in) :: nsubsteps !! Number of substeps used to compute the time per step ! Internals character(len=*), parameter :: nosubstepfmt = '" Total wall time: ", es12.5, "; Interval wall time: ", es12.5 ' @@ -53,9 +54,6 @@ module subroutine walltime_report(self, message, nsubsteps) 'Interval wall time/step: ", es12.5' character(len=STRMAX) :: fmt integer(I8B) :: count_delta_step, count_delta_main, count_now - real(DP) :: wall_main !! Value of total elapsed time at the end of a timed step - real(DP) :: wall_step !! Value of elapsed time since the start of a timed step - real(DP) :: wall_per_substep !! Value of time per substep if (.not.self%main_is_started) then write(*,*) "Wall timer error: The step finish time cannot be calculated because the timer is not started!" @@ -65,18 +63,17 @@ module subroutine walltime_report(self, message, nsubsteps) call system_clock(count_now) count_delta_main = count_now - self%count_start_main count_delta_step = count_now - self%count_start_step - wall_main = count_delta_main / (self%count_rate * 1.0_DP) - wall_step = count_delta_step / (self%count_rate * 1.0_DP) + self%wall_main = count_delta_main / (self%count_rate * 1.0_DP) + self%wall_step = count_delta_step / (self%count_rate * 1.0_DP) if (present(nsubsteps)) then - wall_per_substep = wall_step / nsubsteps + self%wall_per_substep = self%wall_step / nsubsteps fmt = '("' // adjustl(message) // '",' // substepfmt // ')' - write(*,trim(adjustl(fmt))) wall_main, self%wall_step, wall_per_substep + write(unit,trim(adjustl(fmt))) self%wall_main, self%wall_step, self%wall_per_substep else fmt = '("' // adjustl(message) // '",' // nosubstepfmt // ')' - write(*,trim(adjustl(fmt))) wall_main, self%wall_step + write(unit,trim(adjustl(fmt))) self%wall_main, self%wall_step end if - return end subroutine walltime_report @@ -92,6 +89,7 @@ module subroutine walltime_reset(self) self%is_paused = .false. self%wall_step = 0.0_DP + self%wall_per_substep = 0.0_DP return end subroutine walltime_reset @@ -107,6 +105,7 @@ module subroutine walltime_start_main(self) call system_clock(self%count_start_main, self%count_rate, self%count_max) self%main_is_started = .true. + self%wall_main = 0.0_DP return end subroutine walltime_start_main @@ -123,7 +122,6 @@ module subroutine walltime_start(self) ! Internals integer(I8B) :: count_resume, count_delta - if (.not.self%main_is_started) then call self%reset() call self%start_main()