GravSynch.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat May 11 08:56:41 2019
@author: nolte
D. D. Nolte, Introduction to Modern Dynamics: Chaos, Networks, Space and Time, 2nd ed. (Oxford,2019)
"""

# https://www.python-course.eu/networkx.php
# https://networkx.github.io/documentation/stable/tutorial.html
# https://networkx.github.io/documentation/stable/reference/functions.html

import numpy as np
from scipy import integrate
from matplotlib import pyplot as plt
import networkx as nx
from UserFunction import linfit
from pathlib import Path
import time

tstart = time.time()

plt.close('all')

Grav = 6.674*10**(-11)
M0 = 1.989*10**30
M = .9*M0
prM = 'Mo' + str(int(M/M0))
c = 3*10**8
Rs = 2*Grav*M/c**2
#R = 19000 + 14959787
R = 4*Rs
scale = 10**-7

Nfac = 25  # 25
N = 20  # 50

# model_case 1 = complete graph (Kuramoto transition)

facoef = .2
nodecouple = nx.complete_graph(N)

# function: omegout, yout = coupleN(G)
def coupleN(G):
    # function: yd = flow_deriv(x_y)
    def flow_deriv(y, t0):

        yp = np.zeros(shape=(N,))
        for omloop in range(N):
            temp = omega[omloop]
            linksz = G.nodes[omloop]['numlink']
            for cloop in range(linksz):
                cindex = G.nodes[omloop]['link'][cloop]
                k = G.nodes[omloop]['coupling'][cloop]

                temp = temp + k * np.sin(y[cindex] - y[omloop])

            yp[omloop] = temp

        yd = np.zeros(shape=(N,))
        for omloop in range(N):
            yd[omloop] = yp[omloop]

        return yd

    # end of function flow_deriv(x_y)

    mnomega = 1.0

    for nodeloop in range(N):
        omega[nodeloop] = G.nodes[nodeloop]['element']

    x_y_z = omega

    # Settle-down Solve for the trajectories
    tsettle = 100
    t = np.linspace(0, tsettle, tsettle)
    x_t = integrate.odeint(flow_deriv, x_y_z, t)
    x0 = x_t[tsettle - 1, 0:N]

    t = np.linspace(1, 1000, 1000)
    y = integrate.odeint(flow_deriv, x0, t)
    siztmp = np.shape(y)
    sy = siztmp[0]

    # Fit the frequency
    m = np.zeros(shape=(N,))
    w = np.zeros(shape=(N,))
    mtmp = np.zeros(shape=(4,))
    btmp = np.zeros(shape=(4,))
    for omloop in range(N):

        if np.remainder(sy, 4) == 0:
            mtmp[0], btmp[0] = linfit(t[0:sy // 2], y[0:sy // 2, omloop])
            mtmp[1], btmp[1] = linfit(t[sy // 2 + 1:sy], y[sy // 2 + 1:sy, omloop])
            mtmp[2], btmp[2] = linfit(t[sy // 4 + 1:3 * sy // 4], y[sy // 4 + 1:3 * sy // 4, omloop])
            mtmp[3], btmp[3] = linfit(t, y[:, omloop])
        else:
            sytmp = 4 * np.floor(sy / 4)
            mtmp[0], btmp[0] = linfit(t[0:sytmp // 2], y[0:sytmp // 2, omloop])
            mtmp[1], btmp[1] = linfit(t[sytmp // 2 + 1:sytmp], y[sytmp // 2 + 1:sytmp, omloop])
            mtmp[2], btmp[2] = linfit(t[sytmp // 4 + 1:3 * sytmp / 4], y[sytmp // 4 + 1:3 * sytmp // 4, omloop])
            mtmp[3], btmp[3] = linfit(t[0:sytmp], y[0:sytmp, omloop])

        # m[omloop] = np.median(mtmp)
        m[omloop] = np.mean(mtmp)

        w[omloop] = mnomega + m[omloop]

    omegout = m
    yout = y

    return omegout, yout
    # end of function: omegout, yout = coupleN(G)

l = []
grav = []
omegatemp = np.array(())
h = 10

print(Rs/R/scale)
h_dil = h * np.sqrt(1-Rs/R)
lab = "_T" + "&S"
for i in range(N):
    h_dil = (i*h) * np.sqrt(1 - Rs / (R+h_dil))
    #l.append(i*h_dil)
    l.append(h_dil)

    grav.append(2*Grav*M/(c**2*R*scale*(1 + l[i]/R)))
    omegatemp = np.append(omegatemp,(grav[i]))

Nlink = N * (N - 1) // 2
omega = np.zeros(shape=(N,))


meanomega = np.mean(omegatemp)
omega = (omegatemp - meanomega)
print(omegatemp)
print(meanomega)
print(omega)
sto = np.std(omega)

lnk = np.zeros(shape=(N,), dtype=int)
for loop in range(N):
    nodecouple.nodes[loop]['element'] = omega[loop]
    nodecouple.nodes[loop]['link'] = list(nx.neighbors(nodecouple, loop))
    nodecouple.nodes[loop]['numlink'] = np.size(list(nx.neighbors(nodecouple, loop)))
    lnk[loop] = np.size(list(nx.neighbors(nodecouple, loop)))

avgdegree = np.mean(lnk)
mnomega = 1

facval = np.zeros(shape=(Nfac,))
yy = np.zeros(shape=(Nfac, N))
xx = np.zeros(shape=(Nfac,))
for facloop in range(Nfac):
    print(facloop)

    fac = facoef * (16 * facloop / (Nfac)) * (1 / (N - 1)) * sto / mnomega
    for nodeloop in range(N):
        nodecouple.nodes[nodeloop]['coupling'] = np.zeros(shape=(lnk[nodeloop],))
        for linkloop in range(lnk[nodeloop]):
            nodecouple.nodes[nodeloop]['coupling'][linkloop] = fac

    facval[facloop] = fac * avgdegree*scale

    omegout, yout = coupleN(nodecouple)  # Here is the subfunction call for the flow

    for omloop in range(N):
        yy[facloop, omloop] = omegout[omloop]

    xx[facloop] = facval[facloop]/scale

plt.figure(1)
lines = plt.plot(xx, yy)
plt.setp(lines, linewidth=0.5)
plt.xticks(fontsize = 8)
plt.xlabel('Couple Constant g')
plt.yticks(fontsize = 8)
plt.ylabel('Frequency')

plt.title('Kuramoto Synchronization Transition')

filename_out = Path("Plot_N" + str(N) + '_' + prM + lab + ".png")
plt.savefig('GravSynch.png',dpi = 300)
plt.show()


elapsed_time = time.time() - tstart
print('elapsed time = ', format(elapsed_time, '.2f'), 'secs')