From efce9f39f72a9a27d65ce27c2cb6db8d8b321773 Mon Sep 17 00:00:00 2001 From: Daniel Garcia-Sanchez Date: Fri, 5 Apr 2019 13:16:28 +0200 Subject: [PATCH] Add step-62 tutorial --- examples/CMakeLists.txt | 1 + examples/step-62/CMakeLists.txt | 39 + examples/step-62/doc/builds-on | 1 + examples/step-62/doc/intro.dox | 220 +++++ examples/step-62/doc/kind | 1 + examples/step-62/doc/results.dox | 67 ++ examples/step-62/doc/tooltip | 1 + examples/step-62/step-62.cc | 1468 ++++++++++++++++++++++++++++++ 8 files changed, 1798 insertions(+) create mode 100644 examples/step-62/CMakeLists.txt create mode 100644 examples/step-62/doc/builds-on create mode 100644 examples/step-62/doc/intro.dox create mode 100644 examples/step-62/doc/kind create mode 100644 examples/step-62/doc/results.dox create mode 100644 examples/step-62/doc/tooltip create mode 100644 examples/step-62/step-62.cc diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 81995abef7..8c4173b41f 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -31,6 +31,7 @@ IF(DEAL_II_COMPONENT_EXAMPLES) PATTERN "*.cc" PATTERN "*.prm" PATTERN "*.inp" + PATTERN "*.ipynb" PATTERN "step*/CMakeLists.txt" PATTERN "doxygen/CMakeLists.txt" # diff --git a/examples/step-62/CMakeLists.txt b/examples/step-62/CMakeLists.txt new file mode 100644 index 0000000000..28f4ae6531 --- /dev/null +++ b/examples/step-62/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-8 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-62") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) + +FIND_PACKAGE(deal.II 9.1.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-62/doc/builds-on b/examples/step-62/doc/builds-on new file mode 100644 index 0000000000..2bdacbaec3 --- /dev/null +++ b/examples/step-62/doc/builds-on @@ -0,0 +1 @@ +step-8 step-18 step-20 step-40 diff --git a/examples/step-62/doc/intro.dox b/examples/step-62/doc/intro.dox new file mode 100644 index 0000000000..7ac70d5c86 --- /dev/null +++ b/examples/step-62/doc/intro.dox @@ -0,0 +1,220 @@ +
+ +This program was contributed by Daniel Garcia-Sanchez. +
+ + +@note As a prerequisite of this program, you need to have HDF5, complex PETSc, +and the p4est libraries installed. The installation of deal.II +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 superlattice cavity is formed by two +[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 +stacks with contrasting physical properties (sound velocity index) is +repeated $N$ times. +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 + +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. +There are two regimes that depend on the waveguide width: +- Single mode: In this case the width of the structure is much + smaller that the wavelength, therefore the waveguide is single mode. + This case can be solved either with FEM (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. + 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 + 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 +[procedure](https://meep.readthedocs.io/en/latest/Python_Tutorials/Resonant_Modes_and_Transmission_in_a_Waveguide_Cavity/) +that is commonly used in time domain [FDTD](https://en.wikipedia.org/wiki/Finite-difference_time-domain_method) +simulations. A pulse at a certain frequency is generated on the left side of the +structure and the transmitted energy is measured on the right side of the structure. +The simulation is run twice. First, we run the simulation with the phononic +structure and measure the transmitted energy. + +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. + +Phononic superlattice cavity + +The transmission coefficient corresponds to the energy of the first simulation +divided by the calibration energy. +We repeat this procedure for each frequency step. + + +

Elastic equations

+The elastic equations in the time domain are +@f[ +\rho\partial_{tt} u_i - \partial_j (c_{ijkl} \varepsilon_{kl}) = f_i, +\qquad i=0,1,2 +@f] +where the stiffness tensor $c_{ijkl}$ depends on the spacial coordinates and +the strain is given by +@f[ +\varepsilon_{kl} =\frac{1}{2}(\partial_k u_l + \partial_l u_k) +@f] + +[A perfectly matched layer (PML)](https://en.wikipedia.org/wiki/Perfectly_matched_layer) +can be used to truncate the solution at the boundaries. +A PML is a transformation that results in a complex coordinate stretching. +The elastic equations in the frequency domain read as follows +@f{eqnarray*} +\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) +@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. +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 +medium is slowly varying, see +[the adiabatic theorem](https://doi.org/10.1103/PhysRevE.66.066608). +In the code a quadratic turn-on of the PML has been used. +A linear and cubic turn-on is also +[known to work](https://doi.org/10.1364/OE.16.011376). +These equations can be expanded into +@f[ +-\omega^2\rho \xi u_m - \partial_n \left(\frac{\xi}{s_n}c_{mnkl} +\varepsilon_{kl}\right) = f_m +@f] +@f[ +\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 +@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$). +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 ++ \beta_{mnkl}\partial_l u_k\right) = f_m +@f] + +We can multiply by $\varphi_m$ and integrate over the domain $\Omega$ and integrate by parts. +@f{eqnarray*} +-\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 +@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} +@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) +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. + +Each of the simulations (displacement and calibration) is stored in a separate HDF5 group: +@code{.py} +import numpy as np +import h5py +import matplotlib.pyplot as plt +import subprocess + +h5_file = h5py.File('results.h5', 'w') +data = h5_file.create_group('data') +displacement = data.create_group('displacement') +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 + + # 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 + + # Number of points in the probe + group.attrs['nb_probe_points'] = 5 + + # Global refinement + group.attrs['grid_level'] = 1 + + # Cavity + group.attrs['cavity_resonance_frequency'] = 20000000.0 + group.attrs['nb_mirror_pairs'] = 30 + + # Material + group.attrs['poissons_ratio'] = 0.27 + group.attrs['youngs_modulus'] = 270000000000.0 + group.attrs['material_a_rho'] = 3200 + if group == displacement: + group.attrs['material_b_rho'] = 2000 + else: + group.attrs['material_b_rho'] = 3200 + group.attrs['lambda'] = (group.attrs['youngs_modulus'] * group.attrs['poissons_ratio'] / + ((1 + group.attrs['poissons_ratio']) * + (1 - 2 * group.attrs['poissons_ratio']))) + 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 + + # 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_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 + + # 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 + +h5_file.close() +@endcode diff --git a/examples/step-62/doc/kind b/examples/step-62/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-62/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-62/doc/results.dox b/examples/step-62/doc/results.dox new file mode 100644 index 0000000000..43fc5d9f2e --- /dev/null +++ b/examples/step-62/doc/results.dox @@ -0,0 +1,67 @@ +

Results

+ +The results are analyzed in the +[jupyter notebook](https://github.com/dangars/dealii/blob/phononic-cavity/examples/step-62/step-62.ipynb) +with the following code +@code{.py} +h5_file = h5py.File('results.h5', 'r') +data = h5_file['data'] + +# Gaussian function that we use to fit the resonance +def resonance_f(freq, freq_m, quality_factor, max_amplitude): + omega = 2 * constants.pi * freq + omega_m = 2 * constants.pi * freq_m + gamma = omega_m / quality_factor + return max_amplitude * omega_m**2 * gamma**2 / (((omega_m**2 - omega**2)**2 + gamma**2 * omega**2)) + +frequency = data['displacement']['frequency'][...] +# Average the probe points +displacement = np.mean(data['displacement']['displacement'], axis=0) +calibration_displacement = np.mean(data['calibration']['displacement'], axis=0) +reflection_coefficient = displacement / calibration_displacement +reflectivity = (np.abs(np.mean(data['displacement']['displacement'][...]**2, axis=0))/ + np.abs(np.mean(data['calibration']['displacement'][...]**2, axis=0))) + +try: + x_data = frequency + y_data = reflectivity + quality_factor_guess = 1e3 + freq_guess = x_data[np.argmax(y_data)] + amplitude_guess = np.max(y_data) + fit_result, covariance = scipy.optimize.curve_fit(resonance_f, x_data, y_data, + [freq_guess, quality_factor_guess, amplitude_guess]) + freq_m = fit_result[0] + quality_factor = np.abs(fit_result[1]) + max_amplitude = fit_result[2] + 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.ylabel('amplitude^2 (a.u.)') + plt.title('Transmission\n' + 'freq = ' + "%.7g" % (freq_guess / 1e6) + 'MHz Q = ' + "%.6g" % quality_factor) +except: + fig = plt.figure() + plt.plot(frequency / 1e6, reflectivity) + plt.xlabel('frequency (MHz)') + 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.ylabel('phase (rad)') +plt.title('Phase (reflection coefficient)\n') + +plt.show() +h5_file.close() +@endcode + +The micropillar cavity exhibits a mechanical resonance at 20MHz and a quality factor of 5091 + +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. + +Phononic superlattice cavity diff --git a/examples/step-62/doc/tooltip b/examples/step-62/doc/tooltip new file mode 100644 index 0000000000..b57f7fe2aa --- /dev/null +++ b/examples/step-62/doc/tooltip @@ -0,0 +1 @@ +Systems of PDE. Elasticity. Tensors. diff --git a/examples/step-62/step-62.cc b/examples/step-62/step-62.cc new file mode 100644 index 0000000000..65e4d9acfb --- /dev/null +++ b/examples/step-62/step-62.cc @@ -0,0 +1,1468 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2000 - 2018 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + + * + * Author: Daniel Garcia-Sanchez, CNRS, 2019 + */ + +// @sect3{Include files} + +// Most of the include files we need for this program have already been +// discussed in previous programs, in particular in step-40. +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +// The following header provides the Tensor class that we use represent the +// material properties. +#include + + +// The following header is necessary for the HDF5 interface of deal.II. +#include + +// This header is required for the function VectorTools::point_value that we use +// to read 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 +#include + +namespace step62 +{ + using namespace dealii; + + // @sect3{Auxiliary classes and functions} + // The following classes are used to store the parameters of the simulation. + + // @sect4{`RightHandSide` class} + // This class is used to define the force pulse on the left side of the + // structure. + template + class RightHandSide : public Function + { + public: + RightHandSide(HDF5::Group &data); + virtual double value(const Point & p, + const unsigned int component) const; + + 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. + HDF5::Group data; + + // The simulation parameters are stored in `data` as HDF5 attributes. The + // following attributes are defined in the jupyter notebook, stored in + // `data` as HDF5 attributes and then read by the constructor. + const double max_force_amplitude; + const double force_sigma_x; + const double force_sigma_y; + const double max_force_width_x; + const double max_force_width_y; + const Point force_center; + + public: + // In this particular simulation the force has only a $x$ component, + // $F_y=0$. + const unsigned int force_component = 0; + }; + + // @sect4{`PML` class} + // This class is used to define the shape of the PML. + template + class PML : public Function> + { + public: + PML(HDF5::Group &data); + virtual std::complex value(const Point & p, + const unsigned int component) const; + + private: + // HDF5::Group in which all the simulation results will be stored. + HDF5::Group data; + + // The same as before, the following attributes are defined in the jupyter + // notebook, stored in `data` as HDF5 attributes and then read by the + // constructor. + const double pml_coeff; + const int pml_coeff_degree; + const double dimension_x; + const double dimension_y; + const bool pml_x; + const bool pml_y; + const double pml_width_x; + const double pml_width_y; + const double a_coeff_x; + const double a_coeff_y; + }; + + + + // @sect4{`Rho` class} + // This class is used to define the mass density. + template + class Rho : public Function + { + public: + Rho(HDF5::Group &data); + virtual double value(const Point & p, + const unsigned int component = 0) const; + + private: + // HDF5::Group in which all the simulation results will be stored. + HDF5::Group data; + + // The same as before, the following attributes are defined in the jupyter + // notebook, stored in `data` as HDF5 attributes and then read by the + // constructor. + const double lambda; + const double mu; + const double material_a_rho; + const double material_b_rho; + const double cavity_resonance_frequency; + const unsigned int nb_mirror_pairs; + const double dimension_y; + const unsigned int grid_level; + double average_rho_width; + }; + + + + // @sect4{`Parameters` class} + // This class contains all the parameters that will be used in the simulation. + template + class Parameters + { + public: + Parameters(HDF5::Group &data); + + // HDF5::Group in which all the simulation results will be stored. + HDF5::Group data; + + // The same as before, the following attributes are defined in the jupyter + // notebook, stored in `data` as HDF5 attributes and then read by the + // constructor. + const std::string simulation_name; + bool save_vtu_files; + const double start_frequency; + const double stop_frequency; + const unsigned int nb_frequency_points; + const double lambda; + const double mu; + const double dimension_x; + const double dimension_y; + const unsigned int nb_probe_points; + const unsigned int grid_level; + Point probe_start_point; + Point probe_stop_point; + const RightHandSide right_hand_side; + const PML pml; + const Rho rho; + + private: + const double comparison_float_constant = 1e-12; + }; + + + + // @sect4{`PointHistory` 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 values and re-use them at each frequency step. The + // `PointHistory` class has already been used in step-18. + + template + class PointHistory + { + public: + PointHistory(unsigned int dofs_per_cell); + + private: + unsigned int dofs_per_cell; + + public: + // We store the mass and stiffness matrices in the variables + // mass_coefficient and stiffness_coefficient. We store as well the + // right_hand_side and JxW values which are going to be the same for all the + // frequency steps. + FullMatrix> mass_coefficient; + FullMatrix> stiffness_coefficient; + std::vector> right_hand_side; + std::complex JxW; + }; + + + + // @sect4{`get_stiffness_tensor` function} + + // This class 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 + // @f[ + // c_{ijkl} + // = + // \lambda \delta_{ij} \delta_{kl} + + // \mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}). + // @f] + template + SymmetricTensor<4, dim> get_stiffness_tensor(const double lambda, + const double mu) + { + SymmetricTensor<4, dim> stiffness_tensor; + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + stiffness_tensor[i][j][k][l] = + (((i == k) && (j == l) ? mu : 0.0) + + ((i == l) && (j == k) ? mu : 0.0) + + ((i == j) && (k == l) ? lambda : 0.0)); + return stiffness_tensor; + } + + + + // @sect3{`ElasticWave` class} + + // 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 vector. + // - We save the stiffness and mass matrices in `quadrature_point_history` and + // use them for each frequency step. + // - We store the measured energy by the probe for each frequency step in the + // HDF5 file. + template + class ElasticWave + { + public: + ElasticWave(Parameters parameters_); + ~ElasticWave(); + void run(); + + private: + void setup_system(); + void assemble_system(double omega, bool calculate_quadrature_data); + void solve(); + void set_position_vector(); + void store_frequency_step_data(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 function loops over the frequency vector and runs the simulation for + // each frequency step. + void frequency_sweep(); + + // The parameters are stored in this variable. + Parameters parameters; + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + QGauss quadrature_formula; + const unsigned int n_q_points; + + // We store the mass and stiffness matrices in this vector. + std::vector> quadrature_point_history; + + DoFHandler dof_handler; + + FESystem fe; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + AffineConstraints> constraints; + + LinearAlgebraPETSc::MPI::SparseMatrix system_matrix; + LinearAlgebraPETSc::MPI::Vector locally_relevant_solution; + LinearAlgebraPETSc::MPI::Vector system_rhs; + + + // This vector contains the range of frequencies that we are going to + // simulate + std::vector frequency; + + // This vector contains the coordinates $(x,y)$ of the points of the + // measurement probe. + FullMatrix position; + + // HDF5 datasets to store the frequency and position vectors. + HDF5::DataSet frequency_dataset; + HDF5::DataSet position_dataset; + + // HDF5 dataset that stores the values of the energy measured by the proble. + HDF5::DataSet displacement; + + + ConditionalOStream pcout; + TimerOutput computing_timer; + }; + + + + // @sect3{Implementation of the auxiliary classes and functions} + + // @sect4{`RightHandSide` class} + + // The constructor reads all the parameters from the HDF5::Group `data` using + // the HDF5::Group::get_attribute function. + template + RightHandSide::RightHandSide(HDF5::Group &data) + : Function(dim) + , data(data) + , max_force_amplitude(data.get_attribute("max_force_amplitude")) + , force_sigma_x(data.get_attribute("force_sigma_x")) + , force_sigma_y(data.get_attribute("force_sigma_y")) + , max_force_width_x(data.get_attribute("max_force_width_x")) + , max_force_width_y(data.get_attribute("max_force_width_y")) + , force_center(Point(data.get_attribute("force_x_pos"), + 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 + // @f{align*} + // F_x &= + // \left\{ + // \begin{array}{ll} + // a \exp(- (\frac{(x-b_x)^2 }{ 2 \sigma_x^2}+\frac{(y-b_y)^2 }{ 2 + // \sigma_y^2})) + // & \text{if}\, x_\textrm{min} + double RightHandSide::value(const Point & p, + const unsigned int component) const + { + if (component == force_component) + { + if (std::abs(p[0] - force_center[0]) < max_force_width_x / 2 && + std::abs(p[1] - force_center[1]) < max_force_width_y / 2) + { + return max_force_amplitude * + std::exp(-(std::pow(p[0] - force_center[0], 2) / + (2 * std::pow(force_sigma_x, 2)) + + std::pow(p[1] - force_center[1], 2) / + (2 * std::pow(force_sigma_y, 2)))); + } + else + { + return 0; + } + } + else + { + return 0; + } + } + + + + // @sect4{`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. + template + PML::PML(HDF5::Group &data) + : Function>(dim) + , data(data) + , pml_coeff(data.get_attribute("pml_coeff")) + , pml_coeff_degree(data.get_attribute("pml_coeff_degree")) + , dimension_x(data.get_attribute("dimension_x")) + , dimension_y(data.get_attribute("dimension_y")) + , pml_x(data.get_attribute("pml_x")) + , pml_y(data.get_attribute("pml_y")) + , pml_width_x(data.get_attribute("pml_width_x")) + , pml_width_y(data.get_attribute("pml_width_y")) + , a_coeff_x(pml_coeff / std::pow(pml_width_x, pml_coeff_degree)) + , a_coeff_y(pml_coeff / std::pow(pml_width_y, pml_coeff_degree)) + {} + + + + // The PML coefficient for the `x` component takes the form + // $s'_x = a_x x^{\textrm{degree}}$ + template + std::complex PML::value(const Point & p, + const unsigned int component) const + { + double calculated_pml_x_coeff = 0; + double calculated_pml_y_coeff = 0; + + if ((component == 0) && pml_x) + { + const double pml_x_start_position = dimension_x / 2 - pml_width_x; + if (std::abs(p[0]) > pml_x_start_position) + { + const double x_prime = std::abs(p[0]) - pml_x_start_position; + calculated_pml_x_coeff = + a_coeff_x * std::pow(x_prime, pml_coeff_degree); + } + } + + if ((component == 1) && pml_y) + { + const double pml_y_start_position = dimension_y / 2 - pml_width_y; + if (std::abs(p[1]) > pml_y_start_position) + { + const double y_prime = std::abs(p[1]) - pml_y_start_position; + calculated_pml_y_coeff = + a_coeff_y * std::pow(y_prime, pml_coeff_degree); + } + } + + return std::complex(1, + std::max(calculated_pml_x_coeff, + calculated_pml_y_coeff)); + } + + + + // @sect4{`Rho` class} + + // This class is used to define the mass density. As we have explained, before + // a phononic superlattice cavity is formed by two + //[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 + // stacks with contrasting physical properties (sound velocity index) is + // repeated $N$ times. The change of in the velocity will be obtained by + // alternating layers with different density. + template + Rho::Rho(HDF5::Group &data) + : Function(1) + , data(data) + , lambda(data.get_attribute("lambda")) + , mu(data.get_attribute("mu")) + , material_a_rho(data.get_attribute("material_a_rho")) + , material_b_rho(data.get_attribute("material_b_rho")) + , cavity_resonance_frequency( + data.get_attribute("cavity_resonance_frequency")) + , nb_mirror_pairs(data.get_attribute("nb_mirror_pairs")) + , dimension_y(data.get_attribute("dimension_y")) + , grid_level(data.get_attribute("grid_level")) + { + // In order to increase the precision we use + // [subpixel + // smoothing](https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/). + average_rho_width = dimension_y / (std::pow(2.0, grid_level)); + data.set_attribute("average_rho_width", average_rho_width); + } + + + + template + double Rho::value(const Point &p, + const unsigned int /*component*/) const + { + // The speed of sound is defined by + // @f[ + // c = \frac{K_e}{\rho} + // @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 + // 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. + // @f[ + // K_e = 4\mu\frac{\lambda +\mu}{\lamda+2\mu} + // @f] + double elastic_constant; + if (dim == 2) + { + elastic_constant = 4 * mu * (lambda + mu) / (lambda + 2 * mu); + } + else if (dim == 3) + { + elastic_constant = mu * (3 * lambda + 2 * mu) / (lambda + mu); + } + else + { + Assert(false, ExcInternalError()); + } + const double material_a_speed_of_sound = + std::sqrt(elastic_constant / material_a_rho); + const double material_a_wavelength = + material_a_speed_of_sound / cavity_resonance_frequency; + const double material_b_speed_of_sound = + std::sqrt(elastic_constant / material_b_rho); + const double material_b_wavelength = + 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 = + 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; + return (1 - coefficient) * material_a_rho + + coefficient * material_b_rho; + } + } + + // Here we define the + // [subpixel + // smoothing](https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/) + // which improves the precision of the simulation. + for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++) + { + double layer_transition_center = + material_a_wavelength / 2 + + idx * (material_b_wavelength / 4 + material_a_wavelength / 4) + + material_b_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; + return (1 - coefficient) * material_b_rho + + coefficient * material_a_rho; + } + } + + // then the cavity + if (std::abs(p[0]) <= material_a_wavelength / 2) + { + return material_a_rho; + } + + // the material_a layers + for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++) + { + 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; + if (std::abs(p[0]) >= (layer_center - layer_width / 2) && + std::abs(p[0]) <= (layer_center + layer_width / 2)) + { + return material_a_rho; + } + } + + // the material_b layers + for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++) + { + 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; + if (std::abs(p[0]) >= (layer_center - layer_width / 2) && + std::abs(p[0]) <= (layer_center + layer_width / 2)) + { + return material_b_rho; + } + } + + // and finally the default is material_a. + return material_a_rho; + } + + + + // @sect4{`Parameters` class} + + // The constructor reads all the parameters from the HDF5::Group `data` using + // the HDF5::Group::get_attribute function. + template + Parameters::Parameters(HDF5::Group &data) + : data(data) + , simulation_name(data.get_attribute("simulation_name")) + , save_vtu_files(data.get_attribute("save_vtu_files")) + , start_frequency(data.get_attribute("start_frequency")) + , stop_frequency(data.get_attribute("stop_frequency")) + , nb_frequency_points(data.get_attribute("nb_frequency_points")) + , lambda(data.get_attribute("lambda")) + , mu(data.get_attribute("mu")) + , dimension_x(data.get_attribute("dimension_x")) + , dimension_y(data.get_attribute("dimension_y")) + , nb_probe_points(data.get_attribute("nb_probe_points")) + , grid_level(data.get_attribute("grid_level")) + , 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} + + // 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) + : dofs_per_cell(dofs_per_cell) + , mass_coefficient(dofs_per_cell, dofs_per_cell) + , stiffness_coefficient(dofs_per_cell, dofs_per_cell) + , right_hand_side(dofs_per_cell) + {} + + + + // @sect3{Implementation of the `ElasticWave` class} + + // @sect4{Constructors and destructors} + + // 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. + template + ElasticWave::ElasticWave(Parameters parameters_) + : 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) + , frequency(parameters.nb_frequency_points) + , position(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( + "position", + std::vector{parameters.nb_probe_points, dim})) + , displacement( + parameters.data.template create_dataset>( + "displacement", + std::vector{parameters.nb_probe_points, + parameters.nb_frequency_points})) + , pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + , computing_timer(mpi_communicator, + pcout, + TimerOutput::summary, + TimerOutput::wall_times) + {} + + + + template + ElasticWave::~ElasticWave() + { + dof_handler.clear(); + } + + + + // @sect4{ElasticWave::setup_system} + + // There is nothing new in this function, the only difference with step-40 is + // that we don't have to apply boundary conditions because we use the PMLs to + // truncate the domain. + template + void ElasticWave::setup_system() + { + TimerOutput::Scope t(computing_timer, "setup"); + + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern( + dsp, + dof_handler.n_locally_owned_dofs_per_processor(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + + + // @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. + template + void ElasticWave::assemble_system(double omega, + bool calculate_quadrature_data) + { + TimerOutput::Scope t(computing_timer, "assembly"); + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + FullMatrix> cell_matrix(dofs_per_cell, dofs_per_cell); + Vector> cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + // Here we store the value of the right hand side, rho and the PML. + std::vector> rhs_values(n_q_points, Vector(dim)); + std::vector rho_values(n_q_points); + std::vector>> pml_values( + n_q_points, Vector>(dim)); + + // We calculate the stiffness tensor for the $\lambda$ and $\mu$ that has + // 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 = + get_stiffness_tensor(parameters.lambda, parameters.mu); + + // We use the same method of step-20 for vector-valued problems. + const FEValuesExtractors::Vector displacement(0); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + cell_matrix = 0; + cell_rhs = 0; + + // We have to calculate the values of the right hand side, rho and + // the PML only if we are going to calculate the mass and the + // stiffness matrices. Otherwise we can skip this calculation which + // considerably reduces the total calculation time. + if (calculate_quadrature_data) + { + fe_values.reinit(cell); + + parameters.right_hand_side.vector_value_list( + fe_values.get_quadrature_points(), rhs_values); + parameters.rho.value_list(fe_values.get_quadrature_points(), + rho_values); + parameters.pml.vector_value_list( + fe_values.get_quadrature_points(), pml_values); + } + + // 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 + // 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(), + ExcInternalError()); + Assert(local_quadrature_points_data < + &quadrature_point_history.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 = + local_quadrature_points_data[q]; + + // Below we declare the force vector and the parameters of the + // PML $s$ and $\xi$. + Tensor<1, dim> force; + Tensor<1, dim, std::complex> s; + std::complex xi(1, 0); + + // The following block is calculated only in the first frequency + // step. + if (calculate_quadrature_data) + { + // Store the value of `JxW`. + quadrature_data.JxW = fe_values.JxW(q); + + for (unsigned int component = 0; component < dim; + ++component) + { + // Convert vectors to tensors and calculate xi + force[component] = rhs_values[q][component]; + s[component] = pml_values[q][component]; + xi *= s[component]; + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const Tensor<1, dim> phi_i = + fe_values[displacement].value(i, q); + const Tensor<2, dim> grad_phi_i = + fe_values[displacement].gradient(i, q); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const Tensor<1, dim> phi_j = + fe_values[displacement].value(j, q); + const Tensor<2, dim> grad_phi_j = + fe_values[displacement].gradient(j, q); + + // calculate the values of the mass matrix. + quadrature_data.mass_coefficient[i][j] = + rho_values[q] * xi * phi_i * phi_j; + + // Loop over the $mnkl$ indices of the stiffness + // tensor. + std::complex stiffness_coefficient = 0; + for (unsigned int m = 0; m < dim; ++m) + { + for (unsigned int n = 0; n < dim; ++n) + { + for (unsigned int k = 0; k < dim; ++k) + { + for (unsigned int l = 0; l < dim; ++l) + { + // Here we calculate the tensors + // $\alpha_{mnkl}$ and + // $\beta_{mnkl}$ + const std::complex alpha = + xi * + stiffness_tensor[m][n][k][l] / + (2.0 * s[n] * s[k]); + const std::complex beta = + xi * + stiffness_tensor[m][n][k][l] / + (2.0 * s[n] * s[l]); + + // Here we calculate the stiffness + // matrix. Note that the stiffness + // 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) + // which is a + // Tensor@<2,dim@>, + // The matrix $G_{ij}$ + // consists of entries + // @f[ + // G_{ij}= + // \frac{\partial\phi_i}{\partial + // x_j} + // =\partial_j \phi_i + // @f] + // Note the position of the indices + // $i$ and $j$ and the notation that + // we use in this tutorial: + // $\partial_j\phi_i$. As the + // stiffness tensor is not + // symmetric, it is very easy to + // make a mistake. + stiffness_coefficient += + grad_phi_i[m][n] * + (alpha * grad_phi_j[l][k] + + beta * grad_phi_j[k][l]); + } + } + } + } + + // We save the value of the stiffness matrix in + // quadrature_data + quadrature_data.stiffness_coefficient[i][j] = + stiffness_coefficient; + } + + // and the value of the right hand side in + // quadrature_data. + quadrature_data.right_hand_side[i] = + phi_i * force * fe_values.JxW(q); + } + } + + // We loop again over the degrees of freedom of the cells to + // calculate the system matrix. These loops are really quick + // because we have already calculated the stiffness and mass + // matrices, only the value of $\omega$ changes. + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + std::complex matrix_sum = 0; + matrix_sum += -std::pow(omega, 2) * + quadrature_data.mass_coefficient[i][j]; + matrix_sum += + quadrature_data.stiffness_coefficient[i][j]; + cell_matrix(i, j) += matrix_sum * quadrature_data.JxW; + } + cell_rhs(i) += quadrature_data.right_hand_side[i]; + } + } + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + // @sect4{ElasticWave::solve} + + // This is even more simple than in step-40. We use the parallel direct solver + // MUMPS which requires less options than an iterative solver. The drawback is + // that it does not scale very well. It is not straightforward to solve the + // Helmholtz equation with an iterative solver. The shifted Laplacian + // multigrid method is a well known approach to precondition this system, but + // this is beyond the scope of this tutorial. + template + void ElasticWave::solve() + { + TimerOutput::Scope t(computing_timer, "solve"); + LinearAlgebraPETSc::MPI::Vector completely_distributed_solution( + locally_owned_dofs, mpi_communicator); + + SolverControl solver_control; + PETScWrappers::SparseDirectMUMPS solver(solver_control, mpi_communicator); + solver.solve(system_matrix, completely_distributed_solution, system_rhs); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + constraints.distribute(completely_distributed_solution); + locally_relevant_solution = completely_distributed_solution; + } + + // @sect4{ElasticWave::set_position_vector} + + // We use this function to calculate the values of the position vector. + template + void ElasticWave::set_position_vector() + { + Point 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 + (-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]; + if (dim == 3) + { + position[position_idx][2] = p[2]; + } + } + } + + // @sect4{ElasticWave::store_frequency_step_data} + + // This function stores in the HDF5 file the measured energy by the probe. + template + void ElasticWave::store_frequency_step_data(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; + + // 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 + // displacement_data contains the value of the displacement at these points. + std::vector coordinates; + std::vector> displacement_data; + for (unsigned int position_idx = 0; + position_idx < parameters.nb_probe_points; + ++position_idx) + { + Point point; + for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) + { + point(dim_idx) = position[position_idx][dim_idx]; + } + bool point_in_locally_owned_cell; + { + // First we have to find out if the point is in a locally owned cell. + auto mapping = StaticMappingQ1::mapping; + const std::pair::active_cell_iterator, + Point> + cell_point = GridTools::find_active_cell_around_point(mapping, + dof_handler, + point); + + point_in_locally_owned_cell = cell_point.first->is_locally_owned(); + } + if (point_in_locally_owned_cell) + { + // Then we can store the values of the displacement in the points of + // the probe in `displacement_data`. + Vector> tmp_vector(dim); + VectorTools::point_value(dof_handler, + locally_relevant_solution, + point, + tmp_vector); + coordinates.emplace_back(position_idx); + coordinates.emplace_back(frequency_idx); + displacement_data.emplace_back( + tmp_vector(probe_displacement_component)); + } + } + + // We write the displacement data in the HDF5 file. The call + // HDF5::DataSet::write_selection() is MPI collective which means that all + // the processes have to participate. + if (coordinates.size() > 0) + { + displacement.write_selection(displacement_data, coordinates); + } + // Therefore even if the process has no data to write it has to participate + // in the collective call. For this we can use HDF5::DataSet::write_none(). + // Note that we have to specify the data type, in this case + // `std::complex`. + else + { + displacement.write_none>(); + } + + // If the variable of the jupyter notbook `save_vtu_files = 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) + { + std::vector solution_names(1, "displacement_x"); + if (dim >= 2) + { + solution_names.emplace_back("displacement_y"); + } + if (dim == 3) + { + solution_names.emplace_back("displacement_z"); + } + std::vector + interpretation(dim, DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.add_data_vector(dof_handler, + locally_relevant_solution, + solution_names, + interpretation); + Vector subdomain(triangulation.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + subdomain(i) = triangulation.locally_owned_subdomain(); + data_out.add_data_vector(subdomain, "subdomain"); + + std::vector> force( + dim, Vector(triangulation.n_active_cells())); + std::vector> pml( + dim, Vector(triangulation.n_active_cells())); + Vector rho(triangulation.n_active_cells()); + + for (auto cell : triangulation.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + 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()) = + parameters.pml.value(cell->center(), dim_idx).imag(); + } + rho(cell->active_cell_index()) = + parameters.rho.value(cell->center()); + } + // And on the cells that we are not interested in, set the + // respective value to a bogus value in order to make sure that if + // we were somehow wrong about our assumption we would find out by + // looking at the graphical output: + else + { + 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; + } + } + } + + for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx) + { + data_out.add_data_vector(force[dim_idx], + "force_" + std::to_string(dim_idx)); + data_out.add_data_vector(pml[dim_idx], + "pml_" + std::to_string(dim_idx)); + } + data_out.add_data_vector(rho, "rho"); + + data_out.build_patches(); + + unsigned int nb_number_positions; + std::stringstream frequency_idx_stream; + 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; + std::string filename = (parameters.simulation_name + "_" + + frequency_idx_stream.str() + ".vtu"); + data_out.write_vtu_in_parallel(filename.c_str(), mpi_communicator); + } + } + + + + // @sect4{ElasticWave::output_results} + + // This function writes the datasets that have not already been written. + template + void ElasticWave::output_results() + { + // The vectors `frequency` and `position` are the same for all the + // processes. Therefore any of the processes can write the corresponding + // `datasets`. Because the call HDF5::DataSet::write is MPI collective, the + // rest of the processes will have to call HDF5::DataSet::write_none. + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + frequency_dataset.write(frequency); + position_dataset.write(position); + } + else + { + frequency_dataset.write_none(); + position_dataset.write_none(); + } + } + + + + // @sect4{ElasticWave::setup_quadrature_point_history} + + // 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. + template + void ElasticWave::setup_quadrature_point_history() + { + 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); + } + + quadrature_point_history.resize(our_cells * quadrature_formula.size(), + PointHistory(fe.dofs_per_cell)); + unsigned int history_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(); + } + Assert(history_index == quadrature_point_history.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. + 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; + + + + 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; + } + + if (frequency_idx == 0) + { + // Write the simulation parameters only once + parameters.data.set_attribute("active_cells", + triangulation.n_active_cells()); + parameters.data.set_attribute("degrees_of_freedom", + dof_handler.n_dofs()); + } + + // We calculate the frequency and omega values for this particular step. + 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; + + // In the first frequency step we calculate the mass and stiffness + // matrices and the right hand side. In the subsequent frequency steps + // we will use those values. This improves considerably the calculation + // time. + assemble_system(current_loop_omega, + (frequency_idx == 0) ? true : false); + solve(); + + frequency[frequency_idx] = current_loop_frequency; + store_frequency_step_data(frequency_idx); + + computing_timer.print_summary(); + computing_timer.reset(); + pcout << std::endl; + } + } + + + + // @sect4{ElasticWave::run} + + // This function is very similar to the one in step-40. + template + void ElasticWave::run() + { +#ifdef DEBUG + std::cout << "Debug mode" << std::endl; +#else + std::cout << "Release mode" << std::endl; +#endif + + { + Point p1; + p1(0) = -parameters.dimension_x / 2; + p1(1) = -parameters.dimension_y / 2; + if (dim == 3) + { + p1(2) = -parameters.dimension_y / 2; + } + Point p2; + p2(0) = parameters.dimension_x / 2; + p2(1) = parameters.dimension_y / 2; + if (dim == 3) + { + p2(2) = parameters.dimension_y / 2; + } + std::vector divisions(dim); + divisions[0] = int(parameters.dimension_x / parameters.dimension_y); + divisions[1] = 1; + if (dim == 3) + { + divisions[2] = 1; + } + GridGenerator::subdivided_hyper_rectangle(triangulation, + divisions, + p1, + p2); + } + + triangulation.refine_global(parameters.grid_level); + + setup_quadrature_point_history(); + + set_position_vector(); + + frequency_sweep(); + + output_results(); + } +} // namespace step62 + +using namespace dealii; + +// @sect4{The main function} + +// The main function is very similar to the one in step-40. +int main(int argc, char *argv[]) +{ + try + { + const int dim = 2; + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + HDF5::File data_file("results.h5", + HDF5::File::FileAccessMode::open, + MPI_COMM_WORLD); + HDF5::Group data = data_file.open_group("data"); + + { + // Displacement simulation. The parameters are read from the + // displacement HDF5 group and the results are saved in the same HDF5 + // group. + auto displacement = data.open_group("displacement"); + step62::Parameters parameters(displacement); + + step62::ElasticWave elastic_problem(parameters); + elastic_problem.run(); + } + + { + // Calibration simulation. The parameters are read from the displacement + // HDF5 group and the results are saved in the same HDF5 group. + auto calibration = data.open_group("calibration"); + step62::Parameters parameters(calibration); + + step62::ElasticWave elastic_problem(parameters); + elastic_problem.run(); + } + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} -- 2.39.5