Commit 1162d74e authored by Mathieu RASSON's avatar Mathieu RASSON

Info sync

parents 008e9450 8c012b8d
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Crank-Nicholson scheme"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Time dispersion term"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Add time dispersion term as a new diffusion source $D_2$.\n",
"* Add one dimension to tables;\n",
"* Add one parameter to initialisation;\n",
"* Border conditions: null conditions at $t_{min}$ and $t_{max}$;\n",
"* Fill another scheme matrix;\n",
"* Split step into two half steps: fill one dimension for each step inside the $z$ loop, separately."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.sparse\n",
"import scipy.linalg as la\n",
"\n",
"class CrankNicolson:\n",
" \"\"\"A class that solves dE/dz = D*d2E/dr2 + V(r)*dE/dr + f(E)\"\"\"\n",
" \n",
" # Cylindrical grid\n",
" def set_grid(self, r_max, n_r, z_min, z_max, n_z):\n",
"\n",
" self.r_max, self.n_r = r_max, n_r\n",
" self.z_min, self.z_max, self.n_z = z_min, z_max, n_z\n",
" self.r_pts, self.delta_r = np.linspace(0, r_max, n_r, retstep=True, endpoint=False)\n",
" self.z_pts, self.delta_z = np.linspace(z_min, z_max, n_z, retstep=True, endpoint=False)\n",
" \n",
" # Parameters of the scheme\n",
" def set_parameters(self, D, V, f):\n",
" \n",
" # V has to be vectorised\n",
" self.D, self.V, self.f = D, V, f\n",
" \n",
" \n",
" \n",
" # One solving step for r dependency\n",
" def solve(self, E_init):\n",
" \n",
" # Coefficient of matrices\n",
" sig = self.D * self.delta_z / 2. / self.delta_r**2\n",
" nu = -self.V * self.delta_z / 4. / self.delta_r # minus sign for V convention\n",
" \n",
" # Empty solution matrix\n",
" self.E_matrix = np.zeros([self.n_z, self.n_r], dtype=complex)\n",
" \n",
" # Sparse solver\n",
" A = self._fillA_sp(sig, nu, self.n_r)\n",
" B = self._fillB_sp(sig, nu, self.n_r)\n",
" \n",
" # Set boundary conditions \n",
" # Dirichlet at infinity\n",
" A[1,-1] = 1.0\n",
" A[2,-2] = 0.0\n",
" B[-1,-1] = 0.0\n",
" B[-1,-2] = 0.0\n",
" # Neumann at r=0\n",
" A[0,1] = -2*sig\n",
" B[0,1] = 2*sig\n",
" \n",
" # Propagate\n",
" E = E_init\n",
" for n in range(self.n_z):\n",
" self.E_matrix[n,:] = E\n",
" fE = f(E)\n",
" # Non linear term at origin\n",
" if n==0:\n",
" fE_old = fE\n",
" \n",
" # Non linear term with half sum\n",
" E = la.solve_banded((1,1),A, B.dot(E) + self.delta_z * (1.5 * fE - 0.5 * fE_old),\\\n",
" check_finite=False)\n",
" fE_old = fE\n",
" \n",
" \n",
" \n",
" # Deep copy of results for export\n",
" # np.savetxt ?\n",
" def get_E(self):\n",
" \n",
" return self.E_matrix.copy()\n",
" \n",
" # Diagonal packing for banded\n",
" def _fillA_sp(self, sig, nu, n):\n",
" \"\"\"Returns a tridiagonal matrix in compact form ab[1+i-j,j]=a[i,j]\"\"\"\n",
" \n",
" A = np.zeros([3,n], dtype=complex) # A has three diagonals and size n\n",
" # Superdiagonal\n",
" A[0,1:] = -(sig - nu(self.r_pts[:-1]))\n",
" # Diagonal\n",
" A[1] = 1+2*sig\n",
" # Subdiagonal\n",
" A[2,:-1] = -(sig + nu(self.r_pts[1:]))\n",
" \n",
" return A\n",
" \n",
" # Sparse tridiagonal storage\n",
" def _fillB_sp(self, sig, nu, n):\n",
" \"\"\"Returns a tridiagonal sparse matrix in csr-form\"\"\"\n",
" \n",
" _o = np.ones(n, dtype=complex)\n",
" supdiag = (sig - nu(self.r_pts[:-1]))\n",
" diag = (1-2*sig)*_o\n",
" subdiag = (sig + nu(self.r_pts[1:]))\n",
" \n",
" return scipy.sparse.diags([supdiag, diag, subdiag], [1,0,-1], (n,n), format=\"csr\")\n",
" \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Non linear terms"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute $\\mathcal{E} \\mapsto f(\\mathcal{E})$ to add non linear terms. Test one after another, with divergence if only Kerr effect."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}
"""
### Adapted crank Nicholson class.
* Prepare half-step method.
* Prepare non linear term for one of the half-step.
* Use $V$ advection to add first radial derivative.
$V$ is a function of $r$
$$
V(r) = \dfrac{i}{2 k_0} \dfrac{1}{r}
$$
and $V(0) = 0$.
**Warning**: do not use $V(0)=0$ for any calculation.
**Warning**: Test scheme matrices separately.
"""
import numpy as np
import scipy.sparse
import scipy.linalg as la
class CrankNicolson:
"""A class that solves dE/dz = D*d2E/dr2 + V(r)*dE/dr"""
# Cylindrical grid
def set_grid(self, r_max, n_r, z_min, z_max, n_z):
self.r_max, self.n_r = r_max, n_r
self.z_min, self.z_max, self.n_z = z_min, z_max, n_z
self.r_pts, self.delta_r = np.linspace(0, r_max, n_r, retstep=True, endpoint=False)
self.z_pts, self.delta_z = np.linspace(z_min, z_max, n_z, retstep=True, endpoint=False)
# Parameters of the scheme
def set_parameters(self, D, V):
# V has to be vectorised
self.D, self.V = D, V
# One solving step for r dependency
def solve(self, E_init):
# Coefficient of matrices
sig = self.D * self.delta_z / 2. / self.delta_r ** 2
nu = lambda x: - self.delta_z / 4. / self.delta_r * self.V(x) # minus sign for V convention
# Empty solution matrix
self.E_matrix = np.zeros([self.n_z, self.n_r], dtype=complex)
# Sparse solver
A = self._fillA_sp(sig, nu, self.n_r)
B = self._fillB_sp(sig, nu, self.n_r)
# Set boundary conditions
# Dirichlet at infinity
A[1, -1] = 1.0
A[2, -2] = 0.0
B[-1, -1] = 0.0
B[-1, -2] = 0.0
# Neumann at r=0
A[0, 1] = -2 * sig
B[0, 1] = 2 * sig
# Propagate
E = E_init.copy()
for n in range(self.n_z):
self.E_matrix[n, :] = E
E = la.solve_banded((1, 1), A, B.dot(E), check_finite=False)
# Deep copy of results for export
# np.savetxt ?
def get_E(self):
return self.E_matrix.copy()
# Diagonal packing for banded
def _fillA_sp(self, sig, nu, n):
"""Returns a tridiagonal matrix in compact form ab[1+i-j,j]=a[i,j]"""
A = np.zeros([3, n], dtype=complex) # A has three diagonals and size n
# Superdiagonal
A[0, 1:] = -(sig - nu(self.r_pts[:-1]))
# Diagonal
A[1] = 1 + 2 * sig
# Subdiagonal
A[2, :-1] = -(sig + nu(self.r_pts[1:]))
return A
# Sparse tridiagonal storage
def _fillB_sp(self, sig, nu, n):
"""Returns a tridiagonal sparse matrix in csr-form"""
_o = np.ones(n, dtype=complex)
supdiag = (sig - nu(self.r_pts[:-1]))
diag = (1 - 2 * sig) * _o
subdiag = (sig + nu(self.r_pts[1:]))
return scipy.sparse.diags([supdiag, diag, subdiag], [1, 0, -1], (n, n), format="csr")
\ No newline at end of file
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Data Production\n",
"\n",
"We define some constants for the whole problem, as well as some functions :\n",
"\n",
"- wavevector\n",
"- potential\n",
"- initialisation of E"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from CCCN import CrankNicolson\n",
"\n",
"lambd = 0.7e-6\n",
"\n",
"k = 2 * np.pi / lambd\n",
"\n",
"w0 = 1e-3\n",
"Pin = 1\n",
"\n",
"diff_coeff = 1j*1/(2 * k)\n",
"#diff_coeff = 0\n",
"#diff_coeff = 1e-4\n",
"\n",
"print('diff :', diff_coeff)\n",
"\n",
"\n",
"def potential(r):\n",
" try:\n",
" return diff_coeff / r\n",
" except ZeroDivisionError:\n",
" return 0\n",
"\n",
"\n",
"def gaussian(r, r0=0, w0=1, Pin=1):\n",
" return np.sqrt(2*Pin/(np.pi*w0**2)) * np.exp(-(r-r0)** 2/(w0**2))\n",
"\n",
"\n",
"def initial_enveloppe(r_pts, w0, Pin):\n",
" return np.array([gaussian(r_pts[i], w0=w0, Pin=Pin) for i in range(len(r_pts))])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Instanciation of the CN Class"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"crank = CrankNicolson()\n",
"\n",
"crank.set_grid(r_max=1e-2, n_r=100, z_min=0, z_max=10, n_z=200)\n",
"crank.set_parameters(D=diff_coeff, V=potential)\n",
"\n",
"crank.solve(initial_enveloppe(crank.r_pts, w0, Pin))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Save the result for later analysis"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"np.savetxt(\"../CCCN_E.dat\", np.abs(crank.E_matrix))\n",
"np.savetxt(\"../CCCN_r_pts.dat\", crank.r_pts)\n",
"np.savetxt(\"../CCCN_z_pts.dat\", crank.z_pts)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
"""
Solution of filamentation with charge density coupled
"""
import numpy as np
import scipy.sparse
import scipy.linalg as la
class CrankNicolsonFilamentation:
"""A class that solves the entire filamentation"""
# Cylindrical grid with rectangular time
def set_grid(self, r_max, n_r, z_min, z_max, n_z, t_min, t_max, n_t):
self.r_max, self.n_r = r_max, n_r
self.z_min, self.z_max, self.n_z = z_min, z_max, n_z
self.t_min, self.t_max, self.n_t = t_min, t_max, n_t
self.r_pts, self.delta_r = np.linspace(0, r_max, n_r, retstep=True, endpoint=False, dtype=complex)
self.z_pts, self.delta_z = np.linspace(z_min, z_max, n_z, retstep=True, endpoint=False, dtype=complex)
self.t_pts, self.delta_t = np.linspace(t_min, t_max, n_t, retstep=True, endpoint=False, dtype=complex)
# Parameters of the scheme
# Dr for diffraction, Dt for dispersion
def set_parameters(self, Dr, Dt, V, ofi, ava, K, sigma, omega0, tau, f):
# V has to be vectorised
self.Dr, self.Dt, self.V, self.ofi, self.ava, self.K, self.sigma, self.omega0, self.tau, self.f = \
Dr, Dt, V, ofi, ava, K, sigma, omega0, tau, f
# One solving step for r dependency and t dependency
def solve(self, E_init):
# Coefficient of matrices
sig_r = self.Dr * (self.delta_z) / 2. / self.delta_r**2
sig_t = self.Dt * (self.delta_z/2) / 2. / (self.delta_t)**2 # (half t operator)
nu = lambda x : -self.V(x) * (self.delta_z) / 4. / self.delta_r # minus sign for V convention
# 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)
Br = self._fillBr_sp(sig_r, nu, self.n_r)
At = self._fillAt_sp(sig_t, self.n_t) # null advection coefficient
Bt = self._fillBt_sp(sig_t, self.n_t) # null advection coefficient
# Set boundary conditions
# Spatial:
# Dirichlet at infinity
Ar[1,-1] = 1.0
Ar[2,-2] = 0.0
Br[-1,-1] = 0.0
Br[-1,-2] = 0.0
# Neumann at r=0
Ar[0,1] = -2*sig_r
Br[0,1] = 2*sig_r
# Time
# Neumann at t=t_min, t_max
for b in [0,1]:
At[1,-b] = 1.0
At[2*b,1-3*b] = 0.0
Bt[-b,-b] = 0.0
Bt[-b,1-3*b] = 0.0
# Propagate
E = E_init.copy()
E1 = E_init.copy()
E2 = E_init.copy()
rho = np.zeros([self.n_r, self.n_t])
for n in range(self.n_z):
# General step initialisation
self.E_matrix[n,:,:] = E
self.rho_matrix[n,:,:] = rho
# 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[:,:]
#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()
# Half non linear term with half sum
for i in range(self.n_r):
E1[i,:] = la.solve_banded((1,1), At, Bt.dot(E[i,:]) + \
(self.delta_z/2.) * (1.5 * fE1[i,:] - 0.5 * fE_old[i,:]), \
check_finite=False)
# Non linear term update
fE_old = fE1.copy()
# Middle half step
# r solver solver for each t
for i in range(self.n_t):
E2[:,i] = la.solve_banded((1,1), Ar, Br.dot(E1[:,i]), \
check_finite=False)
# 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[:,:]
#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,:]) + \
(self.delta_z/2.) * (1.5 * fE2[i,:] - 0.5 * fE_old[i,:]), \
check_finite=False)
# Non linear term update
fE_old = fE2.copy()
# E has been updated, put into E_matrix at next step
# Density coupling (notations article)
# Vectorised Calculation for all r
# Time integration from null conditions
for k in range(self.n_t-1):
# 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*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
# Deep copy of results for export
# np.savetxt ?
def get_data(self):
return self.E_matrix.copy(), self.rho_matrix.copy()
# Diagonal packing for banded r
def _fillAr_sp(self, sig, nu, n):
"""Returns a tridiagonal matrix in compact form ab[1+i-j,j]=a[i,j]"""
A = np.zeros([3,n], dtype=complex) # A has three diagonals and size n
# Superdiagonal
A[0,1:] = -(sig - nu(self.r_pts[:-1]))
# Diagonal
A[1] = 1+2*sig
# Subdiagonal
A[2,:-1] = -(sig + nu(self.r_pts[1:]))
return A
# Diagonal packing for banded t
def _fillAt_sp(self, sig, n):
"""Returns a tridiagonal matrix in compact form ab[1+i-j,j]=a[i,j]"""
_o = np.ones(n, dtype=complex)
A = np.zeros([3,n], dtype=complex) # A has three diagonals and size n
# Superdiagonal
A[0,1:] = -sig
# Diagonal
A[1] = 1+2*sig
# Subdiagonal
A[2,:-1] = -sig
return A
# Sparse tridiagonal storage for r
def _fillBr_sp(self, sig, nu, n):
"""Returns a tridiagonal sparse matrix in csr-form"""
_o = np.ones(n, dtype=complex)
supdiag = (sig - nu(self.r_pts[:-1]))
diag = (1-2*sig)*_o
subdiag = (sig + nu(self.r_pts[1:]))
return scipy.sparse.diags([supdiag, diag, subdiag], [1,0,-1], (n,n), format="csr")
# Sparse tridiagonal storage for t
def _fillBt_sp(self, sig, n):
"""Returns a tridiagonal sparse matrix in csr-form"""
_o = np.ones(n, dtype=complex)
supdiag = sig*_o[:-1]
diag = (1-2*sig)*_o
subdiag = sig*_o[:-1]
return scipy.sparse.diags([supdiag, diag, subdiag], [1,0,-1], (n,n), format="csr")
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Functions"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def potential(r):\n",
" try:\n",
" return diff_coeff / r\n",
" except ZeroDivisionError:\n",
" return 0\n",
"\n",
"\n",
"def gaussian(r, t, r0=0, w0=1, t0=0, tp=1, Pin=1):\n",
" return np.sqrt(2*Pin/(np.pi*w0**2)) * np.exp(-(r-r0)** 2/(w0**2) - (t-t0)** 2/(tp**2))\n",
"\n",
"\n",
"def initial_enveloppe(r_pts, t_pts, w0, tp, Pin):\n",
" return np.array([[gaussian(r_pts[i], t_pts[j], w0=w0, tp=tp, Pin=Pin) for j in range(len(t_pts))] \\\n",
" for i in range(len(r_pts))])\n",
"\n",
"#----- Non linearities --------------\n",
"\n",
"def kerr(E):\n",
" return 1j * k * n2 * np.abs(E)**2 * E\n",
"\n",
"def answer(t):\n",
" pass\n",
" \n",
"def kerr_delay(E):\n",
" pass\n",
" \n",
"\n",
"def kerr_MPA(E):\n",
" return 1j*k*n2*np.abs(E)**2 * E - beta/2 * np.abs(E)**(2*K-2) * E"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Parameters"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from FCN import CrankNicolsonFilamentation\n",
"\n",
"# Physics\n",
"c = 3e8\n",
"mu = 4*np.pi*1e-7\n",
"\n",
"# Beam\n",
"lambd = 0.7e-6\n",
"w0 = 1e-3\n",
"tp = 85e-15\n",
"#tp = 1\n",
"Pcr = 1.8e9\n",
"#Pin = 5 * Pcr\n",
"Pin = 2*Pcr\n",
"\n",
"# Dispersion\n",
"k = 2 * np.pi / lambd\n",
"k_2 = 2e-28\n",
"\n",
"# Kerr Effect\n",
"n2 = 5.57e-23\n",
"#n2 = 7e-23\n",
"\n",
"# MPA\n",
"K = 7\n",
"beta = 6.5e-104\n",
"\n",
"# Density\n",
"omega0 = k * c\n",
"ofi = 8.84e-71/omega0\n",
"#ofi = 0\n",
"ava = 2.9e-6\n",
"#ava = 0\n",
"sigma = 5.1e-24\n",
"tau = 3.5e-13\n",
"\n",
"diff_coeff = 1j*1/(2 * k)\n",
"#diff_coeff = 0\n",
"#diff_coeff = 1e-4\n",
"disp_coeff = -1j*k_2/2\n",
"#disp_coeff = 1j * 1e-25\n",
"#disp_coeff = 0\n",
"\n",
"print('diff :', diff_coeff)\n",
"print('disp :', disp_coeff)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Instantiation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"delta_z = 3e-3\n",
"z_max = 2.7\n",
"n_z = int(z_max/delta_z)\n",
"\n",
"crank = CrankNicolsonFilamentation()\n",
"\n",