Commit 6eedab29 authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

clean-up of test.c and changes to some of the initial conditions.


Former-commit-id: 36fd1847b38d8480624737d8f82133f8d22900f0
parent 0c94b5ac
###############################################################################
# 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
import sys
import math
from numpy import *
# Reads the HDF5 output of SWIFT and generates a radial density profile
# of the different physical values.
# Input values?
if len(sys.argv) < 3 :
print "Useage: " , sys.argv[0] , " <filename> <nr. bins>"
exit()
# Get the input arguments
fileName = sys.argv[1];
nr_bins = int( sys.argv[2] );
# Open the file
fileName = sys.argv[1];
file = h5py.File( fileName , 'r' )
# Get the space dimensions.
grp = file[ "/Header" ]
boxsize = grp.attrs[ 'BoxSize' ]
boxsize = boxsize[0]
# Get the particle data
grp = file.get( '/PartType0' )
ds = grp.get( 'Coordinates' )
coords = ds[...]
ds = grp.get( 'Velocities' )
v = ds[...]
ds = grp.get( 'Masses' )
m = ds[...]
ds = grp.get( 'SmoothingLength' )
h = ds[...]
ds = grp.get( 'InternalEnergy' )
u = ds[...]
ds = grp.get( 'ParticleIDs' )
ids = ds[...]
ds = grp.get( 'Density' )
rho = ds[...]
# Set the origin to be the middle of the box
origin = [ boxsize / 2 , boxsize / 2 , boxsize / 2 ]
# Get the maximum radius
r_max = boxsize / 2
# Init the bins
nr_parts = coords.shape[0]
bins_v = zeros( nr_bins )
bins_m = zeros( nr_bins )
bins_h = zeros( nr_bins )
bins_u = zeros( nr_bins )
bins_rho = zeros( nr_bins )
bins_count = zeros( nr_bins )
bins_P = zeros( nr_bins )
# Loop over the particles and fill the bins.
for i in range( nr_parts ):
# Get the box index.
r = coords[i] - origin
r = sqrt( r[0]*r[0] + r[1]*r[1] + r[2]*r[2] )
ind = floor( r / r_max * nr_bins )
# Update the bins
if ( ind < nr_bins ):
bins_count[ind] += 1
bins_v[ind] += sqrt( v[i,0]*v[i,0] + v[i,1]*v[i,1] + v[i,2]*v[i,2] )
bins_m[ind] += m[i]
bins_h[ind] += h[i]
bins_u[ind] += u[i]
bins_rho[ind] += rho[i]
bins_P[ind] += (2.0/3)*u[i]*rho[i]
# Loop over the bins and dump them
print "# bucket left right count v m h u rho"
for i in range( nr_bins ):
# Normalize by the bin volume.
r_left = r_max * i / nr_bins
r_right = r_max * (i+1) / nr_bins
vol = 4/3*math.pi*(r_right*r_right*r_right - r_left*r_left*r_left)
ivol = 1.0 / vol
print "%i %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e" % \
( i , r_left , r_right , \
bins_count[i] * ivol , \
bins_v[i] / bins_count[i] , \
bins_m[i] * ivol , \
bins_h[i] / bins_count[i] , \
bins_u[i] / bins_count[i] , \
bins_rho[i] / bins_count[i] ,
bins_P[i] / bins_count[i] )
......@@ -33,7 +33,7 @@ P_0 = 1.e-5 # Background Pressure
E_0 = 1.e2 # Energy of the explosion
gamma = 5./3. # Gas polytropic index
t = 0.2 # Time of the solution
t = 0.15 # Time of the solution
N = 1000 # Number of radial points
R_max = 3. # Maximal radius
......
......@@ -37,36 +37,21 @@ fileName = "sodShock.hdf5"
#---------------------------------------------------
#Read in high density glass
glass1 = h5py.File("../Glass/glass_200000.hdf5")
pos1 = glass1["/PartType1/Coordinates"][:,:]
pos1 = pos1 / 20. # Particles are in [0:0.5, 0:0.5, 0:0.5]
toRemove = array(0)
for i in range(size(pos1)/3-1):
if pos1[i,0] == pos1[i+1,0]:
if pos1[i,1] == pos1[i+1,1]:
if pos1[i,2] == pos1[i+1,2]:
toRemove = append(toRemove, i)
pos1 = delete(pos1, toRemove, 0)
# glass1 = h5py.File("../Glass/glass_200000.hdf5")
glass1 = h5py.File("glass_001.hdf5")
pos1 = glass1["/PartType0/Coordinates"][:,:]
pos1 = pos1 / 2. # Particles are in [0:0.5, 0:0.5, 0:0.5]
#Read in high density glass
glass2 = h5py.File("../Glass/glass_50000.hdf5")
pos2 = glass2["/PartType1/Coordinates"][:,:]
pos2 = pos2 / 20. # Particles are in [0:0.5, 0:0.5, 0:0.5]
toRemove = array(0)
for i in range(size(pos2)/3-1):
if pos2[i,0] == pos2[i+1,0]:
if pos2[i,1] == pos2[i+1,1]:
if pos2[i,2] == pos2[i+1,2]:
toRemove = append(toRemove, i)
pos2 = delete(pos2, toRemove, 0)
# glass2 = h5py.File("../Glass/glass_50000.hdf5")
glass2 = h5py.File("glass_002.hdf5")
pos2 = glass2["/PartType0/Coordinates"][:,:]
pos2 = pos2 / 2. # Particles are in [0:0.5, 0:0.5, 0:0.5]
#Generate high density region
rho1 = size(pos1) / size(pos2)
rho1 = 1.
coord1 = append(pos1, pos1 + [0, 0.5, 0], 0)
coord1 = append(coord1, pos1 + [0, 0, 0.5], 0)
coord1 = append(coord1, pos1 + [0, 0.5, 0.5], 0)
......@@ -77,7 +62,7 @@ u1 = ones((N1, 1)) * P1 / ((gamma - 1.) * rho1)
m1 = ones((N1, 1)) * 0.5 * rho1 / N1
#Generate low density region
rho2 = 1.
rho2 = 0.25
coord2 = append(pos2, pos2 + [0, 0.5, 0], 0)
coord2 = append(coord2, pos2 + [0, 0, 0.5], 0)
coord2 = append(coord2, pos2 + [0, 0.5, 0.5], 0)
......@@ -133,4 +118,7 @@ ds[()] = ids[:]
file.close()
print "particle 31075: " , coords[31075,:]
print "particle 31076: " , coords[31076,:]
###############################################################################
# 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
import sys
import math
from numpy import *
# Reads the HDF5 output of SWIFT and generates a radial density profile
# of the different physical values.
# Input values?
if len(sys.argv) < 3 :
print "Useage: " , sys.argv[0] , " <filename> <nr. bins>"
exit()
# Get the input arguments
fileName = sys.argv[1];
nr_bins = int( sys.argv[2] );
# Open the file
fileName = sys.argv[1];
file = h5py.File( fileName , 'r' )
# Get the space dimensions.
grp = file[ "/Header" ]
boxsize = grp.attrs[ 'BoxSize' ]
boxsize = boxsize[0]
# Get the particle data
grp = file.get( '/PartType0' )
ds = grp.get( 'Coordinates' )
coords = ds[...]
ds = grp.get( 'Velocities' )
v = ds[...]
# ds = grp.get( 'Mass' )
# m = ds[...]
ds = grp.get( 'SmoothingLength' )
h = ds[...]
ds = grp.get( 'InternalEnergy' )
u = ds[...]
ds = grp.get( 'ParticleIDs' )
ids = ds[...]
ds = grp.get( 'Density' )
rho = ds[...]
# Get the maximum radius
r_max = boxsize
# Init the bins
nr_parts = coords.shape[0]
bins_v = zeros( nr_bins )
bins_m = zeros( nr_bins )
bins_h = zeros( nr_bins )
bins_u = zeros( nr_bins )
bins_rho = zeros( nr_bins )
bins_count = zeros( nr_bins )
bins_P = zeros( nr_bins )
# Loop over the particles and fill the bins.
for i in range( nr_parts ):
# Get the box index.
r = coords[i,0]
ind = floor( r / r_max * nr_bins )
# Update the bins
bins_count[ind] += 1
bins_v[ind] += v[i,0] # sqrt( v[i,0]*v[i,0] + v[i,1]*v[i,1] + v[i,2]*v[i,2] )
# bins_m[ind] += m[i]
bins_h[ind] += h[i]
bins_u[ind] += u[i]
bins_rho[ind] += rho[i]
bins_P[ind] += (2.0/3)*u[i]*rho[i]
# Loop over the bins and dump them
print "# bucket left right count v m h u rho"
for i in range( nr_bins ):
# Normalize by the bin volume.
r_left = r_max * i / nr_bins
r_right = r_max * (i+1) / nr_bins
vol = 4/3*math.pi*(r_right*r_right*r_right - r_left*r_left*r_left)
ivol = 1.0 / vol
print "%i %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e" % \
( i , r_left , r_right , \
bins_count[i] * ivol , \
bins_v[i] / bins_count[i] , \
bins_m[i] * ivol , \
bins_h[i] / bins_count[i] , \
bins_u[i] / bins_count[i] , \
bins_rho[i] / bins_count[i] ,
bins_P[i] / bins_count[i] )
......@@ -29,11 +29,11 @@ from numpy import *
# Parameters
rho_L = 4.
P_L = 1.
rho_L = 1
P_L = 1
v_L = 0.
rho_R = 1.
rho_R = 0.25
P_R = 0.1795
v_R = 0.
......@@ -42,8 +42,8 @@ gamma = 5./3. # Polytropic index
t = 0.12 # Time of the evolution
N = 1000 # Number of points
x_min = -0.15
x_max = 0.15
x_min = -0.25
x_max = 0.25
# ---------------------------------------------------------------
......
......@@ -885,7 +885,7 @@ int main ( int argc , char *argv[] ) {
engine_step( &e , 0 );
if(j % 100 == 0)
write_output(&e);
write_output(&e);
/* Dump the first few particles. */
......@@ -913,7 +913,7 @@ int main ( int argc , char *argv[] ) {
for ( k = 0 ; k < s.nr_parts ; k++ )
if ( s.parts[k].dt < p->dt )
p = &s.parts[k];
printf( "main: particle %lli/%i at [ %e %e %e ] has smallest dt=%.3e (h=%.3e,u=%.3e).\n" ,
printf( "main: particle %lli/%i at [ %e %e %e ] has smallest dt=%e (h=%.3e,u=%.3e).\n" ,
p->id , (int)(p - s.parts) , p->x[0] , p->x[1] , p->x[2] , p->dt , p->h , p->u ); */
/* Output. */
......@@ -982,6 +982,9 @@ int main ( int argc , char *argv[] ) {
printf( "main: particle %lli/%i at [ %e %e %e ] (h=%e) has maximum wcount %.3f.\n" ,
p->id , (int)(p - s.parts) , p->x[0] , p->x[1] , p->x[2] , p->h , p->density.wcount );
/* Write final output. */
write_output( &e );
/* Get the average interactions per particle. */
// icount = 0;
// space_map_parts( &s , &map_icount , &icount );
......
......@@ -64,23 +64,28 @@ int main ( int argc , char *argv[] ) {
int k, N = 100;
struct part p1, p2;
float x, w, dwdx, r2, dx[3] = { 0.0f , 0.0f , 0.0f };
float x, w, dwdx, r2, dx[3] = { 0.0f , 0.0f , 0.0f }, gradw[3];
/* Init the particles. */
for ( k = 0 ; k < 3 ; k++ ) {
p1.a[k] = 0.0f; p1.v[k] = 0.0f; p1.x[k] = 0.0;
p2.a[k] = 0.0f; p2.v[k] = 0.0f; p2.x[k] = 0.0;
}
p1.v[0] = 100.0f;
p1.id = 0; p2.id = 1;
p1.density.wcount = 48.0f; p2.density.wcount = 48.0f;
p1.rho = 1.0f; p1.mass = 9.7059e-4; p1.h = 0.222871287 / 2;
p2.rho = 1.0f; p2.mass = 9.7059e-4; p2.h = 0.222871287 / 2;
p1.force.c = 0.0f; p1.force.balsara = 0.0f;
p2.force.c = 0.0f; p2.force.balsara = 0.0f;
p1.force.c = 0.0040824829f; p1.force.balsara = 0.0f;
p2.force.c = 58.8972740361f; p2.force.balsara = 0.0f;
p1.u = 1.e-5 / ((const_gamma - 1.)*p1.rho);
p2.u = 1.e-5 / ((const_gamma - 1.)*p2.rho) + 100.0f / ( 33 * p2.mass );
p1.force.POrho2 = p1.u * ( const_gamma - 1.0f ) / p1.rho;
p2.force.POrho2 = p2.u * ( const_gamma - 1.0f ) / p2.rho;
/* Dump a header. */
printParticle_single( &p1 );
printParticle_single( &p2 );
printf( "# r a_1 udt_1 a_2 udt_2\n" );
/* Loop over the different radii. */
......@@ -91,28 +96,31 @@ int main ( int argc , char *argv[] ) {
r2 = dx[0]*dx[0];
/* Clear the particle fields. */
/* p1.a[0] = 0.0f; p1.force.u_dt = 0.0f;
p2.a[0] = 0.0f; p2.force.u_dt = 0.0f; */
p1.a[0] = 0.0f; p1.force.u_dt = 0.0f;
p2.a[0] = 0.0f; p2.force.u_dt = 0.0f;
/* Interact the particles. */
// runner_iact_force( r2 , dx , p1.h , p2.h , &p1 , &p2 );
runner_iact_force( r2 , dx , p1.h , p2.h , &p1 , &p2 );
/* Clear the particle fields. */
p1.rho = 0.0f; p1.density.wcount = 0.0f;
p2.rho = 0.0f; p2.density.wcount = 0.0f;
/* p1.rho = 0.0f; p1.density.wcount = 0.0f;
p2.rho = 0.0f; p2.density.wcount = 0.0f; */
/* Interact the particles. */
runner_iact_density( r2 , dx , p1.h , p2.h , &p1 , &p2 );
// runner_iact_density( r2 , dx , p1.h , p2.h , &p1 , &p2 );
/* Evaluate just the kernel. */
x = fabsf( dx[0] ) / p1.h;
kernel_deval( x , &w , &dwdx );
gradw[0] = dwdx / (p1.h*p1.h*p1.h*p1.h) * dx[0] / sqrtf( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
gradw[1] = dwdx / (p1.h*p1.h*p1.h*p1.h) * dx[1] / sqrtf( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
gradw[2] = dwdx / (p1.h*p1.h*p1.h*p1.h) * dx[2] / sqrtf( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
/* Output the results. */
printf( "%.3e %.3e %.3e %.3e %.3e %.3e %.3e\n" ,
// -dx[0] , p1.a[0] , p1.force.u_dt , p2.a[0] , p2.force.u_dt ,
-dx[0] , p1.rho , p1.density.wcount , p2.rho , p2.density.wcount ,
w , dwdx );
printf( "%.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e\n" ,
-dx[0] , p1.a[0] , p1.a[1] , p1.a[2] , p1.force.u_dt ,
/// -dx[0] , p1.rho , p1.density.wcount , p2.rho , p2.density.wcount ,
w , dwdx , gradw[0] , gradw[1] , gradw[2] );
} /* loop over radii. */
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment