"[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"
]
}
##
-# CMake script for the step-8 tutorial program:
+# CMake script for the step-62 tutorial program:
##
# Set the name of the project and target:
href="../../readme.html" target="body">README</a> file.
<h1>Introduction</h1>
-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$.
+<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.01.svg" height="200" />
-<img alt="Phononic superlattice cavity" src="https://raw.githubusercontent.com/dangars/dealii/phononic-cavity/examples/step-62/doc/step-62.01.svg?sanitize=true" height="200" />
+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
The simulation is run twice. First, we run the simulation with the phononic
structure and measure the transmitted energy:
-<img alt="Phononic superlattice cavity" src="https://raw.githubusercontent.com/dangars/dealii/phononic-cavity/examples/step-62/doc/step-62.02.svg?sanitize=true" height="200" />
+<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.02.svg" height="200" />
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:
-<img alt="Phononic superlattice cavity" src="https://raw.githubusercontent.com/dangars/dealii/phononic-cavity/examples/step-62/doc/step-62.03.svg?sanitize=true" height="200" />
+<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.03.svg" height="200" />
The transmission coefficient corresponds to the energy of the first simulation
divided by the calibration energy.
\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
\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
-\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}
<h3>Simulation parameters</h3>
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.
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')
# 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
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
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
<h1>Results</h1>
+<h3>Resonance frequency and bandgap</h3>
+
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')
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:
+
+<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.05.png" height="400" />
+<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.06.png" height="400" />
+
+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:
+
+<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.07.png" height="400" />
+
+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.
+
+<h3>Mode profile</h3>
+
+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:
+
+<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.08.png" height="400" />
+
+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}$.
-<img alt="Phononic superlattice cavity" src="https://raw.githubusercontent.com/dangars/dealii/phononic-cavity/examples/step-62/doc/step-62.05.svg?sanitize=true" height="400" />
-<img alt="Phononic superlattice cavity" src="https://raw.githubusercontent.com/dangars/dealii/phononic-cavity/examples/step-62/doc/step-62.06.svg?sanitize=true" height="400" />
+<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.09.png" height="400" />
-To obtain the phononic bandgap around the mechanical resonance, the parameter frequency range can be set to 16 MHz.
+<h3>Experimental applications</h3>
-<img alt="Phononic superlattice cavity" src="https://raw.githubusercontent.com/dangars/dealii/phononic-cavity/examples/step-62/doc/step-62.07.svg?sanitize=true" height="400" />
+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.
-Systems of PDE. Elasticity. Tensors.
+Elastic equation in the frequency domain. Calculate the transmission and resonance frequency of a phononic structure.
/* ---------------------------------------------------------------------
*
- * 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.
*
#include <fstream>
#include <iostream>
-// 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 <deal.II/base/tensor.h>
#include <deal.II/base/hdf5.h>
// 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 <deal.II/numerics/vector_tools.h>
// We need this header for the function GridTools::find_active_cell_around_point
-// that we use in the function ElasticWave<dim>::store_frequency_step_data
+// that we use in the function `ElasticWave<dim>::store_frequency_step_data()`
#include <deal.II/grid/grid_tools.h>
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<hid_t>
- // 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
// 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;
const double dimension_y;
const unsigned int nb_probe_points;
const unsigned int grid_level;
- Point<dim> probe_start_point;
- Point<dim> probe_stop_point;
+ const Point<dim> probe_start_point;
+ const Point<dim> probe_stop_point;
const RightHandSide<dim> right_hand_side;
const PML<dim> pml;
const Rho<dim> rho;
- // @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 <int dim>
- class PointHistory
+ class QuadratureCache
{
public:
- PointHistory(unsigned int dofs_per_cell);
+ QuadratureCache(const unsigned int dofs_per_cell);
private:
unsigned int dofs_per_cell;
FullMatrix<std::complex<double>> mass_coefficient;
FullMatrix<std::complex<double>> stiffness_coefficient;
std::vector<std::complex<double>> right_hand_side;
- std::complex<double> 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}
// =
// 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.
class ElasticWave
{
public:
- ElasticWave(Parameters<dim> parameters_);
- ~ElasticWave();
+ ElasticWave(const Parameters<dim> ¶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.
parallel::distributed::Triangulation<dim> triangulation;
- QGauss<dim> quadrature_formula;
- const unsigned int n_q_points;
+ QGauss<dim> quadrature_formula;
- // We store the mass and stiffness matrices in this vector.
- std::vector<PointHistory<dim>> quadrature_point_history;
+ // We store the mass and stiffness matrices for each cell this vector.
+ std::vector<QuadratureCache<dim>> quadrature_cache;
- DoFHandler<dim> dof_handler;
- FESystem<dim> fe;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
// This vector contains the range of frequencies that we are going to
- // simulate
+ // simulate.
std::vector<double> frequency;
// This vector contains the coordinates $(x,y)$ of the points of the
// measurement probe.
- FullMatrix<double> position;
+ FullMatrix<double> 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;
- // @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 <int dim>
RightHandSide<dim>::RightHandSide(HDF5::Group &data)
: Function<dim>(dim)
data.get_attribute<double>("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\{
// \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}<x<x_\textrm{max}$ and
// $y_\textrm{min} <y<y_\textrm{max}$.
- // @sect4{`PML` class}
+ // @sect4{The `PML` class}
// As before, the constructor reads all the parameters from the HDF5::Group
- // `data` using the HDF5::Group::get_attribute function. As we have discussed,
- // a quadratic turn-on of the PML has been defined in the jupyter notebook. It
- // is possible to use a linear, cubic or another power degree by changing the
- // parameter pml_coeff_degree. The parameters `pml_x` and `pml_y` can be used
- // to turn on and off the `x` and `y` PMLs.
+ // `data` using the HDF5::Group::get_attribute() function. As we have
+ // discussed, a quadratic turn-on of the PML has been defined in the jupyter
+ // notebook. It is possible to use a linear, cubic or another power degree by
+ // changing the parameter `pml_coeff_degree`. The parameters `pml_x` and
+ // `pml_y` can be used to turn on and off the `x` and `y` PMLs.
template <int dim>
PML<dim>::PML(HDF5::Group &data)
: Function<dim, std::complex<double>>(dim)
- // @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 <int dim>
Rho<dim>::Rho(HDF5::Group &data)
// @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]
material_b_speed_of_sound / cavity_resonance_frequency;
// The density $\rho$ takes the following form
- //<img alt="Phononic superlattice cavity"
- // src="https://raw.githubusercontent.com/dangars/dealii/phononic-cavity/examples/step-62/doc/step-62.04.svg?sanitize=true"
+ //
<img alt="Phononic superlattice cavity"
+ // src="https://www.dealii.org/images/steps/developer/step-62.04.svg"
// height="200" />
// 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;
}
// 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;
(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;
}
// 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))
{
// 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))
{
- // @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 <int dim>
Parameters<dim>::Parameters(HDF5::Group &data)
: data(data)
, dimension_y(data.get_attribute<double>("dimension_y"))
, nb_probe_points(data.get_attribute<int>("nb_probe_points"))
, grid_level(data.get_attribute<int>("grid_level"))
+ , probe_start_point(
+ Point<dim>(data.get_attribute<double>("probe_pos_x"),
+ data.get_attribute<double>("probe_pos_y") -
+ data.get_attribute<double>("probe_width_y") / 2))
+ , probe_stop_point(
+ Point<dim>(data.get_attribute<double>("probe_pos_x"),
+ data.get_attribute<double>("probe_pos_y") +
+ data.get_attribute<double>("probe_width_y") / 2))
, right_hand_side(data)
, pml(data)
, rho(data)
- {
- probe_start_point =
- Point<dim>(data.get_attribute<double>("probe_pos_x"),
- data.get_attribute<double>("probe_pos_y") -
- data.get_attribute<double>("probe_width_y") / 2);
- probe_stop_point =
- Point<dim>(data.get_attribute<double>("probe_pos_x"),
- data.get_attribute<double>("probe_pos_y") +
- data.get_attribute<double>("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 <int dim>
- PointHistory<dim>::PointHistory(unsigned int dofs_per_cell)
+ QuadratureCache<dim>::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)
// @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 <int dim>
- ElasticWave<dim>::ElasticWave(Parameters<dim> parameters_)
- : parameters(parameters_)
+ ElasticWave<dim>::ElasticWave(const Parameters<dim> ¶meters)
+ : parameters(parameters)
, mpi_communicator(MPI_COMM_WORLD)
, triangulation(mpi_communicator,
typename Triangulation<dim>::MeshSmoothing(
Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening))
, quadrature_formula(2)
- , n_q_points(quadrature_formula.size())
- , dof_handler(triangulation)
, fe(FE_Q<dim>(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<double>(
"frequency",
std::vector<hsize_t>{parameters.nb_frequency_points}))
- , position_dataset(parameters.data.template create_dataset<double>(
+ , probe_positions_dataset(parameters.data.template create_dataset<double>(
"position",
std::vector<hsize_t>{parameters.nb_probe_points, dim}))
, displacement(
- template <int dim>
- ElasticWave<dim>::~ElasticWave()
- {
- dof_handler.clear();
- }
-
-
-
// @sect4{ElasticWave::setup_system}
// There is nothing new in this function, the only difference with step-40 is
// @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 <int dim>
- void ElasticWave<dim>::assemble_system(double omega,
- bool calculate_quadrature_data)
+ void ElasticWave<dim>::assemble_system(const double omega,
+ const bool calculate_quadrature_data)
{
TimerOutput::Scope t(computing_timer, "assembly");
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<std::complex<double>> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<std::complex<double>> cell_rhs(dofs_per_cell);
std::vector<Vector<std::complex<double>>> pml_values(
n_q_points, Vector<std::complex<double>>(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 =
}
// 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<dim> *local_quadrature_points_data =
- reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
- Assert(local_quadrature_points_data >=
- &quadrature_point_history.front(),
+ QuadratureCache<dim> *local_quadrature_points_data =
+ reinterpret_cast<QuadratureCache<dim> *>(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<dim> &quadrature_data =
+ QuadratureCache<dim> &quadrature_data =
local_quadrature_points_data[q];
// Below we declare the force vector and the parameters of the
// 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
- // <code>Tensor@<2,dim@></code>,
+ // <code>Tensor@<2,dim@></code>.
// The matrix $G_{ij}$
// consists of entries
// @f[
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 <int dim>
- void ElasticWave<dim>::set_position_vector()
+ void ElasticWave<dim>::initialize_probe_positions_vector()
{
- Point<dim> p;
for (unsigned int position_idx = 0;
position_idx < parameters.nb_probe_points;
++position_idx)
// Because of the way the operator + and - are overloaded. To substract
// two points, the following has to be done:
// `Point_b<dim> + (-Point_a<dim>)`
- 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<dim> 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];
}
}
}
// This function stores in the HDF5 file the measured energy by the probe.
template <int dim>
- void ElasticWave<dim>::store_frequency_step_data(unsigned int frequency_idx)
+ void
+ ElasticWave<dim>::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
Point<dim> 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;
{
displacement.write_none<std::complex<double>>();
}
- // 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)
{
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;
}
}
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;
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<double>();
- position_dataset.write_none<double>();
+ probe_positions_dataset.write_none<double>();
}
}
- // @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 <int dim>
- void ElasticWave<dim>::setup_quadrature_point_history()
+ void ElasticWave<dim>::setup_quadrature_cache()
{
- unsigned int our_cells = 0;
- for (typename Triangulation<dim>::active_cell_iterator cell =
- triangulation.begin_active();
- cell != triangulation.end();
- ++cell)
- if (cell->is_locally_owned())
- ++our_cells;
-
triangulation.clear_user_data();
{
- std::vector<PointHistory<dim>> tmp;
- tmp.swap(quadrature_point_history);
+ std::vector<QuadratureCache<dim>> tmp;
+ quadrature_cache.swap(tmp);
}
- quadrature_point_history.resize(our_cells * quadrature_formula.size(),
- PointHistory<dim>(fe.dofs_per_cell));
- unsigned int history_index = 0;
+ quadrature_cache.resize(triangulation.n_locally_owned_active_cells() *
+ quadrature_formula.size(),
+ QuadratureCache<dim>(fe.dofs_per_cell));
+ unsigned int cache_index = 0;
for (typename Triangulation<dim>::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 <int dim>
// 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 <int dim>
void ElasticWave<dim>::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)
}
// 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
void ElasticWave<dim>::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
{
triangulation.refine_global(parameters.grid_level);
- setup_quadrature_point_history();
+ setup_quadrature_cache();
- set_position_vector();
+ initialize_probe_positions_vector();
frequency_sweep();
}
{
- // 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<dim> parameters(calibration);
"# 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",
" 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",
" 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",
" <defs>\n",
" <path d=\"M 0 0 \n",
"L 0 3.5 \n",
- "\" id=\"m617721c33a\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
+ "\" id=\"mb3d45af095\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
" </defs>\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"89.435795\" xlink:href=\"#m617721c33a\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"89.435795\" xlink:href=\"#mb3d45af095\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_1\">\n",
" <g id=\"xtick_2\">\n",
" <g id=\"line2d_2\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"150.308523\" xlink:href=\"#m617721c33a\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"150.308523\" xlink:href=\"#mb3d45af095\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_2\">\n",
" <g id=\"xtick_3\">\n",
" <g id=\"line2d_3\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"211.18125\" xlink:href=\"#m617721c33a\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"211.18125\" xlink:href=\"#mb3d45af095\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_3\">\n",
" <g id=\"xtick_4\">\n",
" <g id=\"line2d_4\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"272.053977\" xlink:href=\"#m617721c33a\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"272.053977\" xlink:href=\"#mb3d45af095\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_4\">\n",
" <g id=\"xtick_5\">\n",
" <g id=\"line2d_5\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"332.926705\" xlink:href=\"#m617721c33a\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"332.926705\" xlink:href=\"#mb3d45af095\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_5\">\n",
" </g>\n",
" </g>\n",
" <g id=\"text_6\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:middle;\" transform=\"rotate(-0, 211.18125, 281.472062)\" x=\"211.18125\" y=\"281.472062\">frequency (MHz)</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:middle;\" transform=\"rotate(-0, 211.18125, 281.472062)\" x=\"211.18125\" y=\"281.472062\">frequency (GHz)</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"matplotlib.axis_2\">\n",
" <defs>\n",
" <path d=\"M 0 0 \n",
"L -3.5 0 \n",
- "\" id=\"m5085cc54b3\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
+ "\" id=\"m626a42a4d9\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
" </defs>\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m5085cc54b3\" y=\"243.323554\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m626a42a4d9\" y=\"243.323682\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_7\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 247.122773)\" x=\"36.78125\" y=\"247.122773\">0.0</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 247.1229)\" x=\"36.78125\" y=\"247.1229\">0.0</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_2\">\n",
" <g id=\"line2d_7\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m5085cc54b3\" y=\"203.429476\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m626a42a4d9\" y=\"203.726948\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_8\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 207.228695)\" x=\"36.78125\" y=\"207.228695\">0.2</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 207.526167)\" x=\"36.78125\" y=\"207.526167\">0.2</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_3\">\n",
" <g id=\"line2d_8\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m5085cc54b3\" y=\"163.535398\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m626a42a4d9\" y=\"164.130215\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_9\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 167.334617)\" x=\"36.78125\" y=\"167.334617\">0.4</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 167.929434)\" x=\"36.78125\" y=\"167.929434\">0.4</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_4\">\n",
" <g id=\"line2d_9\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m5085cc54b3\" y=\"123.64132\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m626a42a4d9\" y=\"124.533482\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_10\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 127.440539)\" x=\"36.78125\" y=\"127.440539\">0.6</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 128.332701)\" x=\"36.78125\" y=\"128.332701\">0.6</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_5\">\n",
" <g id=\"line2d_10\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m5085cc54b3\" y=\"83.747242\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m626a42a4d9\" y=\"84.936749\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_11\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 87.546461)\" x=\"36.78125\" y=\"87.546461\">0.8</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 88.735968)\" x=\"36.78125\" y=\"88.735968\">0.8</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_6\">\n",
" <g id=\"line2d_11\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m5085cc54b3\" y=\"43.853164\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"43.78125\" xlink:href=\"#m626a42a4d9\" y=\"45.340016\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_12\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 47.652383)\" x=\"36.78125\" y=\"47.652383\">1.0</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 36.78125, 49.139235)\" x=\"36.78125\" y=\"49.139235\">1.0</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_13\">\n",
" </g>\n",
" </g>\n",
" <g id=\"line2d_12\">\n",
- " <path clip-path=\"url(#pdfba9970d1)\" d=\"M 58.999432 243.309999 \n",
- "L 170.370587 243.186291 \n",
- "L 187.152542 242.985402 \n",
- "L 194.780703 242.704851 \n",
- "L 199.3576 242.327233 \n",
- "L 202.408864 241.844129 \n",
- "L 204.697313 241.207732 \n",
- "L 206.222945 240.525399 \n",
- "L 207.748577 239.453304 \n",
- "L 208.511394 238.672035 \n",
- "L 209.27421 237.630594 \n",
- "L 210.037026 236.200996 \n",
- "L 210.799842 234.167334 \n",
- "L 211.562658 231.144176 \n",
- "L 212.325474 226.395056 \n",
- "L 213.08829 218.39477 \n",
- "L 213.851106 203.693859 \n",
- "L 214.613923 174.026098 \n",
+ " <path clip-path=\"url(#p020971ab79)\" d=\"M 58.999432 243.309983 \n",
+ "L 170.370587 243.184651 \n",
+ "L 187.152542 242.980594 \n",
+ "L 194.780703 242.69493 \n",
+ "L 199.3576 242.309509 \n",
+ "L 202.408864 241.815212 \n",
+ "L 204.697313 241.162397 \n",
+ "L 206.222945 240.460693 \n",
+ "L 206.985761 239.974925 \n",
+ "L 207.748577 239.355064 \n",
+ "L 208.511394 238.547294 \n",
+ "L 209.27421 237.468163 \n",
+ "L 210.037026 235.982944 \n",
+ "L 210.799842 233.863454 \n",
+ "L 211.562658 230.700476 \n",
+ "L 212.325474 225.708001 \n",
+ "L 213.08829 217.249328 \n",
+ "L 213.851106 201.609491 \n",
+ "L 214.613923 169.958876 \n",
"L 216.139555 45.639136 \n",
- "L 216.902371 89.236544 \n",
- "L 217.665187 159.983972 \n",
- "L 218.428003 196.997182 \n",
- "L 219.190819 214.954687 \n",
- "L 219.953636 224.452124 \n",
- "L 220.716452 229.955602 \n",
- "L 221.479268 233.392639 \n",
- "L 222.242084 235.670014 \n",
- "L 223.0049 237.251643 \n",
- "L 223.767716 238.392511 \n",
- "L 224.530532 239.241402 \n",
- "L 226.056165 240.395366 \n",
- "L 227.581797 241.122327 \n",
- "L 229.107429 241.609167 \n",
- "L 231.395877 242.084662 \n",
- "L 234.447142 242.463115 \n",
- "L 239.024039 242.772681 \n",
- "L 246.6522 243.013231 \n",
- "L 260.38289 243.175159 \n",
- "L 292.421168 243.272418 \n",
- "L 363.363068 243.308318 \n",
- "L 363.363068 243.308318 \n",
+ "L 216.902371 96.886661 \n",
+ "L 217.665187 164.187692 \n",
+ "L 218.428003 198.841602 \n",
+ "L 219.190819 215.822253 \n",
+ "L 219.953636 224.900029 \n",
+ "L 220.716452 230.2053 \n",
+ "L 221.479268 233.540209 \n",
+ "L 222.242084 235.761076 \n",
+ "L 223.0049 237.309603 \n",
+ "L 223.767716 238.430172 \n",
+ "L 224.530532 239.266147 \n",
+ "L 226.056165 240.40596 \n",
+ "L 227.581797 241.126323 \n",
+ "L 229.107429 241.609914 \n",
+ "L 231.395877 242.083283 \n",
+ "L 234.447142 242.460909 \n",
+ "L 239.024039 242.770487 \n",
+ "L 246.6522 243.011581 \n",
+ "L 260.38289 243.174225 \n",
+ "L 292.421168 243.272107 \n",
+ "L 363.363068 243.308298 \n",
+ "L 363.363068 243.308298 \n",
"\" style=\"fill:none;stroke:#1f77b4;stroke-linecap:square;stroke-width:1.5;\"/>\n",
" </g>\n",
" <g id=\"line2d_13\">\n",
- " <path clip-path=\"url(#pdfba9970d1)\" d=\"M 58.999432 243.311864 \n",
- "L 170.370587 243.187576 \n",
- "L 187.152542 242.986329 \n",
- "L 194.780703 242.705432 \n",
- "L 199.3576 242.327459 \n",
- "L 202.408864 241.843987 \n",
- "L 204.697313 241.207191 \n",
- "L 206.222945 240.524494 \n",
- "L 207.748577 239.451915 \n",
- "L 208.511394 238.670339 \n",
- "L 209.27421 237.628534 \n",
- "L 210.037026 236.198498 \n",
- "L 210.799842 234.164308 \n",
- "L 211.562658 231.140511 \n",
- "L 212.325474 226.390637 \n",
- "L 213.08829 218.389545 \n",
- "L 213.851106 203.68813 \n",
- "L 214.613923 174.021596 \n",
- "L 216.139555 45.642454 \n",
- "L 216.902371 89.23038 \n",
- "L 217.665187 159.984195 \n",
- "L 218.428003 197.001686 \n",
- "L 219.190819 214.960534 \n",
- "L 219.953636 224.458192 \n",
- "L 220.716452 229.961522 \n",
- "L 221.479268 233.3983 \n",
- "L 222.242084 235.675398 \n",
- "L 223.0049 237.256766 \n",
- "L 223.767716 238.397397 \n",
- "L 224.530532 239.246076 \n",
- "L 226.056165 240.399685 \n",
- "L 227.581797 241.126364 \n",
- "L 229.107429 241.612978 \n",
- "L 231.395877 242.088208 \n",
- "L 234.447142 242.466401 \n",
- "L 239.024039 242.775699 \n",
- "L 246.6522 243.015975 \n",
- "L 260.38289 243.177647 \n",
- "L 292.421168 243.274681 \n",
- "L 363.363068 243.310526 \n",
- "L 363.363068 243.310526 \n",
+ " <path clip-path=\"url(#p020971ab79)\" d=\"M 58.999432 243.311864 \n",
+ "L 170.370587 243.185936 \n",
+ "L 187.152542 242.981507 \n",
+ "L 194.780703 242.695484 \n",
+ "L 199.3576 242.309692 \n",
+ "L 202.408864 241.815013 \n",
+ "L 204.697313 241.16178 \n",
+ "L 206.222945 240.459696 \n",
+ "L 206.985761 239.973693 \n",
+ "L 207.748577 239.353559 \n",
+ "L 208.511394 238.545467 \n",
+ "L 209.27421 237.465955 \n",
+ "L 210.037026 235.980277 \n",
+ "L 210.799842 233.860233 \n",
+ "L 211.562658 230.696589 \n",
+ "L 212.325474 225.703333 \n",
+ "L 213.08829 217.243853 \n",
+ "L 213.851106 201.60361 \n",
+ "L 214.613923 169.954664 \n",
+ "L 216.139555 45.641554 \n",
+ "L 216.902371 96.880455 \n",
+ "L 217.665187 164.188236 \n",
+ "L 218.428003 198.84626 \n",
+ "L 219.190819 215.828198 \n",
+ "L 219.953636 224.906185 \n",
+ "L 220.716452 230.211307 \n",
+ "L 221.479268 233.545957 \n",
+ "L 222.242084 235.766548 \n",
+ "L 223.0049 237.314812 \n",
+ "L 223.767716 238.435143 \n",
+ "L 224.530532 239.270905 \n",
+ "L 226.056165 240.410358 \n",
+ "L 227.581797 241.130436 \n",
+ "L 229.107429 241.613796 \n",
+ "L 231.395877 242.086894 \n",
+ "L 234.447142 242.464255 \n",
+ "L 239.024039 242.773559 \n",
+ "L 246.6522 243.014372 \n",
+ "L 260.38289 243.176752 \n",
+ "L 292.421168 243.274404 \n",
+ "L 363.363068 243.310534 \n",
+ "L 363.363068 243.310534 \n",
"\" style=\"fill:none;stroke:#ff7f0e;stroke-linecap:square;stroke-width:1.5;\"/>\n",
" </g>\n",
" <g id=\"patch_3\">\n",
" </g>\n",
" <g id=\"text_14\">\n",
" <text style=\"font-family:DejaVu Sans;font-size:12px;font-style:normal;font-weight:400;\" transform=\"translate(171.555 16.318125)\">Transmission</text>\n",
- " <text style=\"font-family:DejaVu Sans;font-size:12px;font-style:normal;font-weight:400;\" transform=\"translate(112.328438 29.7555)\">freq = 20.00815MHz Q = 5091.3</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:12px;font-style:normal;font-weight:400;\" transform=\"translate(109.03875 29.7555)\">freq = 20.00815GHz Q = 5046.91</text>\n",
" </g>\n",
" </g>\n",
" </g>\n",
" <defs>\n",
- " <clipPath id=\"pdfba9970d1\">\n",
+ " <clipPath id=\"p020971ab79\">\n",
" <rect height=\"217.44\" width=\"334.8\" x=\"43.78125\" y=\"35.7555\"/>\n",
" </clipPath>\n",
" </defs>\n",
"</svg>\n"
],
"text/plain": [
- "<matplotlib.figure.Figure at 0x7f2cf4228b38>"
+ "<matplotlib.figure.Figure at 0x7f785c3ecf28>"
]
},
"metadata": {},
" <defs>\n",
" <path d=\"M 0 0 \n",
"L 0 3.5 \n",
- "\" id=\"mb7c8d130cb\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
+ "\" id=\"m71cbb29084\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
" </defs>\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"97.815483\" xlink:href=\"#mb7c8d130cb\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"97.815483\" xlink:href=\"#m71cbb29084\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_1\">\n",
" <g id=\"xtick_2\">\n",
" <g id=\"line2d_2\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"158.68821\" xlink:href=\"#mb7c8d130cb\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"158.68821\" xlink:href=\"#m71cbb29084\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_2\">\n",
" <g id=\"xtick_3\">\n",
" <g id=\"line2d_3\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"219.560938\" xlink:href=\"#mb7c8d130cb\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"219.560938\" xlink:href=\"#m71cbb29084\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_3\">\n",
" <g id=\"xtick_4\">\n",
" <g id=\"line2d_4\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"280.433665\" xlink:href=\"#mb7c8d130cb\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"280.433665\" xlink:href=\"#m71cbb29084\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_4\">\n",
" <g id=\"xtick_5\">\n",
" <g id=\"line2d_5\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"341.306392\" xlink:href=\"#mb7c8d130cb\" y=\"253.1955\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"341.306392\" xlink:href=\"#m71cbb29084\" y=\"253.1955\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_5\">\n",
" </g>\n",
" </g>\n",
" <g id=\"text_6\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:middle;\" transform=\"rotate(-0, 219.560938, 281.472062)\" x=\"219.560938\" y=\"281.472062\">frequency (MHz)</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:middle;\" transform=\"rotate(-0, 219.560938, 281.472062)\" x=\"219.560938\" y=\"281.472062\">frequency (GHz)</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"matplotlib.axis_2\">\n",
" <defs>\n",
" <path d=\"M 0 0 \n",
"L -3.5 0 \n",
- "\" id=\"mb6f7814b25\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
+ "\" id=\"m02a51f3213\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
" </defs>\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#mb6f7814b25\" y=\"222.828373\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#m02a51f3213\" y=\"222.90313\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_7\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 226.627592)\" x=\"45.160938\" y=\"226.627592\">−1.0</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 226.702349)\" x=\"45.160938\" y=\"226.702349\">−1.0</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_2\">\n",
" <g id=\"line2d_7\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#mb6f7814b25\" y=\"186.887394\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#m02a51f3213\" y=\"186.940536\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_8\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 190.686613)\" x=\"45.160938\" y=\"190.686613\">−0.5</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 190.739755)\" x=\"45.160938\" y=\"190.739755\">−0.5</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_3\">\n",
" <g id=\"line2d_8\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#mb6f7814b25\" y=\"150.946415\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#m02a51f3213\" y=\"150.977942\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_9\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 154.745633)\" x=\"45.160938\" y=\"154.745633\">0.0</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 154.777161)\" x=\"45.160938\" y=\"154.777161\">0.0</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_4\">\n",
" <g id=\"line2d_9\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#mb6f7814b25\" y=\"115.005435\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#m02a51f3213\" y=\"115.015348\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_10\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 118.804654)\" x=\"45.160938\" y=\"118.804654\">0.5</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 118.814566)\" x=\"45.160938\" y=\"118.814566\">0.5</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_5\">\n",
" <g id=\"line2d_10\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#mb6f7814b25\" y=\"79.064456\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#m02a51f3213\" y=\"79.052753\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_11\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 82.863675)\" x=\"45.160938\" y=\"82.863675\">1.0</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 82.851972)\" x=\"45.160938\" y=\"82.851972\">1.0</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"ytick_6\">\n",
" <g id=\"line2d_11\">\n",
" <g>\n",
- " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#mb6f7814b25\" y=\"43.123477\"/>\n",
+ " <use style=\"stroke:#000000;stroke-width:0.8;\" x=\"52.160938\" xlink:href=\"#m02a51f3213\" y=\"43.090159\"/>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_12\">\n",
- " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 46.922695)\" x=\"45.160938\" y=\"46.922695\">1.5</text>\n",
+ " <text style=\"font-family:DejaVu Sans;font-size:10px;font-style:normal;font-weight:400;text-anchor:end;\" transform=\"rotate(-0, 45.160938, 46.889378)\" x=\"45.160938\" y=\"46.889378\">1.5</text>\n",
" </g>\n",
" </g>\n",
" <g id=\"text_13\">\n",
" </g>\n",
" </g>\n",
" <g id=\"line2d_12\">\n",
- " <path clip-path=\"url(#p3e468d8608)\" d=\"M 67.379119 166.024468 \n",
- "L 131.455674 202.631496 \n",
- "L 160.442687 218.958165 \n",
- "L 177.224642 228.180039 \n",
- "L 187.904068 233.823839 \n",
- "L 195.532229 237.618911 \n",
- "L 200.871942 240.041934 \n",
- "L 204.686023 241.556175 \n",
- "L 207.737287 242.544312 \n",
- "L 210.025736 243.074145 \n",
- "L 211.551368 243.274447 \n",
- "L 213.077 243.296555 \n",
- "L 214.602633 243.060419 \n",
- "L 215.365449 242.806402 \n",
- "L 216.128265 242.429974 \n",
- "L 216.891081 241.895939 \n",
- "L 217.653897 241.154331 \n",
- "L 218.416713 240.131859 \n",
- "L 219.179529 238.716751 \n",
- "L 219.942346 236.730351 \n",
- "L 220.705162 233.870166 \n",
- "L 221.467978 229.586474 \n",
- "L 222.230794 222.791462 \n",
- "L 222.99361 211.131255 \n",
- "L 223.756426 189.428242 \n",
- "L 225.282058 109.108975 \n",
- "L 226.044875 82.899918 \n",
- "L 226.807691 68.965345 \n",
- "L 227.570507 61.056294 \n",
- "L 228.333323 56.175162 \n",
- "L 229.096139 52.964643 \n",
- "L 229.858955 50.757375 \n",
- "L 230.621771 49.194818 \n",
- "L 231.384588 48.069323 \n",
- "L 232.147404 47.253036 \n",
- "L 232.91022 46.663208 \n",
- "L 233.673036 46.244031 \n",
- "L 235.198668 45.772806 \n",
+ " <path clip-path=\"url(#pbed3d3adf0)\" d=\"M 67.379119 166.088347 \n",
+ "L 130.692858 202.273193 \n",
+ "L 159.679871 218.610148 \n",
+ "L 176.461826 227.843396 \n",
+ "L 187.141252 233.502139 \n",
+ "L 194.769413 237.318877 \n",
+ "L 200.109126 239.771074 \n",
+ "L 203.923207 241.32225 \n",
+ "L 206.974471 242.35993 \n",
+ "L 209.26292 242.949164 \n",
+ "L 211.551368 243.282275 \n",
+ "L 213.077 243.287383 \n",
+ "L 214.602633 243.027369 \n",
+ "L 215.365449 242.757524 \n",
+ "L 216.128265 242.361532 \n",
+ "L 216.891081 241.802844 \n",
+ "L 217.653897 241.029444 \n",
+ "L 218.416713 239.964816 \n",
+ "L 219.179529 238.491884 \n",
+ "L 219.942346 236.422753 \n",
+ "L 220.705162 233.437761 \n",
+ "L 221.467978 228.952754 \n",
+ "L 222.230794 221.805147 \n",
+ "L 222.99361 209.473433 \n",
+ "L 223.756426 186.519787 \n",
+ "L 225.282058 106.374166 \n",
+ "L 226.044875 81.576783 \n",
+ "L 226.807691 68.308942 \n",
+ "L 227.570507 60.700773 \n",
+ "L 228.333323 55.968132 \n",
+ "L 229.096139 52.837706 \n",
+ "L 229.858955 50.676796 \n",
+ "L 230.621771 49.142553 \n",
+ "L 231.384588 48.035069 \n",
+ "L 232.147404 47.230592 \n",
+ "L 232.91022 46.64868 \n",
+ "L 233.673036 46.23489 \n",
+ "L 235.198668 45.769903 \n",
"L 236.7243 45.639136 \n",
- "L 238.249933 45.730986 \n",
- "L 240.538381 46.150976 \n",
- "L 242.82683 46.796841 \n",
- "L 245.878094 47.883863 \n",
- "L 249.692175 49.473546 \n",
- "L 255.031888 51.957515 \n",
- "L 261.897233 55.404039 \n",
- "L 271.813842 60.649188 \n",
- "L 286.307349 68.594483 \n",
- "L 308.429017 81.004387 \n",
- "L 345.044191 101.832791 \n",
- "L 371.742756 117.108544 \n",
- "L 371.742756 117.108544 \n",
+ "L 238.249933 45.732212 \n",
+ "L 240.538381 46.152534 \n",
+ "L 242.82683 46.797943 \n",
+ "L 245.878094 47.883969 \n",
+ "L 249.692175 49.472364 \n",
+ "L 255.031888 51.954866 \n",
+ "L 261.897233 55.400177 \n",
+ "L 271.813842 60.644621 \n",
+ "L 286.307349 68.59029 \n",
+ "L 308.429017 81.002467 \n",
+ "L 345.044191 101.836603 \n",
+ "L 371.742756 117.117192 \n",
+ "L 371.742756 117.117192 \n",
"\" style=\"fill:none;stroke:#1f77b4;stroke-linecap:square;stroke-width:1.5;\"/>\n",
" </g>\n",
" <g id=\"patch_3\">\n",
" </g>\n",
" </g>\n",
" <defs>\n",
- " <clipPath id=\"p3e468d8608\">\n",
+ " <clipPath id=\"pbed3d3adf0\">\n",
" <rect height=\"217.44\" width=\"334.8\" x=\"52.160938\" y=\"35.7555\"/>\n",
" </clipPath>\n",
" </defs>\n",
"</svg>\n"
],
"text/plain": [
- "<matplotlib.figure.Figure at 0x7f2cf60e0630>"
+ "<matplotlib.figure.Figure at 0x7f785eabc5c0>"
]
},
"metadata": {},
" 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": {
projects / dealii.git / commitdiff