diff --git a/examples/SineWavePotential_1D/makeIC.py b/examples/SineWavePotential_1D/makeIC.py
index 321af0714219dfa0c7cbb3d80845881dcbb8416d..afbf1bc0fa47a27677cb9c5645d439432bd9fd9a 100644
--- a/examples/SineWavePotential_1D/makeIC.py
+++ b/examples/SineWavePotential_1D/makeIC.py
@@ -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')
diff --git a/examples/SineWavePotential_1D/plotSolution.py b/examples/SineWavePotential_1D/plotSolution.py
index 65e981e4648fe0fe5d1da6cf3e753fb8a34f0fb4..3bb889aaabd3cdac0274afb09647d0e3aebb02cc 100644
--- a/examples/SineWavePotential_1D/plotSolution.py
+++ b/examples/SineWavePotential_1D/plotSolution.py
@@ -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]))
diff --git a/examples/SineWavePotential_1D/sineWavePotential.yml b/examples/SineWavePotential_1D/sineWavePotential.yml
index 9662841032a12d48870f12de8c1bcfcd579a6b42..e6285785099f10902ea60b21334a0ad26c0593de 100644
--- a/examples/SineWavePotential_1D/sineWavePotential.yml
+++ b/examples/SineWavePotential_1D/sineWavePotential.yml
@@ -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
diff --git a/src/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index d7e376fe667dd3cfeac58d5820053057ab928d3e..2d748d022b40284e399992160fc6c82433bcf9d2 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -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) {