Commit ecfeaf43 authored by Sayah-Sel's avatar Sayah-Sel
Browse files

animation trajectories for several epsilon

parent 870610bc
No preview for this file type
import numpy as np
from numpy import fft as f
import multiprocessing as mp
import time
import finufft
import os
import warnings
os.environ['KMP_DUPLICATE_LIB_OK']='True'
d=2
L=np.pi
N=32
dx=L/N
mu=0
corr_L=["norm1","norm2"]
corr="norm2"
l_cV=0.1*np.pi
l_cA=0.1*np.pi
e_L=[0.01,0.1,0.2,0.3,0.5,0.7,0.8,0.9,0.99]
XI=[]
for k in range(d):
XI.append(f.fftfreq(2*N)*2*N*2*np.pi/(2*L))
xi=np.array(np.meshgrid(*XI,indexing='ij'))
X= np.mgrid[tuple(slice(- L, L, dx) for _ in range(d))].T
Y=X.reshape((-1,d))
I=np.ones_like(Y.sum(axis=-1))
if corr=="norm2":
p=2
elif corr=="norm1":
p=1
def make_V(size=1,d=d,l_c=l_cV,xi=xi,N=N,corr=corr,nb=len(e_L)):
if corr=="norm2":
k2=(xi**2).sum(axis=0)
covVH=(2*np.pi*l_c**2)**(d/2)*np.exp(-l_c**2*k2/4) #(2*N)**d*
elif corr=="norm1":
truc=2*l_c/(1+(l_c*xi)**2)
covVH=2*truc.prod(axis=0) #(2*N)**d*
if size>1:
list_Vh=[]
for k in range(size):
Vh=covVH*(np.random.standard_normal(covVH.shape)+1j*np.random.standard_normal(covVH.shape))
list_Vh.append(Vh)
list_Vh=np.array(list_Vh)
return list_Vh
else:
Vh=covVH*(np.random.standard_normal(covVH.shape)+1j*np.random.standard_normal(covVH.shape))
return Vh
def make_A(size=1,d=d,l_c=l_cA,xi=xi,N=N,corr=corr,nb=len(e_L)):
Ah=np.zeros((d,d)+xi[0].shape,dtype=np.complex128)
if corr=="norm2":
k2=(xi**2).sum(axis=0)
covAH=(2*np.pi*l_c**2)**(d/2)*np.exp(-l_c**2*k2/4) #(2*N)**d*
elif corr=="norm1":
truc=2*l_c/(1+(l_c*xi)**2)
covAH=2*truc.prod(axis=0) #(2*N)**d
if size>1:
list_Ah=[]
for _ in range(size):
for i in range(1,d):
for j in range(i):
Ah[i,j,:]=covAH*(np.random.standard_normal(covAH.shape)+1j*np.random.standard_normal(covAH.shape))
Ah[j,i,:]=-Ah[i,j,:]
list_Ah.append(Ah)
list_Ah=np.array(list_Ah)
return list_Ah
else:
for i in range(1,d):
for j in range(i):
Ah[i,j,:]=covAH*(np.random.standard_normal(covAH.shape)+1j*np.random.standard_normal(covAH.shape))
Ah[j,i,:]=-Ah[i,j,:]
return Ah
#For comparison, we used the same potential A and V, only the force f=e*f_V+(1-e)*f_A changes, so we can really what individual fixed points becomes
Vh=make_V()
Ah=make_A()
eL=np.array([np.ones(samples)*e for e in e_L]).flatten()
plan=finufft.Plan(nufft_type=2,n_modes_or_dim=(2*N,)*d,modeord=1,eps=1e-6,isign=1,n_trans=d,nthreads=1,fftw=64)
Nt=300
Time=3*10**3
dt=3*10**-3
nabVh=np.array([1j*xi[k]*Vh for k in range(d)])
rotAh=np.array([np.array([1j*xi[j]*Ah[i,j,:] for j in range(d)]).sum(axis=0) for i in range(d)])/np.sqrt(d)
Ns=10
start=np.random.uniform(-L,L,size=(Ns,d,1))
#The following code is far more horrible than nearly anything else in the repo,
#because due to periodic b.c the trajectories have to be cut when they reach boundaries
#It is done in three step:
# 1) the trajectories are first cut in sub-array during simulations (because boundary crossing are already detected at this stage so it is faster)
# 2) collections and associated colors are instanciated before plotting: we want all the subparts of a trajectory in the same color
# 3) during plot/animation trajectories are plot up to the current time as a collection of segments
meta_traject=[]
meta_K=[]
max_K=[]
meta_fixe=[]
for e in e_L:
warnings.filterwarnings("ignore")
v=1-e
a=e
fh=-v*nabVh+a*rotAh
#np.random.seed()
T=0.0
threshold=dt*(T+0.1)
#@jit(nopython=True)
def run(x,Time=Time,d=d,dt=dt,T=T,threshold=threshold,plan=plan,save=False): #save = True when we want the full trajectories (as always)
pointss=[]
points=x
fixe=0
force=np.empty((d,1))
x_new=x
K=[]
for k in range(Time):
plan.setpts(*x)
plan.execute(fh,out=force)
#np.random.standard_normal(size=(d,1))*T
x_new=x+dt*force
x=(L+x_new)%(2*L)-L
if len(points[0])>50 and np.std(points[:,-20:-1],axis=1).sum()<threshold:
fixe=1
if save==True:
pointss.append(points.T)
K.append(k)
return pointss,K,fixe
return x.T[0],fixe,k
if np.sqrt((x_new-x)**2).sum()>10**-3:
pointss.append(points.T)
K.append(k)
points=x
else:
points=np.append(points,x,axis=1) #points is treated as a dynamic array
if save==True:
pointss.append(points.T)
K.append(k)
return pointss,K,fixe #Here pointss is a list of arrays...
return x.T[0],fixe,k
def dist(start,d=d,Time=Time,dt=dt,T=T,mu=mu):
traject=[]
K=[]
fixed=[]
for x in start:
traj,k,fixe=run(x,save=True)
traject.append(traj)
K.append(k)
fixed.append(fixe)
return traject,K,fixed
traject,K,fixed=dist(start) #... and trajec is a list of list of arrays...
meta_K.append(K)
meta_traject.append(traject) #... and meta_trajec is a list of list of list of arrays...
meta_fixe.append(fixed)
max_K.append(np.max(np.max(K)))
##
import os
import matplotlib
import matplotlib.pyplot as pl
import matplotlib.animation as anim
import matplotlib.axes as ax
import matplotlib.colors as clr
import matplotlib.collections as mcol
# WARNING: the 2 following lines are workarounds on my personal laptop, they may not work or not be necessary for over configurations
pl.rcParams['animation.ffmpeg_path'] = '/usr/local/bin/ffmpeg'
os.environ['PATH']=os.environ['PATH']+':/usr/local/bin'
matplotlib.use('Agg')
Writer = anim.writers['ffmpeg']
writer = Writer(fps=7, metadata=dict(artist='Me'), bitrate=1800)
fig, axes = pl.subplots()
Zt=[Y.T[k] for k in range(d)]
Z=tuple(Zt)
plan=finufft.Plan(nufft_type=2,n_modes_or_dim=(2*N,)*d,modeord=1,eps=1e-8,isign=1,n_trans=d,fftw=64)
plan.setpts(*Z)
t1=time.time()
nablaV=np.real(plan.execute(nabVh).reshape(X.shape))
rotA=np.real(plan.execute(rotAh).reshape(X.shape))
#the use of a quiver plot is not so usual but very useful for plotting vector fields
e=e_L[0]
force=(-(1-e)*nablaV+e*rotA)
Quiv=axes.quiver(X[:,:,0].flatten(),X[:,:,1].flatten(),force[:,:,0].flatten(),force[:,:,1].flatten(), animated=True,)
axes.set_xlabel('x (µm)')
axes.set_ylabel('u')
axes.set_xlim((-L,L))
axes.set_ylim((-L,L))
t2=time.time()
plotts=[]
color_cycle=axes._get_lines.prop_cycler
for k in range(Ns):
color=next(color_cycle)['color']
line=mcol.LineCollection(segments=[],color=color,animated=True)
axes.add_collection(line)
#plotts.append(line)
fig.suptitle('dynamics with '+r'$l_c=$'+str(l_cV/np.pi)[0:4]+', '+str(Ns)+' trajectories', fontsize=12)
skip=30
TT=list(np.concatenate([np.arange(0,int(max_K[k]/skip)+1) for k in range(len(e_L))]))
#p=-np.ones((1))
# an ad-hoc object, created to store efficiently the indice of the current value of epsilon...
class repere:
def __init__(self,val):
repere.value=val
p_rep=repere(-2)
#This function is properly horrible, but this is because all the tricky parts of the animation are inside
def plot_frame(i):
if i==0:
p=p_rep.value+1
setattr(p_rep,'value',p)
e=e_L[p]
force=(-(1-e)*nablaV+e*rotA)
Quiv.set_UVC(force[:,:,0].flatten(),force[:,:,1].flatten())
else:
p=p_rep.value
j=i*skip
for k in range(Ns):
segS=meta_traject[p][k]
KK=meta_K[p][k]
K_max=KK[-1]
if j<K_max:
seg=[]
q=0
j_prime=KK[q]
while j_prime<j:
seg.append(segS[q])
q+=1
j_prime=KK[q]
if q>0:
j_prime_old=KK[q-1]
else:
j_prime_old=0
seg.append(segS[q][0:j-j_prime_old])
axes.collections[k+1].set_segments(seg)
else:
axes.collections[k+1].set_segments(segS)
ti=dt*i*skip
t=str(ti)
if ti<0.1:
axes.set_title(r'$t=$'+t[0:7]+'s'+', '+r'$\varepsilon=$'+str(e_L[p]))
if ti<0.1:
axes.set_title(r'$t=$'+t[0:7]+'s'+', '+r'$\varepsilon=$'+str(e_L[p]))
elif ti<1:
axes.set_title(r'$t=$'+t[0:5]+'s'+', '+r'$\varepsilon=$'+str(e_L[p]))
elif ti<10:
axes.set_title(r'$t=$'+t[0:3]+'s'+', '+r'$\varepsilon=$'+str(e_L[p]))
else:
k=int(np.log10(ti))+1
axes.set_title('t='+t[0:k]+'s'+', '+r'$\varepsilon=$'+str(e_L[p]))
fig.canvas.draw()
return fig
# Animate the solution
animation=anim.FuncAnimation(fig, plot_frame,frames=TT, interval=100, repeat=False,blit=False)
animation.save("/Users/sayah/Desktop/3A/stage/code/videos/quiver_005_1.mp4", writer=writer)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment