From 969fdb0b44a833ba043b001acc8a5eb33f80ae87 Mon Sep 17 00:00:00 2001 From: Yael-II <154762804+Yael-II@users.noreply.github.com> Date: Tue, 21 Jan 2025 14:44:47 +0100 Subject: [PATCH] new test --- README.md | 10 ++- Source/test_integrators.py | 169 +++++++++++++++++++++++++++++++++++++ test_integrators.sh | 4 + 3 files changed, 180 insertions(+), 3 deletions(-) create mode 100644 Source/test_integrators.py create mode 100755 test_integrators.sh diff --git a/README.md b/README.md index f541910..bbf487c 100644 --- a/README.md +++ b/README.md @@ -44,7 +44,7 @@ The result of the simulations are saved in the Output directory, with the prefix ./area.sh ``` 3. To run the tests, use: - - To test the potentials + - To test the potentials: ```bash ./test_potentials.sh ``` @@ -53,11 +53,15 @@ The result of the simulations are saved in the Output directory, with the prefix ./test_evolution.sh ./test_evolution_chaotic.sh ``` - - To test the generation of particles in this potential with a given energy + - To test the generation of particles in this potential with a given energy: ```bash ./test_initial_E.sh ``` - - To get the running time of both Poincaré sections computations (parallel vs. linear algorithms) + - To test the different integrators we tried: + ```bash + ./test_integrators + ``` + - To get the running time of both Poincaré sections computations (parallel vs. linear algorithms): ```bash ./time_poincare_sections.sh ``` diff --git a/Source/test_integrators.py b/Source/test_integrators.py new file mode 100644 index 0000000..5d68a2f --- /dev/null +++ b/Source/test_integrators.py @@ -0,0 +1,169 @@ +""" +Test: Integrators + +Demonstrating Keplerian 2-body orbits using various integrators, +and comparing accuracy and runtime over a range of step sizes. + +@ Author: Moussouni, Yaël (MSc student) & Bhat, Junaid Ramzan (MSc student) +@ Institution: Université de Strasbourg, CNRS, Observatoire astronomique + de Strasbourg, UMR 7550, F-67000 Strasbourg, France +@ Date: 2025-01-01 + +Licence: +Order and Chaos in a 2D potential +Copyright (C) 2025 Yaël Moussouni (yael.moussouni@etu.unistra.fr) + Bhat, Junaid Ramzan (junaid-ramzan.bhat@etu.unistra.fr) + +test_integrators.py +Copyright (C) 2025 Bhat, Junaid Ramzan (junaid-ramzan.bhat@etu.unistra.fr) + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see https://www.gnu.org/licenses/. +""" + +import numpy as np +import matplotlib.pyplot as plt +import time + +import integrator as intg # integrator module +import kepler_orbit as kep # Analytical & numeric Kepler functions +import initial_conditions as ic # For initial condition + +# ---------------------------- +# 1. Setup & global parameters +# ---------------------------- +t0 = 0.0 +T_final = 8.0 +y0 = ic.one_part(1, 0, 0, 1) # [x, y, vx, vy] + + +h_range = np.logspace(-3.5, -0.1, 25) + +# For plotting lines/colors +methods = [ + ("Euler", intg.euler, 'o--', 'red'), + ("RK2", intg.rk2, 's--', 'royalblue'), + ("RK4", intg.rk4, '^--', 'limegreen') +] +colors = {'Analytical': 'navy', 'Euler': 'red', 'RK2': 'royalblue', 'RK4': 'limegreen'} + +# Compute machine epsilon +eps = 1.0 +while 1.0 + eps/2 > 1.0: + eps /= 2.0 +print(f"Machine epsilon: {eps}") + +# Arrays to store final energy errors & times +err_euler, err_rk2, err_rk4 = [], [], [] +time_euler, time_rk2, time_rk4 = [], [], [] + +# ------------------------------------------ +# 2. Main loop over step sizes h in h_range +# ------------------------------------------ +for h in h_range: + N = int(T_final / h) + + # Analytical solution + t_ana, pos_ana, vel_ana, en_ana = kep.kepler_analytical_orb(y0, h, N) + E_analytical_final = en_ana[-1] + + # Numerical integrators + timing + all_solutions = {} + for (label, method, *_), store_err, store_t in zip( + methods, + [err_euler, err_rk2, err_rk4], + [time_euler, time_rk2, time_rk4] + ): + start_time = time.time() + t_num, y_num = intg.integrator_type(t0, y0, h, N, kep.kepler_orbnum, method) + elapsed = time.time() - start_time + + store_t.append(elapsed) + + # Final energy error + en_num = kep.kepler_enrgnum(y_num) + E_numerical_final = en_num[-1] + store_err.append(abs(E_analytical_final - E_numerical_final)) + + all_solutions[label] = y_num + + # Orbit plot (optional, can comment out if too many figures) + eu_vals = all_solutions["Euler"] + rk2_vals = all_solutions["RK2"] + rk4_vals = all_solutions["RK4"] + + fig, ax = plt.subplots(figsize=(6, 6)) + ax.plot(pos_ana[:, 0], pos_ana[:, 1], '-.', color=colors['Analytical'], + label="Analytical", zorder=4) + ax.plot(eu_vals[:, 0, 0], eu_vals[:, 0, 1], color=colors['Euler'], label="Euler") + ax.plot(rk2_vals[:, 0, 0], rk2_vals[:, 0, 1], color=colors['RK2'], label="RK2") + ax.plot(rk4_vals[:, 0, 0], rk4_vals[:, 0, 1], color=colors['RK4'], label="RK4") + ax.set_title(f"Orbit for dt = {h:.4g}") + ax.set_xlabel("x position(reduced units)") + ax.set_ylabel("y position(reduced units)") + ax.legend(loc="best") + ax.grid(True, which='both', ls='--', alpha=0.5) + plt.tight_layout() + plt.savefig(f"orbit_dt_{h:.4g}.png", dpi=300, bbox_inches='tight') + plt.show() + +# --------------------------------------------------------------- +# 3. Summary Plots: CPU time and final energy error (Log-Log) +# --------------------------------------------------------------- + +# --- Step size vs. CPU Time (Log-Log) --- +fig, ax = plt.subplots(figsize=(8, 6)) +ax.loglog(h_range, time_euler, 'o--', color='red', label="Euler") +ax.loglog(h_range, time_rk2, 's--', color='royalblue', label="RK2") +ax.loglog(h_range, time_rk4, '^--', color='limegreen', label="RK4") + +ax.set_xlabel(r"Step size $(\Delta t)$", fontsize=12) +ax.set_ylabel("CPU Time (s)", fontsize=12) +ax.set_title(r"$\Delta t$ vs. CPU Time ", fontsize=14) +ax.minorticks_on() +ax.grid(True, which="major", linestyle="--", linewidth=0.5, alpha=0.7) +ax.grid(True, which="minor", linestyle=":", linewidth=0.5, alpha=0.5) +ax.legend(loc="best", fontsize=12) +plt.tight_layout() +plt.savefig("dt_vs_cpu_time_loglog.png", dpi=300, bbox_inches='tight') +plt.show() + +# --- Step size vs. Final Energy Error (Log-Log) --- +fig, ax = plt.subplots(figsize=(8, 6)) +ax.loglog(h_range, err_euler, 'o--', color='red', label="Euler") +ax.loglog(h_range, err_rk2, 's--', color='royalblue', label="RK2") +ax.loglog(h_range, err_rk4, '^--', color='limegreen', label="RK4") + +# Machine Epsilon line (horizontal) +ax.axhline(eps, color='darkred', ls='--', label='Machine Epsilon') + +ax.set_xlabel(r"Step size $(\Delta t)$", fontsize=12) +ax.set_ylabel(r"$|E_{\mathrm{analytical}} - E_{\mathrm{numerical}}|$", fontsize=12) +ax.set_title(r"$\Delta t$ vs. Final Energy Error ", fontsize=14) +ax.minorticks_on() +ax.grid(True, which="major", linestyle="--", linewidth=0.5, alpha=0.7) +ax.grid(True, which="minor", linestyle=":", linewidth=0.5, alpha=0.5) + +# Ensure 'Machine Epsilon' is in legend +handles, labels = ax.get_legend_handles_labels() +if 'Machine Epsilon' not in labels: + import matplotlib.lines as mlines + h_me = mlines.Line2D([], [], color='darkred', ls='--', label='Machine Epsilon') + handles.append(h_me) + labels.append('Machine Epsilon') +ax.legend(handles, labels, loc="best", fontsize=12) + +plt.tight_layout() +plt.savefig("timestep_vs_final_energy_error_loglog1.png", dpi=300, bbox_inches='tight') +plt.show() + diff --git a/test_integrators.sh b/test_integrators.sh new file mode 100755 index 0000000..b871563 --- /dev/null +++ b/test_integrators.sh @@ -0,0 +1,4 @@ +#!/usr/bin/env bash + +source activate.sh +venv/bin/python Source/test_integrators.py