Skip to content
Snippets Groups Projects
Commit dde8809a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Modified all IC generation scripts to make sure all IDs are >0

parent 6335b70f
Branches
Tags
2 merge requests!136Master,!116Basic implementation of gparts
...@@ -64,7 +64,7 @@ v = inputFile["/PartType0/Velocities"][:,:] ...@@ -64,7 +64,7 @@ v = inputFile["/PartType0/Velocities"][:,:]
m = inputFile["/PartType0/Masses"][:] m = inputFile["/PartType0/Masses"][:]
h = inputFile["/PartType0/SmoothingLength"][:] h = inputFile["/PartType0/SmoothingLength"][:]
u = inputFile["/PartType0/InternalEnergy"][:] u = inputFile["/PartType0/InternalEnergy"][:]
ids = np.array(range(np.size(u)), dtype='L') ids = np.array(range(np.size(u)), dtype='L') + 1
# Downsample # Downsample
print "Downsampling..." print "Downsampling..."
......
...@@ -111,6 +111,6 @@ ds[()] = h ...@@ -111,6 +111,6 @@ ds[()] = h
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f') ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u ds[()] = u
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L') ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids ds[()] = ids + 1
file.close() file.close()
...@@ -82,7 +82,7 @@ for i in range(L): ...@@ -82,7 +82,7 @@ for i in range(L):
else: else:
P = P + 3. + 4.*log(2.) P = P + 3. + 4.*log(2.)
u[index] = P / ((gamma - 1.)*rho) u[index] = P / ((gamma - 1.)*rho)
ids[index] = partId ids[index] = partId + 1
partId = partId + 1 partId = partId + 1
......
...@@ -103,6 +103,6 @@ ds[()] = h ...@@ -103,6 +103,6 @@ ds[()] = h
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f') ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u ds[()] = u
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L') ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids ds[()] = ids + 1
file.close() file.close()
...@@ -69,7 +69,7 @@ for i in range(L): ...@@ -69,7 +69,7 @@ for i in range(L):
m[index] = mass m[index] = mass
h[index] = 1.1255 * boxSize / L h[index] = 1.1255 * boxSize / L
u[index] = internalEnergy u[index] = internalEnergy
ids[index] = index ids[index] = index + 1
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L: if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L:
u[index] = u[index] + E0 / (33. * mass) u[index] = u[index] + E0 / (33. * mass)
print "Particle " , index , " set to detonate." print "Particle " , index , " set to detonate."
......
...@@ -72,7 +72,7 @@ for i in range(L): ...@@ -72,7 +72,7 @@ for i in range(L):
m[index] = mass m[index] = mass
h[index] = 1.1255 * hbox h[index] = 1.1255 * hbox
u[index] = internalEnergy u[index] = internalEnergy
ids[index] = index ids[index] = index + 1
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox: if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox:
u[index] = u[index] + E0 / (28. * mass) u[index] = u[index] + E0 / (28. * mass)
print "Particle " , index , " set to detonate." print "Particle " , index , " set to detonate."
......
...@@ -86,7 +86,7 @@ h = append(h1, h2,0) ...@@ -86,7 +86,7 @@ h = append(h1, h2,0)
u = append(u1, u2,0) u = append(u1, u2,0)
m = append(m1, m2,0) m = append(m1, m2,0)
ids = zeros(numPart, dtype='L') ids = zeros(numPart, dtype='L')
for i in range(numPart): for i in range(1, numPart+1):
ids[i] = i ids[i] = i
#Final operations #Final operations
......
...@@ -86,7 +86,7 @@ u = zeros(1) ...@@ -86,7 +86,7 @@ u = zeros(1)
ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1)) ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L') ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids ds[()] = ids + 1
x = ids % L; x = ids % L;
y = ((ids - x) / L) % L; y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2; z = (ids - x - L * y) / L**2;
......
...@@ -98,7 +98,7 @@ for n in range(n_iterations): ...@@ -98,7 +98,7 @@ for n in range(n_iterations):
u = zeros(1) u = zeros(1)
ids = linspace(offset, offset+N, N, endpoint=False).reshape((N,1)) ids = linspace(offset, offset+N, N, endpoint=False).reshape((N,1))
ds_id[offset:offset+N] = ids ds_id[offset:offset+N] = ids + 1
x = ids % L; x = ids % L;
y = ((ids - x) / L) % L; y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2; z = (ids - x - L * y) / L**2;
...@@ -130,10 +130,8 @@ u = full((remainder, 1), internalEnergy) ...@@ -130,10 +130,8 @@ u = full((remainder, 1), internalEnergy)
ds_u[offset:offset+remainder] = u ds_u[offset:offset+remainder] = u
u = zeros(1) u = zeros(1)
print "Done", offset+remainder,"/", numPart
ids = linspace(offset, offset+remainder, remainder, endpoint=False).reshape((remainder,1)) ids = linspace(offset, offset+remainder, remainder, endpoint=False).reshape((remainder,1))
ds_id[offset:offset+remainder] = ids ds_id[offset:offset+remainder] = ids + 1
x = ids % L; x = ids % L;
y = ((ids - x) / L) % L; y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2; z = (ids - x - L * y) / L**2;
...@@ -143,6 +141,7 @@ coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L) ...@@ -143,6 +141,7 @@ coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L) coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
ds_x[offset:offset+remainder,:] = coords ds_x[offset:offset+remainder,:] = coords
print "Done", offset+remainder,"/", numPart
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment