Commit 70a4877b authored by Mathieu RASSON's avatar Mathieu RASSON

First filamentation 2 Pcr

parent 5c173a45
......@@ -7,9 +7,6 @@ import numpy as np
import scipy.sparse
import scipy.linalg as la
import matplotlib.pyplot as plt
%matplotlib inline
class CrankNicolsonFilamentation:
"""A class that solves the entire filamentation"""
......@@ -44,6 +41,8 @@ class CrankNicolsonFilamentation:
# Empty solution matrices (z,r,t)
self.E_matrix = np.zeros([self.n_z, self.n_r, self.n_t], dtype=complex)
self.rho_matrix = np.zeros([self.n_z, self.n_r, self.n_t])
#self.nonlin1 = np.zeros([self.n_z, self.n_r, self.n_t], dtype=complex)
#self.nonlin2 = np.zeros([self.n_z, self.n_r, self.n_t], dtype=complex)
# Sparse solver
Ar = self._fillAr_sp(sig_r, nu, self.n_r)
......@@ -85,9 +84,9 @@ class CrankNicolsonFilamentation:
# First fourth step
# t solver for each r with half non linear
fE1 = self.f(E) - self.sigma/2*(1 + 1j*self.omega0*self.tau) * rho[:,:] * E[:,:]
print('Defocus plasma :', np.amax(np.abs(self.sigma/2*(1 + self.omega0*self.tau) * rho[:,:] * E[:,:])**2))
#print('fE1 :', fE1)
#print('max fE1 :', np.amax(fE1))
#self.nonlin1[n] = fE1.copy()
#print('Focus :', self.f(E)[0,self.n_t//2] / E[0,self.n_t//2])
#print('Defocus :', fE1[0,self.n_t//2] / E[0,self.n_t//2])
# Non linear term at origin
if n==0:
fE_old = fE1.copy()
......@@ -109,8 +108,7 @@ class CrankNicolsonFilamentation:
# Second fourth step
# t solver for each r with half non linear
fE2 = self.f(E2) - self.sigma/2*(1 + 1j*self.omega0*self.tau) * rho[:,:] * E2[:,:]
#print('fE2 :', fE2)
#print('max fE2 :', np.amax(fE2))
#self.nonlin2[n] = fE2.copy()
# Half non linear term with half sum (old already defined even at n=0)
for i in range(self.n_r):
E[i,:] = la.solve_banded((1,1), At, Bt.dot(E2[i,:]) + \
......@@ -128,8 +126,8 @@ class CrankNicolsonFilamentation:
# Trapeze integral
a = np.exp(self.ava*self.delta_t * (np.abs(E[:,k+1])**2 + np.abs(E[:,k])**2)*0.5)
# Next t step calculation for next z step
rho[:,k+1] = a * (rho[:,k] + (self.delta_t/2) * self.ofi*np.abs(E[:,k])**(2*K)) \
+ (self.delta_t/2) * self.ofi*np.abs(E[:,k+1])**(2*K)
rho[:,k+1] = a * (rho[:,k] + (self.delta_t/2) * self.ofi*np.abs(E[:,k])**(2*self.K)) \
+ (self.delta_t/2) * self.ofi*np.abs(E[:,k+1])**(2*self.K)
# rho has been updated, put into rho_matrix at next step
......
......@@ -115,9 +115,13 @@
"metadata": {},
"outputs": [],
"source": [
"crank = CrankNicolsonFocus()\n",
"delta_z = 2e-4\n",
"z_max = 4\n",
"n_z = int(z_max/delta_z)\n",
"\n",
"crank.set_grid(r_max=2.8e-3, n_r=100, z_min=0, z_max=10, n_z=100, t_min=-5e-13, t_max=5e-13, n_t=100)\n",
"crank = CrankNicolsonFilamentation()\n",
"\n",
"crank.set_grid(r_max=2e-3, n_r=150, z_min=0, z_max=z_max, n_z=n_z, t_min=-2e-13, t_max=2e-13, n_t=150)\n",
"crank.set_parameters(Dr=diff_coeff, Dt=disp_coeff, V=potential, ofi=ofi, ava=ava ,K=K, sigma=sigma, \\\n",
" omega0=omega0, tau=tau, f=kerr_MPA)\n",
"\n",
......@@ -137,6 +141,9 @@
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"mat = initial_enveloppe(crank.r_pts, crank.t_pts, w0, tp, Pin)\n",
"\n",
"fig, ax = plt.subplots()\n",
......@@ -159,7 +166,8 @@
"metadata": {},
"outputs": [],
"source": [
"E, rho = crank.get_data()"
"E, rho = crank.get_data()\n",
"#E, rho = crank.E_matrix.copy(), crank.rho_matrix.copy()"
]
},
{
......@@ -169,10 +177,12 @@
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"mesh = ax.pcolormesh(np.real(crank.r_pts), np.real(crank.z_pts), np.abs(E[:,:,crank.n_t//2])**2)\n",
"#mesh = ax.pcolormesh(np.real(crank.r_pts), np.real(crank.z_pts), np.abs(E[:,:,crank.n_t//2])**2)\n",
"mesh = ax.pcolormesh(np.abs(E[:,:,crank.n_t//2])**2)\n",
"#mesh = ax.pcolormesh(np.log(np.abs(E[:,:,crank.n_t//2])**2))\n",
"plt.colorbar(mesh, ax=ax)\n",
"#ax.set_ylim(1.29999999,1.30000001)\n",
"#ax.set_xlim(0, 4e-5)\n",
"#ax.set_ylim(2,2.5)\n",
"ax.set_xlim(0, 20)\n",
"#plt.savefig('kerr', dpi=300)\n",
"fig.show()"
]
......@@ -183,10 +193,48 @@
"metadata": {},
"outputs": [],
"source": [
"listemax = np.array([np.amax(np.abs(E[i,:,:])**2) for i in range(crank.n_z)])\n",
"#listemax = np.array([np.amax(np.abs(E[i,:,:])**2) for i in range(crank.n_z)])\n",
"listemax = np.array([np.amax(np.abs(E[i,2,:])**2) for i in range(crank.n_z)])\n",
"\n",
"maxi = crank.n_z\n",
"\n",
"plt.figure()\n",
"plt.plot(np.real(crank.z_pts), listemax)\n",
"plt.plot(np.real(crank.z_pts[:maxi]), listemax[:maxi]/Pcr)\n",
"plt.xlabel('Propagation distance (m)')\n",
"plt.ylabel('Squared field amplitude ($V^2 \\cdot m^{-2}$)')\n",
"plt.savefig('filamentation_2P.png', dpi=300)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('Kerr: ', 1/(n2*k*3e17))\n",
"print('Plasma: ', 1/(sigma/2*omega0*tau*5e22))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The two caracteristic distances are around 5 mm. Therefore $\\Delta z$ has to be chosen around 0.5 mm and 0.1 mm is even better."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"listemax = np.array([np.amax(rho[i,:,:]) for i in range(crank.n_z)])\n",
"\n",
"maxi = crank.n_z\n",
"\n",
"plt.figure()\n",
"plt.plot(np.real(crank.z_pts[:maxi]), listemax[:maxi])\n",
"plt.show()"
]
},
......@@ -197,13 +245,31 @@
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"mesh = ax.pcolormesh(np.real(crank.r_pts), np.real(crank.z_pts), np.abs(rho[:,:,-1])**2)\n",
"mesh = ax.pcolormesh(np.real(crank.r_pts), np.real(crank.z_pts), rho[:,:,-1])\n",
"plt.colorbar(mesh, ax=ax)\n",
"#ax.set_ylim(1.29999999,1.30000001)\n",
"#ax.set_xlim(0, 4e-5)\n",
"#plt.savefig('kerr', dpi=300)\n",
"fig.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.plot(np.abs(E[4500,:,crank.n_t//2])**2)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
......
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