Commit 2bbbda1d by Antoine RAVETTA

### verify Marburger's law

parent fda2201b
 """ # 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