Skip to content
Snippets Groups Projects
Commit f8035b34 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Removed EoS constraint in HLLC Riemann solver. Added missing parameter for...

Removed EoS constraint in HLLC Riemann solver. Added missing parameter for SineWavePotential_1D. Played around with ICs and plot script of 1D sine wave potential.
parent c904a715
Branches
Tags
1 merge request!588Gizmo mfm clean
......@@ -31,12 +31,12 @@ cs2 = 2.*uconst/3.
A = 10.
fileName = "sineWavePotential.hdf5"
numPart = 100
numPart = 1000
boxSize = 1.
coords = np.zeros((numPart, 3))
v = np.zeros((numPart, 3))
m = np.zeros(numPart) + 1.
m = np.zeros(numPart) + 1000. / numPart
h = np.zeros(numPart) + 2./numPart
u = np.zeros(numPart) + uconst
ids = np.arange(numPart, dtype = 'L')
......
......@@ -23,8 +23,13 @@
import numpy as np
import h5py
import sys
import matplotlib
matplotlib.use("Agg")
import pylab as pl
pl.rcParams["figure.figsize"] = (12, 10)
pl.rcParams["text.usetex"] = True
# these should be the same as in makeIC.py
uconst = 20.2615290634
cs2 = 2.*uconst/3.
......@@ -39,15 +44,20 @@ fileName = sys.argv[1]
file = h5py.File(fileName, 'r')
coords = np.array(file["/PartType0/Coordinates"])
rho = np.array(file["/PartType0/Density"])
P = np.array(file["/PartType0/Pressure"])
u = np.array(file["/PartType0/InternalEnergy"])
agrav = np.array(file["/PartType0/GravAcceleration"])
m = np.array(file["/PartType0/Masses"])
vs = np.array(file["/PartType0/Velocities"])
ids = np.array(file["/PartType0/ParticleIDs"])
x = np.linspace(0., 1., 1000)
rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
P = cs2*rho
a = A * np.sin(2. * np.pi * x)
apart = A * np.sin(2. * np.pi * coords[:,0])
tkin = -0.5 * np.dot(apart, coords[:,0])
print tkin, 0.5 * np.dot(m, vs[:,0]**2)
ids_reverse = np.argsort(ids)
......@@ -65,13 +75,38 @@ for i in range(len(P)):
corr = 1.
idxp1 = ids_reverse[ip1]
idxm1 = ids_reverse[im1]
gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0])
gradP[i] = (P[idxp1] - P[idxm1])/(coords[idxp1,0] - coords[idxm1,0])
fig, ax = pl.subplots(2, 2, sharex = True)
ax[0][0].plot(coords[:,0], rho, ".", label = "gizmo-mfm")
ax[0][0].plot(x, rho_x, "-", label = "stable solution")
ax[0][0].set_ylabel("$\\rho{}$ (code units)")
ax[0][0].legend(loc = "best")
ax[0][1].plot(x, a, "-", label = "$\\nabla{}\\Phi{}$ external")
ax[0][1].plot(coords[:,0], gradP/rho, ".",
label = "$\\nabla{}P/\\rho{}$ gizmo-mfm")
ax[0][1].set_ylabel("force (code units)")
ax[0][1].legend(loc = "best")
ax[1][0].axhline(y = uconst, label = "isothermal $u$", color = "k",
linestyle = "--")
ax[1][0].plot(coords[:,0], u, ".", label = "gizmo-mfm")
ax[1][0].set_ylabel("$u$ (code units)")
ax[1][0].set_xlabel("$x$ (code units)")
ax[1][0].legend(loc = "best")
#ax[1][1].plot(coords[:,0], m, "y.")
#ax[1][1].set_ylabel("$m_i$ (code units)")
#ax[1][1].set_xlabel("$x$ (code units)")
fig, ax = pl.subplots(2, 2)
ax[1][1].axhline(y = 0., label = "target", color = "k",
linestyle = "--")
ax[1][1].plot(coords[:,0], vs[:,0], ".", label = "gizmo-mfm")
ax[1][1].set_ylabel("$v_x$ (code units)")
ax[1][1].set_xlabel("$x$ (code units)")
ax[1][1].legend(loc = "best")
ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
ax[0][0].plot(x, rho_x, "g-")
ax[0][1].plot(coords[:,0], gradP/rho, "b.")
ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
ax[1][1].plot(coords[:,0], m, "y.")
pl.tight_layout()
pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))
......@@ -10,7 +10,7 @@ InternalUnitSystem:
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 10. # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_min: 1e-8 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
......@@ -36,3 +36,6 @@ InitialConditions:
SineWavePotential:
amplitude: 10.
timestep_limit: 1.
EoS:
isothermal_internal_energy: 20.2615290634
......@@ -33,11 +33,6 @@
#include "riemann_checks.h"
#include "riemann_vacuum.h"
#ifndef EOS_IDEAL_GAS
#error \
"The HLLC Riemann solver currently only supports and ideal gas equation of state. Either select this equation of state, or try using another Riemann solver!"
#endif
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
const float *WL, const float *WR, const float *n, const float *vij,
float *totflux) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment