...

Commits (4)
 ... ... @@ -476,21 +476,20 @@ "source": [ "# Energy evolution\n", "\n", "ener = np.zeros(crank.n_z)\n", "for i in range(crank.n_z):\n", " # Filled at each z\n", " for j in range(crank.n_r):\n", " # Term integrated along t at r added against r value\n", " ener[i] += crank.delta_t * np.sum(np.abs(E[i,j,:])**2) * np.real(crank.r_pts[j])\n", "# r integration step\n", "ener *= crank.delta_r\n", "ener = np.array([crank.delta_r*crank.delta_t*np.sum(np.abs(E[i,:,:])**2) for i in range(crank.n_z)])\n", "\n", "plt.figure()\n", "plt.plot(np.real(crank.z_pts), ener)\n", "plt.ylim(0,2e-14)\n", "plt.plot(ener)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, ... ... @@ -718,9 +717,9 @@ "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*k_2/2\n", "#disp_coeff = 1j * 1e-25\n", "#disp_coeff = 0\n", "disp_coeff = 0\n", "\n", "print('diff :', diff_coeff)\n", "print('disp :', disp_coeff)\n", ... ... @@ -847,29 +846,6 @@ "fig.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Energy evolution\n", "\n", "ener = np.zeros(crank.n_z)\n", "for i in range(crank.n_z):\n", " # Filled at each z\n", " for j in range(crank.n_r):\n", " # Term integrated along t at r added against r value\n", " ener[i] += crank.delta_t * np.sum(np.abs(E[i,j,:])**2) * np.real(crank.r_pts[j])\n", "# r integration step\n", "ener *= crank.delta_r\n", "\n", "plt.figure()\n", "plt.plot(np.real(crank.z_pts), ener)\n", "plt.ylim(0,2e-14)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, ... ... @@ -1022,18 +998,10 @@ "source": [ "# Energy evolution\n", "\n", "ener = np.zeros(crank.n_z)\n", "for i in range(crank.n_z):\n", " # Filled at each z\n", " for j in range(crank.n_r):\n", " # Term integrated along t at r added against r value\n", " ener[i] += crank.delta_t * np.sum(np.abs(E[i,j,:])**2) * np.real(crank.r_pts[j])\n", "# r integration step\n", "ener *= crank.delta_r\n", "ener = np.array([crank.delta_r*crank.delta_t*np.sum(np.abs(E[i,:,:])**2) for i in range(crank.n_z)])\n", "\n", "plt.figure()\n", "plt.plot(np.real(crank.z_pts), ener)\n", "plt.ylim(0,2e-5)\n", "plt.plot(ener)\n", "plt.show()" ] }, ... ...
 { "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": 1, "metadata": {}, "outputs": [ { "ename": "ModuleNotFoundError", "evalue": "No module named 'CrankNicolson'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mCrankNicolson\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mlambd\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.7e-6\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'CrankNicolson'" ] } ], "source": [ "import numpy as np\n", "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": 4, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'CrankNicolson' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mcrank\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mCrankNicolson\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mcrank\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_grid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr_max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1e-2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_r\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz_min\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz_max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_z\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m200\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mcrank\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_parameters\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mD\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdiff_coeff\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mV\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpotential\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'CrankNicolson' is not defined" ] } ], "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(\"../CN_cylindric_complex_E.dat\", np.abs(crank.E_matrix))\n", "np.savetxt(\"../CN_cylindric_complex_r_pts.dat\", crank.r_pts)\n", "np.savetxt(\"../CN_cylindric_complex_z_pts.dat\", crank.z_pts)" ] } ], "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 }
 { "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": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "diff : 5.570423008216336e-08j\n" ] } ], "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": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/antoine/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:20: RuntimeWarning: divide by zero encountered in true_divide\n", "/home/antoine/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:20: RuntimeWarning: invalid value encountered in true_divide\n" ] } ], "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": 7, "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": [] }, { "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 }
 """ ### Adapted class Term added with new table dimension: $$D_t \dfrac{\partial^2 \mathcal{E}}{\partial t^2}$$ Matrix similar to TP4 without advection. 3D matrix format: $$E(z_i, r_j, t_k)$$ **Warning**: Non linear term treated with **space half step**. Initial conditions at $z=0$ in **2D**. Half step $\Delta z/2$ in the scheme. """ import numpy as np import scipy.sparse import scipy.linalg as la class CrankNicolsonDisp: """A class that solves dE/dz = D1*d2E/dr2 + V(r)*dE/dr + D2*d2E/dt2 + f(E)""" # 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): # V has to be vectorised self.Dr, self.Dt, self.V = Dr, Dt, V # 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. / self.delta_t**2 nu = lambda x : -self.V(x) * (self.delta_z) / 4. / self.delta_r # minus sign for V convention # Empty solution matrix (z,r,t) self.E_matrix = 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() E_intermediate = E_init.copy() for n in range(self.n_z): # General step initialisation self.E_matrix[n,:,:] = E # First half step # r operator solver for each t for i in range(self.n_t): E_intermediate[:,i] = la.solve_banded((1,1), Ar, Br.dot(E[:,i]), \ check_finite=False) # Second half step # t operator solver for each r for i in range(self.n_r): E[i,:] = la.solve_banded((1,1), At, Bt.dot(E_intermediate[i,:]), \ check_finite=False) # E has been updated, put into E_matrix at next step # Deep copy of results for export # np.savetxt ? def get_E(self): return self.E_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") \ No newline at end of file
 { "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "diff : 5.570423008216336e-08j\n", "disp : -1e-28j\n" ] } ], "source": [ "import numpy as np\n", "from HSCN import CrankNicolsonDisp\n", "\n", "lambd = 0.7e-6\n", "\n", "k = 2 * np.pi / lambd\n", "k_2 = 2e-28\n", "\n", "w0 = 1e-3\n", "tp = 85e-15\n", "#tp = 1\n", "Pin = 1\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)\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, 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))])" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/antoine/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:27: RuntimeWarning: divide by zero encountered in true_divide\n", "/home/antoine/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:27: RuntimeWarning: invalid value encountered in true_divide\n" ] } ], "source": [ "crank = CrankNicolsonDisp()\n", "\n", "crank.set_grid(r_max=1e-2, 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.set_parameters(Dr=diff_coeff, Dt=disp_coeff, V=potential)\n", "\n", "crank.solve(initial_enveloppe(crank.r_pts, crank.t_pts, w0, tp, Pin))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "np.save(\"../HSCN_E.npy\", np.abs(crank.E_matrix))\n", "np.save(\"../HSCN_r_pts.npy\", crank.r_pts)\n", "np.save(\"../HSCN_z_pts.npy\", crank.z_pts)\n", "np.save(\"../HSCN_t_pts.npy\", crank.t_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 }
 """ # Non linear terms Compute $\mathcal{E} \mapsto f(\mathcal{E})$ to add non linear terms. Test one after another, with divergence if only Kerr effect. Add non linear integrations in a symmetric manner: * Fourth $z$ step and half $t$ step with half non linear terms; * Half $z$ step and total $r$ step; * Fourth $z$ step and half $t$ step with half non linear terms. **Warning**: do not forget $\Delta z/4$, $\Delta z/2$ and $\Delta t/2$ factors to obtain the right scale. """ import numpy as np import scipy.sparse import scipy.linalg as la import matplotlib.pyplot as plt class CrankNicolsonFocus: """A class that solves dE/dz = D1*d2E/dr2 + V(r)*dE/dr + D2*d2E/dt2 + f(E)""" # 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, f): # V has to be vectorised self.Dr, self.Dt, self.V, self.f = Dr, Dt, V, 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 matrix (z,r,t) self.E_matrix = 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() for n in range(self.n_z): # General step initialisation self.E_matrix[n,:,:] = E # First fourth step # t solver for each r with half non linear fE1 = self.f(E) #print('fE1 :', fE1) #print('max fE1 :', np.amax(fE1)) # 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) #print('fE2 :', fE2) #print('max fE2 :', np.amax(fE2)) # 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 # Deep copy of results for export # np.savetxt ? def get_E(self): return self.E_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") \ No newline at end of file
This diff is collapsed.