Commit 4b35d04d authored by Antoine RAVETTA's avatar Antoine RAVETTA

reordering of the files

parent 48aefbc9
This diff is collapsed.
{
"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<ipython-input-1-bc37d7a89e1a>\u001b[0m in \u001b[0;36m<module>\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<ipython-input-4-4030111dc1b1>\u001b[0m in \u001b[0;36m<module>\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": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"ERROR:root:Internal Python error in the inspect module.\n",
"Below is the traceback from this internal error.\n",
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Traceback (most recent call last):\n",
" File \"/home/antoine/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py\", line 3325, in run_code\n",
" exec(code_obj, self.user_global_ns, self.user_ns)\n",
" File \"<ipython-input-9-de921bce0e9e>\", line 2, in <module>\n",
" from ..CN_classesCCCN import CrankNicolson\n",
"ValueError: attempted relative import beyond top-level package\n",
"\n",
"During handling of the above exception, another exception occurred:\n",
"\n",
"Traceback (most recent call last):\n",
" File \"/home/antoine/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py\", line 2039, in showtraceback\n",
" stb = value._render_traceback_()\n",
"AttributeError: 'ValueError' object has no attribute '_render_traceback_'\n",
"\n",
"During handling of the above exception, another exception occurred:\n",
"\n",
"Traceback (most recent call last):\n",
" File \"/home/antoine/anaconda3/lib/python3.7/site-packages/IPython/core/ultratb.py\", line 1101, in get_records\n",
" return _fixed_getinnerframes(etb, number_of_lines_of_context, tb_offset)\n",
" File \"/home/antoine/anaconda3/lib/python3.7/site-packages/IPython/core/ultratb.py\", line 319, in wrapped\n",
" return f(*args, **kwargs)\n",
" File \"/home/antoine/anaconda3/lib/python3.7/site-packages/IPython/core/ultratb.py\", line 353, in _fixed_getinnerframes\n",
" records = fix_frame_records_filenames(inspect.getinnerframes(etb, context))\n",
" File \"/home/antoine/anaconda3/lib/python3.7/inspect.py\", line 1502, in getinnerframes\n",
" frameinfo = (tb.tb_frame,) + getframeinfo(tb, context)\n",
" File \"/home/antoine/anaconda3/lib/python3.7/inspect.py\", line 1460, in getframeinfo\n",
" filename = getsourcefile(frame) or getfile(frame)\n",
" File \"/home/antoine/anaconda3/lib/python3.7/inspect.py\", line 696, in getsourcefile\n",
" if getattr(getmodule(object, filename), '__loader__', None) is not None:\n",
" File \"/home/antoine/anaconda3/lib/python3.7/inspect.py\", line 725, in getmodule\n",
" file = getabsfile(object, _filename)\n",
" File \"/home/antoine/anaconda3/lib/python3.7/inspect.py\", line 709, in getabsfile\n",
" return os.path.normcase(os.path.abspath(_filename))\n",
" File \"/home/antoine/anaconda3/lib/python3.7/posixpath.py\", line 383, in abspath\n",
" cwd = os.getcwd()\n",
"FileNotFoundError: [Errno 2] No such file or directory\n"
]
},
{
"ename": "ValueError",
"evalue": "attempted relative import beyond top-level package",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m"
]
}
],
"source": [
"import numpy as np\n",
"from ..CN_classesCCCN 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": 3,
"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",
"/home/antoine/Desktop/PHY571/codes/CN_classes/CCCN.py:42: RuntimeWarning: invalid value encountered in multiply\n",
" nu = lambda x: - self.delta_z / 4. / self.delta_r * self.V(x) # minus sign for V convention\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": 4,
"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)"
]
},
{
"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
}
This diff is collapsed.
This diff is collapsed.
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