From af0e98cb7da9fa5e44debf97013d34d41694807c Mon Sep 17 00:00:00 2001
From: Rainer Weinberger <rainer.weinberger@cfa.harvard.edu>
Date: Tue, 19 Feb 2019 17:07:11 -0500
Subject: [PATCH] Added example shocktube_1d
---
examples/shocktube_1d/Config.sh | 30 +++
examples/shocktube_1d/Riemann.py | 179 ++++++++++++++++++
examples/shocktube_1d/check.py | 302 +++++++++++++++++++++++++++++++
examples/shocktube_1d/create.py | 108 +++++++++++
examples/shocktube_1d/param.txt | 102 +++++++++++
test.sh | 1 +
6 files changed, 722 insertions(+)
create mode 100644 examples/shocktube_1d/Config.sh
create mode 100644 examples/shocktube_1d/Riemann.py
create mode 100644 examples/shocktube_1d/check.py
create mode 100644 examples/shocktube_1d/create.py
create mode 100644 examples/shocktube_1d/param.txt
diff --git a/examples/shocktube_1d/Config.sh b/examples/shocktube_1d/Config.sh
new file mode 100644
index 0000000..1f9c471
--- /dev/null
+++ b/examples/shocktube_1d/Config.sh
@@ -0,0 +1,30 @@
+#!/bin/bash # this line only there to enable syntax highlighting in this file
+
+## examples/shocktube_1d/Config.sh
+## config file for 1d shocktube probelm
+
+
+#--------------------------------------- Basic operation mode of code
+ONEDIMS # 1d simulation
+REFLECTIVE_X=2 # X Boundary; 1: Reflective, 2: Inflow/Outflow; not active: periodic
+GAMMA=1.4 # Adiabatic index of gas; 5/3 if not set
+
+#--------------------------------------- Mesh motion and regularization
+REGULARIZE_MESH_CM_DRIFT # Mesh regularization; Move mesh generating point towards center of mass to make cells rounder.
+REGULARIZE_MESH_CM_DRIFT_USE_SOUNDSPEED # Limit mesh regularization speed by local sound speed
+REGULARIZE_MESH_FACE_ANGLE # Use maximum face angle as roundness criterion in mesh regularization
+
+#--------------------------------------- Time integration options
+FORCE_EQUAL_TIMESTEPS # variable but global timestep
+
+#---------------------------------------- Single/Double Precision
+DOUBLEPRECISION=1 # Mode of double precision: not defined: single; 1: full double precision 2: mixed, 3: mixed, fewer single precisions; unless short of memory, use 1.
+INPUT_IN_DOUBLEPRECISION # initial conditions are in double precision
+OUTPUT_IN_DOUBLEPRECISION # snapshot files will be written in double precision
+
+#--------------------------------------- Output/Input options
+HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
+
+#--------------------------------------- Testing and Debugging options
+DEBUG # enables core-dumps
+
diff --git a/examples/shocktube_1d/Riemann.py b/examples/shocktube_1d/Riemann.py
new file mode 100644
index 0000000..e7baa27
--- /dev/null
+++ b/examples/shocktube_1d/Riemann.py
@@ -0,0 +1,179 @@
+""" @package examples/shocktube_1d/Riemann.py
+exact solution to Riemann problem
+
+created by Rainer Weinberger, last modified: 19.02.2019
+"""
+
+import numpy as np
+
+
+def f_K(pres, W_K, gamma):
+ """
+ Flux functions for left or right state of a Riemann problem;
+ Auxiliary function for NewtonRaphson.
+
+ \param[in] pres float pressure of central state
+ \param[in] W_K array, shape (3,), left/right hand side state of Riemann problem
+ \param[in] gamma float, adiabatic index of gas
+ """
+ integy_K = W_K[2] / ( gamma - 1.0 ) / W_K[0]
+ a_K = np.sqrt( gamma * W_K[2] / W_K[0] )
+ A = 2.0 / ( gamma + 1.0 ) / W_K[0]
+ B = ( gamma - 1.0 ) / ( gamma + 1.0 ) * W_K[2]
+
+ if pres > W_K[2] : ##shock
+ return ( pres - W_K[2] ) * np.sqrt( A / ( pres + B ) )
+ else: ## rarefaction
+ exponent = ( gamma - 1.0) / 2.0 / gamma
+ return ( 2.0 * a_K / ( gamma - 1.0 ) ) * ( ( pres / W_K[2] )**exponent - 1.0 )
+
+def f_K_prime(pres, W_K, gamma_K):
+ """
+ Derivative of flux functions for left or right state of a Riemann problem;
+ Auxiliary function for NewtonRaphson.
+
+ \param[in] pres float pressure of central state
+ \param[in] W_K array, shape (3,), left/right hand side state of Riemann problem
+ \param[in] gamma float, adiabatic index of gas
+ """
+ integy_K = W_K[2] / ( gamma_K - 1.0 ) / W_K[0]
+ a_K = np.sqrt( gamma_K * W_K[2] / W_K[0] )
+ A = 2.0 / ( gamma_K + 1.0 ) / W_K[0]
+ B = ( gamma_K - 1.0 ) / ( gamma_K + 1.0 ) * W_K[2]
+
+ if pres > W_K[2] : ##shock
+ return np.sqrt( A / ( B + pres ) ) * ( 1.0 - ( pres - W_K[2] ) / ( 2.0 * ( B + pres ) ) )
+ else: ## rarefaction
+ exponent = -( gamma_K + 1.0 ) / 2.0 / gamma_K
+ return ( pres / W_K[2] )**exponent / W_K[0] / a_K
+
+def NewtonRaphson(pres, W_L, W_R, gamma):
+ """
+ NewtonRaphson iteration to determine central pressure of a Riemann problem
+
+ \param[in] W_L array, shape (3,) left hand side density, velocity and pressure
+ \param[in] W_R array, shape (3,) right hand side density, velocity and pressure
+ \param[in] gamma float, adiabatic index of gas
+
+ \return pres pressure of central region
+ """
+ change = 1.0
+ iter = 0
+ while np.abs(change) > 1e-6:
+ iter += 1
+ ## root finding for function
+ function = f_K(pres, W_L, gamma) + f_K(pres, W_R, gamma) + W_R[1] - W_L[1]
+ function_prime = f_K_prime(pres, W_L, gamma) + f_K_prime(pres, W_R, gamma)
+ pres_new = pres - function / function_prime
+ change = 2.0 * ( pres_new - pres ) / ( pres_new + pres )
+ pres = pres_new
+ if(iter > 100):
+ print("ERROR: NewtonRaphson: maximum number of iterations reached! Could not converge!")
+ break
+ return pres
+
+def rarefaction_fan_left(W_K, gamma, xx, time):
+ gm1 = (gamma - 1.0)
+ gp1 = (gamma + 1.0)
+ gamma_ratio = gm1 / gp1
+ a = np.sqrt( gamma * W_K[2] / W_K[0] )
+ rho = W_K[0] * (2.0 / gp1 + gamma_ratio / a * ( W_K[1] - xx / time ) )**( 2.0 / gm1 )## wrong?
+ vel = 2.0 / gp1 * ( a + gm1 / 2.0 * W_K[1] + xx / time )
+ pres = W_K[2] * ( 2.0 / gp1 + gamma_ratio / a * ( W_K[1] - xx / time ) )**( 2.0 * gamma / gm1 )
+ return np.array([rho, vel, pres]).T
+
+def rarefaction_fan_right(W_K, gamma, xx, time):
+ gm1 = (gamma - 1.0)
+ gp1 = (gamma + 1.0)
+ gamma_ratio = gm1 / gp1
+ a = np.sqrt( gamma * W_K[2] / W_K[0] )
+ rho = W_K[0] * (2.0 / gp1 - gamma_ratio / a * ( W_K[1] - xx / time ) )**( 2.0 / gm1 )
+ vel = 2.0 / gp1 * ( -a + gm1 / 2.0 * W_K[1] + xx / time )
+ pres = W_K[2] * ( 2.0 / gp1 - gamma_ratio / a * ( W_K[1] - xx / time ) )**( 2.0 * gamma / gm1 )
+ return np.array([rho, vel, pres]).T
+
+def RiemannProblem(xx, x0, W_L, W_R, gamma, time):
+ """
+ RiemannProblem: Samples the solution of a Riemann problem at a given time at positions xx
+
+ \param[in] xx, array, shape (n,), 1d coordinate of grid
+ \param[in] x0, float initial coordinate of discontinuety
+ \param[in] W_L, array, shape(3,), density, velocity and pressure of left state
+ \param[in] W_R, array, shape(3,), density, velocity and pressure of right state
+ \param[in] gamma, float, adiabatic index of gas
+ \param[in] time, time since initial conditions
+
+ \return xx, W, x_shock; xx: as input; W: array, shape(n,3) densiy, velocity and pressure at pos xx
+ x_shock: postition of charactaeristic features
+ """
+
+ ## some auxiliary quantitites
+ a_L = np.sqrt( gamma * W_L[2] / W_L[0] ) ## sound speeds
+ a_R = np.sqrt( gamma * W_R[2] / W_R[0] )
+ gm1 = gamma - 1.0 ## gamma +- 1
+ gp1 = gamma + 1.0
+ gamma_ratio = gm1 / gp1
+ PosOfCharacteristics = np.zeros(5) ## 0: left most to 4: right most; 2 is rarefaction; if 0 and 1 (3 and 4) have the same value, there is a shock, otherwise rarefaction wave
+
+ ## calculate central pressure; initial guess: arithmetic mean
+ pstar = 0.5 * ( W_L[2] + W_R[2] )
+ pstar = NewtonRaphson( pstar, W_L, W_R, gamma)
+
+ ## calculate central velocity and sound speed
+ ustar = 0.5 * ( W_L[1] + W_R[1] ) + 0.5 * ( f_K(pstar, W_R, gamma) - f_K(pstar, W_L, gamma) )
+ astar_L = a_L * ( pstar / W_L[2] )**( (gamma - 1.0) / (2.0 * gamma) )
+ astar_R = a_R * ( pstar / W_R[2] )**( (gamma - 1.0) / (2.0 * gamma) )
+
+ ## decide on shock or rarefaction on left and right; get velocities of features
+ PosOfCharacteristics[2] = ustar
+ if pstar > W_L[2]: ## shock on left
+ PosOfCharacteristics[0] = W_L[1] - a_L * np.sqrt( gp1 / 2.0 / gamma * pstar / W_L[2] + gm1 / 2.0 / gamma )
+ PosOfCharacteristics[1] = PosOfCharacteristics[0]
+ leftShock = True
+ else: ## rarefaction wave on left
+ PosOfCharacteristics[0] = W_L[1] - a_L ## outer
+ PosOfCharacteristics[1] = ustar - astar_L ## inner
+ if pstar > W_R[2]: ## shock on right
+ PosOfCharacteristics[4] = W_R[1] + a_R * np.sqrt( gp1 / 2.0 / gamma * pstar / W_R[2] + gm1 / 2.0 / gamma )
+ PosOfCharacteristics[3] = PosOfCharacteristics[4]
+ rightShock = True
+ else: ## rarefaction wave on right
+ PosOfCharacteristics[4] = W_R[1] + a_R ## outer
+ PosOfCharacteristics[3] = ustar + astar_R ## inner
+
+ ## self-similar, i.e. pos = velocity * time
+ PosOfCharacteristics[:] *= time
+
+ ## select different regions
+ i_L, = np.where(xx-x0 < PosOfCharacteristics[0])
+ i_star_L, = np.where( (xx-x0 >= PosOfCharacteristics[1]) & (xx-x0 < PosOfCharacteristics[2]) )
+ i_star_R, = np.where( (xx-x0 >= PosOfCharacteristics[2]) & (xx-x0 < PosOfCharacteristics[3]) )
+ i_R, = np.where(xx-x0 >= PosOfCharacteristics[4])
+
+ ## sample solution
+ W = np.zeros( [len(xx), 3], dtype=np.float64 )
+ ## left most and right most state: same as initial conditions
+ W[i_L,:] = W_L[:]
+ W[i_R, :] = W_R[:]
+
+ if PosOfCharacteristics[0] != PosOfCharacteristics[1]: ## rarefaction wave on left
+ i_fan_L, = np.where( (xx-x0 >= PosOfCharacteristics[0]) & (xx-x0 < PosOfCharacteristics[1]) )
+ W[i_fan_L,:] = rarefaction_fan_left(W_L, gamma, xx[i_fan_L]-x0, time)
+ W[i_star_L, 0] = (pstar / W_L[2])**(1./gamma) * W_L[0]
+ else: ## shock on left
+ W[i_star_L, 0] = W_L[0] * ( ( pstar / W_L[2] + gamma_ratio ) / ( gamma_ratio * pstar / W_L[2] + 1.0 ) )
+
+ W[i_star_L, 1] = ustar
+ W[i_star_L, 2] = pstar
+
+ if PosOfCharacteristics[3] != PosOfCharacteristics[4]: ## rarefaction wave on right
+ i_fan_R, = np.where( (xx-x0 >= PosOfCharacteristics[3]) & (xx-x0 < PosOfCharacteristics[4]) )
+ W[i_fan_R, :] = rarefaction_fan_right(W_R, gamma, xx[i_fan_R]-x0, time)
+ W[i_star_R, 0] = (pstar / W_R[2])**(1./gamma) * W_R[0]
+ else: ## shock on right
+ W[i_star_R, 0] = W_R[0] * ( ( pstar / W_R[2] + gamma_ratio ) / ( gamma_ratio * pstar / W_R[2] + 1.0 ) )
+
+ W[i_star_R, 1] = ustar
+ W[i_star_R, 2] = pstar
+
+ return xx, W, PosOfCharacteristics
diff --git a/examples/shocktube_1d/check.py b/examples/shocktube_1d/check.py
new file mode 100644
index 0000000..64436b6
--- /dev/null
+++ b/examples/shocktube_1d/check.py
@@ -0,0 +1,302 @@
+""" @package ./examples/shocktube_1d/check
+Code that checks results of 1d shocktube problem
+
+created by Rainer Weinberger, last modified: 19.02.2019
+"""
+
+""" load libraries """
+import sys # needed for exit codes
+import numpy as np # scientific computing package
+import h5py # hdf5 format
+import matplotlib.pyplot as plt ## needs to be active for plotting!
+
+from Riemann import * ## Riemann-solver
+
+createFigures = True
+forceExitOnError=False ## exits immediately when tolerance is exceeded
+
+""" check functiions """
+def CheckL1Error(Pos, W, W_L, W_R, gamma, position_0, time):
+ """
+ Compare the hydrodynamical quantities of the simulation with the exact
+ solution of the Riemann problem, calculate the L1 error and check whether
+ avarge L1 error is acceptable
+
+ \param[in] Pos: shape (n,) array with 1d positions
+ \param[in] W: shape (n,3) array with density, velocity and pressure of each cell
+ \param[in] W_L: shape (3,) array with left initial condition state (density, velocity and pressure)
+ \param[in] W_R: shape (3,) array with right initial condition state (density, velocity and pressure)
+ \param[in] gamma: adiabatic index of gas
+ \param[in] position_0: initial position of discontinuity
+ \param[in] time: time of snapshot data
+
+ \return status: zero if everything is within tolerances, otherwise 1
+ """
+ xx, W_exact, PosOfCharacteristics = RiemannProblem(Pos, position_0, W_L, W_R, gamma, time)
+
+ norm = np.abs(W_exact)
+ norm[:,1] = 1 ## use absolute error in velocity component
+
+ ## calculate L1 norm
+ delta = np.abs(W_exact - W) / norm
+ L1 = np.average(delta, axis=0)
+
+ ## tolarance value; found empirically, fist order convergence!
+ L1MaxAllowed = 2.0 / np.float(Pos.shape[0])
+
+ if np.any(L1 > L1MaxAllowed):
+ print("CheckL1Error: ERROR: L1 error too large: %g %g %g; tolerance %g"% (L1[0], L1[1], L1[2], L1MaxAllowed) )
+ return 1
+ else:
+ print("CheckL1Error: L1 error fine: %g %g %g; tolerance %g"% (L1[0], L1[1], L1[2], L1MaxAllowed) )
+ return 0
+
+
+def CheckTotalVariation(Pos, W, W_L, W_R, gamma, position_0, time):
+ """
+ Compare the total variation in simulation quantities with the total
+ variation in the analytic solution of the Riemann problem
+
+ \param[in] Pos: shape (n,) array with 1d positions
+ \param[in] W: shape (n,3) array with density, velocity and pressure of each cell
+ \param[in] W_L: shape (3,) array with left initial condition state (density, velocity and pressure)
+ \param[in] W_R: shape (3,) array with right initial condition state (density, velocity and pressure)
+ \param[in] gamma: adiabatic index of gas
+ \param[in] position_0: initial position of discontinuity
+ \param[in] time: time of snapshot data
+
+ \return status: zero if everything is within tolerances, otherwise 1
+ """
+ xx, W_exact, PosOfCharacteristics = RiemannProblem(Pos, position_0, W_L, W_R, gamma, time)
+
+ TotalVariationSim = np.zeros(3)
+ TotalVariationExact = np.zeros(3)
+
+ i_sorted = np.argsort(Pos) ## sorted by position
+ dW = W[i_sorted[1:],:] - W[i_sorted[:-1],:] ## difference of neighbouring cells
+ dW_exact = W_exact[i_sorted[1:],:] - W_exact[i_sorted[:-1],:]
+
+ for i in np.arange(3):
+ i1_pos, = np.where(dW[:,i] >= 0)
+ i1_neg, = np.where(dW[:,i] < 0)
+ TotalVariationSim[i] = np.sum(dW[i1_pos, i], axis=0) - np.sum(dW[i1_neg, i])
+ TotalVariationExact[i] = np.sum(dW_exact[i1_pos, i], axis=0) - np.sum(dW_exact[i1_neg, i])
+
+ MaxRatioAllowed = 1.01
+ if np.any( TotalVariationSim / TotalVariationExact > MaxRatioAllowed):
+ print("CheckTotalVariation: ERROR: TotalVariation Sim/Exact: %g %g %g, tolerance: %g" % (TotalVariationSim[0] / TotalVariationExact[0], TotalVariationSim[1] / TotalVariationExact[1], TotalVariationSim[2] / TotalVariationExact[2], MaxRatioAllowed) )
+ return 1
+ else:
+ print("CheckTotalVariation: TotalVariation Sim/Exact fine: %g %g %g, tolerance: %g" % (TotalVariationSim[0] / TotalVariationExact[0], TotalVariationSim[1] / TotalVariationExact[1], TotalVariationSim[2] / TotalVariationExact[2], MaxRatioAllowed) )
+ return 0
+
+
+def CheckWidthOfDiscontinuities(Pos, W, W_L, W_R, gamma, position_0, time):
+ """
+ Measure the width of the fluid discontinuities in simulation quantities
+ to assess the numerical diffusivity of the scheme
+
+ \param[in] Pos: shape (n,) array with 1d positions
+ \param[in] W: shape (n,3) array with density, velocity and pressure of each cell
+ \param[in] W_L: shape (3,) array with left initial condition state (density, velocity and pressure)
+ \param[in] W_R: shape (3,) array with right initial condition state (density, velocity and pressure)
+ \param[in] gamma: adiabatic index of gas
+ \param[in] position_0: initial position of discontinuity
+ \param[in] time: time of snapshot data
+
+ \return status: zero if everything is within tolerances, otherwise 1
+ """
+ xx, W_exact, PosOfCharacteristics = RiemannProblem(Pos, position_0, W_L, W_R, gamma, time)
+
+ for i, pos_char in enumerate(PosOfCharacteristics):
+ ReturnFlag = 0
+ ## PosOfCharacteristics will give different values in index 0 and 1 (3 and 4)
+ ## if it is a rarefaction; in this case, don't check, otherwise do check
+ if i == 1 or i == 4:
+ continue
+ if i == 0 or i == 3:
+ if PosOfCharacteristics[i] != PosOfCharacteristics[i+1]:
+ continue
+
+ ## get hydrodynamic states left and right of jump
+ xx = np.array( [10.0+pos_char-1e-3, 10.0+pos_char+1e-3] )
+ xx, W_exact, dummy = RiemannProblem(xx, 10.0, W_L, W_R, gamma, time)
+
+ ## get threshold values to measure how many cells are within this interval
+ jump = W_exact[1,:] - W_exact[0,:]
+ percentiles = np.array([W_exact[0,:] + 0.05 * jump, W_exact[0,:] + 0.95 * jump])
+ percentile_05 = np.min( percentiles, axis=0 )
+ percentile_95 = np.max( percentiles, axis=0 )
+ i_sorted = np.argsort( np.abs(Pos-pos_char-position_0) )
+
+ i_low = np.full(3, i_sorted[0], dtype=np.int32)
+ i_high = np.full(3, i_sorted[0], dtype=np.int32)
+
+ ## check by how many cells 5th to 95th percentile of jump are sampled
+ for j in np.arange(3):
+ while W[i_low[j],j] > percentile_05[j] and W[i_low[j],j] < percentile_95[j]:
+ i_low[j] -= 1
+ while W[i_low[j],j] > percentile_05[j] and W[i_low[j],j] < percentile_95[j]:
+ i_high[j] += 1
+
+ ## sufficient for exact Riemann solver
+ MaxNumerOfCells = 4
+
+ if(i == 2):
+ print("CheckWidthOfDiscontinuities: density jump at contact discontinuity resolved by %d cells (5th to 95th precentile), tolerance: %d" \
+ % (i_high[0]-i_low[0], MaxNumerOfCells) )
+ else:
+ print("CheckWidthOfDiscontinuities: density, velocity and pressure jump at shock resolved by %d, %d and %d cells (5th to 95th precentile), tolerance: %d" \
+ % (i_high[0]-i_low[0], i_high[1]-i_low[1], i_high[2]-i_low[2], MaxNumerOfCells) )
+ if np.any(i_high-i_low > MaxNumerOfCells):
+ print("CheckWidthOfDiscontinuities: ERROR: discontinuity too wide!")
+ ReturnFlag += 1
+
+ return ReturnFlag
+
+def PlotSimulationData(Pos, W, W_L, W_R, gamma, position_0, time, simulation_directory):
+ """
+ Plot density, velocity, specific internal energy and pressure of
+ Simulation and exact solution on top
+
+ \param[in] Pos: shape (n,) array with 1d positions
+ \param[in] W: shape (n,3) array with density, velocity and pressure of each cell
+ \param[in] W_L: shape (3,) array with left initial condition state (density, velocity and pressure)
+ \param[in] W_R: shape (3,) array with right initial condition state (density, velocity and pressure)
+ \param[in] gamma: adiabatic index of gas
+ \param[in] position_0: initial position of discontinuity
+ \param[in] time: time of snapshot data
+ \param[in] simulation_directory: path to simulation
+
+ \return status: zero if everything went fine
+ """
+ xx = np.linspace(Pos.min(),Pos.max(),1000)
+ xx, W_exact, dummy = RiemannProblem(xx, position_0, W_L, W_R, gamma, time)
+
+ min = np.zeros(4)
+ max = np.zeros(4)
+ min[0] = np.min( W_exact[:,0] ) - 0.1
+ max[0] = np.max( W_exact[:,0] ) + 0.1
+ min[1] = np.min( W_exact[:,1] ) - 0.1
+ max[1] = np.max( W_exact[:,1] ) + 0.1
+ min[2] = np.min( W_exact[:,2] / W_exact[:,0] / (gamma - 1.0) ) - 0.1
+ max[2] = np.max( W_exact[:,2] / W_exact[:,0] / (gamma - 1.0) ) + 0.1
+ min[3] = np.min( W_exact[:,2] ) - 0.1
+ max[3] = np.max( W_exact[:,2] ) + 0.1
+
+ ## plot:
+ fig, ax = plt.subplots(4, sharex=True, figsize=np.array([6.9,6.0]) )
+ fig.subplots_adjust(left = 0.13, bottom = 0.09,right = 0.98, top = 0.98)
+
+ nskip = 0
+ i_show = np.arange(nskip, len(density)-nskip)
+ ax[0].plot(Pos, W[:,0], ls="", marker='o')
+ ax[0].plot(xx, W_exact[:,0], color="k")
+ ax[1].plot(Pos, W[:,1], ls="", marker='o')
+ ax[1].plot(xx, W_exact[:,1], color="k")
+ ax[2].plot(Pos, W[:,2]/W[:,0]/(gamma - 1.0), ls="", marker='o')
+ ax[2].plot(xx, W_exact[:,2]/W_exact[:,0]/(gamma - 1.0), color="k")
+ ax[3].plot(Pos, W[:,2], ls="", marker='o')
+ ax[3].plot(xx, W_exact[:,2], color="k")
+
+ ax[3].set_xlim([0,20])
+ for i_plot in np.arange(4):
+ ax[i_plot].set_ylim([ min[i_plot], max[i_plot] ])
+
+ ## set labels
+ ax[3].set_xlabel(r"pos")
+ ax[0].set_ylabel(r"density")
+ ax[1].set_ylabel(r"velocity")
+ ax[2].set_ylabel(r"spec. int. energy")
+ ax[3].set_ylabel(r"pressure")
+
+ print(simulation_directory+"/output/figure_%03d.pdf" % (i_file) )
+ fig.savefig(simulation_directory+"/output/figure_%03d.pdf" % (i_file))
+ plt.close(fig)
+
+ return 0
+
+simulation_directory = str(sys.argv[1])
+print("examples/wave_1d/check.py: checking simulation output in directory " + simulation_directory)
+
+Dtype = np.float64 # double precision: np.float64, for single use np.float32
+## open initial conditiions to get parameters
+directory = simulation_directory+""
+filename = "IC.hdf5"
+try:
+ data = h5py.File(directory+filename, "r")
+except:
+ sys.exit(1)
+IC_position = np.array(data["PartType0"]["Coordinates"], dtype = np.float64)
+IC_mass = np.array(data["PartType0"]["Masses"], dtype = np.float64)
+IC_velocity = np.array(data["PartType0"]["Velocities"], dtype = np.float64)
+IC_internalEnergy = np.array(data["PartType0"]["InternalEnergy"], dtype = np.float64)
+NumberOfCells = np.int32(data["Header"].attrs["NumPart_Total"][0])
+DeltaMaxAllowed = 2.0 / np.float(NumberOfCells)
+
+## get parameters of Riemann problem from ICs
+gamma = 1.4 ## needs to be identical to Config.sh =(5./3.) if not specified there!
+dx = IC_position[1,0]-IC_position[0,0] ## assuming a uniform grid
+W_L = np.array([IC_mass[0]/dx, IC_velocity[0,0], (gamma - 1.0)*IC_internalEnergy[0]*IC_mass[0]/dx])
+W_R = np.array([IC_mass[-1]/dx, IC_velocity[-1,0], (gamma - 1.0)*IC_internalEnergy[-1]*IC_mass[-1]/dx])
+i_0, = np.where( (IC_mass[:]/dx == W_R[0]) & (IC_velocity[:,0] == W_R[1]) & ((gamma - 1.0)*IC_internalEnergy[-1]*IC_mass[-1]/dx == W_R[2]) )
+position_0 = 0.5 * (IC_position[i_0[0]-1,0]+IC_position[i_0[0],0]) ## discontinuity at interface position
+
+""" loop over all output files """
+i_file = -1
+ReturnFlag = 0
+while True:
+ i_file += 1
+
+ ## read in data
+ directory = simulation_directory+"/output/"
+ filename = "snap_%03d.hdf5" % (i_file)
+ print("try to open "+directory+filename)
+ ## open hdf5 file
+ try:
+ data = h5py.File(directory+filename, "r")
+ except:
+ break
+ print("analyzing "+ filename)
+
+ ## get data from snapshot
+ time = np.float( data["Header"].attrs["Time"] )
+ position = np.array(data["PartType0"]["Coordinates"], dtype = np.float64)
+ density = np.array(data["PartType0"]["Density"], dtype = np.float64)
+ vel = np.array(data["PartType0"]["Velocities"], dtype = np.float64)
+ internalEnergy = np.array(data["PartType0"]["InternalEnergy"], dtype = np.float64)
+ ## convert to more useful data structure
+ W = np.array([density, vel[:,0], (gamma-1.0)*internalEnergy*density], dtype=np.float64).T ## shape: (n,3)
+
+
+ """ plot data if you want """
+ if createFigures:
+ ReturnFlag += PlotSimulationData(position[:,0], W, W_L, W_R, gamma, position_0, time, simulation_directory)
+ if ReturnFlag > 0 and forceExitOnError:
+ print('ERROR: something went wrong in plot')
+ sys.exit(ReturnFlag)
+
+ """ perform checks """
+ ReturnFlag += CheckL1Error(position[:,0], W, W_L, W_R, gamma, position_0, time)
+ if ReturnFlag > 0 and forceExitOnError:
+ print('ERROR: exceeding tolerance!')
+ sys.exit(ReturnFlag)
+
+ ReturnFlag += CheckTotalVariation(position[:,0], W, W_L, W_R, gamma, position_0, time)
+ if ReturnFlag > 0 and forceExitOnError:
+ print('ERROR: exceeding tolerance!')
+ sys.exit(ReturnFlag)
+
+ ReturnFlag += CheckWidthOfDiscontinuities(position[:,0], W, W_L, W_R, gamma, position_0, time)
+ if ReturnFlag > 0 and forceExitOnError:
+ print('ERROR: exceeding tolerance!')
+ sys.exit(ReturnFlag)
+
+if ReturnFlag == 0:
+ print("check.py: success!")
+else:
+ print("check.py: failed!")
+
+sys.exit(ReturnFlag)
+
+
diff --git a/examples/shocktube_1d/create.py b/examples/shocktube_1d/create.py
new file mode 100644
index 0000000..2854575
--- /dev/null
+++ b/examples/shocktube_1d/create.py
@@ -0,0 +1,108 @@
+""" @package examples/shocktube_1d/create.py
+Code that creates 1d shocktube problem initial conditions
+supposed to be as siple as possible and self-contained
+
+created by Rainer Weinberger, 19.02.2019
+"""
+
+""" load libraries """
+import sys # system specific calls
+import numpy as np # scientific computing package
+import h5py # hdf5 format
+
+simulation_directory = str(sys.argv[1])
+print("examples/shocktube_1d/create.py: creating ICs in directory " + simulation_directory)
+
+FloatType = np.float64 # double precision: np.float64, for single use np.float32
+IntType = np.int32
+
+Boxsize = FloatType(20.0)
+NumberOfCells = IntType(128)
+
+""" initial condition parameters """
+FilePath = simulation_directory + '/IC.hdf5'
+
+## Riemann problem
+position_0 = FloatType(10.0) # initial position of discontinuity
+density_0 = FloatType(1.0)
+density_1 = FloatType(0.125)
+velocity_0 = FloatType(0.0)
+velocity_1 = FloatType(0.0)
+pressure_0 = FloatType(1.0)
+pressure_1 = FloatType(0.1)
+
+gamma = FloatType(1.4) ## note: this has to be consistent with the parameter settings for Arepo!
+gamma_minus1 = gamma - FloatType(1.0)
+uthermal_0 = pressure_0 / gamma_minus1 / density_0
+uthermal_1 = pressure_1 / gamma_minus1 / density_1
+
+""" set up grid: uniform 1d grid """
+## spacing
+dx = Boxsize / FloatType(NumberOfCells)
+## position of first and last cell
+pos_first, pos_last = FloatType(0.5) * dx, Boxsize - FloatType(0.5) * dx
+
+## set up grid
+Pos = np.zeros([NumberOfCells, 3], dtype=FloatType)
+Pos[:,0] = np.linspace(pos_first, pos_last, NumberOfCells, dtype=FloatType)
+Volume = np.full(NumberOfCells, dx, dtype=FloatType)
+
+""" set up hydrodynamical quantitites """
+## left state
+Density = np.full(Pos.shape[0], density_0, dtype=FloatType)
+Velocity = np.zeros(Pos.shape, dtype=FloatType)
+Velocity[:,0] = velocity_0
+Uthermal = np.full(Pos.shape[0], uthermal_0, dtype=FloatType)
+
+## right state
+i_right, = np.where( Pos[:,0] > position_0)
+Density[i_right] = density_1
+Velocity[i_right,0] = velocity_1
+Uthermal[i_right] = uthermal_1
+
+## mass instead of density needed for input
+Mass = Density * Volume
+
+""" write *.hdf5 file; minimum number of fields required by Arepo """
+IC = h5py.File(FilePath, 'w') # open/create file
+
+## create hdf5 groups
+header = IC.create_group("Header") # create header group
+part0 = IC.create_group("PartType0") # create particle group for gas cells
+
+## write header entries
+NumPart = np.array([NumberOfCells, 0, 0, 0, 0, 0], dtype=IntType)
+header.attrs.create("NumPart_ThisFile", NumPart)
+header.attrs.create("NumPart_Total", NumPart)
+header.attrs.create("NumPart_Total_HighWord", np.zeros(6, dtype=IntType) )
+header.attrs.create("MassTable", np.zeros(6, dtype=IntType) )
+header.attrs.create("Time", 0.0)
+header.attrs.create("Redshift", 0.0)
+header.attrs.create("BoxSize", Boxsize)
+header.attrs.create("NumFilesPerSnapshot", 1)
+header.attrs.create("Omega0", 0.0)
+header.attrs.create("OmegaB", 0.0)
+header.attrs.create("OmegaLambda", 0.0)
+header.attrs.create("HubbleParam", 1.0)
+header.attrs.create("Flag_Sfr", 0)
+header.attrs.create("Flag_Cooling", 0)
+header.attrs.create("Flag_StellarAge", 0)
+header.attrs.create("Flag_Metals", 0)
+header.attrs.create("Flag_Feedback", 0)
+if Pos.dtype == np.float64:
+ header.attrs.create("Flag_DoublePrecision", 1)
+else:
+ header.attrs.create("Flag_DoublePrecision", 0)
+
+## write cell data
+part0.create_dataset("ParticleIDs", data=np.arange(1, NumberOfCells+1) )
+part0.create_dataset("Coordinates", data=Pos)
+part0.create_dataset("Masses", data=Mass)
+part0.create_dataset("Velocities", data=Velocity)
+part0.create_dataset("InternalEnergy", data=Uthermal)
+
+## close file
+IC.close()
+
+""" normal exit """
+sys.exit(0)
diff --git a/examples/shocktube_1d/param.txt b/examples/shocktube_1d/param.txt
new file mode 100644
index 0000000..60a410a
--- /dev/null
+++ b/examples/shocktube_1d/param.txt
@@ -0,0 +1,102 @@
+%% examples/shocktube_1d/param.txt
+% parameter file for 1d shocktube problem
+
+%% system options
+MaxMemSize 2500
+CpuTimeBetRestartFile 9000
+TimeLimitCPU 90000
+
+%% initial conditions
+InitCondFile ./run/examples/shocktube_1d/IC
+ICFormat 3
+
+%% output options
+OutputDir ./run/examples/shocktube_1d/output/
+SnapshotFileBase snap
+SnapFormat 3
+NumFilesPerSnapshot 1
+NumFilesWrittenInParallel 1
+
+%% resubmit opitions
+ResubmitOn 0
+ResubmitCommand my-scriptfile
+OutputListFilename ol
+OutputListOn 0
+
+%% simulation mode
+CoolingOn 0
+StarformationOn 0
+PeriodicBoundariesOn 1
+ComovingIntegrationOn 0
+
+%% Cosmological parameters
+Omega0 0.0
+OmegaBaryon 0.0
+OmegaLambda 0.0
+HubbleParam 1.0
+
+%% Simulation parameters
+BoxSize 20.0
+TimeOfFirstSnapshot 0.0
+TimeBetStatistics 0.1
+TimeBegin 0.0
+TimeMax 5.0
+TimeBetSnapshot 1.0
+
+%% Units
+UnitVelocity_in_cm_per_s 1.0
+UnitLength_in_cm 1.0
+UnitMass_in_g 1.0
+GravityConstantInternal 0.0
+
+%% settings for gravity
+ErrTolIntAccuracy 0.1
+ErrTolTheta 0.1
+ErrTolForceAcc 0.1
+
+%% timestepping
+MaxSizeTimestep 0.1
+MinSizeTimestep 1e-5
+TypeOfTimestepCriterion 0
+
+%% moving mesh
+CellShapingSpeed 0.5
+CellMaxAngleFactor 2.25
+TypeOfOpeningCriterion 1
+
+%% hydrodynamics
+CourantFac 0.3
+LimitUBelowThisDensity 0.0
+LimitUBelowCertainDensityToThisValue 0.0
+DesNumNgb 64
+MaxNumNgbDeviation 2
+InitGasTemp 0.0
+MinGasTemp 0.0
+MinEgySpec 0.0
+MinimumDensityOnStartUp 0.0
+
+%% domain decomposition
+MultipleDomains 2
+TopNodeFactor 4
+ActivePartFracForNewDomainDecomp 0.01
+
+%% gravitational softening
+GasSoftFactor 0.01
+SofteningComovingType0 0.1
+SofteningComovingType1 0.1
+SofteningComovingType2 0.1
+SofteningComovingType3 0.1
+SofteningComovingType4 0.1
+SofteningComovingType5 0.1
+SofteningMaxPhysType0 0.1
+SofteningMaxPhysType1 0.1
+SofteningMaxPhysType2 0.1
+SofteningMaxPhysType3 0.1
+SofteningMaxPhysType4 0.1
+SofteningMaxPhysType5 0.1
+SofteningTypeOfPartType0 0
+SofteningTypeOfPartType1 0
+SofteningTypeOfPartType2 0
+SofteningTypeOfPartType3 0
+SofteningTypeOfPartType4 0
+SofteningTypeOfPartType5 0
diff --git a/test.sh b/test.sh
index 81dacc9..7314dc6 100644
--- a/test.sh
+++ b/test.sh
@@ -13,6 +13,7 @@ NUMBER_OF_TASKS=1
TESTS=""
## available 1d test cases
TESTS+="wave_1d "
+TESTS+="shocktube_1d "
## loop over all tests
for TEST in $TESTS
--
GitLab