Added example shocktube_1d

examples/shocktube_1d/Config.sh

+#!/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
+

examples/shocktube_1d/Riemann.py

+""" @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

examples/shocktube_1d/check.py

+""" @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)
+
+

examples/shocktube_1d/create.py

+""" @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