Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Merge branch 'encounter_storage' into debug
Browse files Browse the repository at this point in the history
  • Loading branch information
daminton committed Dec 8, 2022
2 parents ea48c5f + ca87dd1 commit 569b06d
Show file tree
Hide file tree
Showing 25 changed files with 739 additions and 635 deletions.
39 changes: 29 additions & 10 deletions examples/Fragmentation/Fragmentation_Movie.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@

import swiftest
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from pathlib import Path
Expand Down Expand Up @@ -75,21 +76,37 @@
for k,v in body_Gmass.items():
body_radius[k] = [((Gmass/GU)/(4./3.*np.pi*density))**(1./3.) for Gmass in v]


# ----------------------------------------------------------------------------------------------------------------------
# Define the animation class that will generate the movies of the fragmentation outcomes
# ----------------------------------------------------------------------------------------------------------------------
figsize = (4,4)


def encounter_combiner(sim):
"""
Combines simulation data with encounter data to produce a dataset that contains the position,
mass, radius, etc. of both. It will interpolate over empty time values to fill in gaps.
"""

# Only keep a minimal subset of necessary data from the simulation and encounter datasets
keep_vars = ['rh','Gmass','radius']
data = sim.data[keep_vars]
enc = sim.enc[keep_vars].load()

# Remove any encounter data at the same time steps that appear in the data to prevent duplicates
t_not_duplicate = ~enc['time'].isin(data['time'])
enc = enc.where(t_not_duplicate,drop=True)

# The following will combine the two datasets along the time dimension, sort the time dimension, and then fill in any time gaps with interpolation
ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time")

return ds

class AnimatedScatter(object):
"""An animated scatter plot using matplotlib.animations.FuncAnimation."""

def __init__(self, sim, animfile, title, style, nskip=1):
if style == "disruption_headon" or style == "hitandrun":
tgood = sim.enc.where(sim.enc['Gmass'] > 0.8 * body_Gmass[style][0]).time
else:
tgood = sim.enc.where(sim.enc['Gmass'] > 0.01 * body_Gmass[style][1]).time

self.ds = sim.enc[['rh','vh','Gmass','radius']].sel(time=tgood).load().interpolate_na(dim="time")
self.ds = encounter_combiner(sim)

nframes = int(self.ds['time'].size)
self.sim = sim
Expand All @@ -101,6 +118,7 @@ def __init__(self, sim, animfile, title, style, nskip=1):
'Central body': 'xkcd:almost black'}

# Set up the figure and axes...
self.figsize = (4,4)
self.fig, self.ax = self.setup_plot()

# Then setup FuncAnimation.
Expand All @@ -110,7 +128,7 @@ def __init__(self, sim, animfile, title, style, nskip=1):
print(f"Finished writing {animfile}")

def setup_plot(self):
fig = plt.figure(figsize=figsize, dpi=300)
fig = plt.figure(figsize=self.figsize, dpi=300)
plt.tight_layout(pad=0)

# Calculate the distance along the y-axis between the colliding bodies at the start of the simulation.
Expand All @@ -120,7 +138,7 @@ def setup_plot(self):

scale_frame = abs(rhy1) + abs(rhy2)
ax = plt.Axes(fig, [0.1, 0.1, 0.8, 0.8])
self.ax_pt_size = figsize[0] * 0.8 * 72 / (2 * scale_frame)
self.ax_pt_size = self.figsize[0] * 0.8 * 72 / (2 * scale_frame)
ax.set_xlim(-scale_frame, scale_frame)
ax.set_ylim(-scale_frame, scale_frame)
ax.set_xticks([])
Expand Down Expand Up @@ -185,6 +203,7 @@ def data_stream(self, frame=0):
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, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
sim.run(dt=1e-4, tstop=4.0e-3, istep_out=1, dump_cadence=0)
sim.run(dt=1e-4, tstop=1.0e-3, istep_out=1, dump_cadence=0)

print("Generating animation")
anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1)
2 changes: 2 additions & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ SET(FAST_MATH_FILES
${SRC}/encounter/encounter_check.f90
${SRC}/encounter/encounter_setup.f90
${SRC}/encounter/encounter_util.f90
${SRC}/encounter/encounter_io.f90
${SRC}/fraggle/fraggle_generate.f90
${SRC}/fraggle/fraggle_io.f90
${SRC}/fraggle/fraggle_placeholder.f90
Expand Down Expand Up @@ -74,6 +75,7 @@ SET(FAST_MATH_FILES
${SRC}/util/util_dealloc.f90
${SRC}/util/util_exit.f90
${SRC}/util/util_fill.f90
${SRC}/util/util_final.f90
${SRC}/util/util_flatten.f90
${SRC}/util/util_get_energy_momentum.f90
${SRC}/util/util_index_array.f90
Expand Down
207 changes: 207 additions & 0 deletions src/encounter/encounter_io.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,207 @@
!! 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.

submodule (encounter_classes) s_encounter_io
use swiftest
contains


module subroutine encounter_io_dump(self, param)
!! author: David A. Minton
!!
!! Dumps the time history of an encounter to file.
implicit none
! Arguments
class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: i

! Most of this is just temporary test code just to get something working. Eventually this should get cleaned up.

do i = 1, self%nframes
if (allocated(self%frame(i)%item)) then
select type(snapshot => self%frame(i)%item)
class is (encounter_snapshot)
param%ioutput = self%tslot(i)
call snapshot%write_frame(self%nc,param)
end select
else
exit
end if
end do


return
end subroutine encounter_io_dump


module subroutine encounter_io_initialize(self, param)
!! author: David A. Minton
!!
!! Initialize a NetCDF encounter file system. This is a simplified version of the main simulation output NetCDF file, but with fewer variables.
use, intrinsic :: ieee_arithmetic
use netcdf
implicit none
! Arguments
class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: nvar, varid, vartype
real(DP) :: dfill
real(SP) :: sfill
logical :: fileExists
character(len=STRMAX) :: errmsg
integer(I4B) :: ndims

associate(nc => self)
dfill = ieee_value(dfill, IEEE_QUIET_NAN)
sfill = ieee_value(sfill, IEEE_QUIET_NAN)

select case (param%out_type)
case("NETCDF_FLOAT")
self%out_type = NF90_FLOAT
case("NETCDF_DOUBLE")
self%out_type = NF90_DOUBLE
end select

! Check if the file exists, and if it does, delete it
inquire(file=nc%enc_file, exist=fileExists)
if (fileExists) then
open(unit=LUN, file=nc%enc_file, status="old", err=667, iomsg=errmsg)
close(unit=LUN, status="delete")
end if

call check( nf90_create(nc%enc_file, NF90_NETCDF4, nc%id), "encounter_io_initialize nf90_create" )

! Dimensions
call check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_initialize nf90_def_dim time_dimid" ) ! Simulation time dimension
call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM , nc%space_dimid), "encounter_io_initialize nf90_def_dim space_dimid" ) ! 3D space dimension
call check( nf90_def_dim(nc%id, nc%id_dimname, param%maxid, nc%id_dimid), "encounter_io_initialize nf90_def_dim id_dimid" ) ! dimension to store particle id numbers
call check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "encounter_io_initialize nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays)

! Dimension coordinates
call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "encounter_io_initialize nf90_def_var time_varid" )
call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "encounter_io_initialize nf90_def_var space_varid" )
call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "encounter_io_initialize nf90_def_var id_varid" )

! Variables
call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%name_varid), "encounter_io_initialize nf90_def_var name_varid" )
call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "encounter_io_initialize nf90_def_var ptype_varid" )
call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rh_varid), "encounter_io_initialize nf90_def_var rh_varid" )
call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%vh_varid), "encounter_io_initialize nf90_def_var vh_varid" )
call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_initialize nf90_def_var Gmass_varid" )
call check( nf90_def_var(nc%id, nc%level_varname, NF90_INT, [nc%id_dimid, nc%time_dimid], nc%level_varid), "encounter_io_initialize nf90_def_var level_varid" )
if (param%lclose) then
call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_initialize nf90_def_var radius_varid" )
end if
if (param%lrotation) then
call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%Ip_varid), "encounter_io_initialize nf90_def_var Ip_varid" )
call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rot_varid), "encounter_io_initialize nf90_def_var rot_varid" )
end if

call check( nf90_inquire(nc%id, nVariables=nvar), "encounter_io_initialize nf90_inquire nVariables" )
do varid = 1, nvar
call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize nf90_inquire_variable" )
select case(vartype)
case(NF90_INT)
call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT" )
case(NF90_FLOAT)
call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT" )
case(NF90_DOUBLE)
call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE" )
case(NF90_CHAR)
call check( nf90_def_var_fill(nc%id, varid, 0, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR" )
end select
end do

! Take the file out of define mode
call check( nf90_enddef(nc%id), "encounter_io_initialize nf90_enddef" )

! Add in the space dimension coordinates
call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize nf90_put_var space" )
end associate

return

667 continue
write(*,*) "Error creating encounter output file. " // trim(adjustl(errmsg))
call util_exit(FAILURE)
end subroutine encounter_io_initialize


module subroutine encounter_io_write_frame(self, nc, param)
!! author: David A. Minton
!!
!! Write a frame of output of an encounter trajectory.
use netcdf
implicit none
! Arguments
class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure
class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: i, tslot, idslot, old_mode, npl, ntp
character(len=NAMELEN) :: charstring

tslot = param%ioutput

call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame nf90_set_fill" )

call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid" )

associate(pl => self%pl, tp => self%tp)
npl = pl%nbody
do i = 1, npl
idslot = pl%id(i)
call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var pl id_varid" )
call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rh_varid" )
call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl vh_varid" )
call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl Gmass_varid" )
select type(pl)
class is (symba_pl)
call check( nf90_put_var(nc%id, nc%level_varid, pl%levelg(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl level_varid" )
end select

if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl radius_varid" )

if (param%lrotation) then
call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl Ip_varid" )
call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rotx_varid" )
end if

charstring = trim(adjustl(pl%info(i)%name))
call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var pl name_varid" )
charstring = trim(adjustl(pl%info(i)%particle_type))
call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var pl particle_type_varid" )
end do

ntp = tp%nbody
do i = 1, ntp
idslot = tp%id(i)
call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var tp id_varid" )
call check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp rh_varid" )
call check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp vh_varid" )

charstring = trim(adjustl(tp%info(i)%name))
call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var tp name_varid" )
charstring = trim(adjustl(tp%info(i)%particle_type))
call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_frame nf90_put_var tp particle_type_varid" )
end do
end associate

call check( nf90_set_fill(nc%id, old_mode, old_mode) )

return
end subroutine encounter_io_write_frame




end submodule s_encounter_io
30 changes: 30 additions & 0 deletions src/encounter/encounter_util.f90
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,36 @@ module subroutine encounter_util_final_list(self)
end subroutine encounter_util_final_list


module subroutine encounter_util_final_snapshot(self)
!! author: David A. Minton
!!
!! Deallocates allocatable arrays in an encounter snapshot
implicit none
! Arguments
type(encounter_snapshot), intent(inout) :: self !! Encounter storage object

if (allocated(self%pl)) deallocate(self%pl)
if (allocated(self%tp)) deallocate(self%tp)
self%t = 0.0_DP

return
end subroutine encounter_util_final_snapshot


module subroutine encounter_util_final_storage(self)
!! author: David A. Minton
!!
!! Deallocates allocatable arrays in an encounter snapshot
implicit none
! Arguments
type(encounter_storage(*)), intent(inout) :: self !! Encounter storage object

call util_final_storage(self%swiftest_storage)

return
end subroutine encounter_util_final_storage


module subroutine encounter_util_resize_list(self, nnew)
!! author: David A. Minton
!!
Expand Down
Loading

0 comments on commit 569b06d

Please sign in to comment.