Add files to G2-gassphere

examples/G2-gassphere/create_initial_conditions.py

+#!/usr/bin/env python3
+"""
+Code that creates intial conditions for the Evrard collapse
+
+"""
+# load libraries
+import numpy as np  # load numpy
+import h5py    # hdf5 format 
+
+""" initial condition parameters """
+FilePath = 'IC.hdf5'
+
+FloatType = np.float32  # double precision: np.float64, for single use np.float32
+IntType = np.int32
+
+Grid= 14 #size of the grid used to generate particle distribution
+Mtot= 1.0 #total mass of the sphere
+
+x = np.zeros((Grid,Grid,Grid))
+y = np.zeros((Grid,Grid,Grid))
+z = np.zeros((Grid,Grid,Grid))
+vx = np.zeros((Grid,Grid,Grid))
+vy = np.zeros((Grid,Grid,Grid))
+vz = np.zeros((Grid,Grid,Grid))
+m = np.zeros((Grid,Grid,Grid))
+u = np.zeros((Grid,Grid,Grid))
+
+xx = ((np.arange(Grid)-(Grid/2-0.5))/(Grid/2.0))
+
+for i in range(Grid):
+	for j in range(Grid):
+		x[:,i,j] = xx[:]
+		y[i,:,j] = xx[:]
+		z[i,j,:] = xx[:]
+
+x = x.flatten()
+y = y.flatten()
+z = z.flatten()
+
+
+#stretching initial conditons to get 1/r density distribution
+r = np.sqrt(x**2+y**2+z**2)**0.5
+
+x=x*r
+y=y*r
+z=z*r
+
+rad = np.sqrt(x**2+y**2+z**2)
+
+j = np.argwhere(rad < 1.0)
+
+
+
+number_particles = j.size
+
+particle_mass = Mtot / number_particles
+
+print("We use "+str(number_particles)+ " particles")
+
+
+x = x[j]
+y = y[j]
+z = z[j]
+
+
+Pos = np.zeros((number_particles,3), dtype=FloatType)
+Vel = np.zeros((number_particles,3), dtype=FloatType)
+Uthermal = np.zeros((number_particles,3), dtype=FloatType)
+Mass = np.zeros((number_particles,3), dtype=FloatType)
+ids = np.arange(number_particles)
+
+Pos[:,0] = x[:,0]
+Pos[:,1] = y[:,0]
+Pos[:,2] = z[:,0]
+
+
+
+
+Mass[:] = particle_mass
+
+Uthermal[:] = 0.05
+
+#write intial conditions file
+
+IC = h5py.File('IC.hdf5', 'w')
+
+## create hdf5 groups
+header = IC.create_group("Header")
+part0 = IC.create_group("PartType0")
+
+## header entries
+NumPart = np.array([number_particles], dtype=IntType)
+header.attrs.create("NumPart_ThisFile", NumPart)
+header.attrs.create("NumPart_Total", NumPart)
+header.attrs.create("NumPart_Total_HighWord", np.zeros(1, dtype=IntType) )
+header.attrs.create("MassTable", np.zeros(1, dtype=IntType) )
+header.attrs.create("Time", 0.0)
+header.attrs.create("Redshift", 0.0)
+header.attrs.create("BoxSize", 0)
+header.attrs.create("NumFilesPerSnapshot", 1)
+header.attrs.create("Omega0", 0.0)
+header.attrs.create("OmegaB", 0.0)
+header.attrs.create("OmegaLambda", 0.0)
+header.attrs.create("HubbleParam", 1.0)
+header.attrs.create("Flag_Sfr", 0)
+header.attrs.create("Flag_Cooling", 0)
+header.attrs.create("Flag_StellarAge", 0)
+header.attrs.create("Flag_Metals", 0)
+header.attrs.create("Flag_Feedback", 0)
+header.attrs.create("Flag_Entropy_ICs", 0)
+if Pos.dtype == np.float64:
+    header.attrs.create("Flag_DoublePrecision", 1)
+else:
+    header.attrs.create("Flag_DoublePrecision", 0)
+
+## copy datasets
+part0.create_dataset("ParticleIDs", data=ids )
+part0.create_dataset("Coordinates", data=Pos)
+
+part0.create_dataset("Masses", data=Mass)
+
+part0.create_dataset("Velocities", data=Vel)
+part0.create_dataset("InternalEnergy", data=Uthermal)
+
+IC.close()
\ No newline at end of file

examples/G2-gassphere/plot_energy_gassphere.py

+#!/usr/bin/env python3
+"""
+Code that plots the energy evolution for the Evrard collapse
+
+"""
+# load libraries
+import numpy as np
+import matplotlib.pyplot as plt
+import csv
+
+
+time = np.zeros(61)
+thermal_energy = np.zeros(61)
+potential_energy = np.zeros(61)
+kinetic_energy = np.zeros(61)
+
+with open('output/energy.txt') as csvfile:
+    readCSV = csv.reader(csvfile)
+    line = 0
+    for row in readCSV:
+        values = row[0].split()
+        if len(values) != 8:
+            break
+        time[line] = values[0]
+        thermal_energy[line] = values[1]
+        potential_energy[line] = values[2]
+        kinetic_energy[line] = values[3]
+
+
+        line = line + 1
+
+total_energy = thermal_energy + potential_energy + kinetic_energy
+
+plt.plot(time, thermal_energy,label="Thermal energy")
+plt.plot(time, potential_energy,label="Potential energy")
+plt.plot(time, kinetic_energy,label="Kinetic energy")
+plt.plot(time, total_energy,label="Total energy")
+plt.plot(time, (total_energy-total_energy[0])/total_energy[0])
+plt.ylabel('Energy')
+plt.xlabel('Time')
+plt.legend()
+plt.show()
\ No newline at end of file