From: Daniel Garcia-Sanchez Date: Thu, 25 Apr 2019 18:05:53 +0000 (+0200) Subject: Address Wolfgang's comments X-Git-Tag: v9.1.0-rc1~148^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=31e1658c52c2aef8fa81b960be6af9e7da814835;p=dealii.git Address Wolfgang's comments --- diff --git a/contrib/python-bindings/notebooks/index.ipynb b/contrib/python-bindings/notebooks/index.ipynb index 7a06cdc3c2..4611e00ecd 100644 --- a/contrib/python-bindings/notebooks/index.ipynb +++ b/contrib/python-bindings/notebooks/index.ipynb @@ -11,7 +11,7 @@ "[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", + "[step-62](https://github.com/dealii/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/examples/step-62/CMakeLists.txt b/examples/step-62/CMakeLists.txt index 28f4ae6531..5ebaaa6b42 100644 --- a/examples/step-62/CMakeLists.txt +++ b/examples/step-62/CMakeLists.txt @@ -1,5 +1,5 @@ ## -# CMake script for the step-8 tutorial program: +# CMake script for the step-62 tutorial program: ## # Set the name of the project and target: diff --git a/examples/step-62/doc/intro.dox b/examples/step-62/doc/intro.dox index 3483a4ee09..67bcad9988 100644 --- a/examples/step-62/doc/intro.dox +++ b/examples/step-62/doc/intro.dox @@ -10,39 +10,57 @@ together with these additional libraries is described in the README file.

Introduction

-In this tutorial we calculate the -[energy gap](https://en.wikipedia.org/wiki/Band_gap) and the -mechanical resonance of a -[phononic superlattice cavity](https://doi.org/10.1103/PhysRevA.94.033813). - +A phononic crystal is a periodic nanostructure that modifies the motion of +mechanical vibrations or [phonons](https://en.wikipedia.org/wiki/Phonon). +Phononic structures can be used to disperse, route and confine mechanical vibrations. +These structures have potential applications in +[quantum information](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.86.1391) +and have been used to study +[macroscopic quantum phenomena](https://science.sciencemag.org/content/358/6360/203). +Phononic crystals are usually fabricated in +[cleanrooms](https://en.wikipedia.org/wiki/Cleanroom). -A phononic superlattice cavity is formed by two -[Distributed Reflector](https://en.wikipedia.org/wiki/Band_gap), +In this tutorial we show how to a design a +[phononic superlattice cavity](https://doi.org/10.1103/PhysRevA.94.033813) +which is a particular type of phononic crystal that can be used to confine +mechanical vibrations. A phononic superlattice cavity is formed by two +[Distributed Bragg Reflector](https://en.wikipedia.org/wiki/Distributed_Bragg_reflector), mirrors and a $\lambda/2$ cavity where $\lambda$ is the acoustic wavelength. Acoustic DBRs are periodic structures where a set of bilayer stacks with contrasting physical properties (sound velocity index) is repeated $N$ times. +Superlattice cavities are usually grown on a +[Gallium Arsenide](https://en.wikipedia.org/wiki/Gallium_arsenide) +wafer by +[Molecular Beam Epitaxy](https://en.wikipedia.org/wiki/Molecular-beam_epitaxy). +The bilayers correspond to GaAs/AlAs mirror pairs. As shown below, the thickness of the mirror layers (brown and green) is $\lambda/4$ and the thickness of the cavity (blue) is $\lambda/2$. +Phononic superlattice cavity -Phononic superlattice cavity +In this tutorial we calculate the +[band gap](https://en.wikipedia.org/wiki/Band_gap) and the +mechanical resonance of a phononic superlattice cavity but the code presented here +can be easily used to design and calculate other types of +[phononic crystals](https://science.sciencemag.org/content/358/6360/203). The device is a waveguide in which the wave goes from left to right. -The simulations of this tutorial are done in 2D; -although because we use templates it is very easy to convert this program to 3D. +The simulations of this tutorial are done in 2D, but the code is dimension +independent and can be easily used with 3D simulations. +The waveguide width is equal to the $y$ dimension of the domain and the +waveguide length is equal to the $x$ dimension of the domain. There are two regimes that depend on the waveguide width: - Single mode: In this case the width of the structure is much - smaller than the wavelength, therefore the waveguide is single mode. + smaller than the wavelength. This case can be solved either with FEM (the approach that we take here) or with a simple semi-analytical [1D transfer matrix formalism](https://en.wikipedia.org/wiki/Transfer_matrix). -- Multimode: In this case the width of the structure is larger than the - wavelength, therefore the waveguide is multimode. +- Multimode: In this case the width of the structure is larger than the wavelength. This case can be solved using FEM or with a [scattering matrix formalism](https://doi.org/10.1103/PhysRevA.94.033813). Although we do not study this case in this tutorial, it is very easy to reach the multimode - case by increasing the parameter waveguide width (`dimension_y` in the jupyter notebook). + regime by increasing the parameter waveguide width (`dimension_y` in the jupyter notebook). The simulations of this tutorial are performed in the frequency domain. To calculate the transmission spectrum, we use a @@ -53,12 +71,12 @@ structure and the transmitted energy is measured on the right side of the struct The simulation is run twice. First, we run the simulation with the phononic structure and measure the transmitted energy: -Phononic superlattice cavity +Phononic superlattice cavity Then, we run the simulation without the phononic structure and measure the transmitted -energy; we use the simulation without the structure for the calibration: +energy. We use the simulation without the structure for the calibration: -Phononic superlattice cavity +Phononic superlattice cavity The transmission coefficient corresponds to the energy of the first simulation divided by the calibration energy. @@ -92,11 +110,13 @@ The elastic equations in the frequency domain then read as follows \nabla\cdot(\boldsymbol{\bar\sigma} \xi \boldsymbol{\Lambda})&=&-\omega^2\rho\xi\mathbf{\bar u}\\ \boldsymbol{\bar \sigma} &=&\mathbf{C}\boldsymbol{\bar\varepsilon}\\ \boldsymbol{\bar\varepsilon}&=&\frac{1}{2}[(\nabla\mathbf{\bar{u}}\boldsymbol{\Lambda}+\boldsymbol{\Lambda}^\mathrm{T}(\nabla\mathbf{\bar{u}})^\mathrm{T})]\\ -\xi &=&s_0\cdot s_1\cdot s_2\\ -\boldsymbol{\Lambda} &=& \operatorname{diag}(1/s_0,1/s_1,1/s_2) +\xi &=&\prod_i^\textrm{dim}s_i\\ +\boldsymbol{\Lambda} &=& \operatorname{diag}(1/s_0,1/s_1,1/s_2)\qquad\textrm{for 3D}\\ +\boldsymbol{\Lambda} &=& \operatorname{diag}(1/s_0,1/s_1)\qquad\textrm{for 2D} @f} where the coefficients $s_i = 1+is_i'(x,y,z)$ account for the absorption. -The imaginary par of $s_i$ is equal to zero outside of the PML. +There are 3 $s_i$ coefficients in 3D and 2 in 2D. +The imaginary par of $s_i$ is equal to zero outside the PML. The PMLs are reflectionless only for the exact wave equations. When the set of equations is discretized the PML is no longer reflectionless. The reflections can be made arbitrarily small as long as the @@ -114,13 +134,17 @@ These equations can be expanded into \varepsilon_{kl} =\frac{1}{2}\left(\frac{1}{s_k}\partial_k u_l + \frac{1}{s_l}\partial_l u_k\right) @f] -which can be written as +where summation over repeated indices (here $n$, as well as $k$ and $l$) is as always implied. +Note that the strain is no longer symmetric after applying the complex coordinate +stretching of the PML. +This set of equations can be written as @f[ -\omega^2\rho \xi u_m - \partial_n \left(\frac{\xi c_{mnkl}}{2s_n s_k} \partial_k u_l + \frac{\xi c_{mnkl}}{2s_n s_l} \partial_l u_k\right) = f_m @f] -Note that the stress tensor is not symmetric inside the PML ($s_j\neq 0$). +The same as the strain, the stress tensor is not symmetric inside the PML ($s_j\neq 0$). +Indeed the fields inside the PML are not physical. It is useful to introduce the tensors $\alpha_{mnkl}$ and $\beta_{mnkl}$. @f[ -\omega^2\rho \xi u_m - \partial_n \left(\alpha_{mnkl}\partial_k u_l @@ -132,16 +156,19 @@ We can multiply by $\varphi_m$ and integrate over the domain $\Omega$ and integr -\omega^2\int_\Omega\rho\xi\varphi_m u_m + \int_\Omega\partial_n\varphi_m \left(\frac{\xi c_{mnkl}}{2s_n s_k} \partial_k u_l + \frac{\xi c_{mnkl}}{2s_n s_l} \partial_l u_k\right) = \int_\Omega\varphi_m f_m @f} - -Then the linear system becomes +It is this set of equations we want to solve for a set of frequencies $\omega$ in order to compute the +transmission coefficient as function of frequency. +The linear system becomes @f{eqnarray*} --\omega^2\int_\Omega\rho \xi\varphi_m^i \varphi_m^j + \int_\Omega\partial_n\varphi_m^i \left(\frac{\xi c_{mnkl}}{2s_n s_k} \partial_k \varphi_l^j -+ \frac{\xi c_{mnkl}}{2s_n s_l} \partial_l \varphi_k^j\right) = A_{ij} +AU&=&F\\ +A_{ij} &=& -\omega^2\int_\Omega\rho \xi\varphi_m^i \varphi_m^j + \int_\Omega\partial_n\varphi_m^i \left(\frac{\xi c_{mnkl}}{2s_n s_k} \partial_k \varphi_l^j ++ \frac{\xi c_{mnkl}}{2s_n s_l} \partial_l \varphi_k^j\right)\\ +F_i &=& \int_\Omega\varphi_m^i f_m @f}

Simulation parameters

In this tutorial we use a python -[jupyter notebook](https://github.com/dangars/dealii/blob/phononic-cavity/examples/step-62/step-62.ipynb) +[jupyter notebook](https://github.com/dealii/dealii/blob/phononic-cavity/examples/step-62/step-62.ipynb) to set up the parameters and run the simulation. First we create a HDF5 file where we store the parameters and the results of the simulation. @@ -152,6 +179,11 @@ import numpy as np import h5py import matplotlib.pyplot as plt import subprocess +import scipy.constants as constants +import scipy.optimize + +# This considerably reduces the size of the svg data +plt.rcParams['svg.fonttype'] = 'none' h5_file = h5py.File('results.h5', 'w') data = h5_file.create_group('data') @@ -161,13 +193,15 @@ calibration = data.create_group('calibration') # Set the parameters for group in [displacement, calibration]: # Dimensions of the domain - group.attrs['dimension_x'] = 0.02 - group.attrs['dimension_y'] = 2e-5 + # The waveguide length is equal to dimension_x + group.attrs['dimension_x'] = 2e-5 + # The waveguide width is equal to dimension_y + group.attrs['dimension_y'] = 2e-8 # Position of the probe that we use to measure the flux - group.attrs['probe_pos_x'] = 0.008 - group.attrs['probe_pos_y'] = 0 - group.attrs['probe_width_y'] = 2e-05 + group.attrs['probe_pos_x'] = 8e-6 + group.attrs['probe_pos_y'] = 0 + group.attrs['probe_width_y'] = 2e-08 # Number of points in the probe group.attrs['nb_probe_points'] = 5 @@ -176,8 +210,8 @@ for group in [displacement, calibration]: group.attrs['grid_level'] = 1 # Cavity - group.attrs['cavity_resonance_frequency'] = 20000000.0 - group.attrs['nb_mirror_pairs'] = 30 + group.attrs['cavity_resonance_frequency'] = 20e9 + group.attrs['nb_mirror_pairs'] = 15 # Material group.attrs['poissons_ratio'] = 0.27 @@ -193,35 +227,35 @@ for group in [displacement, calibration]: group.attrs['mu']= (group.attrs['youngs_modulus'] / (2 * (1 + group.attrs['poissons_ratio']))) # Force - group.attrs['max_force_amplitude'] = 1e20 - group.attrs['force_sigma_x'] = 1 - group.attrs['force_sigma_y'] = 1 - group.attrs['max_force_width_x'] = 0.0003 - group.attrs['max_force_width_y'] = 0.001 - group.attrs['force_x_pos'] = -0.008 - group.attrs['force_y_pos'] = 0 + group.attrs['max_force_amplitude'] = 1e26 + group.attrs['force_sigma_x'] = 1e-7 + group.attrs['force_sigma_y'] = 1 + group.attrs['max_force_width_x'] = 3e-7 + group.attrs['max_force_width_y'] = 2e-8 + group.attrs['force_x_pos'] = -8e-6 + group.attrs['force_y_pos'] = 0 # PML - group.attrs['pml_x'] = True - group.attrs['pml_y'] = False - group.attrs['pml_width_x'] = 0.0018 - group.attrs['pml_width_y'] = 0.0005 - group.attrs['pml_coeff'] = 1.6 + group.attrs['pml_x'] = True + group.attrs['pml_y'] = False + group.attrs['pml_width_x'] = 1.8e-6 + group.attrs['pml_width_y'] = 5e-7 + group.attrs['pml_coeff'] = 1.6 group.attrs['pml_coeff_degree'] = 2 # Frequency sweep - group.attrs['center_frequency'] = 19990180.0 - group.attrs['frequency_range'] = 6000000.0 - group.attrs['start_frequency'] = group.attrs['center_frequency'] - group.attrs['frequency_range'] / 2 - group.attrs['stop_frequency'] = group.attrs['center_frequency'] + group.attrs['frequency_range'] / 2 - group.attrs['nb_frequency_points'] = 10 + group.attrs['center_frequency'] = 20e9 + group.attrs['frequency_range'] = 0.5e9 + group.attrs['start_frequency'] = group.attrs['center_frequency'] - group.attrs['frequency_range'] / 2 + group.attrs['stop_frequency'] = group.attrs['center_frequency'] + group.attrs['frequency_range'] / 2 + group.attrs['nb_frequency_points'] = 400 # Other parameters if group == displacement: group.attrs['simulation_name'] = 'phononic_cavity_displacement' else: group.attrs['simulation_name'] = 'phononic_cavity_calibration' - group.attrs['save_vtu_files'] = True + group.attrs['save_vtu_files'] = False h5_file.close() @endcode diff --git a/examples/step-62/doc/results.dox b/examples/step-62/doc/results.dox index 43fc5d9f2e..04cbc1af6c 100644 --- a/examples/step-62/doc/results.dox +++ b/examples/step-62/doc/results.dox @@ -1,7 +1,9 @@

Results

+

Resonance frequency and bandgap

+ The results are analyzed in the -[jupyter notebook](https://github.com/dangars/dealii/blob/phononic-cavity/examples/step-62/step-62.ipynb) +[jupyter notebook](https://github.com/dealii/dealii/blob/phononic-cavity/examples/step-62/step-62.ipynb) with the following code @code{.py} h5_file = h5py.File('results.h5', 'r') @@ -36,32 +38,97 @@ try: y_data_fit = resonance_f(x_data, freq_m, quality_factor, max_amplitude) fig = plt.figure() - plt.plot(frequency / 1e6, reflectivity, frequency / 1e6, y_data_fit) - plt.xlabel('frequency (MHz)') + plt.plot(frequency / 1e9, reflectivity, frequency / 1e9, y_data_fit) + plt.xlabel('frequency (GHz)') plt.ylabel('amplitude^2 (a.u.)') - plt.title('Transmission\n' + 'freq = ' + "%.7g" % (freq_guess / 1e6) + 'MHz Q = ' + "%.6g" % quality_factor) + plt.title('Transmission\n' + 'freq = ' + "%.7g" % (freq_guess / 1e9) + 'GHz Q = ' + "%.6g" % quality_factor) except: fig = plt.figure() - plt.plot(frequency / 1e6, reflectivity) - plt.xlabel('frequency (MHz)') + plt.plot(frequency / 1e9, reflectivity) + plt.xlabel('frequency (GHz)') plt.ylabel('amplitude^2 (a.u.)') plt.title('Transmission') fig = plt.figure() -plt.plot(frequency / 1e6, np.angle(reflection_coefficient)) -plt.xlabel('frequency (MHz)') +plt.plot(frequency / 1e9, np.angle(reflection_coefficient)) +plt.xlabel('frequency (GHz)') plt.ylabel('phase (rad)') -plt.title('Phase (reflection coefficient)\n') +plt.title('Phase (transmission coefficient)\n') plt.show() h5_file.close() @endcode -The micropillar cavity exhibits a mechanical resonance at 20MHz and a quality factor of 5091 +A phononic cavity is characterized by the +[resonance frequency](https://en.wikipedia.org/wiki/Resonance) and the +[the quality factor](https://en.wikipedia.org/wiki/Q_factor). +The quality factor is equal to the ratio between the stored energy in the resonator and the energy +dissipated energy per cycle, which is approximately equivalent to the ratio between the +resonance frequency and the +[full width at half maximum (FWHM)](https://en.wikipedia.org/wiki/Full_width_at_half_maximum). +The FWHM is equal to the bandwidth over which the power of vibration is greater than half the +power at the resonant frequency. +@f[ +Q = \frac{f_r}{\Delta f} = \frac{\omega_r}{\Delta \omega} = +2 \pi \times \frac{\text{energy stored}}{\text{energy dissipated per cycle}} +@f] + +The square of the amplitude of the mechanical resonance $a^2$ as a function of the frequency +has a gaussian shape +@f[ +a^2 = a_\textrm{max}^2\frac{\omega^2\Gamma^2}{(\omega_r^2-\omega^2)^2+\Gamma^2\omega^2} +@f] +where $f_r = \frac{\omega_r}{2\pi}$ is the resonance frequency and $\Gamma=\frac{\omega_r}{Q}$ is the dissipation rate. +We used the previous equation in the jupyter notebook to fit the mechanical resonance. + +Given the values we have chosen for the parameters, one could estimate the resonance frequency +analytically. Indeed, this is then confirmed by what we get in this program: +the phononic superlattice cavity exhibits a mechanical resonance at 20GHz and a quality factor of 5046. +The following images show the transmission amplitude and phase as a function of frequency in the +vicinity of the resonance frequency: + +Phononic superlattice cavity +Phononic superlattice cavity + +The images above suggest that the periodic structure has its intended effect: It really only lets waves of a very +specific frequency pass through, whereas all other waves are reflected. This is of course precisely what one builds +these sorts of devices for. +But it is not quite this easy. In practice, there is really only a "band gap", i.e., the device blocks waves other than +the desired one at 20GHz only within a certain frequency range. Indeed, to find out how large this "gap" is within +which waves are blocked, we can extend the frequency range to 16 GHz through the appropriate parameters in the +input file. We then obtain the following image: + +Phononic superlattice cavity + +What this image suggests is that in the range of around 18 to around 22 GHz, really only the waves with a frequency +of 20 GHz are allowed to pass through, but beyond this range, there are plenty of other frequencies that can pass +through the device. + +

Mode profile

+ +We can inspect the mode profile with Paraview or Visit. +As we have discussed, at resonance all the mechanical +energy is transmitted and the amplitude of motion is amplified inside the cavity. +It can be observed that the PMLs are quite effective to truncate the solution. +The following image shows the mode profile at resonance: + +Phononic superlattice cavity + +On the other hand, out of resonance all the mechanical energy is +reflected. The following image shows the profile at 19.75 GHz. +Note the interference between the force pulse and the reflected wave +at the position $x=-8\mu\textrm{m}$. -Phononic superlattice cavity -Phononic superlattice cavity +Phononic superlattice cavity -To obtain the phononic bandgap around the mechanical resonance, the parameter frequency range can be set to 16 MHz. +

Experimental applications

-Phononic superlattice cavity +Phononic superlattice cavities find application in +[quantum optomechanics](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.86.1391). +Here we have presented the simulation of a 2D superlattice cavity, +but this code can be used as well to simulate "real world" 3D devices such as +[micropillar superlattice cavities](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.060101), +which are promising candidates to study macroscopic quantum phenomena. +The 20GHz mode of a micropillar superlattice cavity is essentially a mechanical harmonic oscillator that is very well isolated +from the environment. If the device is cooled down to 20mK in a dilution fridge, the mode would then become a +macroscopic quantum harmonic oscillator. diff --git a/examples/step-62/doc/tooltip b/examples/step-62/doc/tooltip index b57f7fe2aa..c606b58750 100644 --- a/examples/step-62/doc/tooltip +++ b/examples/step-62/doc/tooltip @@ -1 +1 @@ -Systems of PDE. Elasticity. Tensors. +Elastic equation in the frequency domain. Calculate the transmission and resonance frequency of a phononic structure. diff --git a/examples/step-62/step-62.cc b/examples/step-62/step-62.cc index e67a041d1a..78f5d13a7a 100644 --- a/examples/step-62/step-62.cc +++ b/examples/step-62/step-62.cc @@ -1,6 +1,6 @@ /* --------------------------------------------------------------------- * - * Copyright (C) 2000 - 2018 by the deal.II authors + * Copyright (C) 2018 - 2019 by the deal.II authors * * This file is part of the deal.II library. * @@ -53,7 +53,7 @@ #include #include -// The following header provides the Tensor class that we use represent the +// The following header provides the Tensor class that we use to represent the // material properties. #include @@ -62,12 +62,11 @@ #include // This header is required for the function VectorTools::point_value that we use -// to read the result of the simulation. - +// to evaluate the result of the simulation. #include // We need this header for the function GridTools::find_active_cell_around_point -// that we use in the function ElasticWave::store_frequency_step_data +// that we use in the function `ElasticWave::store_frequency_step_data()` #include namespace step62 @@ -89,12 +88,11 @@ namespace step62 const unsigned int component) const override; private: - // `data` is the HDF5::Group in which all the simulation results will be - // stored. Note that this variable points to the same HDF5::Group of - // `RightHandSide::data`, `PML::data` and `Parameters::data`. When a - // HDF5::Group is copied, it will point to the same HDF5 Group; this is - // achieved with the protected std::shared_ptr - // HDF5::Group::hdf5_reference. + // The variable `data` is the HDF5::Group in which all the simulation + // results will be stored. Note that the variables `RightHandSide::data`, + // `PML::data`, `Rho::data` and `Parameters::data` point to the same group + // of the HDF5 file. When a HDF5::Group is copied, it will point to the same + // group of the HDF5 file. HDF5::Group data; // The simulation parameters are stored in `data` as HDF5 attributes. The @@ -190,7 +188,7 @@ namespace step62 // notebook, stored in `data` as HDF5 attributes and then read by the // constructor. const std::string simulation_name; - bool save_vtu_files; + const bool save_vtu_files; const double start_frequency; const double stop_frequency; const unsigned int nb_frequency_points; @@ -200,8 +198,8 @@ namespace step62 const double dimension_y; const unsigned int nb_probe_points; const unsigned int grid_level; - Point probe_start_point; - Point probe_stop_point; + const Point probe_start_point; + const Point probe_stop_point; const RightHandSide right_hand_side; const PML pml; const Rho rho; @@ -212,17 +210,19 @@ namespace step62 - // @sect4{`PointHistory` class} + // @sect4{The `QuadratureCache` class} // The calculation of the mass and stiffness matrices is very expensive. These // matrices are the same for all the frequency steps. The right hand side // vector is also the same for all the frequency steps. We use this class to - // store these objects and re-use them at each frequency step. The - // `PointHistory` class has already been used in step-18. + // store these objects and re-use them at each frequency step. Note that here + // we don't store the assembled mass and stiffness matrices and right hand + // sides, but instead the data for a single cell. `QuadratureCache` class is + // very similar to the `PointHistory` class that has been used in step-18. template - class PointHistory + class QuadratureCache { public: - PointHistory(unsigned int dofs_per_cell); + QuadratureCache(const unsigned int dofs_per_cell); private: unsigned int dofs_per_cell; @@ -235,19 +235,19 @@ namespace step62 FullMatrix> mass_coefficient; FullMatrix> stiffness_coefficient; std::vector> right_hand_side; - std::complex JxW; + double JxW; }; // @sect4{The `get_stiffness_tensor()` function} - // This class returns the stiffness tensor of the material. For the sake of + // This function returns the stiffness tensor of the material. For the sake of // simplicity we consider the stiffness to be isotropic and homogeneous; only - // the density $\rho$ depends on the position. As we have previously done in - // step-8. The stiffness coefficients $c_{ijkl}$ can be expressed in function - // of the two coefficients $\lambda$ and $\mu$. The coefficient tensor reduces - // to + // the density $\rho$ depends on the position. As we have previously shown in + // step-8, if the stiffness is isotropic and homogeneous, the stiffness + // coefficients $c_{ijkl}$ can be expressed as a function of the two + // coefficients $\lambda$ and $\mu$. The coefficient tensor reduces to // @f[ // c_{ijkl} // = @@ -277,7 +277,7 @@ namespace step62 // Next let's declare the main class of this program. Its structure is very // similar to the step-40 tutorial program. The main differences are: // - The sweep over the frequency values. - // - We save the stiffness and mass matrices in `quadrature_point_history` and + // - We save the stiffness and mass matrices in `quadrature_cache` and // use them for each frequency step. // - We store the measured energy by the probe for each frequency step in the // HDF5 file. @@ -285,21 +285,21 @@ namespace step62 class ElasticWave { public: - ElasticWave(Parameters parameters_); - ~ElasticWave(); + ElasticWave(const Parameters ¶meters); void run(); private: void setup_system(); - void assemble_system(double omega, bool calculate_quadrature_data); + void assemble_system(const double omega, + const bool calculate_quadrature_data); void solve(); - void set_position_vector(); - void store_frequency_step_data(unsigned int frequency_idx); + void initialize_probe_positions_vector(); + void store_frequency_step_data(const unsigned int frequency_idx); void output_results(); - // This is called before every time step to set up a pristine state for the - // history variables. - void setup_quadrature_point_history(); + // This is called before every frequency step to set up a pristine state + // for the cache variables. + void setup_quadrature_cache(); // This function loops over the frequency vector and runs the simulation for // each frequency step. @@ -312,15 +312,14 @@ namespace step62 parallel::distributed::Triangulation triangulation; - QGauss quadrature_formula; - const unsigned int n_q_points; + QGauss quadrature_formula; - // We store the mass and stiffness matrices in this vector. - std::vector> quadrature_point_history; + // We store the mass and stiffness matrices for each cell this vector. + std::vector> quadrature_cache; - DoFHandler dof_handler; - FESystem fe; + FESystem fe; + DoFHandler dof_handler; IndexSet locally_owned_dofs; IndexSet locally_relevant_dofs; @@ -333,16 +332,16 @@ namespace step62 // This vector contains the range of frequencies that we are going to - // simulate + // simulate. std::vector frequency; // This vector contains the coordinates $(x,y)$ of the points of the // measurement probe. - FullMatrix position; + FullMatrix probe_positions; - // HDF5 datasets to store the frequency and position vectors. + // HDF5 datasets to store the frequency and `probe_positions` vectors. HDF5::DataSet frequency_dataset; - HDF5::DataSet position_dataset; + HDF5::DataSet probe_positions_dataset; // HDF5 dataset that stores the values of the energy measured by the proble. HDF5::DataSet displacement; @@ -354,12 +353,12 @@ namespace step62 - // @sect3{Implementation of the auxiliary classes and functions} + // @sect3{Implementation of the auxiliary classes} - // @sect4{`RightHandSide` class} + // @sect4{The `RightHandSide` class} // The constructor reads all the parameters from the HDF5::Group `data` using - // the HDF5::Group::get_attribute function. + // the HDF5::Group::get_attribute() function. template RightHandSide::RightHandSide(HDF5::Group &data) : Function(dim) @@ -373,8 +372,8 @@ namespace step62 data.get_attribute("force_y_pos"))) {} - // This function defines the spacial shape of the force vector pulse which - // takes the form of a gaussian function + // This function defines the spatial shape of the force vector pulse which + // takes the form of a Gaussian function // @f{align*} // F_x &= // \left\{ @@ -389,7 +388,7 @@ namespace step62 // \right.\\ // F_y &= 0 // @f} - // where a is the maximum amplitude that takes the force and $\sigma_x$ and + // where $a$ is the maximum amplitude that takes the force and $\sigma_x$ and // $\sigma_y$ are the standard deviations for the $x$ and $y$ components. Note // that the pulse has been cropped to $x_\textrm{min} PML::PML(HDF5::Group &data) : Function>(dim) @@ -485,15 +484,15 @@ namespace step62 - // @sect4{`Rho` class} + // @sect4{The `Rho` class} - // This class is used to define the mass density. As we have explained, before + // This class is used to define the mass density. As we have explaine before, // a phononic superlattice cavity is formed by two - //[Distributed Reflector](https://en.wikipedia.org/wiki/Band_gap), + // [Distributed Reflector](https://en.wikipedia.org/wiki/Band_gap), // mirrors and a $\lambda/2$ cavity where $\lambda$ is the acoustic - // wavelength. Acoustic DBRs are periodic structures where a set of bilayer + // wavelength. Acoustic DBRs are periodic structures where a set of bilayer // stacks with contrasting physical properties (sound velocity index) is - // repeated $N$ times. The change of in the velocity will be obtained by + // repeated $N$ times. The change of in the wave velocity is generated by // alternating layers with different density. template Rho::Rho(HDF5::Group &data) @@ -528,12 +527,12 @@ namespace step62 // @f] // where $K_e$ is the effective elastic constant and $\rho$ the density. // Here we consider the case in which the waveguide width is much smaller - // than the wavelength. In this case it can be shown that for a two + // than the wavelength. In this case it can be shown that for the two // dimensional case // @f[ // K_e = 4\mu\frac{\lambda +\mu}{\lamda+2\mu} // @f] - // and for a three dimensional case $K_e$ is equal to the Young's modulus. + // and for the three dimensional case $K_e$ is equal to the Young's modulus. // @f[ // K_e = 4\mu\frac{\lambda +\mu}{\lamda+2\mu} // @f] @@ -560,23 +559,24 @@ namespace step62 material_b_speed_of_sound / cavity_resonance_frequency; // The density $\rho$ takes the following form - //Phononic superlattice cavity // where the brown color represents material_a and the green color // represents material_b. for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++) { - double layer_transition_center = + const double layer_transition_center = material_a_wavelength / 2 + idx * (material_b_wavelength / 4 + material_a_wavelength / 4); if (std::abs(p[0]) >= (layer_transition_center - average_rho_width / 2) && std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2)) { - double coefficient = (std::abs(p[0]) - (layer_transition_center - - average_rho_width / 2)) / - average_rho_width; + const double coefficient = + (std::abs(p[0]) - + (layer_transition_center - average_rho_width / 2)) / + average_rho_width; return (1 - coefficient) * material_a_rho + coefficient * material_b_rho; } @@ -588,7 +588,7 @@ namespace step62 // which improves the precision of the simulation. for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++) { - double layer_transition_center = + const double layer_transition_center = material_a_wavelength / 2 + idx * (material_b_wavelength / 4 + material_a_wavelength / 4) + material_b_wavelength / 4; @@ -596,9 +596,10 @@ namespace step62 (layer_transition_center - average_rho_width / 2) && std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2)) { - double coefficient = (std::abs(p[0]) - (layer_transition_center - - average_rho_width / 2)) / - average_rho_width; + const double coefficient = + (std::abs(p[0]) - + (layer_transition_center - average_rho_width / 2)) / + average_rho_width; return (1 - coefficient) * material_b_rho + coefficient * material_a_rho; } @@ -613,11 +614,11 @@ namespace step62 // the material_a layers for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++) { - double layer_center = + const double layer_center = material_a_wavelength / 2 + idx * (material_b_wavelength / 4 + material_a_wavelength / 4) + material_b_wavelength / 4 + material_a_wavelength / 8; - double layer_width = material_a_wavelength / 4; + const double layer_width = material_a_wavelength / 4; if (std::abs(p[0]) >= (layer_center - layer_width / 2) && std::abs(p[0]) <= (layer_center + layer_width / 2)) { @@ -628,11 +629,11 @@ namespace step62 // the material_b layers for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++) { - double layer_center = + const double layer_center = material_a_wavelength / 2 + idx * (material_b_wavelength / 4 + material_a_wavelength / 4) + material_b_wavelength / 8; - double layer_width = material_b_wavelength / 4; + const double layer_width = material_b_wavelength / 4; if (std::abs(p[0]) >= (layer_center - layer_width / 2) && std::abs(p[0]) <= (layer_center + layer_width / 2)) { @@ -646,10 +647,10 @@ namespace step62 - // @sect4{`Parameters` class} + // @sect4{The `Parameters` class} // The constructor reads all the parameters from the HDF5::Group `data` using - // the HDF5::Group::get_attribute function. + // the HDF5::Group::get_attribute() function. template Parameters::Parameters(HDF5::Group &data) : data(data) @@ -664,28 +665,27 @@ namespace step62 , dimension_y(data.get_attribute("dimension_y")) , nb_probe_points(data.get_attribute("nb_probe_points")) , grid_level(data.get_attribute("grid_level")) + , probe_start_point( + Point(data.get_attribute("probe_pos_x"), + data.get_attribute("probe_pos_y") - + data.get_attribute("probe_width_y") / 2)) + , probe_stop_point( + Point(data.get_attribute("probe_pos_x"), + data.get_attribute("probe_pos_y") + + data.get_attribute("probe_width_y") / 2)) , right_hand_side(data) , pml(data) , rho(data) - { - probe_start_point = - Point(data.get_attribute("probe_pos_x"), - data.get_attribute("probe_pos_y") - - data.get_attribute("probe_width_y") / 2); - probe_stop_point = - Point(data.get_attribute("probe_pos_x"), - data.get_attribute("probe_pos_y") + - data.get_attribute("probe_width_y") / 2); - } + {} - // @sect4{`PointHistory` class} + // @sect4{The `QuadratureCache` class} // We need to reserve enough space for the mass and stiffness matrices and the // right hand side vector. template - PointHistory::PointHistory(unsigned int dofs_per_cell) + QuadratureCache::QuadratureCache(const unsigned int dofs_per_cell) : dofs_per_cell(dofs_per_cell) , mass_coefficient(dofs_per_cell, dofs_per_cell) , stiffness_coefficient(dofs_per_cell, dofs_per_cell) @@ -696,31 +696,30 @@ namespace step62 // @sect3{Implementation of the `ElasticWave` class} - // @sect4{Constructors and destructors} + // @sect4{Constructor} // This is very similar to the constructor of step-40. In addition we create // the HDF5 datasets `frequency_dataset`, `position_dataset` and - // `displacement`. Note the use of the `template` for the creation of the HDF5 - // datasets. It is a C++ requirement to use the `template` keyword in order to - // treat `create_dataset` as a dependent template name. + // `displacement`. Note the use of the `template` keyword for the creation of + // the HDF5 datasets. It is a C++ requirement to use the `template` keyword in + // order to treat `create_dataset` as a dependent template name. template - ElasticWave::ElasticWave(Parameters parameters_) - : parameters(parameters_) + ElasticWave::ElasticWave(const Parameters ¶meters) + : parameters(parameters) , mpi_communicator(MPI_COMM_WORLD) , triangulation(mpi_communicator, typename Triangulation::MeshSmoothing( Triangulation::smoothing_on_refinement | Triangulation::smoothing_on_coarsening)) , quadrature_formula(2) - , n_q_points(quadrature_formula.size()) - , dof_handler(triangulation) , fe(FE_Q(1), dim) + , dof_handler(triangulation) , frequency(parameters.nb_frequency_points) - , position(parameters.nb_probe_points, dim) + , probe_positions(parameters.nb_probe_points, dim) , frequency_dataset(parameters.data.template create_dataset( "frequency", std::vector{parameters.nb_frequency_points})) - , position_dataset(parameters.data.template create_dataset( + , probe_positions_dataset(parameters.data.template create_dataset( "position", std::vector{parameters.nb_probe_points, dim})) , displacement( @@ -738,14 +737,6 @@ namespace step62 - template - ElasticWave::~ElasticWave() - { - dof_handler.clear(); - } - - - // @sect4{ElasticWave::setup_system} // There is nothing new in this function, the only difference with step-40 is @@ -792,14 +783,14 @@ namespace step62 // @sect4{ElasticWave::assemble_system} - // This very similar to step-40. Although there are notable differences. We - // assembly the system for each frequency/omega step. In the first step we set - // `calculate_quadrature_data = True` and we calculate the mass and stiffness - // matrices and the right hand side vector. In the subsequent steps we will - // use that data to accelerate the calculation. + // This function is also very similar to step-40, though there are notable + // differences. We assemble the system for each frequency/omega step. In the + // first step we set `calculate_quadrature_data = True` and we calculate the + // mass and stiffness matrices and the right hand side vector. In the + // subsequent steps we will use that data to accelerate the calculation. template - void ElasticWave::assemble_system(double omega, - bool calculate_quadrature_data) + void ElasticWave::assemble_system(const double omega, + const bool calculate_quadrature_data) { TimerOutput::Scope t(computing_timer, "assembly"); @@ -808,6 +799,7 @@ namespace step62 update_values | update_gradients | update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); FullMatrix> cell_matrix(dofs_per_cell, dofs_per_cell); Vector> cell_rhs(dofs_per_cell); @@ -820,7 +812,7 @@ namespace step62 std::vector>> pml_values( n_q_points, Vector>(dim)); - // We calculate the stiffness tensor for the $\lambda$ and $\mu$ that has + // We calculate the stiffness tensor for the $\lambda$ and $\mu$ that have // been defined in the jupyter notebook. Note that contrary to $\rho$ the // stiffness is constant among for the whole domain. const SymmetricTensor<4, dim> stiffness_tensor = @@ -853,23 +845,21 @@ namespace step62 } // We have done this in step-18. Get a pointer to the quadrature - // point history data local to the present cell, and, as a defensive + // cache data local to the present cell, and, as a defensive // measure, make sure that this pointer is within the bounds of the // global array: - PointHistory *local_quadrature_points_data = - reinterpret_cast *>(cell->user_pointer()); - Assert(local_quadrature_points_data >= - &quadrature_point_history.front(), + QuadratureCache *local_quadrature_points_data = + reinterpret_cast *>(cell->user_pointer()); + Assert(local_quadrature_points_data >= &quadrature_cache.front(), ExcInternalError()); - Assert(local_quadrature_points_data < - &quadrature_point_history.back(), + Assert(local_quadrature_points_data <= &quadrature_cache.back(), ExcInternalError()); for (unsigned int q = 0; q < n_q_points; ++q) { // The quadrature_data variable is used to store the mass and // stiffness matrices, the right hand side vector and the value // of `JxW`. - PointHistory &quadrature_data = + QuadratureCache &quadrature_data = local_quadrature_points_data[q]; // Below we declare the force vector and the parameters of the @@ -939,9 +929,9 @@ namespace step62 // matrix is not symmetric because // of the PMLs. We use the gradient // function (see the - // [documentation](https://www.dealii.org/current/doxygen/deal.II/group__vector__valued.html) + // [documentation](https://www.dealii.org/current/doxygen/deal.II/group__vector__valued.html)) // which is a - // Tensor@<2,dim@>, + // Tensor@<2,dim@>. // The matrix $G_{ij}$ // consists of entries // @f[ @@ -1035,13 +1025,12 @@ namespace step62 locally_relevant_solution = completely_distributed_solution; } - // @sect4{ElasticWave::set_position_vector} + // @sect4{ElasticWave::initialize_position_vector} // We use this function to calculate the values of the position vector. template - void ElasticWave::set_position_vector() + void ElasticWave::initialize_probe_positions_vector() { - Point p; for (unsigned int position_idx = 0; position_idx < parameters.nb_probe_points; ++position_idx) @@ -1049,14 +1038,15 @@ namespace step62 // Because of the way the operator + and - are overloaded. To substract // two points, the following has to be done: // `Point_b + (-Point_a)` - p = (position_idx / ((double)(parameters.nb_probe_points - 1))) * - (parameters.probe_stop_point + (-parameters.probe_start_point)) + - parameters.probe_start_point; - position[position_idx][0] = p[0]; - position[position_idx][1] = p[1]; + const Point p = + (position_idx / ((double)(parameters.nb_probe_points - 1))) * + (parameters.probe_stop_point + (-parameters.probe_start_point)) + + parameters.probe_start_point; + probe_positions[position_idx][0] = p[0]; + probe_positions[position_idx][1] = p[1]; if (dim == 3) { - position[position_idx][2] = p[2]; + probe_positions[position_idx][2] = p[2]; } } } @@ -1065,13 +1055,14 @@ namespace step62 // This function stores in the HDF5 file the measured energy by the probe. template - void ElasticWave::store_frequency_step_data(unsigned int frequency_idx) + void + ElasticWave::store_frequency_step_data(const unsigned int frequency_idx) { TimerOutput::Scope t(computing_timer, "store_frequency_step_data"); // We store the displacement in the $x$ direction; the displacement in the // $y$ direction is negligible. - const int probe_displacement_component = 0; + const unsigned int probe_displacement_component = 0; // The vector coordinates contains the coordinates in the HDF5 file of the // points of the probe that are located in locally owned cells. The vector @@ -1085,7 +1076,7 @@ namespace step62 Point point; for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) { - point(dim_idx) = position[position_idx][dim_idx]; + point[dim_idx] = probe_positions[position_idx][dim_idx]; } bool point_in_locally_owned_cell; { @@ -1131,7 +1122,7 @@ namespace step62 displacement.write_none>(); } - // If the variable of the jupyter notbook `save_vtu_files = True` then all + // If the variable `save_vtu_files` in the input file equals `True` then all // the data will be saved as vtu. The procedure to write `vtu` files has // been described in step-40. if (parameters.save_vtu_files) @@ -1186,10 +1177,10 @@ namespace step62 { for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) { - force[dim_idx](cell->active_cell_index()) = - parameters.right_hand_side.value(cell->center(), dim_idx); - pml[dim_idx](cell->active_cell_index()) = -1e+20; + force[dim_idx](cell->active_cell_index()) = -1e+20; + pml[dim_idx](cell->active_cell_index()) = -1e+20; } + rho(cell->active_cell_index()) = -1e+20; } } @@ -1204,9 +1195,8 @@ namespace step62 data_out.build_patches(); - unsigned int nb_number_positions; - std::stringstream frequency_idx_stream; - nb_number_positions = + std::stringstream frequency_idx_stream; + const unsigned int nb_number_positions = ((unsigned int)std::log10(parameters.nb_frequency_points)) + 1; frequency_idx_stream << std::setw(nb_number_positions) << std::setfill('0') << frequency_idx; @@ -1231,83 +1221,75 @@ namespace step62 if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) { frequency_dataset.write(frequency); - position_dataset.write(position); + probe_positions_dataset.write(probe_positions); } else { frequency_dataset.write_none(); - position_dataset.write_none(); + probe_positions_dataset.write_none(); } } - // @sect4{ElasticWave::setup_quadrature_point_history} + // @sect4{ElasticWave::setup_quadrature_cache} // We use this function at the beginning of our computations to set up initial - // values of the history variables. This function has been described in - // step-18. There are no differences with the function of step-18. + // values of the cache variables. This function has been described in step-18. + // There are no differences with the function of step-18. template - void ElasticWave::setup_quadrature_point_history() + void ElasticWave::setup_quadrature_cache() { - unsigned int our_cells = 0; - for (typename Triangulation::active_cell_iterator cell = - triangulation.begin_active(); - cell != triangulation.end(); - ++cell) - if (cell->is_locally_owned()) - ++our_cells; - triangulation.clear_user_data(); { - std::vector> tmp; - tmp.swap(quadrature_point_history); + std::vector> tmp; + quadrature_cache.swap(tmp); } - quadrature_point_history.resize(our_cells * quadrature_formula.size(), - PointHistory(fe.dofs_per_cell)); - unsigned int history_index = 0; + quadrature_cache.resize(triangulation.n_locally_owned_active_cells() * + quadrature_formula.size(), + QuadratureCache(fe.dofs_per_cell)); + unsigned int cache_index = 0; for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) if (cell->is_locally_owned()) { - cell->set_user_pointer(&quadrature_point_history[history_index]); - history_index += quadrature_formula.size(); + cell->set_user_pointer(&quadrature_cache[cache_index]); + cache_index += quadrature_formula.size(); } - Assert(history_index == quadrature_point_history.size(), - ExcInternalError()); + Assert(cache_index == quadrature_cache.size(), ExcInternalError()); } // @sect4{ElasticWave::frequency_sweep} - template // For clarity we divide the function `run` of step-40 into the functions // `run` and `frequency_sweep`. In the function `frequency_sweep` we place the // iteration over the frequency vector. + template void ElasticWave::frequency_sweep() { for (unsigned int frequency_idx = 0; frequency_idx < parameters.nb_frequency_points; ++frequency_idx) { - std::cout << parameters.simulation_name + " frequency idx: " - << frequency_idx << '/' << parameters.nb_frequency_points - 1 - << std::endl; + pcout << parameters.simulation_name + " frequency idx: " + << frequency_idx << '/' << parameters.nb_frequency_points - 1 + << std::endl; setup_system(); if (frequency_idx == 0) { - std::cout << " Number of active cells : " - << triangulation.n_active_cells() << std::endl; - std::cout << " Number of degrees of freedom : " - << dof_handler.n_dofs() << std::endl; + pcout << " Number of active cells : " + << triangulation.n_active_cells() << std::endl; + pcout << " Number of degrees of freedom : " + << dof_handler.n_dofs() << std::endl; } if (frequency_idx == 0) @@ -1320,12 +1302,13 @@ namespace step62 } // We calculate the frequency and omega values for this particular step. - double current_loop_frequency = + const double current_loop_frequency = (parameters.start_frequency + frequency_idx * (parameters.stop_frequency - parameters.start_frequency) / (parameters.nb_frequency_points - 1)); - double current_loop_omega = 2 * numbers::PI * current_loop_frequency; + const double current_loop_omega = + 2 * numbers::PI * current_loop_frequency; // In the first frequency step we calculate the mass and stiffness // matrices and the right hand side. In the subsequent frequency steps @@ -1353,9 +1336,9 @@ namespace step62 void ElasticWave::run() { #ifdef DEBUG - std::cout << "Debug mode" << std::endl; + pcout << "Debug mode" << std::endl; #else - std::cout << "Release mode" << std::endl; + pcout << "Release mode" << std::endl; #endif { @@ -1388,9 +1371,9 @@ namespace step62 triangulation.refine_global(parameters.grid_level); - setup_quadrature_point_history(); + setup_quadrature_cache(); - set_position_vector(); + initialize_probe_positions_vector(); frequency_sweep(); @@ -1429,7 +1412,7 @@ int main(int argc, char *argv[]) } { - // Calibration simulation. The parameters are read from the displacement + // Calibration simulation. The parameters are read from the calibration // HDF5 group and the results are saved in the same HDF5 group. auto calibration = data.open_group("calibration"); step62::Parameters parameters(calibration); diff --git a/examples/step-62/step-62.ipynb b/examples/step-62/step-62.ipynb index 0e318a3d59..2639aeb0d2 100644 --- a/examples/step-62/step-62.ipynb +++ b/examples/step-62/step-62.ipynb @@ -42,13 +42,15 @@ "# 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", + " # The waveguide length is equal to dimension_x\n", + " group.attrs['dimension_x'] = 2e-5\n", + " # The waveguide width is equal to dimension_y\n", + " group.attrs['dimension_y'] = 2e-8\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", + " group.attrs['probe_pos_x'] = 8e-6\n", + " group.attrs['probe_pos_y'] = 0\n", + " group.attrs['probe_width_y'] = 2e-08\n", " \n", " # Number of points in the probe\n", " group.attrs['nb_probe_points'] = 5\n", @@ -57,8 +59,8 @@ " group.attrs['grid_level'] = 1\n", "\n", " # Cavity\n", - " group.attrs['cavity_resonance_frequency'] = 20e6\n", - " group.attrs['nb_mirror_pairs'] = 15\n", + " group.attrs['cavity_resonance_frequency'] = 20e9\n", + " group.attrs['nb_mirror_pairs'] = 15\n", "\n", " # Material\n", " group.attrs['poissons_ratio'] = 0.27\n", @@ -74,27 +76,27 @@ " 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", + " group.attrs['max_force_amplitude'] = 1e26\n", + " group.attrs['force_sigma_x'] = 1e-7\n", + " group.attrs['force_sigma_y'] = 1\n", + " group.attrs['max_force_width_x'] = 3e-7\n", + " group.attrs['max_force_width_y'] = 2e-8\n", + " group.attrs['force_x_pos'] = -8e-6\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_x'] = True\n", + " group.attrs['pml_y'] = False\n", + " group.attrs['pml_width_x'] = 1.8e-6\n", + " group.attrs['pml_width_y'] = 5e-7\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['center_frequency'] = 20e9\n", + " group.attrs['frequency_range'] = 0.5e9\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", @@ -188,10 +190,10 @@ " \n", " \n", + "\" id=\"mb3d45af095\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -201,7 +203,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -211,7 +213,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -221,7 +223,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -231,7 +233,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -239,7 +241,7 @@ " \n", " \n", " \n", - " frequency (MHz)\n", + " frequency (GHz)\n", " \n", " \n", " \n", @@ -248,64 +250,64 @@ " \n", " \n", + "\" id=\"m626a42a4d9\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 0.0\n", + " 0.0\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 0.2\n", + " 0.2\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 0.4\n", + " 0.4\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 0.6\n", + " 0.6\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 0.8\n", + " 0.8\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 1.0\n", + " 1.0\n", " \n", " \n", " \n", @@ -313,91 +315,93 @@ " \n", " \n", " \n", - " \n", " \n", " \n", - " \n", " \n", " \n", @@ -422,19 +426,19 @@ " \n", " \n", " Transmission\n", - " freq = 20.00815MHz Q = 5091.3\n", + " freq = 20.00815GHz Q = 5046.91\n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -477,10 +481,10 @@ " \n", " \n", + "\" id=\"m71cbb29084\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -490,7 +494,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -500,7 +504,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -510,7 +514,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -520,7 +524,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -528,7 +532,7 @@ " \n", " \n", " \n", - " frequency (MHz)\n", + " frequency (GHz)\n", " \n", " \n", " \n", @@ -537,64 +541,64 @@ " \n", " \n", + "\" id=\"m02a51f3213\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " −1.0\n", + " −1.0\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " −0.5\n", + " −0.5\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 0.0\n", + " 0.0\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 0.5\n", + " 0.5\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 1.0\n", + " 1.0\n", " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", - " 1.5\n", + " 1.5\n", " \n", " \n", " \n", @@ -602,58 +606,58 @@ " \n", " \n", " \n", - " \n", " \n", " \n", @@ -683,14 +687,14 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -730,35 +734,26 @@ " 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.plot(frequency / 1e9, reflectivity, frequency / 1e9, y_data_fit)\n", + " plt.xlabel('frequency (GHz)')\n", " plt.ylabel('amplitude^2 (a.u.)')\n", - " plt.title('Transmission\\n' + 'freq = ' + \"%.7g\" % (freq_guess / 1e6) + 'MHz Q = ' + \"%.6g\" % quality_factor)\n", + " plt.title('Transmission\\n' + 'freq = ' + \"%.7g\" % (freq_guess / 1e9) + 'GHz Q = ' + \"%.6g\" % quality_factor)\n", "except:\n", " fig = plt.figure()\n", - " plt.plot(frequency / 1e6, reflectivity)\n", - " plt.xlabel('frequency (MHz)')\n", + " plt.plot(frequency / 1e9, reflectivity)\n", + " plt.xlabel('frequency (GHz)')\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.plot(frequency / 1e9, np.angle(reflection_coefficient))\n", + "plt.xlabel('frequency (GHz)')\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": {