Commit e471eb6c authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Fixed Sod shock ICs.

parent 528e6cf8
......@@ -43,14 +43,14 @@ vol = boxSize[0] * boxSize[1] * boxSize[2]
glass1 = h5py.File("glass_001.hdf5")
pos1 = glass1["/PartType0/Coordinates"][:,:]
pos1 = pos1 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
glass_h1 = glass1["/PartType0/SmoothingLength"][:] / factor
#Read in high density glass
# glass2 = h5py.File("../Glass/glass_50000.hdf5")
glass2 = h5py.File("glass_002.hdf5")
pos2 = glass2["/PartType0/Coordinates"][:,:]
pos2 = pos2 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
glass_h2 = glass2["/PartType0/SmoothingLength"][:] / factor
#Generate high density region
rho1 = 1.
......@@ -61,9 +61,10 @@ coord1 = append(coord1, coord1 + [0.25, 0, 0], 0)
# coord1 = append(coord1, pos1 + [0, 0.5, 0.5], 0)
N1 = size(coord1)/3
v1 = zeros((N1, 3))
h1 = ones(N1) * 2.251 * 0.5 * vol / (size(pos1)/3)**(1./3.)
u1 = ones(N1) * P1 / ((gamma - 1.) * rho1)
m1 = ones(N1) * vol * 0.5 * rho1 / N1
h1 = append(glass_h1, glass_h1, 0)
h1 = append(h1, h1, 0)
#Generate low density region
rho2 = 0.25
......@@ -74,9 +75,10 @@ coord2 = append(coord2, coord2 + [0.25, 0, 0], 0)
# coord2 = append(coord2, pos2 + [0, 0.5, 0.5], 0)
N2 = size(coord2)/3
v2 = zeros((N2, 3))
h2 = ones(N2) * 2.251 * 0.5 * vol / (size(pos2)/3)**(1./3.)
u2 = ones(N2) * P2 / ((gamma - 1.) * rho2)
m2 = ones(N2) * vol * 0.5 * rho2 / N2
h2 = append(glass_h2, glass_h2, 0)
h2 = append(h2, h2, 0)
#Merge arrays
numPart = N1 + N2
......@@ -89,8 +91,8 @@ ids = zeros(numPart, dtype='L')
for i in range(1, numPart+1):
ids[i-1] = i
#Final operations
h /= 2
#Final operation since we come from Gadget-2 cubic spline ICs
h /= 1.825752
#File
file = h5py.File(fileName, 'w')
......
......@@ -60,7 +60,7 @@ void printParticle(struct part *parts, struct xpart *xparts, long long int id,
/* Look for the particle. */
for (size_t i = 0; i < N; i++)
if (parts[i].id == id) {
printf("## Particle[%zd]:\n id=%lld", i, parts[i].id);
printf("## Particle[%zd]:\n id=%lld ", i, parts[i].id);
hydro_debug_particle(&parts[i], &xparts[i]);
found = 1;
break;
......@@ -76,12 +76,12 @@ void printgParticle(struct gpart *gparts, long long int id, size_t N) {
/* Look for the particle. */
for (size_t i = 0; i < N; i++)
if (gparts[i].id == -id) {
printf("## gParticle[%zd] (DM) :\n id=%lld", i, -gparts[i].id);
printf("## gParticle[%zd] (DM) :\n id=%lld ", i, -gparts[i].id);
gravity_debug_particle(&gparts[i]);
found = 1;
break;
} else if (gparts[i].id > 0 && gparts[i].part->id == id) {
printf("## gParticle[%zd] (hydro) :\n id=%lld", i, gparts[i].id);
printf("## gParticle[%zd] (hydro) :\n id=%lld ", i, gparts[i].id);
gravity_debug_particle(&gparts[i]);
found = 1;
break;
......
......@@ -23,13 +23,13 @@ __attribute__((always_inline))
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, "
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
"wcount=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
"dS/dt=%.3e, c=%.3e\n"
"divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e] \n "
"v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho,
p->h, (int)p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho,
p->force.pressure, p->entropy, p->entropy_dt, p->force.soundspeed,
p->div_v, p->force.curl_v, p->density.rot_v[0], p->density.rot_v[1],
p->density.rot_v[2], p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end);
......
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