From: Daniel Garcia-Sanchez Date: Fri, 5 Apr 2019 11:16:54 +0000 (+0200) Subject: Add Jupyter notebook X-Git-Tag: v9.1.0-rc1~148^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0378a543123709026454f3b3e6c59974333718b5;p=dealii.git Add Jupyter notebook --- diff --git a/contrib/python-bindings/notebooks/index.ipynb b/contrib/python-bindings/notebooks/index.ipynb index aac4262709..7a06cdc3c2 100644 --- a/contrib/python-bindings/notebooks/index.ipynb +++ b/contrib/python-bindings/notebooks/index.ipynb @@ -8,7 +8,11 @@ "This page contains a link to every notebook as well as a short explanation on what the notebooks do.\n", "\n", "## tutorial-1\n", - "[tutorial-1](https://github.com/dealii/dealii/tree/master/contrib/python-bindings/notebooks/tutorial-1.ipynb) shows how to create a **Triangulation** using the python binding. We also show how to refine the grid, merge two **Triangulations**, and finally how to output the **Triangulation** so that it can be loaded in a C++ code." + "[tutorial-1](https://github.com/dealii/dealii/tree/master/contrib/python-bindings/notebooks/tutorial-1.ipynb) shows how to create a **Triangulation** using the python binding. We also show how to refine the grid, merge two **Triangulations**, and finally how to output the **Triangulation** so that it can be loaded in a C++ code.\n", + "\n", + "## step-62\n", + "[step-62](https://github.com/dangars/dealii/tree/phononic-cavity/examples/step-62/step-62.ipynb) shows how to to calculate the [energy band gap](https://en.wikipedia.org/wiki/Band_gap) and the\n", + "mechanical resonance of a [micropillar superlattice cavity](https://doi.org/10.1103/PhysRevA.94.033813).\n" ] } ], diff --git a/contrib/python-bindings/notebooks/step-62.ipynb b/contrib/python-bindings/notebooks/step-62.ipynb new file mode 120000 index 0000000000..cdeb6e9202 --- /dev/null +++ b/contrib/python-bindings/notebooks/step-62.ipynb @@ -0,0 +1 @@ +../../../examples/step-62/step-62.ipynb \ No newline at end of file diff --git a/examples/step-62/step-62.ipynb b/examples/step-62/step-62.ipynb new file mode 100644 index 0000000000..0e318a3d59 --- /dev/null +++ b/examples/step-62/step-62.ipynb @@ -0,0 +1,785 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Step-62" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Set parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%config InlineBackend.figure_formats = ['svg']\n", + "%matplotlib inline\n", + "import numpy as np\n", + "import h5py\n", + "import matplotlib.pyplot as plt\n", + "import subprocess\n", + "import scipy.constants as constants\n", + "import scipy.optimize\n", + "\n", + "# This considerably reduces the size of the svg data\n", + "plt.rcParams['svg.fonttype'] = 'none'\n", + "\n", + "h5_file = h5py.File('results.h5', 'w')\n", + "data = h5_file.create_group('data')\n", + "displacement = data.create_group('displacement')\n", + "calibration = data.create_group('calibration')\n", + "\n", + "# Set the parameters\n", + "for group in [displacement, calibration]:\n", + " # Dimensions of the domain\n", + " group.attrs['dimension_x'] = 0.02\n", + " group.attrs['dimension_y'] = 2e-5\n", + " \n", + " # Position of the probe that we use to measure the flux\n", + " group.attrs['probe_pos_x'] = 0.008\n", + " group.attrs['probe_pos_y'] = 0\n", + " group.attrs['probe_width_y'] = 2e-05\n", + " \n", + " # Number of points in the probe\n", + " group.attrs['nb_probe_points'] = 5\n", + " \n", + " # Global refinement\n", + " group.attrs['grid_level'] = 1\n", + "\n", + " # Cavity\n", + " group.attrs['cavity_resonance_frequency'] = 20e6\n", + " group.attrs['nb_mirror_pairs'] = 15\n", + "\n", + " # Material\n", + " group.attrs['poissons_ratio'] = 0.27\n", + " group.attrs['youngs_modulus'] = 270000000000.0\n", + " group.attrs['material_a_rho'] = 3200\n", + " if group == displacement:\n", + " group.attrs['material_b_rho'] = 2000\n", + " else:\n", + " group.attrs['material_b_rho'] = 3200 \n", + " group.attrs['lambda'] = (group.attrs['youngs_modulus'] * group.attrs['poissons_ratio'] /\n", + " ((1 + group.attrs['poissons_ratio']) *\n", + " (1 - 2 * group.attrs['poissons_ratio'])))\n", + " group.attrs['mu']= (group.attrs['youngs_modulus'] / (2 * (1 + group.attrs['poissons_ratio'])))\n", + "\n", + " # Force\n", + " group.attrs['max_force_amplitude'] = 1e20\n", + " group.attrs['force_sigma_x'] = 1\n", + " group.attrs['force_sigma_y'] = 1\n", + " group.attrs['max_force_width_x'] = 0.0003\n", + " group.attrs['max_force_width_y'] = 0.001\n", + " group.attrs['force_x_pos'] = -0.008\n", + " group.attrs['force_y_pos'] = 0\n", + "\n", + " # PML\n", + " group.attrs['pml_x'] = True\n", + " group.attrs['pml_y'] = False\n", + " group.attrs['pml_width_x'] = 0.0018\n", + " group.attrs['pml_width_y'] = 0.0005\n", + " group.attrs['pml_coeff'] = 1.6\n", + " group.attrs['pml_coeff_degree'] = 2\n", + "\n", + " # Frequency sweep\n", + " group.attrs['center_frequency'] = 20e6\n", + " group.attrs['frequency_range'] = 0.5e6\n", + " group.attrs['start_frequency'] = group.attrs['center_frequency'] - group.attrs['frequency_range'] / 2\n", + " group.attrs['stop_frequency'] = group.attrs['center_frequency'] + group.attrs['frequency_range'] / 2\n", + " group.attrs['nb_frequency_points'] = 400\n", + "\n", + " # Other parameters\n", + " if group == displacement:\n", + " group.attrs['simulation_name'] = 'phononic_cavity_displacement'\n", + " else:\n", + " group.attrs['simulation_name'] = 'phononic_cavity_calibration'\n", + " group.attrs['save_vtu_files'] = False\n", + " \n", + "h5_file.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run simulation" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "CompletedProcess(args=['mpiexec', './step-62'], returncode=0)" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Now you can run your simulation locally or in a cluster\n", + "subprocess.run(['mpiexec','./step-62'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Analyze data" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 19.8\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 19.9\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 20.0\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 20.1\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 20.2\n", + " \n", + " \n", + " \n", + " frequency (MHz)\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 0.0\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 0.2\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 0.4\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 0.6\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 0.8\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 1.0\n", + " \n", + " \n", + " \n", + " amplitude^2 (a.u.)\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " Transmission\n", + " freq = 20.00815MHz Q = 5091.3\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 19.8\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 19.9\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 20.0\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 20.1\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 20.2\n", + " \n", + " \n", + " \n", + " frequency (MHz)\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " −1.0\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " −0.5\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 0.0\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 0.5\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 1.0\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " 1.5\n", + " \n", + " \n", + " \n", + " phase (rad)\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " Phase (transmission coefficient)\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "h5_file = h5py.File('results.h5', 'r')\n", + "data = h5_file['data']\n", + "\n", + "# Gaussian function that we use to fit the resonance\n", + "def resonance_f(freq, freq_m, quality_factor, max_amplitude):\n", + " omega = 2 * constants.pi * freq\n", + " omega_m = 2 * constants.pi * freq_m\n", + " gamma = omega_m / quality_factor\n", + " return max_amplitude * omega_m**2 * gamma**2 / (((omega_m**2 - omega**2)**2 + gamma**2 * omega**2))\n", + "\n", + "frequency = data['displacement']['frequency'][...]\n", + "# Average the probe points\n", + "displacement = np.mean(data['displacement']['displacement'], axis=0)\n", + "calibration_displacement = np.mean(data['calibration']['displacement'], axis=0)\n", + "reflection_coefficient = displacement / calibration_displacement\n", + "reflectivity = (np.abs(np.mean(data['displacement']['displacement'][...]**2, axis=0))/\n", + " np.abs(np.mean(data['calibration']['displacement'][...]**2, axis=0)))\n", + "\n", + "try:\n", + " x_data = frequency\n", + " y_data = reflectivity\n", + " quality_factor_guess = 1e3\n", + " freq_guess = x_data[np.argmax(y_data)]\n", + " amplitude_guess = np.max(y_data)\n", + " fit_result, covariance = scipy.optimize.curve_fit(resonance_f, x_data, y_data,\n", + " [freq_guess, quality_factor_guess, amplitude_guess])\n", + " freq_m = fit_result[0]\n", + " quality_factor = np.abs(fit_result[1])\n", + " max_amplitude = fit_result[2]\n", + " y_data_fit = resonance_f(x_data, freq_m, quality_factor, max_amplitude)\n", + "\n", + " fig = plt.figure()\n", + " plt.plot(frequency / 1e6, reflectivity, frequency / 1e6, y_data_fit)\n", + " plt.xlabel('frequency (MHz)')\n", + " plt.ylabel('amplitude^2 (a.u.)')\n", + " plt.title('Transmission\\n' + 'freq = ' + \"%.7g\" % (freq_guess / 1e6) + 'MHz Q = ' + \"%.6g\" % quality_factor)\n", + "except:\n", + " fig = plt.figure()\n", + " plt.plot(frequency / 1e6, reflectivity)\n", + " plt.xlabel('frequency (MHz)')\n", + " plt.ylabel('amplitude^2 (a.u.)')\n", + " plt.title('Transmission')\n", + "\n", + "fig = plt.figure()\n", + "plt.plot(frequency / 1e6, np.angle(reflection_coefficient))\n", + "plt.xlabel('frequency (MHz)')\n", + "plt.ylabel('phase (rad)')\n", + "plt.title('Phase (transmission coefficient)\\n')\n", + "\n", + "plt.show()\n", + "h5_file.close()" + ] + }, + { + "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.5.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}