added new example for coldplasmavlasov, removed examples from gitignore

analysis/check_distribution_moments.py

@@ -0,0 +1,165 @@
+import os
+import numpy as np
+
+BASE = (
+    "thesis_postprocess_validation/post_processing/"
+    "kinetic_data/hot_elec/distribution_function/v1_v3_density"
+)
+
+TIME_PATH = (
+    "thesis_postprocess_validation/post_processing/t_grid.npy"
+)
+
+t = np.load(TIME_PATH)
+
+v1 = np.load(os.path.join(BASE, "grid_v1.npy"))
+v3 = np.load(os.path.join(BASE, "grid_v3.npy"))
+F = np.load(os.path.join(BASE, "f_binned.npy"))
+
+dv1 = np.mean(np.diff(v1))
+dv3 = np.mean(np.diff(v3))
+
+
+def moments(Fi):
+
+    norm = np.sum(Fi) * dv1 * dv3
+
+    mean1 = (
+        np.sum(v1[:, None] * Fi)
+        * dv1 * dv3 / norm
+    )
+
+    mean3 = (
+        np.sum(v3[None, :] * Fi)
+        * dv1 * dv3 / norm
+    )
+
+    var1 = (
+        np.sum(
+            (v1[:, None] - mean1) ** 2 * Fi
+        )
+        * dv1 * dv3 / norm
+    )
+
+    var3 = (
+        np.sum(
+            (v3[None, :] - mean3) ** 2 * Fi
+        )
+        * dv1 * dv3 / norm
+    )
+
+    covariance = (
+        np.sum(
+            (v1[:, None] - mean1)
+            * (v3[None, :] - mean3)
+            * Fi
+        )
+        * dv1 * dv3 / norm
+    )
+
+    return (
+        norm,
+        mean1,
+        mean3,
+        np.sqrt(var1),
+        np.sqrt(var3),
+        covariance,
+    )
+
+
+results = np.array([moments(Fi) for Fi in F])
+
+norm = results[:, 0]
+mean1 = results[:, 1]
+mean3 = results[:, 2]
+sigma1 = results[:, 3]
+sigma3 = results[:, 4]
+covariance = results[:, 5]
+
+anisotropy = sigma1**2 / sigma3**2
+
+relative_change = np.array(
+    [
+        np.linalg.norm(Fi - F[0])
+        / np.linalg.norm(F[0])
+        for Fi in F
+    ]
+)
+
+marginal_v1_initial = np.sum(F[0], axis=1) * dv3
+marginal_v1_final = np.sum(F[-1], axis=1) * dv3
+
+marginal_v3_initial = np.sum(F[0], axis=0) * dv1
+marginal_v3_final = np.sum(F[-1], axis=0) * dv1
+
+v1_change = (
+    np.linalg.norm(
+        marginal_v1_final - marginal_v1_initial
+    )
+    / np.linalg.norm(marginal_v1_initial)
+)
+
+v3_change = (
+    np.linalg.norm(
+        marginal_v3_final - marginal_v3_initial
+    )
+    / np.linalg.norm(marginal_v3_initial)
+)
+
+print("============================================================")
+print("DISTRIBUTION MOMENT ANALYSIS")
+print("============================================================")
+
+print("\nInitial values:")
+print("integral =", norm[0])
+print("sigma_v1 =", sigma1[0])
+print("sigma_v3 =", sigma3[0])
+print("anisotropy =", anisotropy[0])
+print("covariance =", covariance[0])
+
+print("\nFinal values:")
+print("integral =", norm[-1])
+print("sigma_v1 =", sigma1[-1])
+print("sigma_v3 =", sigma3[-1])
+print("anisotropy =", anisotropy[-1])
+print("covariance =", covariance[-1])
+
+print("\nMarginal changes:")
+print("v1 relative L2 change =", v1_change)
+print("v3 relative L2 change =", v3_change)
+
+print("\nFull 2D distribution change:")
+print("initial -> final =", relative_change[-1])
+print("maximum =", np.max(relative_change))
+print(
+    "time of maximum change =",
+    t[np.argmax(relative_change)],
+)
+
+print("\nRanges over complete run:")
+
+print(
+    "sigma_v1:",
+    np.min(sigma1),
+    np.max(sigma1),
+)
+
+print(
+    "sigma_v3:",
+    np.min(sigma3),
+    np.max(sigma3),
+)
+
+print(
+    "anisotropy:",
+    np.min(anisotropy),
+    np.max(anisotropy),
+)
+
+print(
+    "covariance:",
+    np.min(covariance),
+    np.max(covariance),
+)
+
+print("\nAnalysis complete.")

analysis/plot_e3_density.py

@@ -0,0 +1,81 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+BASE = (
+    "thesis_postprocess_validation/post_processing/"
+    "kinetic_data/hot_elec/distribution_function/e3_density"
+)
+
+TIME_PATH = (
+    "thesis_postprocess_validation/post_processing/t_grid.npy"
+)
+
+OUT = "validation_plots"
+os.makedirs(OUT, exist_ok=True)
+
+t = np.load(TIME_PATH)
+
+e3 = np.load(os.path.join(BASE, "grid_e3.npy"))
+F = np.load(os.path.join(BASE, "f_binned.npy"))
+DF = np.load(os.path.join(BASE, "delta_f_binned.npy"))
+
+de3 = np.mean(np.diff(e3))
+
+integrals = np.sum(F, axis=1) * de3
+
+print("============================================================")
+print("E3 SPATIAL DENSITY VALIDATION")
+print("============================================================")
+
+print("grid shape =", e3.shape)
+print("distribution shape =", F.shape)
+
+print("\ne3 range:")
+print(np.min(e3), np.max(e3))
+
+print("\nIntegral:")
+print("initial =", integrals[0])
+print("final =", integrals[-1])
+print("minimum =", np.min(integrals))
+print("maximum =", np.max(integrals))
+
+relative_variation = np.max(
+    np.abs(integrals - integrals[0])
+) / abs(integrals[0])
+
+print("maximum relative variation =", relative_variation)
+
+print(
+    "\nDelta-f identically zero =",
+    np.all(DF == 0),
+)
+
+indices = [0, len(t) // 2, len(t) - 1]
+
+plt.figure(figsize=(10, 6))
+
+for i in indices:
+    plt.plot(
+        e3,
+        F[i],
+        label=f"t = {t[i]:.3f}",
+    )
+
+plt.xlabel(r"$e_3$")
+plt.ylabel("f")
+plt.title("Hot-Electron Spatial Density Distribution")
+plt.legend()
+plt.grid()
+
+plt.tight_layout()
+
+path = os.path.join(
+    OUT,
+    "e3_density_selected_times.png",
+)
+
+plt.savefig(path, dpi=200)
+plt.close()
+
+print("\nSaved:", path)

analysis/plot_energy_diagnostics.py

@@ -0,0 +1,105 @@
+import os
+import h5py
+import numpy as np
+import matplotlib.pyplot as plt
+
+DATA = (
+    "thesis_postprocess_validation/data/data_proc0.hdf5"
+)
+
+OUT = "validation_plots"
+os.makedirs(OUT, exist_ok=True)
+
+with h5py.File(DATA, "r") as f:
+
+    t = f["time/value"][:]
+
+    en_B = f["scalar/en_B"][:]
+    en_E = f["scalar/en_E"][:]
+    en_J = f["scalar/en_J"][:]
+    en_f = f["scalar/en_f"][:]
+    en_tot = f["scalar/en_tot"][:]
+
+relative_energy_error = (
+    (en_tot - en_tot[0])
+    / en_tot[0]
+)
+
+print("============================================================")
+print("ENERGY DIAGNOSTICS")
+print("============================================================")
+
+print("number of states =", len(t))
+print("time range =", t[0], t[-1])
+
+print("\nInitial energies:")
+print("E_B =", en_B[0])
+print("E_E =", en_E[0])
+print("E_J =", en_J[0])
+print("E_f =", en_f[0])
+print("E_total =", en_tot[0])
+
+print("\nFinal energies:")
+print("E_B =", en_B[-1])
+print("E_E =", en_E[-1])
+print("E_J =", en_J[-1])
+print("E_f =", en_f[-1])
+print("E_total =", en_tot[-1])
+
+print("\nEnergy conservation:")
+print(
+    "final relative change =",
+    relative_energy_error[-1],
+)
+
+print(
+    "maximum absolute relative change =",
+    np.max(np.abs(relative_energy_error)),
+)
+
+plt.figure(figsize=(10, 6))
+
+plt.plot(t, en_B, label=r"$E_B$")
+plt.plot(t, en_E, label=r"$E_E$")
+plt.plot(t, en_J, label=r"$E_J$")
+plt.plot(t, en_f, label=r"$E_f$")
+
+plt.xlabel("Time")
+plt.ylabel("Energy")
+plt.title("Energy Evolution — Validation Run")
+plt.legend()
+plt.grid()
+
+plt.tight_layout()
+
+path = os.path.join(
+    OUT,
+    "energy_components.png",
+)
+
+plt.savefig(path, dpi=200)
+plt.close()
+
+print("\nSaved:", path)
+
+
+plt.figure(figsize=(10, 6))
+
+plt.plot(t, relative_energy_error)
+
+plt.xlabel("Time")
+plt.ylabel("Relative total-energy change")
+plt.title("Total Energy Conservation")
+plt.grid()
+
+plt.tight_layout()
+
+path = os.path.join(
+    OUT,
+    "total_energy_error.png",
+)
+
+plt.savefig(path, dpi=200)
+plt.close()
+
+print("Saved:", path)