makeIC_fcc.py 4 KB
Newer Older
Pedro Gonnet's avatar
Pedro Gonnet committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
###############################################################################
 # This file is part of SWIFT.
 # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
 #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 # 
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 # 
 # This program 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 Lesser General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 # 
 ##############################################################################

import h5py
import random
from numpy import *

# Generates a swift IC file for the Sedov blast test in a periodic cubic box

# Parameters
periodic= 1      # 1 For periodic box
29
boxSize = 5.
Pedro Gonnet's avatar
Pedro Gonnet committed
30
31
32
33
L = 64            # Number of particles boxes along one axis
rho = 1.          # Density
P = 1.e-5         # Pressure
E0= 1.e2          # Energy of the explosion
34
pert = 0.025
Pedro Gonnet's avatar
Pedro Gonnet committed
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
gamma = 5./3.     # Gas adiabatic index
fileName = "sedov.hdf5" 


#---------------------------------------------------
numPart = 4*(L**3) 
mass = boxSize**3 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
off = array( [ [ 0.0 , 0.0 , 0.0 ] , [ 0.0 , 0.5 , 0.5 ] , [ 0.5 , 0.0 , 0.5 ] , [ 0.5 , 0.5 , 0.0 ] ] );
hbox = boxSize / L

# if L%2 == 0:
#     print "Number of particles along each dimension must be odd."
#     exit()

#Generate particles
coords = zeros((numPart, 3))
v      = zeros((numPart, 3))
53
54
55
56
m      = zeros(numPart)
h      = zeros(numPart)
u      = zeros(numPart)
ids    = zeros(numPart, dtype='L')
Pedro Gonnet's avatar
Pedro Gonnet committed
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

for i in range(L):
    for j in range(L):
        for k in range(L):
            x = (i + 0.25) * hbox
            y = (j + 0.25) * hbox
            z = (k + 0.25) * hbox
            for ell in range(4): 
                index = 4*(i*L*L + j*L + k) + ell
                coords[index,0] = x + off[ell,0] * hbox
                coords[index,1] = y + off[ell,1] * hbox
                coords[index,2] = z + off[ell,2] * hbox
                v[index,0] = 0.
                v[index,1] = 0.
                v[index,2] = 0.
                m[index] = mass
                h[index] = 2.251 / 2 * hbox
                u[index] = internalEnergy
                ids[index] = index
                if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox:
                    u[index] = u[index] + E0 / (28. * mass)
                    print "Particle " , index , " set to detonate."
                coords[index,0] += random.random() * pert * hbox
                coords[index,1] += random.random() * pert * hbox
                coords[index,2] += random.random() * pert * hbox


#--------------------------------------------------

#File
file = h5py.File(fileName, 'w')

# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]

#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic

#Particle group
grp = file.create_group("/PartType0")
106
107
108
109
110
111
grp.create_dataset('Coordinates', data=coords, dtype='d')
grp.create_dataset('Velocities', data=v, dtype='f')
grp.create_dataset('Masses', data=m, dtype='f')
grp.create_dataset('SmoothingLength', data=h, dtype='f')
grp.create_dataset('InternalEnergy', data=u, dtype='f')
grp.create_dataset('ParticleIDs', data=ids, dtype='L')
Pedro Gonnet's avatar
Pedro Gonnet committed
112
113

file.close()