From 3614020ede36e96f02cebcae32db942e085b81ab Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 17 Dec 2021 15:17:34 +0100 Subject: [PATCH] switch back to unit time in seconds, prepare for parallel computing --- lib/LeapFrog.py | 4 +-- lib/plots.py | 11 ++++---- lib/units.py | 2 +- main.py | 4 +-- main_concurrent.py | 68 ++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 79 insertions(+), 10 deletions(-) create mode 100755 main_concurrent.py diff --git a/lib/LeapFrog.py b/lib/LeapFrog.py index 5cfd77e..ef3a388 100644 --- a/lib/LeapFrog.py +++ b/lib/LeapFrog.py @@ -58,9 +58,9 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals if display and j % 10 == 0: # display progression if len(dyn_syst.bodylist) == 1: - d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt)) + d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt/yr)) else: - d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt)) + d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt/yr)) if display: d.close() if not savename is None: diff --git a/lib/plots.py b/lib/plots.py index c645018..d3029d5 100755 --- a/lib/plots.py +++ b/lib/plots.py @@ -110,12 +110,12 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""): bodies = "" for body in dyn_syst.bodylist: bodies += str(body)+" ; " - title1, title2 = "Relative difference of the {0:s} ","for a system composed of {0:s}\n integrated with {1:s} for a duration of {2:.2f} years ".format(bodies, integrator, duration) + title1, title2 = "Relative difference of the {0:s} ","for a system composed of {0:s}\n integrated with {1:s} for a duration of {2:.2f} years ".format(bodies, integrator, duration/yr) fig1 = plt.figure(figsize=(15,7)) ax1 = fig1.add_subplot(111) for i in range(len(E)): - ax1.plot(np.arange(E[i].shape[0]-1)*step[i], np.abs((E[i][1:]-E[i][0])/E[i][0]), label="step of {0:.2e}yr".format(step[i])) + ax1.plot(np.arange(E[i].shape[0]-1)*step[i]/yr, np.abs((E[i][1:]-E[i][0])/E[i][0]), label="step of {0:.2e}s".format(step[i])) ax1.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta E_m}{E_m(t=0)}\right|$", yscale='log') ax1.legend() fig1.suptitle(title1.format("mechanical energy")+title2) @@ -126,7 +126,7 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""): for i in range(len(L)): dL = ((L[i]-L[i][0])/L[i][0]) dL[np.isnan(dL)] = 0. - ax2.plot(np.arange(L[i].shape[0]-1)*step[i], np.abs(np.sum(dL[1:],axis=1)), label="step of {0:.2e}yr".format(step[i])) + ax2.plot(np.arange(L[i].shape[0]-1)*step[i]/yr, np.abs(np.sum(dL[1:],axis=1)), label="step of {0:.2e}s".format(step[i])) ax2.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$",yscale='log') ax2.legend() fig2.suptitle(title1.format("kinetic moment")+title2) @@ -134,8 +134,9 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""): fig3 = plt.figure(figsize=(15,7)) ax3 = fig3.add_subplot(111) - ax3.plot(np.arange(sma.shape[0])*step[i], sma, label="a (semi major axis)") - ax3.plot(np.arange(ecc.shape[0])*step[i], ecc, label="e (eccentricity)") + for i in range(len(L)): + ax3.plot(np.arange(sma[i].shape[0])*step[i]/yr, sma[i], label="a (step of {0:.2e}s)".format(step[i])) + ax3.plot(np.arange(ecc[i].shape[0])*step[i]/yr, ecc[i], label="e (step of {0:.2e}s)".format(step[i])) ax3.set(xlabel=r"$t \, [yr]$", ylabel=r"$a \, [au] \, or \, e$") ax3.legend() fig3.suptitle("Semi major axis and eccentricity "+title2) diff --git a/lib/units.py b/lib/units.py index db33ba9..e9516e5 100644 --- a/lib/units.py +++ b/lib/units.py @@ -8,4 +8,4 @@ globals()['G'] = 6.67e-11 #Gravitational constant in SI units globals()['Ms'] = 2e30 #Solar mass in kg globals()['au'] = 1.5e11 #Astronomical unit in m globals()['yr'] = 3.15576e7 #year in seconds -globals()['Ga'] = G*Ms*yr**2/au**3 #Gravitational constant dimensionless \ No newline at end of file +globals()['Ga'] = G*Ms/au**3 #Gravitational constant dimensionless \ No newline at end of file diff --git a/main.py b/main.py index 4fb0ef7..ac213e3 100755 --- a/main.py +++ b/main.py @@ -27,9 +27,9 @@ def main(): v = np.array([v1, v2, v3],dtype=np.longdouble) #integration parameters - duration, step = 100*yr/yr, np.array([1./(365.25*2.), 1./(365.25*1.), 5./(365.25*1.)],dtype=np.longdouble)*yr/yr #integration time and step in years + duration, step = 100*yr, np.array([60.],dtype=np.longdouble) #integration time and step in seconds step = np.sort(step)[::-1] - integrator = "hermite" + integrator = "leapfrog" n_bodies = 3 display = False savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator) diff --git a/main_concurrent.py b/main_concurrent.py new file mode 100755 index 0000000..067673d --- /dev/null +++ b/main_concurrent.py @@ -0,0 +1,68 @@ +#!/usr/bin/python +# -*- coding:utf-8 -*- +from sys import exit as sysexit +from copy import deepcopy +import numpy as np +import matplotlib.pyplot as plt +from concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor +from lib.objects import Body, System +from lib.LeapFrog import leapfrog +from lib.hermite import hermite +from lib.plots import display_parameters +from lib.units import * + +def main(): + #initialisation + m = np.array([1., 1., 1e-1],dtype=np.longdouble)*Ms/Ms # Masses in Solar mass + a = np.array([1., 1., 5.],dtype=np.longdouble)*au/au # Semi-major axis in astronomical units + e = np.array([0., 0., 0.],dtype=np.longdouble) # Eccentricity + psi = np.array([0., 0., 0.],dtype=np.longdouble)*np.pi/180. # Inclination of the orbital plane in degrees + + x1 = np.array([0., -1., 0.],dtype=np.longdouble)*a[0]*(1.+e[0]) + x2 = np.array([0., 1., 0.],dtype=np.longdouble)*a[1]*(1.+e[1]) + x3 = np.array([np.cos(psi[2]), 0., np.sin(psi[2])],dtype=np.longdouble)*a[2]*(1.+e[2]) + q = np.array([x1, x2, x3],dtype=np.longdouble) + + v1 = np.array([np.sqrt(Ga*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble) + v2 = np.array([-np.sqrt(Ga*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble) + v3 = np.array([0., np.sqrt(Ga*(m[0]+m[1])*(2./np.sqrt(np.sum(q[2]**2))-1./a[2])), 0.],dtype=np.longdouble) + v = np.array([v1, v2, v3],dtype=np.longdouble) + + #integration parameters + duration, step = 100*yr, np.array([1./(365.25*24.), 12./(365.25*24.), 24./(365.25*24.)],dtype=np.longdouble)*yr #integration time and step in seconds + step = np.sort(step)[::-1] + integrator = "leapfrog" + n_bodies = 3 + display = False + savename = "{0:d}bodies_conc_{1:s}".format(n_bodies, integrator) + + bodies, bodysyst = [],[] + for j in range(n_bodies): + bodies.append(Body(m[j], q[j], v[j])) + bin_syst = System(bodies[0:2]) + dyn_syst = System(bodies, main=True) + bodysyst = [[deepcopy(bin_syst), deepcopy(dyn_syst)] for _ in range(n_bodies)] + #simulation start + exe = ProcessPoolExecutor() + future_ELae = [] + for i,step0 in enumerate(step): + if i != 0: + display = False + if integrator.lower() in ['leapfrog', 'frogleap', 'frog']: + future_ELae.append(exe.submit(leapfrog, bodysyst[i][1], bodysyst[i][0], duration, step0, recover_param=True, display=display, savename=savename)) + elif integrator.lower() in ['hermite','herm']: + future_ELae.append(exe.submit(hermite, bodysyst[i][1], bodysyst[i][0], duration, step0, recover_param=True, display=display, savename=savename)) + + E, L, sma, ecc = [], [], [], [] + for future in future_ELae: + E0, L0, sma0, ecc0 = future.result() + E.append(E0) + L.append(L0) + sma.append(sma0) + ecc.append(ecc0) + parameters = [duration, step, dyn_syst, integrator] + display_parameters(E, L, sma, ecc, parameters=parameters, savename=savename) + return 0 + +if __name__ == '__main__': + sysexit(main())