1019 * The speed of sound is defined by
1021 * c = \frac{K_e}{\rho}
1023 * where @f$K_e@f$ is the effective elastic constant and @f$\rho@f$ the density.
1024 * Here we consider the
case in which the waveguide width is much smaller
1025 * than the wavelength. In
this case it can be shown that
for the two
1028 * K_e = 4\mu\frac{\lambda +\mu}{\lambda+2\mu}
1030 * and
for the three dimensional
case @f$K_e@f$ is
equal to the Young
's modulus.
1032 * K_e = \mu\frac{3\lambda +2\mu}{\lambda+\mu}
1036 * double elastic_constant;
1039 * elastic_constant = 4 * mu * (lambda + mu) / (lambda + 2 * mu);
1041 * else if (dim == 3)
1043 * elastic_constant = mu * (3 * lambda + 2 * mu) / (lambda + mu);
1047 * Assert(false, ExcInternalError());
1049 * const double material_a_speed_of_sound =
1050 * std::sqrt(elastic_constant / material_a_rho);
1051 * const double material_a_wavelength =
1052 * material_a_speed_of_sound / cavity_resonance_frequency;
1053 * const double material_b_speed_of_sound =
1054 * std::sqrt(elastic_constant / material_b_rho);
1055 * const double material_b_wavelength =
1056 * material_b_speed_of_sound / cavity_resonance_frequency;
1060 * The density @f$\rho@f$ takes the following form
1061 * <img alt="Phononic superlattice cavity"
1062 * src="https://www.dealii.org/images/steps/developer/step-62.04.svg"
1064 * where the brown color represents material_a and the green color
1065 * represents material_b.
1068 * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
1070 * const double layer_transition_center =
1071 * material_a_wavelength / 2 +
1072 * idx * (material_b_wavelength / 4 + material_a_wavelength / 4);
1073 * if (std::abs(p[0]) >=
1074 * (layer_transition_center - average_rho_width / 2) &&
1075 * std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1077 * const double coefficient =
1079 * (layer_transition_center - average_rho_width / 2)) /
1080 * average_rho_width;
1081 * return (1 - coefficient) * material_a_rho +
1082 * coefficient * material_b_rho;
1088 * Here we define the
1090 * smoothing](https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/)
1091 * which improves the precision of the simulation.
1094 * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
1096 * const double layer_transition_center =
1097 * material_a_wavelength / 2 +
1098 * idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1099 * material_b_wavelength / 4;
1100 * if (std::abs(p[0]) >=
1101 * (layer_transition_center - average_rho_width / 2) &&
1102 * std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1104 * const double coefficient =
1106 * (layer_transition_center - average_rho_width / 2)) /
1107 * average_rho_width;
1108 * return (1 - coefficient) * material_b_rho +
1109 * coefficient * material_a_rho;
1118 * if (std::abs(p[0]) <= material_a_wavelength / 2)
1120 * return material_a_rho;
1125 * the material_a layers
1128 * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
1130 * const double layer_center =
1131 * material_a_wavelength / 2 +
1132 * idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1133 * material_b_wavelength / 4 + material_a_wavelength / 8;
1134 * const double layer_width = material_a_wavelength / 4;
1135 * if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1136 * std::abs(p[0]) <= (layer_center + layer_width / 2))
1138 * return material_a_rho;
1144 * the material_b layers
1147 * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
1149 * const double layer_center =
1150 * material_a_wavelength / 2 +
1151 * idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1152 * material_b_wavelength / 8;
1153 * const double layer_width = material_b_wavelength / 4;
1154 * if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1155 * std::abs(p[0]) <= (layer_center + layer_width / 2))
1157 * return material_b_rho;
1163 * and finally the default is material_a.
1166 * return material_a_rho;
1174 * <a name="TheParametersclassimplementation"></a>
1175 * <h4>The `Parameters` class implementation</h4>
1179 * The constructor reads all the parameters from the HDF5::Group `data` using
1180 * the HDF5::Group::get_attribute() function.
1183 * template <int dim>
1184 * Parameters<dim>::Parameters(HDF5::Group &data)
1186 * , simulation_name(data.get_attribute<std::string>("simulation_name"))
1187 * , save_vtu_files(data.get_attribute<bool>("save_vtu_files"))
1188 * , start_frequency(data.get_attribute<double>("start_frequency"))
1189 * , stop_frequency(data.get_attribute<double>("stop_frequency"))
1190 * , nb_frequency_points(data.get_attribute<int>("nb_frequency_points"))
1191 * , lambda(data.get_attribute<double>("lambda"))
1192 * , mu(data.get_attribute<double>("mu"))
1193 * , dimension_x(data.get_attribute<double>("dimension_x"))
1194 * , dimension_y(data.get_attribute<double>("dimension_y"))
1195 * , nb_probe_points(data.get_attribute<int>("nb_probe_points"))
1196 * , grid_level(data.get_attribute<int>("grid_level"))
1197 * , probe_start_point(data.get_attribute<double>("probe_pos_x"),
1198 * data.get_attribute<double>("probe_pos_y") -
1199 * data.get_attribute<double>("probe_width_y") / 2)
1200 * , probe_stop_point(data.get_attribute<double>("probe_pos_x"),
1201 * data.get_attribute<double>("probe_pos_y") +
1202 * data.get_attribute<double>("probe_width_y") / 2)
1203 * , right_hand_side(data)
1213 * <a name="TheQuadratureCacheclassimplementation"></a>
1214 * <h4>The `QuadratureCache` class implementation</h4>
1218 * We need to reserve enough space for the mass and stiffness matrices and the
1219 * right hand side vector.
1222 * template <int dim>
1223 * QuadratureCache<dim>::QuadratureCache(const unsigned int dofs_per_cell)
1224 * : dofs_per_cell(dofs_per_cell)
1225 * , mass_coefficient(dofs_per_cell, dofs_per_cell)
1226 * , stiffness_coefficient(dofs_per_cell, dofs_per_cell)
1227 * , right_hand_side(dofs_per_cell)
1235 * <a name="ImplementationoftheElasticWaveclass"></a>
1236 * <h3>Implementation of the `ElasticWave` class</h3>
1241 * <a name="Constructor"></a>
1242 * <h4>Constructor</h4>
1246 * This is very similar to the constructor of @ref step_40 "step-40". In addition we create
1247 * the HDF5 datasets `frequency_dataset`, `position_dataset` and
1248 * `displacement`. Note the use of the `template` keyword for the creation of
1249 * the HDF5 datasets. It is a C++ requirement to use the `template` keyword in
1250 * order to treat `create_dataset` as a dependent template name.
1253 * template <int dim>
1254 * ElasticWave<dim>::ElasticWave(const Parameters<dim> ¶meters)
1255 * : parameters(parameters)
1256 * , mpi_communicator(MPI_COMM_WORLD)
1257 * , triangulation(mpi_communicator,
1258 * typename Triangulation<dim>::MeshSmoothing(
1259 * Triangulation<dim>::smoothing_on_refinement |
1260 * Triangulation<dim>::smoothing_on_coarsening))
1261 * , quadrature_formula(2)
1262 * , fe(FE_Q<dim>(1), dim)
1263 * , dof_handler(triangulation)
1264 * , frequency(parameters.nb_frequency_points)
1265 * , probe_positions(parameters.nb_probe_points, dim)
1266 * , frequency_dataset(parameters.data.template create_dataset<double>(
1268 * std::vector<hsize_t>{parameters.nb_frequency_points}))
1269 * , probe_positions_dataset(parameters.data.template create_dataset<double>(
1271 * std::vector<hsize_t>{parameters.nb_probe_points, dim}))
1273 * parameters.data.template create_dataset<std::complex<double>>(
1275 * std::vector<hsize_t>{parameters.nb_probe_points,
1276 * parameters.nb_frequency_points}))
1277 * , pcout(std::cout,
1278 * (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
1279 * , computing_timer(mpi_communicator,
1281 * TimerOutput::summary,
1282 * TimerOutput::wall_times)
1290 * <a name="ElasticWavesetup_system"></a>
1291 * <h4>ElasticWave::setup_system</h4>
1295 * There is nothing new in this function, the only difference with @ref step_40 "step-40" is
1296 * that we don't have to
apply boundary conditions because we use the PMLs to
1297 * truncate the domain.
1300 *
template <
int dim>
1301 *
void ElasticWave<dim>::setup_system()
1305 * dof_handler.distribute_dofs(fe);
1307 * locally_owned_dofs = dof_handler.locally_owned_dofs();
1310 * locally_relevant_solution.reinit(locally_owned_dofs,
1311 * locally_relevant_dofs,
1312 * mpi_communicator);
1314 * system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1316 * constraints.clear();
1317 * constraints.reinit(locally_relevant_dofs);
1320 * constraints.close();
1326 * locally_owned_dofs,
1328 * locally_relevant_dofs);
1330 * system_matrix.reinit(locally_owned_dofs,
1331 * locally_owned_dofs,
1333 * mpi_communicator);
1341 * <a name=
"ElasticWaveassemble_system"></a>
1342 * <h4>ElasticWave::assemble_system</h4>
1346 * This
function is also very similar to @ref step_40
"step-40", though there are notable
1347 * differences. We
assemble the system
for each frequency/omega step. In the
1348 *
first step we
set `calculate_quadrature_data = True` and we calculate the
1349 * mass and stiffness matrices and the right hand side vector. In the
1350 * subsequent steps we will use that data to accelerate the calculation.
1353 *
template <
int dim>
1354 *
void ElasticWave<dim>::assemble_system(
const double omega,
1355 *
const bool calculate_quadrature_data)
1360 * quadrature_formula,
1363 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
1364 *
const unsigned int n_q_points = quadrature_formula.size();
1369 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1373 * Here we store the
value of the right hand side, rho and the PML.
1376 * std::vector<Vector<double>> rhs_values(n_q_points,
Vector<double>(dim));
1377 * std::vector<double> rho_values(n_q_points);
1378 * std::vector<Vector<std::complex<double>>> pml_values(
1379 * n_q_points,
Vector<std::complex<double>>(dim));
1383 * We calculate the stiffness tensor
for the @f$\lambda@f$ and @f$\mu@f$ that have
1384 * been defined in the jupyter notebook. Note that contrary to @f$\rho@f$ the
1385 * stiffness is constant among
for the whole domain.
1389 * get_stiffness_tensor<dim>(parameters.lambda, parameters.mu);
1393 * We use the same method of @ref step_20
"step-20" for vector-valued problems.
1398 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1399 *
if (cell->is_locally_owned())
1406 * We have to calculate the values of the right hand side, rho and
1407 * the PML only
if we are going to calculate the mass and the
1408 * stiffness matrices. Otherwise we can skip
this calculation which
1409 * considerably reduces the total calculation time.
1412 *
if (calculate_quadrature_data)
1414 * fe_values.reinit(cell);
1416 * parameters.right_hand_side.vector_value_list(
1417 * fe_values.get_quadrature_points(), rhs_values);
1418 * parameters.rho.value_list(fe_values.get_quadrature_points(),
1420 * parameters.pml.vector_value_list(
1421 * fe_values.get_quadrature_points(), pml_values);
1426 * We have done
this in @ref step_18
"step-18". Get a pointer to the quadrature
1427 * cache data local to the present cell, and, as a defensive
1428 * measure, make sure that
this pointer is within the bounds of the
1432 * QuadratureCache<dim> *local_quadrature_points_data =
1433 *
reinterpret_cast<QuadratureCache<dim> *
>(cell->user_pointer());
1434 *
Assert(local_quadrature_points_data >= &quadrature_cache.front(),
1436 *
Assert(local_quadrature_points_data <= &quadrature_cache.back(),
1438 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1442 * The quadrature_data variable is used to store the mass and
1443 * stiffness matrices, the right hand side vector and the
value
1447 * QuadratureCache<dim> &quadrature_data =
1448 * local_quadrature_points_data[q];
1452 * Below we declare the force vector and the parameters of the
1453 * PML @f$s@f$ and @f$\xi@f$.
1458 * std::complex<double> xi(1, 0);
1462 * The following block is calculated only in the
first frequency
1466 *
if (calculate_quadrature_data)
1470 * Store the
value of `JxW`.
1473 * quadrature_data.JxW = fe_values.JxW(q);
1475 *
for (
unsigned int component = 0; component < dim; ++component)
1479 * Convert vectors to tensors and calculate xi
1482 * force[component] = rhs_values[q][component];
1483 * s[component] = pml_values[q][component];
1484 * xi *= s[component];
1486 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1489 * fe_values[displacement].value(i, q);
1491 * fe_values[displacement].gradient(i, q);
1493 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1496 * fe_values[displacement].value(j, q);
1498 * fe_values[displacement].gradient(j, q);
1502 * calculate the values of the mass
matrix.
1505 * quadrature_data.mass_coefficient[i][j] =
1506 * rho_values[q] * xi * phi_i * phi_j;
1510 * Loop over the @f$mnkl@f$ indices of the stiffness
1514 * std::complex<double> stiffness_coefficient = 0;
1515 *
for (
unsigned int m = 0; m < dim; ++m)
1516 *
for (
unsigned int n = 0; n < dim; ++n)
1517 *
for (
unsigned int k = 0; k < dim; ++k)
1518 *
for (
unsigned int l = 0;
l < dim; ++
l)
1522 * Here we calculate the tensors
1523 * @f$\alpha_{mnkl}@f$ and @f$\beta_{mnkl}@f$
1526 *
const std::complex<double> alpha =
1527 * xi * stiffness_tensor[m][n][k][
l] /
1528 * (2.0 * s[n] * s[k]);
1529 *
const std::complex<double> beta =
1530 * xi * stiffness_tensor[m][n][k][
l] /
1531 * (2.0 * s[n] * s[
l]);
1535 * Here we calculate the stiffness
matrix.
1536 * Note that the stiffness
matrix is not
1537 *
symmetric because of the PMLs. We use the
1538 * gradient
function (see the
1539 * [documentation](https:
1540 * which is a <code>
Tensor@<2,dim@></code>.
1541 * The
matrix @f$G_{ij}@f$ consists of entries
1544 * \frac{\partial\phi_i}{\partial x_j}
1545 * =\partial_j \phi_i
1547 * Note the position of the indices @f$i@f$ and
1548 * @f$j@f$ and the notation that we use in
this
1549 * tutorial: @f$\partial_j\phi_i@f$. As the
1550 * stiffness tensor is not
symmetric, it is
1551 * very easy to make a mistake.
1554 * stiffness_coefficient +=
1555 * grad_phi_i[m][n] *
1556 * (alpha * grad_phi_j[
l][k] +
1557 * beta * grad_phi_j[k][
l]);
1566 * quadrature_data.stiffness_coefficient[i][j] =
1567 * stiffness_coefficient;
1572 * and the
value of the right hand side in
1576 * quadrature_data.right_hand_side[i] =
1577 * phi_i * force * fe_values.JxW(q);
1583 * We
loop again over the degrees of freedom of the cells to
1584 * calculate the system
matrix. These loops are really quick
1585 * because we have already calculated the stiffness and mass
1586 * matrices, only the
value of @f$\omega@f$ changes.
1589 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1591 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1593 * std::complex<double> matrix_sum = 0;
1594 * matrix_sum += -
std::pow(omega, 2) *
1595 * quadrature_data.mass_coefficient[i][j];
1596 * matrix_sum += quadrature_data.stiffness_coefficient[i][j];
1597 *
cell_matrix(i, j) += matrix_sum * quadrature_data.JxW;
1599 * cell_rhs(i) += quadrature_data.right_hand_side[i];
1602 * cell->get_dof_indices(local_dof_indices);
1603 * constraints.distribute_local_to_global(
cell_matrix,
1605 * local_dof_indices,
1617 * <a name=
"ElasticWavesolve"></a>
1618 * <h4>ElasticWave::solve</h4>
1622 * This is even more simple than in @ref step_40
"step-40". We use the
parallel direct solver
1623 * MUMPS which requires less options than an iterative solver. The drawback is
1624 * that it does not
scale very well. It is not straightforward to solve the
1625 * Helmholtz equation with an iterative solver. The shifted Laplacian
1626 * multigrid method is a well known approach to precondition
this system, but
1627 *
this is beyond the scope of
this tutorial.
1630 *
template <
int dim>
1631 *
void ElasticWave<dim>::solve()
1635 * locally_owned_dofs, mpi_communicator);
1639 * solver.solve(system_matrix, completely_distributed_solution, system_rhs);
1641 * pcout <<
" Solved in " << solver_control.
last_step() <<
" iterations."
1643 * constraints.distribute(completely_distributed_solution);
1644 * locally_relevant_solution = completely_distributed_solution;
1650 * <a name=
"ElasticWaveinitialize_position_vector"></a>
1651 * <h4>ElasticWave::initialize_position_vector</h4>
1655 * We use
this function to calculate the values of the position vector.
1658 *
template <
int dim>
1659 *
void ElasticWave<dim>::initialize_probe_positions_vector()
1661 *
for (
unsigned int position_idx = 0;
1662 * position_idx < parameters.nb_probe_points;
1667 * Because of the way the
operator + and - are overloaded to subtract
1668 * two points, the following has to be done:
1669 * `Point_b<dim> + (-Point_a<dim>)`
1673 * (position_idx / ((
double)(parameters.nb_probe_points - 1))) *
1674 * (parameters.probe_stop_point + (-parameters.probe_start_point)) +
1675 * parameters.probe_start_point;
1676 * probe_positions[position_idx][0] = p[0];
1677 * probe_positions[position_idx][1] = p[1];
1680 * probe_positions[position_idx][2] = p[2];
1688 * <a name=
"ElasticWavestore_frequency_step_data"></a>
1689 * <h4>ElasticWave::store_frequency_step_data</h4>
1693 * This
function stores in the
HDF5 file the measured energy by the probe.
1696 *
template <
int dim>
1698 * ElasticWave<dim>::store_frequency_step_data(
const unsigned int frequency_idx)
1704 * We store the displacement in the @f$x@f$ direction; the displacement in the
1705 * @f$y@f$ direction is negligible.
1708 *
const unsigned int probe_displacement_component = 0;
1712 * The vector coordinates contains the coordinates in the
HDF5 file of the
1713 * points of the probe that are located in locally owned cells. The vector
1714 * displacement_data contains the
value of the displacement at these points.
1717 * std::vector<hsize_t> coordinates;
1718 * std::vector<std::complex<double>> displacement_data;
1719 *
for (
unsigned int position_idx = 0;
1720 * position_idx < parameters.nb_probe_points;
1724 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1726 *
point[dim_idx] = probe_positions[position_idx][dim_idx];
1728 *
bool point_in_locally_owned_cell;
1732 * First we have to find out
if the
point is in a locally owned cell.
1736 *
const std::pair<typename DoFHandler<dim>::active_cell_iterator,
1742 * point_in_locally_owned_cell = cell_point.first->is_locally_owned();
1744 *
if (point_in_locally_owned_cell)
1748 * Then we can store the values of the displacement in the points of
1749 * the probe in `displacement_data`.
1754 * locally_relevant_solution,
1757 * coordinates.emplace_back(position_idx);
1758 * coordinates.emplace_back(frequency_idx);
1759 * displacement_data.emplace_back(
1760 * tmp_vector(probe_displacement_component));
1766 * We write the displacement data in the
HDF5 file. The
call
1768 * the processes have to participate.
1771 * if (coordinates.size() > 0)
1773 * displacement.write_selection(displacement_data, coordinates);
1777 * Therefore even
if the process has no data to write it has to participate
1779 * Note that we have to specify the data type, in
this case
1780 * `std::complex<double>`.
1785 * displacement.write_none<std::complex<double>>();
1790 * If the variable `save_vtu_files` in the input file equals `True` then all
1791 * the data will be saved as
vtu. The procedure to write `
vtu` files has
1792 * been described in @ref step_40
"step-40".
1795 *
if (parameters.save_vtu_files)
1797 * std::vector<std::string> solution_names(dim,
"displacement");
1798 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1804 * locally_relevant_solution,
1808 *
for (
unsigned int i = 0; i < subdomain.size(); ++i)
1812 * std::vector<Vector<double>> force(
1814 * std::vector<Vector<double>> pml(
1820 *
if (cell->is_locally_owned())
1822 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1824 * force[dim_idx](cell->active_cell_index()) =
1825 * parameters.right_hand_side.value(cell->center(), dim_idx);
1826 * pml[dim_idx](cell->active_cell_index()) =
1827 * parameters.pml.value(cell->center(), dim_idx).imag();
1829 * rho(cell->active_cell_index()) =
1830 * parameters.rho.value(cell->center());
1834 * And on the cells that we are not interested in,
set the
1835 * respective
value to a bogus
value in order to make sure that
if
1836 * we were somehow wrong about our assumption we would find out by
1837 * looking at the graphical output:
1842 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1844 * force[dim_idx](cell->active_cell_index()) = -1
e+20;
1845 * pml[dim_idx](cell->active_cell_index()) = -1
e+20;
1847 * rho(cell->active_cell_index()) = -1
e+20;
1851 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1862 * std::stringstream frequency_idx_stream;
1863 *
const unsigned int nb_number_positions =
1864 * ((
unsigned int)
std::log10(parameters.nb_frequency_points)) + 1;
1865 * frequency_idx_stream << std::setw(nb_number_positions)
1866 * << std::setfill(
'0') << frequency_idx;
1867 * std::string filename = (parameters.simulation_name +
"_" +
1868 * frequency_idx_stream.str() +
".vtu");
1878 * <a name=
"ElasticWaveoutput_results"></a>
1879 * <h4>ElasticWave::output_results</h4>
1883 * This
function writes the datasets that have not already been written.
1886 *
template <
int dim>
1887 *
void ElasticWave<dim>::output_results()
1891 * The vectors `frequency` and `position` are the same
for all the
1892 * processes. Therefore any of the processes can write the corresponding
1899 * frequency_dataset.write(frequency);
1900 * probe_positions_dataset.write(probe_positions);
1904 * frequency_dataset.write_none<
double>();
1905 * probe_positions_dataset.write_none<
double>();
1914 * <a name=
"ElasticWavesetup_quadrature_cache"></a>
1915 * <h4>ElasticWave::setup_quadrature_cache</h4>
1919 * We use
this function at the beginning of our computations to
set up
initial
1920 * values of the cache variables. This
function has been described in @ref step_18
"step-18".
1921 * There are no differences with the
function of @ref step_18
"step-18".
1924 *
template <
int dim>
1925 *
void ElasticWave<dim>::setup_quadrature_cache()
1930 * std::vector<QuadratureCache<dim>> tmp;
1931 * quadrature_cache.swap(tmp);
1934 * quadrature_cache.resize(
triangulation.n_locally_owned_active_cells() *
1935 * quadrature_formula.size(),
1936 * QuadratureCache<dim>(fe.dofs_per_cell));
1937 *
unsigned int cache_index = 0;
1938 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1939 *
if (cell->is_locally_owned())
1941 * cell->set_user_pointer(&quadrature_cache[cache_index]);
1942 * cache_index += quadrature_formula.size();
1952 * <a name=
"ElasticWavefrequency_sweep"></a>
1953 * <h4>ElasticWave::frequency_sweep</h4>
1957 * For clarity we divide the
function `
run` of @ref step_40
"step-40" into the
functions
1958 * `
run` and `frequency_sweep`. In the
function `frequency_sweep` we place the
1959 * iteration over the frequency vector.
1962 *
template <
int dim>
1963 *
void ElasticWave<dim>::frequency_sweep()
1965 *
for (
unsigned int frequency_idx = 0;
1966 * frequency_idx < parameters.nb_frequency_points;
1969 * pcout << parameters.simulation_name +
" frequency idx: "
1970 * << frequency_idx <<
'/' << parameters.nb_frequency_points - 1
1976 *
if (frequency_idx == 0)
1978 * pcout <<
" Number of active cells : "
1980 * pcout <<
" Number of degrees of freedom : "
1981 * << dof_handler.n_dofs() << std::endl;
1984 *
if (frequency_idx == 0)
1988 * Write the simulation parameters only once
1991 * parameters.data.set_attribute(
"active_cells",
1993 * parameters.data.set_attribute(
"degrees_of_freedom",
1994 * dof_handler.n_dofs());
1999 * We calculate the frequency and omega values
for this particular step.
2002 *
const double current_loop_frequency =
2003 * (parameters.start_frequency +
2005 * (parameters.stop_frequency - parameters.start_frequency) /
2006 * (parameters.nb_frequency_points - 1));
2007 *
const double current_loop_omega =
2012 * In the
first frequency step we calculate the mass and stiffness
2013 * matrices and the right hand side. In the subsequent frequency steps
2014 * we will use those values. This improves considerably the calculation
2018 * assemble_system(current_loop_omega,
2019 * (frequency_idx == 0) ?
true :
false);
2022 * frequency[frequency_idx] = current_loop_frequency;
2023 * store_frequency_step_data(frequency_idx);
2025 * computing_timer.print_summary();
2026 * computing_timer.reset();
2027 * pcout << std::endl;
2036 * <a name=
"ElasticWaverun"></a>
2041 * This
function is very similar to the
one in @ref step_40
"step-40".
2044 *
template <
int dim>
2048 * pcout <<
"Debug mode" << std::endl;
2050 * pcout <<
"Release mode" << std::endl;
2055 * p1(0) = -parameters.dimension_x / 2;
2056 * p1(1) = -parameters.dimension_y / 2;
2059 * p1(2) = -parameters.dimension_y / 2;
2062 * p2(0) = parameters.dimension_x / 2;
2063 * p2(1) = parameters.dimension_y / 2;
2066 * p2(2) = parameters.dimension_y / 2;
2068 * std::vector<unsigned int> divisions(dim);
2069 * divisions[0] =
int(parameters.dimension_x / parameters.dimension_y);
2083 * setup_quadrature_cache();
2085 * initialize_probe_positions_vector();
2087 * frequency_sweep();
2098 * <a name=
"Themainfunction"></a>
2099 * <h4>The main
function</h4>
2103 * The main
function is very similar to the
one in @ref step_40
"step-40".
2106 *
int main(
int argc,
char *argv[])
2110 *
using namespace dealii;
2111 *
const unsigned int dim = 2;
2123 * Displacement simulation. The parameters are read from the
2124 * displacement
HDF5 group and the results are saved in the same
HDF5
2128 *
auto displacement = data.
open_group(
"displacement");
2129 * step62::Parameters<dim> parameters(displacement);
2131 * step62::ElasticWave<dim> elastic_problem(parameters);
2132 * elastic_problem.run();
2138 * Calibration simulation. The parameters are read from the calibration
2139 *
HDF5 group and the results are saved in the same
HDF5 group.
2142 *
auto calibration = data.
open_group(
"calibration");
2143 * step62::Parameters<dim> parameters(calibration);
2145 * step62::ElasticWave<dim> elastic_problem(parameters);
2146 * elastic_problem.run();
2149 *
catch (std::exception &exc)
2151 * std::cerr << std::endl
2153 * <<
"----------------------------------------------------"
2155 * std::cerr <<
"Exception on processing: " << std::endl
2156 * << exc.what() << std::endl
2157 * <<
"Aborting!" << std::endl
2158 * <<
"----------------------------------------------------"
2165 * std::cerr << std::endl
2167 * <<
"----------------------------------------------------"
2169 * std::cerr <<
"Unknown exception!" << std::endl
2170 * <<
"Aborting!" << std::endl
2171 * <<
"----------------------------------------------------"
2179 <a name=
"Results"></a><h1>Results</h1>
2182 <a name=
"Resonancefrequencyandbandgap"></a><h3>Resonance frequency and bandgap</h3>
2185 The results are analyzed in the
2186 [jupyter notebook](https:
2187 with the following code
2189 h5_file = h5py.File(
'results.h5',
'r')
2190 data = h5_file[
'data']
2192 # Gaussian function that we use to fit the resonance
2193 def resonance_f(freq, freq_m, quality_factor, max_amplitude):
2194 omega = 2 * constants.pi * freq
2195 omega_m = 2 * constants.pi * freq_m
2196 gamma = omega_m / quality_factor
2197 return max_amplitude * omega_m**2 *
gamma**2 / (((omega_m**2 - omega**2)**2 +
gamma**2 * omega**2))
2199 frequency = data[
'displacement'][
'frequency'][...]
2200 # Average the probe points
2201 displacement = np.
mean(data[
'displacement'][
'displacement'], axis=0)
2202 calibration_displacement = np.
mean(data[
'calibration'][
'displacement'], axis=0)
2203 reflection_coefficient = displacement / calibration_displacement
2204 reflectivity = (np.
abs(np.
mean(data[
'displacement'][
'displacement'][...]**2, axis=0))/
2205 np.
abs(np.
mean(data[
'calibration'][
'displacement'][...]**2, axis=0)))
2209 y_data = reflectivity
2210 quality_factor_guess = 1e3
2211 freq_guess = x_data[np.argmax(y_data)]
2212 amplitude_guess = np.
max(y_data)
2213 fit_result, covariance = scipy.optimize.curve_fit(resonance_f, x_data, y_data,
2214 [freq_guess, quality_factor_guess, amplitude_guess])
2215 freq_m = fit_result[0]
2216 quality_factor = np.
abs(fit_result[1])
2217 max_amplitude = fit_result[2]
2218 y_data_fit = resonance_f(x_data, freq_m, quality_factor, max_amplitude)
2221 plt.plot(frequency / 1e9, reflectivity, frequency / 1e9, y_data_fit)
2222 plt.xlabel(
'frequency (GHz)')
2223 plt.ylabel(
'amplitude^2 (a.u.)')
2224 plt.title(
'Transmission\n' +
'freq = ' +
"%.7g" % (freq_guess / 1e9) +
'GHz Q = ' +
"%.6g" % quality_factor)
2227 plt.plot(frequency / 1e9, reflectivity)
2228 plt.xlabel(
'frequency (GHz)')
2229 plt.ylabel(
'amplitude^2 (a.u.)')
2230 plt.title(
'Transmission')
2233 plt.plot(frequency / 1e9, np.
angle(reflection_coefficient))
2234 plt.xlabel(
'frequency (GHz)')
2235 plt.ylabel(
'phase (rad)')
2236 plt.title(
'Phase (transmission coefficient)\n')
2242 A phononic cavity is characterized by the
2243 [resonance frequency](https:
2244 [the quality factor](https:
2245 The quality factor is
equal to the ratio between the stored energy in the resonator and the energy
2246 dissipated energy per cycle, which is approximately equivalent to the ratio between the
2247 resonance frequency and the
2248 [full width at half maximum (FWHM)](https:
2249 The FWHM is
equal to the bandwidth over which the power of vibration is greater than half the
2250 power at the resonant frequency.
2252 Q = \frac{f_r}{\Delta f} = \frac{\omega_r}{\Delta \omega} =
2253 2 \pi \times \frac{\text{energy stored}}{\text{energy dissipated per cycle}}
2256 The square of the amplitude of the mechanical resonance @f$a^2@f$ as a
function of the frequency
2257 has a gaussian shape
2259 a^2 = a_\textrm{
max}^2\frac{\omega^2\Gamma^2}{(\omega_r^2-\omega^2)^2+\Gamma^2\omega^2}
2261 where @f$f_r = \frac{\omega_r}{2\pi}@f$ is the resonance frequency and @f$\Gamma=\frac{\omega_r}{Q}@f$ is the dissipation rate.
2262 We used the previous equation in the jupyter notebook to fit the mechanical resonance.
2264 Given the values we have chosen
for the parameters,
one could estimate the resonance frequency
2265 analytically. Indeed,
this is then confirmed by what we get in
this program:
2266 the phononic superlattice cavity exhibits a mechanical resonance at 20GHz and a quality factor of 5046.
2267 The following images show the transmission amplitude and phase as a
function of frequency in the
2268 vicinity of the resonance frequency:
2270 <img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.05.png" height=
"400" />
2271 <img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.06.png" height=
"400" />
2273 The images above suggest that the periodic structure has its intended effect: It really only lets waves of a very
2274 specific frequency pass through, whereas all other waves are reflected. This is of course precisely what
one builds
2275 these sorts of devices
for.
2276 But it is not quite
this easy. In practice, there is really only a
"band gap", i.e., the device blocks waves other than
2277 the desired
one at 20GHz only within a certain frequency range. Indeed, to find out how large
this "gap" is within
2278 which waves are blocked, we can extend the frequency range to 16 GHz through the appropriate parameters in the
2279 input file. We then obtain the following image:
2281 <img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.07.png" height=
"400" />
2283 What
this image suggests is that in the range of around 18 to around 22 GHz, really only the waves with a frequency
2284 of 20 GHz are allowed to pass through, but beyond
this range, there are plenty of other frequencies that can pass
2287 <a name=
"Modeprofile"></a><h3>Mode profile</h3>
2290 We can inspect the mode profile with Paraview or VisIt.
2291 As we have discussed, at resonance all the mechanical
2292 energy is transmitted and the amplitude of motion is amplified inside the cavity.
2293 It can be observed that the PMLs are quite effective to truncate the solution.
2294 The following image shows the mode profile at resonance:
2296 <img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.08.png" height=
"400" />
2298 On the other hand, out of resonance all the mechanical energy is
2299 reflected. The following image shows the profile at 19.75 GHz.
2300 Note the interference between the force pulse and the reflected wave
2301 at the position @f$x=-8\mu\textrm{m}@f$.
2303 <img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.09.png" height=
"400" />
2305 <a name=
"Experimentalapplications"></a><h3>Experimental applications</h3>
2308 Phononic superlattice cavities find application in
2309 [quantum optomechanics](https:
2310 Here we have presented the simulation of a 2D superlattice cavity,
2311 but
this code can be used as well to simulate
"real world" 3D devices such as
2312 [micropillar superlattice cavities](https:
2313 which are promising candidates to study macroscopic quantum phenomena.
2314 The 20GHz mode of a micropillar superlattice cavity is essentially a mechanical harmonic oscillator that is very well isolated
2315 from the environment. If the device is cooled down to 20mK in a dilution fridge, the mode would then become a
2316 macroscopic quantum harmonic oscillator.
2319 <a name=
"PlainProg"></a>
2320 <h1> The plain program</h1>
2321 @include
"step-62.cc"