Peano 4
Loading...
Searching...
No Matches
plots-time-integration-test.py
Go to the documentation of this file.
1import numpy as np
2import glob, os
3
4import matplotlib
5import matplotlib.pyplot as plt
6from matplotlib import gridspec
7from matplotlib.colors import LogNorm
8from matplotlib import mathtext
9import math
10from scipy import stats
11from tkinter import Tcl
13
14#Plot Setup
15font = {'size':12}
16plt.rc('font', **font)
17plt.rc('font', family='serif')
18plt.rc('text', usetex=False)
19
20# DO NOT TOUCH ANYTHING AFTER THIS ---------------------------------------------
21
22# Import snapshots
23snapshot_dir = ''
24files = glob.glob(snapshot_dir + "particles-*.vtu")
25# Sort the snapshots
26files = Tcl().call('lsort', '-dict', files)
27
28# Simulation parameters
29E_ini = 1.11362
30time_step_sim = 1e-3
31snapshot_rate = 10
32SUN_X = 0.5
33SUN_Y = 0.55
34
35# Initialize arrays
36E_dif_all = np.zeros( len(files) )
37t_all = np.zeros( len(files) )
38
39# Read snapshots and gather data
40for N_snapshot, file in enumerate(files, start=0):
41
42 print("N_snapshot: ", N_snapshot)
43 time = N_snapshot * snapshot_rate * time_step_sim
44 t_all[N_snapshot] = time
45
46 print('Reading snapshot: ' + file )
47 print("time of simulation = ", time)
48 print("********************************************")
49
50 reader = ParticleVTUReader(vtufile=file)
51 partData = reader.load()
52 partData.show_attribute_list()
53
54 print("********************************************")
55
56 # Read coordinates
57 pos = partData.x
58 x = pos[:,0]
59 y = pos[:,1]
60
61 # Read energy
62 E_tot = partData.energyTot
63 E_tot = np.sum(E_tot)
64 E_dif = (E_tot - E_ini) / E_ini
65 E_dif_all[N_snapshot] = E_dif
66
67 print("Measure energy conservation:")
68 print("Initial E = ", E_ini)
69 print("Total E = ", E_tot)
70 print("Relative Energy difference = ", E_dif )
71 print("********************************************")
72
73print("Done reading all snapshots. Total= ", len(files) )
74
75# Plot radial profile
76fig = plt.figure( figsize=(6, 4) )
77
78plt.plot(t_all, E_dif_all, c='r', lw=0.8, zorder=1)
79plt.xlabel(r'$t$')
80plt.ylabel(r'$\Delta E/E_{\rm ini}$')
81#plt.margins(0.05, 0.1)
82plt.tick_params(direction='in', top=True, right=True, left=True, bottom=True)
83plt.axhline(0, c='gray', ls=':', lw=1.2, zorder=-49 )
84plt.xlim(0, t_all.max())
85plt.ylim(-0.005, 0.005 )
86
87plt.tight_layout()
88
89#plt.show()
90fig.savefig('energy_conservation_profile.pdf')
91fig.savefig('energy_conservation_profile.png')
92print("saved energy_conservation_profile")
93#plt.close()
94
95# Check energy conservation
96E_tolerance = 0.05
97
98E_min = E_dif_all.min()
99E_max = E_dif_all.max()
100E_mean= np.mean(E_dif_all)
101
102assert abs(E_min) < E_tolerance, f"Energy not conserved. E_min: {E_min}"
103assert abs(E_mean)< E_tolerance, f"Energy not conserved. E_mean: {E_min}"
104assert abs(E_max) < E_tolerance, f"Energy not conserved. E_max: {E_max}"
105
106# @TODO test if there is any long-term drift, e.g. of the mean.
107
108print('Test passed!')
109
Reads in particle fields from a single .vtu file, or all .vtu at a given snapshot time specified in t...