https://hal.archives-ouvertes.fr/hal-03464411
Raw File
visualization.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys
from decimal import Decimal
import numpy as np
import h5py
import shutil

import matplotlib
#matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from matplotlib import rcParams
from matplotlib import scale as mscale
from matplotlib import transforms as mtransforms
from matplotlib.ticker import StrMethodFormatter
from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size

from mpl_toolkits.mplot3d import Axes3D

from scipy.interpolate import interp2d, interpn

import pyvista as pv
from pyvista import examples


plt.rcParams.update({'figure.max_open_warning': 0})


asinh_scale = 1e2
class AsinhScale(mscale.ScaleBase):
	"""
    Scales data using asinh scale

    The scale function:
      asinh(y) = ln(y+sqrt(y*y+1))

    The inverse scale function:
      sinh(y) = (exp(x)-exp(-x))/2
	"""

	# The scale class must have a member ``name`` that defines the string used
	# to select the scale.  For example, ``gca().set_yscale("asinh")`` would
	# be used to select this scale.
	name = 'asinh'

	def __init__(self, axis, **kwargs):
		#super().__init__(axis)
		mscale.ScaleBase.__init__(self)
		self.thresh = None #thresh
		"""
        Any keyword arguments passed to "set_xscale" and "set_yscale" will
        be passed along to the scale's constructor.
		"""

	def get_transform(self):
		return self.AsinhTransform(self.thresh)
		"""
		Override this method to return a new instance that does the
		actual transformation of the data.
		"""

	def set_default_locators_and_formatters(self, axis):
		pass
		

	class AsinhTransform(mtransforms.Transform):
        # There are two value members that must be defined.
        # ``input_dims`` and ``output_dims`` specify number of input
        # dimensions and output dimensions to the transformation.
        # These are used by the transformation framework to do some
        # error checking and prevent incompatible transformations from
        # being connected together.  When defining transforms for a
        # scale, which are, by definition, separable and have only one
        # dimension, these members should always be set to 1.
		input_dims = output_dims = 1

		def __init__(self, thresh):
			mtransforms.Transform.__init__(self)
			self.thresh = thresh

		def transform_non_affine(self, a):
			"""
            This transform takes an Nx1 ``numpy`` array and returns a
            transformed copy.  Since the range of the Mercator scale
            is limited by the user-specified threshold, the input
            array must be masked to contain only valid values.
            ``matplotlib`` will handle masked arrays and remove the
            out-of-range data from the plot.  Importantly, the
            ``transform`` method *must* return an array that is the
            same shape as the input array, since these values need to
            remain synchronized with values in the other dimension.
			"""
			return np.arcsinh(a*asinh_scale)

		def inverted(self):
			"""
			Override this method so matplotlib knows how to get the
			inverse transform for this transform.
			"""
			return AsinhScale.SinhTransform(self.thresh)

	class SinhTransform(mtransforms.Transform):
		input_dims = output_dims = 1

		def __init__(self, thresh):
			mtransforms.Transform.__init__(self)
			self.thresh = thresh

		def transform_non_affine(self, a):
			return np.sinh(a)/asinh_scale

		def inverted(self):
			return AsinhScale.AsinhTransform(self.thresh)




def read_memprofile(namefile):
	f = open(namefile, 'r')
	lines = f.readlines()

	use_mpi = (("mpirun" in lines[0]) and ("-np" in lines[0]))
	for i in range(len(lines)):
		line = lines[i]
		if (i == 0):
			resident_memory_global = []
			time_global = []
			resident_memory_children = []
			time_children = []
			if use_mpi:
				nchildren = int(line.split(" -np ")[1].split(" ")[0])
				for R in range(nchildren):
					resident_memory_children.append([])
					time_children.append([])
		else:
			data = line.split(" ")
			if ("MEM" in line):
				resident_memory_global.append(float(data[1]))
				time_global.append(float(data[2]))
			elif ("CHLD" in line):
				resident_memory_children[int(data[1])].append(float(data[2]))
				time_children[int(data[1])].append(float(data[3]))

	tinit = time_global[0]
	for i in range(len(time_global)):
		time_global[i] -= tinit
	if use_mpi:
		for R in range(nchildren):
			tinit = time_children[R][0]
			for i in range(len(time_children[R])):
				time_children[R][i] -= tinit

	f.close()

	return [resident_memory_global, time_global, resident_memory_children, time_children]





def plot_relative_error(file_abserror, file_norm, title, outputpng, yscale="linear"):
	f = open(file_abserror, "r")
	lines = f.readlines()
	time = []
	abserror = []
	for line in lines:
		data = line.split(" ")
		new_data = []
		for entry in data:
			if not(entry == ""):
				new_data.append(entry)
		time.append(float(new_data[0]))
		abserror.append(float(new_data[-1]))
	f.close()
	
	f = open(file_norm, "r")
	lines = f.readlines()
	norm = []
	for line in lines:
		data = line.split(" ")
		new_data = []
		for entry in data:
			if not(entry == ""):
				new_data.append(entry)
		norm.append(float(new_data[-1]))
	f.close()
	
	relerror = []
	for i in range(len(abserror)):
		if (norm[i] > 1e-13):
			relerror.append(abserror[i]/norm[i])
		else:
			relerror.append(abserror[i])
	
	
	fs = 30
	plt.rc('xtick',labelsize=fs)
	plt.rc('ytick',labelsize=fs)
	matplotlib.rcParams['font.family'] = 'monospace'

	fig = plt.figure(figsize=(24,24.))
	ax = plt.gca()
	
	im = ax.plot(time, relerror, linewidth=4)
	
	ax.set_yscale(yscale)
	plt.xlim(min(time),max(time))
	
	plt.xlabel("time", fontsize=fs)
	plt.ylabel("Relative error", fontsize=fs)
	plt.title(title, fontsize=fs)
	ttl = ax.title
	ttl.set_position([.5, 1.1])
	ax.grid()
	
	plt.draw()
	if (not(outputpng == '')):
		fig.savefig(outputpng, bbox_extra_artists=(ttl,), bbox_inches='tight')
	
	fig.clf()
	plt.close()


def plot_absolute_error(file_abserror, title, outputpng, yscale="linear"):
	f = open(file_abserror, "r")
	lines = f.readlines()
	time = []
	abserror = []
	for line in lines:
		data = line.split(" ")
		new_data = []
		for entry in data:
			if not(entry == ""):
				new_data.append(entry)
		time.append(float(new_data[0]))
		abserror.append(float(new_data[-1]))
	f.close()
	
	fs = 30
	plt.rc('xtick',labelsize=fs)
	plt.rc('ytick',labelsize=fs)
	matplotlib.rcParams['font.family'] = 'monospace'

	fig = plt.figure(figsize=(24,24.))
	ax = plt.gca()
	
	im = ax.plot(time, abserror, linewidth=4)
	
	ax.set_yscale(yscale)
	plt.xlim(min(time),max(time))
	
	plt.xlabel("time", fontsize=fs)
	plt.ylabel("Absolute error", fontsize=fs)
	plt.title(title, fontsize=fs)
	ttl = ax.title
	ttl.set_position([.5, 1.1])
	ax.grid()
	
	plt.draw()
	if (not(outputpng == '')):
		fig.savefig(outputpng, bbox_extra_artists=(ttl,), bbox_inches='tight')
	
	fig.clf()
	plt.close()


def plot_diagnostic(file_diagnostic, title, ytitle, outputpng, yscale="linear"):
	f = open(file_diagnostic, "r")
	lines = f.readlines()
	time = []
	diag = []
	for line in lines:
		data = line.split(" ")
		new_data = []
		for entry in data:
			if not(entry == ""):
				new_data.append(entry)
		time.append(float(new_data[0]))
		diag.append(float(new_data[-1]))
	f.close()
	
	fs = 30
	plt.rc('xtick',labelsize=fs)
	plt.rc('ytick',labelsize=fs)
	matplotlib.rcParams['font.family'] = 'monospace'

	fig = plt.figure(figsize=(24,24.))
	ax = plt.gca()
	
	im = ax.plot(time, diag, linewidth=4)
	
	ax.set_yscale(yscale)
	plt.xlim(min(time),max(time))

	minval = min(diag)
	maxval = max(diag)
	
	plt.xlabel("time", fontsize=fs)
	plt.ylabel(ytitle, fontsize=fs)
	plt.title(title, fontsize=fs)
	ttl = ax.title
	ttl.set_position([.5, 1.1])
	ax.grid()
	
	plt.draw()
	if (not(outputpng == '')):
		fig.savefig(outputpng, bbox_extra_artists=(ttl,), bbox_inches='tight')
	
	fig.clf()
	plt.close()
	
def plot_diagnostics(files_diagnostics, title, ytitle, outputpng, labels="", yscale="linear"):
	fs = 30
	plt.rc('xtick',labelsize=fs)
	plt.rc('ytick',labelsize=fs)
	matplotlib.rcParams['font.family'] = 'monospace'

	fig = plt.figure(figsize=(24,24.))
	ax = plt.gca()
	i = 0
	for namefile in files_diagnostics:
		f = open(namefile, "r")
		lines = f.readlines()
		time = []
		diag = []
		for line in lines:
			data = line.split(" ")
			new_data = []
			for entry in data:
				if not(entry == ""):
					new_data.append(entry)
			time.append(float(new_data[0]))
			diag.append(float(new_data[-1]))
		f.close()
		"""
		if (yscale == "log"):
			md = min(diag)
			if (md < 0.):
				for k in range(0,len(diag)):
					diag[k] = abs(diag[k])
		if (yscale == "asinh"):
			C = (2*4*np.arctan(1))**6
			#C = 1e6
			for k in range(0,len(diag)):
				diag[k] = diag[k]* C
		"""
		if (labels == ""):
			im = ax.plot(time, diag, linewidth=4)
		else:
			im = ax.plot(time, diag, linewidth=4, label=labels[i])

		ax.set_yscale(yscale)
		plt.xlim(min(time),max(time))
		
		if (i == 0):
			minval = min(diag)
			maxval = max(diag)
		else:
			minval = min([minval, min(diag)])
			maxval = max([maxval, max(diag)])
		i = i+1
	
	plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,e}'))

	
	plt.xlabel("time", fontsize=fs)
	plt.ylabel(ytitle, fontsize=fs)
	plt.title(title, fontsize=fs)
	ttl = ax.title
	ttl.set_position([.5, 1.1])
	leg = plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left", fontsize=fs)
	ax.grid()
	
	plt.draw()
	if (not(outputpng == '')):
		fig.savefig(outputpng, bbox_extra_artists=(ttl,leg,), bbox_inches='tight')
	
	fig.clf()
	plt.close()


def plot_curves(x_list, y_list, xlabel, ylabel, title, outputpng, labels="", markers="", linestyles="", linecolors="", linewidth=2, yscale="linear", xlim="", ylim="", pngheightpx=600, pngwidthpx=800):
	fs = 8
	mydpi = 96

	#fig = plt.figure(figsize=(fs,fs))
	fig = plt.figure(figsize=(pngwidthpx/mydpi, pngheightpx/mydpi), dpi=mydpi)
	
	
	#fs = 30
	plt.rc('xtick',labelsize=fs*1.5)
	plt.rc('ytick',labelsize=fs*1.5)
	matplotlib.rcParams['font.family'] = 'monospace'
	if (markers == ""):
		themarkers = ['o', '.', ',', 'x', '+', 'v', '^', '<', '>', 's', 'd']
	else:
		themarkers = markers
	if (linestyles == ""):
		thelinestyles = ["solid"]*len(x_list)
	else:
		thelinestyles = linestyles
	
	if (linecolors == ""):
		thelinecolors = ["blue", "green", "red", "magenta", "cyan", "yellow", "black"]*((len(x_list)//7)+1)
	else:
		thelinecolors = linecolors

	#fig = plt.figure(figsize=(fs,fs/2))
	ax = plt.gca()

	i = 0
	for i in range(len(x_list)):
		x = x_list[i]
		y = y_list[i]
		if (labels != ""):
			label = labels[i]
		if (markers != ""):
			marker = themarkers[i]
		if (linestyles != ""):
			linestyle = thelinestyles[i]
		if (linecolors != ""):
			color = thelinecolors[i]


		if (labels == ""):
			if (markers == ""):
				if (linestyles == ""):
					if (linecolors == ""):
						im = ax.plot(x, y, linewidth=linewidth)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color)
				else:
					if (linecolors == ""):
						im = ax.plot(x, y, linewidth=linewidth, linestyle=linestyle)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle)
			else:
				if (linestyles == ""):
					if (linecolors == ""):
						im = ax.plot(x, y, linewidth=linewidth, marker=marker, markersize=fs)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, marker=marker, markersize=fs)
				else:
					if (linecolors == ""):
						im = ax.plot(x, y, linewidth=linewidth, linestyle=linestyle, marker=marker, markersize=fs)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, marker=marker, markersize=fs)
		else:
			if (markers == ""):
				if (linestyles == ""):
					if (linecolors == ""):
						im = ax.plot(x, y, linewidth=linewidth, label=label)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, label=label)
				else:
					if (linecolors == ""):
						im = ax.plot(x, y, linewidth=linewidth, linestyle=linestyle, label=label)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, label=label)
			else:
				if (linestyles == ""):
					if (linecolors == ""):
						im = ax.plot(x, y, linewidth=linewidth, marker=marker, markersize=fs, label=label)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, marker=marker, markersize=fs, label=label)
				else:
					if (linecolors == ""):
						im = ax.plot(x, y, linewidth=4, linestyle=linestyle, marker=marker, markersize=fs, label=label)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, marker=marker, markersize=fs, label=label)
		
		ax.set_yscale(yscale)

		if (xlim == ""):
			plt.xlim(min(x),max(x))
		else:
			plt.xlim(min(xlim), max(xlim))

		if (i == 0):
			minval = min(y)
			maxval = max(y)
		else:
			minval = min([minval, min(y)])
			maxval = max([maxval, max(y)])
		i = i+1
		
		if (ylim != ""):
			minval = min(ylim)
			maxval = max(ylim)

		#plt.ylim(minval-abs(minval)*0.05, maxval+abs(maxval)*0.05)
	
	ax.yaxis.set_major_formatter(StrMethodFormatter('{x:+.2e}'))

	
	xl = plt.xlabel(xlabel, fontsize=fs*1.5)
	yl = plt.ylabel(ylabel, fontsize=fs*1.5)
	plt.title(title, fontsize=fs*1.5)
	ttl = ax.title
	ttl.set_position([.5, 1.1])
	
	plt.margins(x=0.5,y=0.5)
	
	#leg = plt.legend(bbox_to_anchor=(1.05, 1.), handlelength=5, loc="upper left", fontsize=fs)
	#leg = plt.legend(bbox_to_anchor=(0.5,-0.1), handlelength=5, loc="upper center", fontsize=fs)
	#leg = plt.legend(handlelength=5, loc="upper right", fontsize=fs)
	if (labels != ""):
		leg = plt.legend(bbox_to_anchor=(0.7,0.5), handlelength=5, loc="center left", fontsize=fs*1.5)
	ax.grid()
	
	plt.subplots_adjust(left=0.15)
	plt.subplots_adjust(right=0.95)
	plt.draw()
	if (not(outputpng == '')):
		if (labels != ""):
			fig.savefig(outputpng, bbox_extra_artists=(ttl,leg,), dpi=mydpi, pad_inches=0.2)
			#fig.savefig(outputpng, bbox_extra_artists=(ttl,leg,), dpi=mydpi, bbox_inches='tight', pad_inches=0.2)
		else:
			fig.savefig(outputpng, bbox_extra_artists=(ttl,), dpi=mydpi)
			#fig.savefig(outputpng, bbox_extra_artists=(ttl,), bbox_inches='tight', pad_inches=0.2)
	
	fig.clf()
	plt.close()
	
	return [minval, maxval]



def exists_in_h5(inputh5, namefield):
	h5file = h5py.File(inputh5, 'r')

	type_H5Group = type(h5file['/'])
	subdirs = namefield.split('/')
	field = subdirs[-1]
	dic = h5file['/']
	type_H5Group = type(dic)

	ans = False
	
	for subdir in subdirs:
		if (subdir == ""):
			pass
		elif (subdir in dic.keys()):
			dic = dic[subdir]
			ans = True
		else:
			ans = False
	
	h5file.close()

	return ans



def extract_from_h5(inputh5, namefield):
	h5file = h5py.File(inputh5, 'r')

	type_H5Group = type(h5file['/'])
	subdirs = namefield.split('/')
	field = subdirs[-1]
	dic = h5file['/']
	type_H5Group = type(dic)
	
	for subdir in subdirs:
		if (subdir == ""):
			pass
		elif (subdir in dic.keys()):
			dic = dic[subdir]
		else:
			print("Error: "+namefield+" cannot be found in file "+inputh5)
			exit()
	if (not(type(dic) is h5py._hl.dataset.Dataset)):
		print("Error: "+namefield+" is not a dataset in file "+inputh5)
		exit()
	
	data = dic[:].copy()
	
	h5file.close()
	
	return data
	



def plot_pseudocolor_2d(data, title, xlim=[0.,1.], ylim=[0.,1.], xlabel="x", ylabel="y", minmax=None, outputpng="", pngheightpx=600, pngwidthpx=800):
	#global fs
	fs = 8
	mydpi = 96

	#fig = plt.figure(figsize=(fs,fs))
	fig = plt.figure(figsize=(pngwidthpx/mydpi, pngheightpx/mydpi), dpi=mydpi)
	plt.rc('xtick',labelsize=fs*1.5)
	plt.rc('ytick',labelsize=fs*1.5)

	im = plt.imshow(data, interpolation="bilinear", extent=xlim[:]+ylim[:], cmap="coolwarm", aspect="equal")
	ax = im.axes
	ax.set_xlabel(xlabel, fontsize=fs*1.5)
	ax.set_ylabel(ylabel, fontsize=fs*1.5)
	ttl = ax.set_title(title, fontsize=fs*1.5)
	ttl.set_position([0.5, 1.1])
	if (minmax != None):
		im.set_clim(minmax)
	
	plt.colorbar(format="%+.2e")

	plt.draw()
	if (not(outputpng == "")):
		#fig.savefig(outputpng, bbox_inches='tight')
		fig.savefig(outputpng, dpi=mydpi)
	



def plot_surface_2d(data, title, zlabel="", xlim=[0.,1.], ylim=[0.,1.], xlabel="x", ylabel="y", minmax=None, outputpng=""):
	global fs
	
	fig = plt.figure(figsize=(fs,fs))
	ax = fig.gca(projection='3d')
	plt.rc('xtick',labelsize=fs*1.5)
	plt.rc('ytick',labelsize=fs*1.5)
	
	(Nx,Ny) = data.shape
	x = np.linspace(xlim[0], xlim[1], Nx)
	y = np.linspace(ylim[0], ylim[1], Ny)
	[X,Y] = np.meshgrid(x,y)
	surf = ax.plot_surface(X, Y, data, cmap="coolwarm", linewidth=0, antialiased=False)
	
	ax.set_xlim3d(min(xlim), max(xlim))
	ax.set_ylim3d(min(ylim), max(ylim))
	if (minmax != None):
		surf.set_clim(minmax)
		ax.set_zlim3d(min(minmax), max(minmax))
	
	"""
	for meth in dir(ax):
		if ("set" in meth):
			print(meth)
	exit()
	"""
	
	ax.set_xlabel(xlabel, fontsize=fs*1.5, labelpad=fs*3)
	ax.set_ylabel(ylabel, fontsize=fs*1.5, labelpad=fs*3)
	"""
	if (zlabel == ""):
		thezlabel = title
		
	else:
		thezlabel = zlabel
	ax.set_zlabel(thezlabel)
	"""
	
	ax.set_title(title, fontsize=fs*1.5)
	ax.tick_params(labelsize=fs*1.5, pad=fs)
	
	ax.view_init(elev=15)
	ax.zaxis._axinfo['juggled'] = (1,1,1)
	
	plt.tight_layout()
	
	
	fig.colorbar(surf, shrink=0.5, aspect=5)

	plt.draw()
	if (not(outputpng == "")):
		fig.savefig(outputpng, bbox_inches='tight')
	
	fig.clf()
	plt.close(fig)

	

def plot_slice_curve(inputh5, p1, p2, projection, linewidth=2, ylabel="", title="", outputpng="", pngheightpx=600, pngwidthpx=800):
	f = extract_from_h5(inputh5, "/distribution/f")
	ndim = len(f.shape)
	nvx = extract_from_h5(inputh5, "/distribution/nvx")[0]
	nvy = extract_from_h5(inputh5, "/distribution/nvy")[0]
	if (ndim == 3):
		nvz = extract_from_h5(inputh5, "/distribution/nvz")[0]
	type_omega = extract_from_h5(inputh5, "/distribution/type_omega")[0]
	L = extract_from_h5(inputh5, "/distribution/L")[0]
	if (type_omega == 1):
		omega = 1.
		u0 = [0.,0.]
		if (ndim == 3):
			u0.append(0.)
	else:
		omega = extract_from_h5(inputh5, "/distribution/omega")[0]
		u0x = extract_from_h5(inputh5, "/distribution/omega/Average_ux")[0]
		u0y = extract_from_h5(inputh5, "/distribution/omega/Average_uy")[0]
		u0 = [u0x, u0y]
		if (ndim == 3):
			u0z = extract_from_h5(inputh5, "/distribution/omega/Average_uz")[0]
			u0.appen(u0z)

	vx = np.linspace(u0[0]-L/omega, u0[0]+L/omega, nvx)
	vy = np.linspace(u0[1]-L/omega, u0[1]+L/omega, nvy)
	if (ndim == 3):
		vz = np.linspace(u0[2]-L/omega, u0[2]+L/omega, nvz)

	ninterp = max(f.shape)
	p1p2_vx = np.linspace(p1[0], p2[0], ninterp)
	p1p2_vy = np.linspace(p1[1], p2[1], ninterp)
	if (ndim == 3):
		p1p2_vz = np.linspace(p1[2], p2[2], ninterp)
	
	if (ndim == 2):
		Pf = interp2d(vx,vy,f)
		slice = np.diag(Pf(p1p2_vx, p1p2_vy))
	elif (ndim == 3):
		points = [[p1p2_vx[i], p1p2_vy[i], p1p2_vz[i]] for i in range(ninterp)]
		slice = interpn((vx,vy,vz), f, points)

	if (outputpng == ""):
		namepng = "test.png"
	else:
		namepng = outputpng
	
	if (projection == "vx"):
		plot_curves([vx], [slice], "vx", ylabel, title, namepng, linecolors=["black"], linewidth=linewidth, pngheightpx=pngheightpx, pngwidthpx=pngwidthpx)
	elif (projection == "vy"):
		plot_curves([vy], [slice], "vy", ylabel, title, namepng, linecolors=["black"], linewidth=linewidth, pngheightpx=pngheightpx, pngwidthpx=pngwidthpx)
	elif (projection == "vz"):
		plot_curves([vz], [slice], "vz", ylabel, title, namepng, linecolors=["black"], linewidth=linewidth, pngheightpx=pngheightpx, pngwidthpx=pngwidthpx)
	else:
		print("Error while using plot_slice_curve: the projection argument must be set to 'x' or 'y'")
		exit()


def plot_slice_pseudocolor(inputh5, projection, porigin, outputpng="", pngheightpx=600, pngwidthpx=800):
	f = extract_from_h5(inputh5, "/distribution/f")
	nvx = extract_from_h5(inputh5, "/distribution/nvx")[0]
	nvy = extract_from_h5(inputh5, "/distribution/nvy")[0]
	nvz = extract_from_h5(inputh5, "/distribution/nvz")[0]
	type_omega = extract_from_h5(inputh5, "/distribution/type_omega")[0]
	L = extract_from_h5(inputh5, "/distribution/L")[0]
	if (type_omega == 1):
		omega = 1.
		u0 = [0.,0.,0.]
	else:
		omega = extract_from_h5(inputh5, "/distribution/omega")[0]
		u0x = extract_from_h5(inputh5, "/distribution/omega/Average_ux")[0]
		u0y = extract_from_h5(inputh5, "/distribution/omega/Average_uy")[0]
		u0z = extract_from_h5(inputh5, "/distribution/omega/Average_uz")[0]
		u0 = [u0x, u0y, u0z]
	
	x = np.linspace(u0[0]-L/omega, u0[0]+L/omega, nvx)
	y = np.linspace(u0[1]-L/omega, u0[1]+L/omega, nvy)
	z = np.linspace(u0[2]-L/omega, u0[2]+L/omega, nvz)
	xlim = [u0[0]-L/omega, u0[0]+L/omega]
	ylim = [u0[1]-L/omega, u0[1]+L/omega]
	zlim = [u0[2]-L/omega, u0[2]+L/omega]

	if (projection == "vxvy"):
		points = []
		for xi in x:
			for yj in y:
				points.append([xi,yj,porigin[2]])
		
		Pf = interpn((x,y,z), f, points)
		Pf = Pf.reshape(nvx,nvy)
		title = "f(vx,vy,pz) with pz = "+str(porigin[2]).format("%.3e")
		plot_pseudocolor_2d(Pf, title, xlim=xlim, ylim=ylim, xlabel="vx", ylabel="vy", minmax=None, outputpng=outputpng, pngheightpx=pngheightpx, pngwidthpx=pngwidthpx)

	elif (projection == "vxvz"):
		points = []
		for xi in x:
			for zk in z:
				points.append([xi,porigin[1],zk])
		
		Pf = interpn((x,y,z), f, points)
		Pf = Pf.reshape(nvx,nvz)
		title = "f(vx,py,vz) with py = "+str(porigin[1]).format("%.3e")
		plot_pseudocolor_2d(Pf, title, xlim=xlim, ylim=zlim, xlabel="vx", ylabel="vz", minmax=None, outputpng=outputpng, pngheightpx=pngheightpx, pngwidthpx=pngwidthpx)
	elif (projection == "vyvz"):
		points = []
		for yj in y:
			for zk in z:
				points.append([porigin[0],yj,zk])

		Pf = interpn((x,y,z), f, points)
		Pf = Pf.reshape(nvy,nvz)
		title = "f(px,vy,vz) with px = "+str(porigin[0]).format("%.3e")
		plot_pseudocolor_2d(Pf, title, xlim=ylim, ylim=zlim, xlabel="vy", ylabel="vz", minmax=None, outputpng=outputpng, pngheightpx=pngheightpx, pngwidthpx=pngwidthpx)



def encode_video(inputpngs, tmppath, outputmovie, slowinglevel=1., forcesave=False):
	if (not(os.path.isdir(tmppath))):
		os.mkdir(tmppath)
		print("Directory "+tmppath+" created")
		creation = True
	else:
		for filetest in os.listdir(tmppath):
			if (((filetest[0:6] == 'frame_') and (filetest[-4:] == '.png')) and not(forcesave)):
				print("Problem with directory "+tmppath+" : it already exists and contains some files named frame_*.png (may induce conflicts). Please choose another directory")
				exit()
		print("Directory "+tmppath+" already exists and is usable for encoding")
		creation = False

	print(str(len(inputpngs)) + ' frames have to be encoded with FFmpeg')

	nframes = 0
	for inputpng in inputpngs:
		shutil.copyfile(inputpng, tmppath+"/frame_"+str(nframes)+".png")
		os.system("cp "+inputpng+" "+tmppath+"/frame_"+str(nframes)+".png")
		nframes = nframes+1

	#cmd = 'ffmpeg -i '+tmppath+'/frame_%d.png -vf scale=640:-1 -filter:v "setpts='+str(slowinglevel)+'*PTS" '+outputmovie+' -loglevel 8'
	cmd = 'ffmpeg -i '+tmppath+'/frame_%d.png -vf "fps=25" '+outputmovie+' -loglevel 8'

	if (forcesave):
		cmd = cmd+' -y'
	print(cmd)
	os.system(cmd)

	for j in range(nframes):
		os.remove(tmppath+"/frame_"+str(j)+".png")
	if (creation):
		os.rmdir(tmppath)

	print('File '+outputmovie+' created')
	
	


def plot_isosurfaces_3d(data, isovalues, title="", xlim=[0.,1.], ylim=[0.,1.], zlim=[0., 1.], xlabel="", ylabel="", zlabel="", zoom_factor=1., opacity=None, outputpng="", pngheight=600, pngwidth=800):
	fs = 16
	
	(nx,ny,nz) = data.shape

	"""
	x, y, z = np.meshgrid(np.linspace(min(xlim), max(xlim), nx), np.linspace(min(ylim), max(ylim), ny), np.linspace(min(zlim), max(zlim), nz))
	
	grid = pv.StructuredGrid(x, y, z)
	grid["f"] = data.flatten(order="F")
	contours = grid.contour(isovalues)
	"""

	grid = pv.UniformGrid()
	grid.dimensions = (nx,ny,nz)
	grid.origin = (min(xlim), min(ylim), min(zlim))
	grid.spacing = ((max(xlim)-min(xlim))/nx, (max(ylim)-min(ylim))/ny, (max(zlim)-min(zlim))/nz)
	grid.point_arrays["f"] = data.flatten(order="F")
	contours = grid.contour(isovalues)

	pv.set_plot_theme('document')
	p = pv.Plotter(off_screen=True)
	#p = pv.Plotter()
	p.add_text(title, font_size=0.5*fs, position='upper_edge')
	if (opacity == None):
		op = [1.]*len(isovalues)
	else:
		if (len(opacity) == len(isovalues)):
			op = opacity
		else:
			print("Error: isovalues and opacity arrays must have same length")
			exit()
	p.add_mesh(contours, opacity=[0.5], show_scalar_bar=False)
	p.add_scalar_bar(title='', vertical=True,  label_font_size=int(fs), position_x=0., position_y=0.1, height=0.8, width=0.1, use_opacity=True)
	bounds = [np.min(xlim), np.max(xlim), np.min(ylim), np.max(ylim), np.min(zlim), np.max(zlim)]
	p.show_bounds(location="outer", grid=True, bounds=bounds, xlabel=xlabel, ylabel=ylabel, zlabel=zlabel, font_size=fs)
	p.camera.zoom(zoom_factor)
	p.show(auto_close=True, screenshot=outputpng, window_size=[pngwidth, pngheight])
	#p.show()


def plot_volume_3d(data, title="", xlim=[0.,1.], ylim=[0.,1.], zlim=[0., 1.], xlabel="", ylabel="", zlabel="", zoom_factor=1., outputpng="", pngheight=600, pngwidth=800):
	fs = 16

	(nx,ny,nz) = data.shape
	#x, y, z = np.meshgrid(np.linspace(min(xlim), max(xlim), nx), np.linspace(min(ylim), max(ylim), ny), np.linspace(min(zlim), max(zlim), nz))

	pv.set_plot_theme('document')

	grid = pv.UniformGrid()
	grid.dimensions = (nx,ny,nz)
	grid.origin = (min(xlim), min(ylim), min(zlim))
	grid.spacing = ((max(xlim)-min(xlim))/nx, (max(ylim)-min(ylim))/ny, (max(zlim)-min(zlim))/nz)
	grid.point_arrays["f"] = data.flatten(order="F")
	
	#p = pv.Plotter()
	p = pv.Plotter(off_screen=True)
	p.add_text(title, font_size=0.5*fs, position='upper_edge')
	p.add_volume(grid, show_scalar_bar=False, opacity="linear")
	p.add_scalar_bar(title='', vertical=True, label_font_size=int(fs), position_x=0., position_y=0.1, height=0.8, width=0.1, use_opacity=False)
	bounds = [np.min(xlim), np.max(xlim), np.min(ylim), np.max(ylim), np.min(zlim), np.max(zlim)]
	p.show_bounds(location="outer", grid=True, xlabel=xlabel, ylabel=ylabel, zlabel=zlabel, font_size=fs)
	p.camera.zoom(zoom_factor)

	#p.show()
	p.show(auto_close=True, screenshot=outputpng, window_size=[pngwidth, pngheight])
	
	

	"""
	grid = pv.StructuredGrid(x, y, z)
	#grid.point_arrays["f"] = data.flatten()
	#grid.origin = (min(xlim), min(ylim), min(zlim))

	#x = np.arange(min(xlim), max(xlim), nx)
	#y = np.arange(min(ylim), max(ylim), ny)
	#z = np.arange(min(zlim), max(zlim), nz)
	#grid = pv.StructuredGrid(x, y, z)
	print(grid)
	#grid["f"] = data.flatten()

	pv.set_plot_theme('document')
	#p = pv.Plotter(off_screen=True)
	p = pv.Plotter()
	p.add_text(title, font_size=0.5*fs, position='upper_edge')
	bounds = [np.min(xlim), np.max(xlim), np.min(ylim), np.max(ylim), np.min(zlim), np.max(zlim)]
	p.show_bounds(location="outer", grid=True, bounds=bounds, xlabel=xlabel, ylabel=ylabel, zlabel=zlabel, font_size=fs)

	p.add_volume(data, show_scalar_bar=False)
	p.add_scalar_bar(title='', vertical=True, label_font_size=int(fs), position_x=0., position_y=0.1, height=0.8, width=0.1, use_opacity=True)
	#p.show_bounds(location="outer", grid=True, bounds=bounds, xlabel=xlabel, ylabel=ylabel, zlabel=zlabel, font_size=fs)

	p.camera.zoom(zoom_factor)

	#p.show(auto_close=True, screenshot=outputpng, window_size=[pngwidth, pngheight])
	p.show()
	"""




"""

def plot_curves_subplot(ax, x_list, y_list, xlabel, ylabel, title, labels=None, markers=None, linestyles=None, linecolors=None, yscale="linear", xlim=None, ylim=None, linewidth=2):
	global fs

	plt.rc('xtick',labelsize=fs)
	plt.rc('ytick',labelsize=fs)
	matplotlib.rcParams['font.family'] = 'monospace'
	if (markers == None):
		themarkers = ['o', '.', ',', 'x', '+', 'v', '^', '<', '>', 's', 'd']
	else:
		themarkers = markers
	if (linestyles == None):
		thelinestyles = ["solid"]*len(x_list)
	else:
		thelinestyles = linestyles
	
	if (linecolors == None):
		thelinecolors = ["blue", "green", "red", "magenta", "cyan", "yellow", "black"]*((len(x_list)//7)+1)
	else:
		thelinecolors = linecolors

	i = 0
	for i in range(len(x_list)):
		x = x_list[i]
		y = y_list[i]
		if (labels != None):
			label = labels[i]
		if (markers != None):
			marker = themarkers[i]
		if (linestyles != None):
			linestyle = thelinestyles[i]
		if (linecolors != None):
			color = thelinecolors[i]


		if (labels == None):
			if (markers == None):
				if (linestyles == None):
					if (linecolors == None):
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color)
				else:
					if (linecolors == None):
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle)
			else:
				if (linestyles == None):
					if (linecolors == None):
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, marker=marker, markersize=fs)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, marker=marker, markersize=fs)
				else:
					if (linecolors == None):
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, marker=marker, markersize=fs)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, marker=marker, markersize=fs)
		else:
			if (markers == None):
				if (linestyles == None):
					if (linecolors == None):
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, label=label)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, label=label)
				else:
					if (linecolors == None):
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, label=label)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, label=label)
			else:
				if (linestyles == None):
					if (linecolors == None):
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, marker=marker, markersize=fs, label=label)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, marker=marker, markersize=fs, label=label)
				else:
					if (linecolors == None):
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, marker=marker, markersize=fs, label=label)
					else:
						im = ax.plot(x, y, linewidth=linewidth, color=color, linestyle=linestyle, marker=marker, markersize=fs, label=label)
		
		ax.set_yscale(yscale)
		if (xlim == None):
			ax.set_xlim(min(x), max(x))
		else:
			ax.set_xlim(min(xlim), max(xlim))
		
		if (i == 0):
			minval = min(y)
			maxval = max(y)
		else:
			minval = min([minval, min(y)])
			maxval = max([maxval, max(y)])
		i = i+1
		
		if (ylim == None):
			ax.set_ylim(minval-abs(minval)*0.5, maxval+abs(maxval)*0.5)
		else:
			minval = min(ylim)
			maxval = max(ylim)
	
	ax.yaxis.set_major_formatter(StrMethodFormatter('{x:+.2e}'))

	ax.set_xlabel(xlabel, fontsize=fs)
	ax.set_ylabel(ylabel, fontsize=fs)
	ttl = ax.set_title(title, fontsize=fs)
	ttl.set_position([0.5, 1.05])

	ax.autoscale(enable=False)
	
	#leg = plt.legend(bbox_to_anchor=(1.05, 1.), handlelength=5, loc="upper left", fontsize=fs)
	#leg = plt.legend(bbox_to_anchor=(0.5,-0.1), handlelength=5, loc="upper center", fontsize=fs)
	#leg = plt.legend(handlelength=5, loc="upper right", fontsize=fs)
	
	
	#leg = ax.legend(bbox_to_anchor=(1,0.5), handlelength=5, loc="center left", fontsize=fs)
	ax.grid()
	
	return im














def pseudocolor_2d_subplot(ax, data, title, xlim=[0.,1.], ylim=[0.,1.], minmax=None):
	global fs
	
	#plt.rc('xtick',labelsize=fs*1.5)
	#plt.rc('ytick',labelsize=fs*1.5)
	im = ax.imshow(data, interpolation="bilinear", extent=xlim[:]+ylim[:], cmap="coolwarm", aspect="equal")
	ax.set_xlabel("x", fontsize=fs)
	ax.set_ylabel("y", fontsize=fs)
	ttl = ax.set_title(title, fontsize=fs)
	ttl.set_position([0.5, 1.05])
	if (minmax != None):
		set_clim(minmax)
	
	return im


def plot_slice_curve_subplot(ax, data, p1, p2, projection, xlim, ylim, ylabel="", title=""):

	nvx = len(data)
	nvy = len(data[0])
	x = np.linspace(min(xlim), max(xlim), nvx)
	y = np.linspace(min(ylim), max(ylim), nvy)
	Pf = interp2d(x, y, data)

	ninterp = max([nvx,nvy])
	p1p2_x = np.linspace(p1[0], p2[0], ninterp)
	p1p2_y = np.linspace(p1[1], p2[1], ninterp)

	slice = np.diag(Pf(p1p2_x, p1p2_y))
	
	if (projection == "x"):
		im = plot_curves_subplot(ax, [x], [slice], "x", ylabel, title, linecolors=["black"], labels=[ylabel])
	elif (projection == "y"):
		im = plot_curves_subplot(ax, [y], [slice], "y", ylabel, title, linecolors=["black"], labels=[ylabel])
	else:
		print("Error while using plot_slice_curve: the projection argument must be set to 'x' or 'y'")
		exit()
	
	return im

	

def multiple_plots(dics):
	global fs

	nrows = 0
	ncols = 0
	if (len(dics) > 0):
		for dic in dics:
			nrows = max([nrows, dic["subplot"][0]])
			ncols = max([ncols, dic["subplot"][1]])
	plt.isinteractive()
	fig, axes = plt.subplots(nrows+1, ncols+1, squeeze=False)
	
	single_plot = ((nrows == 0) and (ncols == 0))

	ims = []
	for dic in dics:
		irow = dic["subplot"][0]
		icol = dic["subplot"][1]
		if single_plot:
			ax = axes
		else:
			ax = axes[irow, icol]
		
		if (dic["type_plot"] == "pseudocolor_2d"):
			im = pseudocolor_2d_subplot(ax, dic["data_2d"], dic["title"], xlim=dic["xlim"], ylim=dic["ylim"], minmax=dic["minmax"])
			#divider = make_axes_locatable(ax)
			#cax = divider.append_axes("right", size="5%", pad=0.2)
			#fig.colorbar(im, ax=ax, format="%+.2e", pad=0.005)
			
		
		elif (dic["type_plot"] == "slice"):
			im = plot_slice_curve_subplot(ax, dic["data_2d"], dic2["segment_proj"][0], dic2["segment_proj"][1], dic["xlabel"], dic["xlim"], dic["ylim"], ylabel=dic["ylabel"], title=dic["title"])
			#leg = im.legend(bbox_to_anchor=(1,0.5), handlelength=5, loc="center left", fontsize=fs)
		
		ims.append(im)
	
	asp = np.diff(axes[0,0].get_xlim())[0] / np.diff(axes[0,0].get_ylim())[0]
	asp /= np.abs(np.diff(axes[1,0].get_xlim())[0] / np.diff(axes[1,0].get_ylim())[0])
	axes[0,0].set_aspect(asp)
	
	plt.tight_layout()

	plt.show()

"""

"""
#inputh5 = "./results/diags_ball_000.h5"
inputh5 = "./results/relaxation/fermions_2d/ball_indicator/hbarmax0.1_L8.0_norescale_residual_RK2_N32/diags_0000.h5"
distrib = extract_from_h5(inputh5, "/distribution/f")
dic1 = {}
dic1["subplot"] = [0,0]
dic1["type_plot"] = "pseudocolor_2d"
dic1["data_2d"] = distrib
dic1["title"] = "distribution"
dic1["xlim"] = [-8.,8.]
dic1["ylim"] = [-8.,8.]
dic1["minmax"] = None

dic2 = {}
dic2["subplot"] = [0,1]
dic2["type_plot"] = "slice"
dic2["segment_proj"] = [[-8.,0.], [8., 0.]]
dic2["xlabel"] = "x"
dic2["ylabel"] = "f(x,0)"
dic2["data_2d"] = distrib
dic2["title"] = "pouet"
dic2["xlim"] = [-8.,8.]
dic2["ylim"] = [-8.,8.]
dic2["minmax"] = None

dics = [dic1, dic2]

fs = 12

multiple_plots(dics)

exit()

"""



# Now that the Scale class has been defined, it must be registered so
# that ``matplotlib`` can find it.
mscale.register_scale(AsinhScale)


fs =16

namefile = "/home/mouton/KINEBEC/code/results/relaxation/fermions_2d/ball_indicator/hbarmax0.1_L6.0_norescale_residual_RK2_convol16-2-2_N64_memprofile.dat"
[mglobal, tglobal, mchild, tchild] = read_memprofile(namefile)

plot_curves(tchild, mchild, "time", "res. mem.", "Evolution of resident memory", "./test.png", labels=["proc 0", "proc 1"], markers="", linestyles="", linecolors="", linewidth=2, yscale="linear", xlim="", ylim="", pngheightpx=600, pngwidthpx=800)



"""
resultsdir = "./results/simu_GDTC"
nfiles = 800
ndigits = 3

outputmovie = resultsdir+"/color/simu_GDTC_color.mp4"
tmppath = resultsdir+"/color/tmp"
inputpngs = []
for i in range(nfiles+1):
	inputpngs.append(resultsdir+"/color/diags_"+str(i).rjust(ndigits, '0')+"_color.png")
encode_video(inputpngs, tmppath, outputmovie)

outputmovie = resultsdir+"/color_clim/simu_GDTC_color_clim.mp4"
tmppath = resultsdir+"/color_clim/tmp"
inputpngs = []
for i in range(nfiles+1):
	inputpngs.append(resultsdir+"/color_clim/diags_"+str(i).rjust(ndigits, '0')+"_color_clim.png")
encode_video(inputpngs, tmppath, outputmovie)

outputmovie = resultsdir+"/surface/simu_GDTC_surface.mp4"
tmppath = resultsdir+"/surface/tmp"
inputpngs = []
for i in range(nfiles+1):
	inputpngs.append(resultsdir+"/surface/diags_"+str(i).rjust(ndigits, '0')+"_surface.png")
encode_video(inputpngs, tmppath, outputmovie)

outputmovie = resultsdir+"/surface_clim/simu_GDTC_surface_clim.mp4"
tmppath = resultsdir+"/surface_clim/tmp"
inputpngs = []
for i in range(nfiles+1):
	inputpngs.append(resultsdir+"/surface_clim/diags_"+str(i).rjust(ndigits, '0')+"_surface_clim.png")
encode_video(inputpngs, tmppath, outputmovie)

"""

"""
file_diagnostic = resultsdir+"/diags_rel_entropy.txt"
title = "Relative entropy"
ytitle = "Relative entropy"
outputpng = resultsdir+"/diags_rel_entropy.png"
plot_diagnostic(file_diagnostic, title, ytitle, outputpng, yscale="log")


minfile = resultsdir+"/diags_minf.txt"
f = open(minfile, "r")
lines = f.readlines()
minf = float(lines[0].split(" ")[-1])
for line in lines:
	minf = min([minf, float(line.split(" ")[-1])])
f.close()

maxfile = resultsdir+"/diags_maxf.txt"
f = open(maxfile, "r")
lines = f.readlines()
maxf = float(lines[0].split(" ")[-1])
for line in lines:
	maxf = max([maxf, float(line.split(" ")[-1])])
f.close()

for i in range(nfiles+1):
	inputh5 = resultsdir+"/diags_"+str(i).rjust(ndigits, '0')+".h5"
	time = extract_from_h5(inputh5, "/distribution/Time")[0]
	distrib = extract_from_h5(inputh5, "/distribution/f")
	title = "Distribution at time t = {:.5e}".format(time)
	
	L = extract_from_h5(inputh5, "/distribution/L")[0]
	type_omega = extract_from_h5(inputh5, "/distribution/type_omega")[0]
	if (type_omega != 1):
		omega = extract_from_h5(inputh5, "/distribution/omega")[0]
		ux = extract_from_h5(inputh5, "/distribution/Average_ux")[0]
		uy = extract_from_h5(inputh5, "/distribution/Average_uy")[0]
	else:
		omega = 1.
		ux = 0.
		uy = 0.
	xlim = [ux-L/omega, ux+L/omega]
	ylim = [uy-L/omega, uy+L/omega]
	
	outputpng = resultsdir+"/diags_"+str(i).rjust(ndigits, '0')+"_surface_clim.png"
	plot_surface_2d(distrib, title, xlim=xlim, ylim=ylim, minmax=[minf,maxf], outputpng=outputpng)
	outputpng = resultsdir+"/diags_"+str(i).rjust(ndigits, '0')+"_color_clim.png"
	plot_pseudocolor_2d(distrib, title, xlim=xlim, ylim=ylim, minmax=[minf,maxf], outputpng=outputpng)
	outputpng = resultsdir+"/diags_"+str(i).rjust(ndigits, '0')+"_surface.png"
	plot_surface_2d(distrib, title, xlim=xlim, ylim=ylim, outputpng=outputpng)
	outputpng = resultsdir+"/diags_"+str(i).rjust(ndigits, '0')+"_color.png"
	plot_pseudocolor_2d(distrib, title, xlim=xlim, ylim=ylim, outputpng=outputpng)

"""

"""
time = []
mass = []
ux = []
uy = []
sxx = []
sxy = []
syx = []
syy = []
T = []
ekin = []
eint = []
entropy = []
rel_entropy = []
mom4 = []
mom6 = []
min_f = []
max_f = []
L1 = []
L2 = []
Linf = []
L1L1 = []
L1L2 = []
L1Linf = []
L2L1 = []
L2L2 = []
L2Linf = []
LinfL1 = []
LinfL2 = []
LinfLinf = []
L1relax = []
L2relax = []
Linfrelax = []
L1L1relax = []
L1L2relax = []
L1Linfrelax = []
L2L1relax = []
L2L2relax = []
L2Linfrelax = []
LinfL1relax = []
LinfL2relax = []
LinfLinfrelax = []
for i in range(nfiles+1):
	inputh5 = resultsdir+"/diags_"+str(i).rjust(ndigits, '0')+".h5"
	time.append(extract_from_h5(inputh5, "/distribution/Time")[0])
	
	mass.append(extract_from_h5(inputh5, "/distribution/Mass")[0])
	ux.append(extract_from_h5(inputh5, "/distribution/Average_ux")[0])
	uy.append(extract_from_h5(inputh5, "/distribution/Average_uy")[0])
	sxx.append(extract_from_h5(inputh5, "/distribution/Stress_energy_xx")[0])
	sxy.append(extract_from_h5(inputh5, "/distribution/Stress_energy_xy")[0])
	syx.append(extract_from_h5(inputh5, "/distribution/Stress_energy_yx")[0])
	syy.append(extract_from_h5(inputh5, "/distribution/Stress_energy_yy")[0])
	T.append(extract_from_h5(inputh5, "/distribution/Temperature")[0])
	ekin.append(extract_from_h5(inputh5, "/distribution/Kinetic_energy")[0])
	eint.append(extract_from_h5(inputh5, "/distribution/Internal_energy")[0])
	entropy.append(extract_from_h5(inputh5, "/distribution/Entropy")[0])
	rel_entropy.append(extract_from_h5(inputh5, "/distribution/Relative_entropy")[0])
	mom4.append(extract_from_h5(inputh5, "/distribution/4th_order_momentum")[0])
	mom6.append(extract_from_h5(inputh5, "/distribution/6th_order_momentum")[0])
	
	min_f.append(extract_from_h5(inputh5, "/distribution/Min_f")[0])
	max_f.append(extract_from_h5(inputh5, "/distribution/Max_f")[0])
	
	L1.append(extract_from_h5(inputh5, "/distribution/L1_norm_space")[0])
	L2.append(extract_from_h5(inputh5, "/distribution/L2_norm_space")[0])
	Linf.append(extract_from_h5(inputh5, "/distribution/Linf_norm_space")[0])
	
	L1L1.append(extract_from_h5(inputh5, "/distribution/L1-L1_norm_time-space")[0])
	L1L2.append(extract_from_h5(inputh5, "/distribution/L1-L2_norm_time-space")[0])
	L1Linf.append(extract_from_h5(inputh5, "/distribution/L1-Linf_norm_time-space")[0])
	L2L1.append(extract_from_h5(inputh5, "/distribution/L2-L1_norm_time-space")[0])
	L2L2.append(extract_from_h5(inputh5, "/distribution/L2-L2_norm_time-space")[0])
	L2Linf.append(extract_from_h5(inputh5, "/distribution/L2-Linf_norm_time-space")[0])
	LinfL1.append(extract_from_h5(inputh5, "/distribution/Linf-L1_norm_time-space")[0])
	LinfL2.append(extract_from_h5(inputh5, "/distribution/Linf-L2_norm_time-space")[0])
	LinfLinf.append(extract_from_h5(inputh5, "/distribution/Linf-Linf_norm_time-space")[0])
	
	L1relax.append(extract_from_h5(inputh5, "/distribution/L1_norm_space_relaxation")[0])
	L2relax.append(extract_from_h5(inputh5, "/distribution/L2_norm_space_relaxation")[0])
	Linfrelax.append(extract_from_h5(inputh5, "/distribution/Linf_norm_space_relaxation")[0])
	
	L1L1relax.append(extract_from_h5(inputh5, "/distribution/L1-L1_norm_time-space_relaxation")[0])
	L1L2relax.append(extract_from_h5(inputh5, "/distribution/L1-L2_norm_time-space_relaxation")[0])
	L1Linfrelax.append(extract_from_h5(inputh5, "/distribution/L1-Linf_norm_time-space_relaxation")[0])
	L2L1relax.append(extract_from_h5(inputh5, "/distribution/L2-L1_norm_time-space_relaxation")[0])
	L2L2relax.append(extract_from_h5(inputh5, "/distribution/L2-L2_norm_time-space_relaxation")[0])
	L2Linfrelax.append(extract_from_h5(inputh5, "/distribution/L2-Linf_norm_time-space_relaxation")[0])
	LinfL1relax.append(extract_from_h5(inputh5, "/distribution/Linf-L1_norm_time-space_relaxation")[0])
	LinfL2relax.append(extract_from_h5(inputh5, "/distribution/Linf-L2_norm_time-space_relaxation")[0])
	LinfLinfrelax.append(extract_from_h5(inputh5, "/distribution/Linf-Linf_norm_time-space_relaxation")[0])
"""
	
	
back to top