Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-62.h
Go to the documentation of this file.
1) const
1018 *   {
1019 * @endcode
1020 *
1021 * The speed of sound is defined by
1022 * @f[
1023 * c = \frac{K_e}{\rho}
1024 * @f]
1025 * where @f$K_e@f$ is the effective elastic constant and @f$\rho@f$ the density.
1026 * Here we consider the case in which the waveguide width is much smaller
1027 * than the wavelength. In this case it can be shown that for the two
1028 * dimensional case
1029 * @f[
1030 * K_e = 4\mu\frac{\lambda +\mu}{\lambda+2\mu}
1031 * @f]
1032 * and for the three dimensional case @f$K_e@f$ is equal to the Young's modulus.
1033 * @f[
1034 * K_e = \mu\frac{3\lambda +2\mu}{\lambda+\mu}
1035 * @f]
1036 *
1037 * @code
1038 *   double elastic_constant;
1039 *   if (dim == 2)
1040 *   {
1041 *   elastic_constant = 4 * mu * (lambda + mu) / (lambda + 2 * mu);
1042 *   }
1043 *   else if (dim == 3)
1044 *   {
1045 *   elastic_constant = mu * (3 * lambda + 2 * mu) / (lambda + mu);
1046 *   }
1047 *   else
1048 *   {
1049 *   Assert(false, ExcInternalError());
1050 *   }
1051 *   const double material_a_speed_of_sound =
1052 *   std::sqrt(elastic_constant / material_a_rho);
1053 *   const double material_a_wavelength =
1054 *   material_a_speed_of_sound / cavity_resonance_frequency;
1055 *   const double material_b_speed_of_sound =
1056 *   std::sqrt(elastic_constant / material_b_rho);
1057 *   const double material_b_wavelength =
1058 *   material_b_speed_of_sound / cavity_resonance_frequency;
1059 *  
1060 * @endcode
1061 *
1062 * The density @f$\rho@f$ takes the following form
1063 * <img alt="Phononic superlattice cavity"
1064 * src="https://www.dealii.org/images/steps/developer/step-62.04.svg"
1065 * height="200" />
1066 * where the brown color represents material_a and the green color
1067 * represents material_b.
1068 *
1069 * @code
1070 *   for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1071 *   {
1072 *   const double layer_transition_center =
1073 *   material_a_wavelength / 2 +
1074 *   idx * (material_b_wavelength / 4 + material_a_wavelength / 4);
1075 *   if (std::abs(p[0]) >=
1076 *   (layer_transition_center - average_rho_width / 2) &&
1077 *   std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1078 *   {
1079 *   const double coefficient =
1080 *   (std::abs(p[0]) -
1081 *   (layer_transition_center - average_rho_width / 2)) /
1082 *   average_rho_width;
1083 *   return (1 - coefficient) * material_a_rho +
1084 *   coefficient * material_b_rho;
1085 *   }
1086 *   }
1087 *  
1088 * @endcode
1089 *
1090 * Here we define the
1091 * [subpixel
1092 * smoothing](https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/)
1093 * which improves the precision of the simulation.
1094 *
1095 * @code
1096 *   for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1097 *   {
1098 *   const double layer_transition_center =
1099 *   material_a_wavelength / 2 +
1100 *   idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1101 *   material_b_wavelength / 4;
1102 *   if (std::abs(p[0]) >=
1103 *   (layer_transition_center - average_rho_width / 2) &&
1104 *   std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1105 *   {
1106 *   const double coefficient =
1107 *   (std::abs(p[0]) -
1108 *   (layer_transition_center - average_rho_width / 2)) /
1109 *   average_rho_width;
1110 *   return (1 - coefficient) * material_b_rho +
1111 *   coefficient * material_a_rho;
1112 *   }
1113 *   }
1114 *  
1115 * @endcode
1116 *
1117 * then the cavity
1118 *
1119 * @code
1120 *   if (std::abs(p[0]) <= material_a_wavelength / 2)
1121 *   {
1122 *   return material_a_rho;
1123 *   }
1124 *  
1125 * @endcode
1126 *
1127 * the material_a layers
1128 *
1129 * @code
1130 *   for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1131 *   {
1132 *   const double layer_center =
1133 *   material_a_wavelength / 2 +
1134 *   idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1135 *   material_b_wavelength / 4 + material_a_wavelength / 8;
1136 *   const double layer_width = material_a_wavelength / 4;
1137 *   if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1138 *   std::abs(p[0]) <= (layer_center + layer_width / 2))
1139 *   {
1140 *   return material_a_rho;
1141 *   }
1142 *   }
1143 *  
1144 * @endcode
1145 *
1146 * the material_b layers
1147 *
1148 * @code
1149 *   for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1150 *   {
1151 *   const double layer_center =
1152 *   material_a_wavelength / 2 +
1153 *   idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1154 *   material_b_wavelength / 8;
1155 *   const double layer_width = material_b_wavelength / 4;
1156 *   if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1157 *   std::abs(p[0]) <= (layer_center + layer_width / 2))
1158 *   {
1159 *   return material_b_rho;
1160 *   }
1161 *   }
1162 *  
1163 * @endcode
1164 *
1165 * and finally the default is material_a.
1166 *
1167 * @code
1168 *   return material_a_rho;
1169 *   }
1170 *  
1171 *  
1172 *  
1173 * @endcode
1174 *
1175 *
1176 * <a name="TheParametersclassimplementation"></a>
1177 * <h4>The `Parameters` class implementation</h4>
1178 *
1179
1180 *
1181 * The constructor reads all the parameters from the HDF5::Group `data` using
1182 * the HDF5::Group::get_attribute() function.
1183 *
1184 * @code
1185 *   template <int dim>
1186 *   Parameters<dim>::Parameters(HDF5::Group &data)
1187 *   : data(data)
1188 *   , simulation_name(data.get_attribute<std::string>("simulation_name"))
1189 *   , save_vtu_files(data.get_attribute<bool>("save_vtu_files"))
1190 *   , start_frequency(data.get_attribute<double>("start_frequency"))
1191 *   , stop_frequency(data.get_attribute<double>("stop_frequency"))
1192 *   , nb_frequency_points(data.get_attribute<int>("nb_frequency_points"))
1193 *   , lambda(data.get_attribute<double>("lambda"))
1194 *   , mu(data.get_attribute<double>("mu"))
1195 *   , dimension_x(data.get_attribute<double>("dimension_x"))
1196 *   , dimension_y(data.get_attribute<double>("dimension_y"))
1197 *   , nb_probe_points(data.get_attribute<int>("nb_probe_points"))
1198 *   , grid_level(data.get_attribute<int>("grid_level"))
1199 *   , probe_start_point(data.get_attribute<double>("probe_pos_x"),
1200 *   data.get_attribute<double>("probe_pos_y") -
1201 *   data.get_attribute<double>("probe_width_y") / 2)
1202 *   , probe_stop_point(data.get_attribute<double>("probe_pos_x"),
1203 *   data.get_attribute<double>("probe_pos_y") +
1204 *   data.get_attribute<double>("probe_width_y") / 2)
1205 *   , right_hand_side(data)
1206 *   , pml(data)
1207 *   , rho(data)
1208 *   {}
1209 *  
1210 *  
1211 *  
1212 * @endcode
1213 *
1214 *
1215 * <a name="TheQuadratureCacheclassimplementation"></a>
1216 * <h4>The `QuadratureCache` class implementation</h4>
1217 *
1218
1219 *
1220 * We need to reserve enough space for the mass and stiffness matrices and the
1221 * right hand side vector.
1222 *
1223 * @code
1224 *   template <int dim>
1225 *   QuadratureCache<dim>::QuadratureCache(const unsigned int dofs_per_cell)
1226 *   : dofs_per_cell(dofs_per_cell)
1227 *   , mass_coefficient(dofs_per_cell, dofs_per_cell)
1228 *   , stiffness_coefficient(dofs_per_cell, dofs_per_cell)
1229 *   , right_hand_side(dofs_per_cell)
1230 *   {}
1231 *  
1232 *  
1233 *  
1234 * @endcode
1235 *
1236 *
1237 * <a name="ImplementationoftheElasticWaveclass"></a>
1238 * <h3>Implementation of the `ElasticWave` class</h3>
1239 *
1240
1241 *
1242 *
1243 * <a name="Constructor"></a>
1244 * <h4>Constructor</h4>
1245 *
1246
1247 *
1248 * This is very similar to the constructor of @ref step_40 "step-40". In addition we create
1249 * the HDF5 datasets `frequency_dataset`, `position_dataset` and
1250 * `displacement`. Note the use of the `template` keyword for the creation of
1251 * the HDF5 datasets. It is a C++ requirement to use the `template` keyword in
1252 * order to treat `create_dataset` as a dependent template name.
1253 *
1254 * @code
1255 *   template <int dim>
1256 *   ElasticWave<dim>::ElasticWave(const Parameters<dim> &parameters)
1257 *   : parameters(parameters)
1258 *   , mpi_communicator(MPI_COMM_WORLD)
1259 *   , triangulation(mpi_communicator,
1260 *   typename Triangulation<dim>::MeshSmoothing(
1261 *   Triangulation<dim>::smoothing_on_refinement |
1262 *   Triangulation<dim>::smoothing_on_coarsening))
1263 *   , quadrature_formula(2)
1264 *   , fe(FE_Q<dim>(1), dim)
1265 *   , dof_handler(triangulation)
1266 *   , frequency(parameters.nb_frequency_points)
1267 *   , probe_positions(parameters.nb_probe_points, dim)
1268 *   , frequency_dataset(parameters.data.template create_dataset<double>(
1269 *   "frequency",
1270 *   std::vector<hsize_t>{parameters.nb_frequency_points}))
1271 *   , probe_positions_dataset(parameters.data.template create_dataset<double>(
1272 *   "position",
1273 *   std::vector<hsize_t>{parameters.nb_probe_points, dim}))
1274 *   , displacement(
1275 *   parameters.data.template create_dataset<std::complex<double>>(
1276 *   "displacement",
1277 *   std::vector<hsize_t>{parameters.nb_probe_points,
1278 *   parameters.nb_frequency_points}))
1279 *   , pcout(std::cout,
1280 *   (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
1281 *   , computing_timer(mpi_communicator,
1282 *   pcout,
1283 *   TimerOutput::never,
1284 *   TimerOutput::wall_times)
1285 *   {}
1286 *  
1287 *  
1288 *  
1289 * @endcode
1290 *
1291 *
1292 * <a name="ElasticWavesetup_system"></a>
1293 * <h4>ElasticWave::setup_system</h4>
1294 *
1295
1296 *
1297 * There is nothing new in this function, the only difference with @ref step_40 "step-40" is
1298 * that we don't have to apply boundary conditions because we use the PMLs to
1299 * truncate the domain.
1300 *
1301 * @code
1302 *   template <int dim>
1303 *   void ElasticWave<dim>::setup_system()
1304 *   {
1305 *   TimerOutput::Scope t(computing_timer, "setup");
1306 *  
1307 *   dof_handler.distribute_dofs(fe);
1308 *  
1309 *   locally_owned_dofs = dof_handler.locally_owned_dofs();
1310 *   locally_relevant_dofs =
1312 *  
1313 *   locally_relevant_solution.reinit(locally_owned_dofs,
1314 *   locally_relevant_dofs,
1315 *   mpi_communicator);
1316 *  
1317 *   system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1318 *  
1319 *   constraints.clear();
1320 *   constraints.reinit(locally_relevant_dofs);
1321 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1322 *  
1323 *   constraints.close();
1324 *  
1325 *   DynamicSparsityPattern dsp(locally_relevant_dofs);
1326 *  
1327 *   DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
1329 *   locally_owned_dofs,
1330 *   mpi_communicator,
1331 *   locally_relevant_dofs);
1332 *  
1333 *   system_matrix.reinit(locally_owned_dofs,
1334 *   locally_owned_dofs,
1335 *   dsp,
1336 *   mpi_communicator);
1337 *   }
1338 *  
1339 *  
1340 *  
1341 * @endcode
1342 *
1343 *
1344 * <a name="ElasticWaveassemble_system"></a>
1345 * <h4>ElasticWave::assemble_system</h4>
1346 *
1347
1348 *
1349 * This function is also very similar to @ref step_40 "step-40", though there are notable
1350 * differences. We assemble the system for each frequency/omega step. In the
1351 * first step we set `calculate_quadrature_data = True` and we calculate the
1352 * mass and stiffness matrices and the right hand side vector. In the
1353 * subsequent steps we will use that data to accelerate the calculation.
1354 *
1355 * @code
1356 *   template <int dim>
1357 *   void ElasticWave<dim>::assemble_system(const double omega,
1358 *   const bool calculate_quadrature_data)
1359 *   {
1360 *   TimerOutput::Scope t(computing_timer, "assembly");
1361 *  
1362 *   FEValues<dim> fe_values(fe,
1363 *   quadrature_formula,
1366 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1367 *   const unsigned int n_q_points = quadrature_formula.size();
1368 *  
1369 *   FullMatrix<std::complex<double>> cell_matrix(dofs_per_cell, dofs_per_cell);
1370 *   Vector<std::complex<double>> cell_rhs(dofs_per_cell);
1371 *  
1372 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1373 *  
1374 * @endcode
1375 *
1376 * Here we store the value of the right hand side, rho and the PML.
1377 *
1378 * @code
1379 *   std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim));
1380 *   std::vector<double> rho_values(n_q_points);
1381 *   std::vector<Vector<std::complex<double>>> pml_values(
1382 *   n_q_points, Vector<std::complex<double>>(dim));
1383 *  
1384 * @endcode
1385 *
1386 * We calculate the stiffness tensor for the @f$\lambda@f$ and @f$\mu@f$ that have
1387 * been defined in the jupyter notebook. Note that contrary to @f$\rho@f$ the
1388 * stiffness is constant among for the whole domain.
1389 *
1390 * @code
1391 *   const SymmetricTensor<4, dim> stiffness_tensor =
1392 *   get_stiffness_tensor<dim>(parameters.lambda, parameters.mu);
1393 *  
1394 * @endcode
1395 *
1396 * We use the same method of @ref step_20 "step-20" for vector-valued problems.
1397 *
1398 * @code
1399 *   const FEValuesExtractors::Vector displacement(0);
1400 *  
1401 *   for (const auto &cell : dof_handler.active_cell_iterators())
1402 *   if (cell->is_locally_owned())
1403 *   {
1404 *   cell_matrix = 0;
1405 *   cell_rhs = 0;
1406 *  
1407 * @endcode
1408 *
1409 * We have to calculate the values of the right hand side, rho and
1410 * the PML only if we are going to calculate the mass and the
1411 * stiffness matrices. Otherwise we can skip this calculation which
1412 * considerably reduces the total calculation time.
1413 *
1414 * @code
1415 *   if (calculate_quadrature_data)
1416 *   {
1417 *   fe_values.reinit(cell);
1418 *  
1419 *   parameters.right_hand_side.vector_value_list(
1420 *   fe_values.get_quadrature_points(), rhs_values);
1421 *   parameters.rho.value_list(fe_values.get_quadrature_points(),
1422 *   rho_values);
1423 *   parameters.pml.vector_value_list(
1424 *   fe_values.get_quadrature_points(), pml_values);
1425 *   }
1426 *  
1427 * @endcode
1428 *
1429 * We have done this in @ref step_18 "step-18". Get a pointer to the quadrature
1430 * cache data local to the present cell, and, as a defensive
1431 * measure, make sure that this pointer is within the bounds of the
1432 * global array:
1433 *
1434 * @code
1435 *   QuadratureCache<dim> *local_quadrature_points_data =
1436 *   reinterpret_cast<QuadratureCache<dim> *>(cell->user_pointer());
1437 *   Assert(local_quadrature_points_data >= &quadrature_cache.front(),
1438 *   ExcInternalError());
1439 *   Assert(local_quadrature_points_data <= &quadrature_cache.back(),
1440 *   ExcInternalError());
1441 *   for (unsigned int q = 0; q < n_q_points; ++q)
1442 *   {
1443 * @endcode
1444 *
1445 * The quadrature_data variable is used to store the mass and
1446 * stiffness matrices, the right hand side vector and the value
1447 * of `JxW`.
1448 *
1449 * @code
1450 *   QuadratureCache<dim> &quadrature_data =
1451 *   local_quadrature_points_data[q];
1452 *  
1453 * @endcode
1454 *
1455 * Below we declare the force vector and the parameters of the
1456 * PML @f$s@f$ and @f$\xi@f$.
1457 *
1458 * @code
1459 *   Tensor<1, dim> force;
1461 *   std::complex<double> xi(1, 0);
1462 *  
1463 * @endcode
1464 *
1465 * The following block is calculated only in the first frequency
1466 * step.
1467 *
1468 * @code
1469 *   if (calculate_quadrature_data)
1470 *   {
1471 * @endcode
1472 *
1473 * Store the value of `JxW`.
1474 *
1475 * @code
1476 *   quadrature_data.JxW = fe_values.JxW(q);
1477 *  
1478 *   for (unsigned int component = 0; component < dim; ++component)
1479 *   {
1480 * @endcode
1481 *
1482 * Convert vectors to tensors and calculate xi
1483 *
1484 * @code
1485 *   force[component] = rhs_values[q][component];
1486 *   s[component] = pml_values[q][component];
1487 *   xi *= s[component];
1488 *   }
1489 *  
1490 * @endcode
1491 *
1492 * Here we calculate the @f$\alpha_{mnkl}@f$ and @f$\beta_{mnkl}@f$
1493 * tensors.
1494 *
1495 * @code
1498 *   for (unsigned int m = 0; m < dim; ++m)
1499 *   for (unsigned int n = 0; n < dim; ++n)
1500 *   for (unsigned int k = 0; k < dim; ++k)
1501 *   for (unsigned int l = 0; l < dim; ++l)
1502 *   {
1503 *   alpha[m][n][k][l] = xi *
1504 *   stiffness_tensor[m][n][k][l] /
1505 *   (2.0 * s[n] * s[k]);
1506 *   beta[m][n][k][l] = xi *
1507 *   stiffness_tensor[m][n][k][l] /
1508 *   (2.0 * s[n] * s[l]);
1509 *   }
1510 *  
1511 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1512 *   {
1513 *   const Tensor<1, dim> phi_i =
1514 *   fe_values[displacement].value(i, q);
1515 *   const Tensor<2, dim> grad_phi_i =
1516 *   fe_values[displacement].gradient(i, q);
1517 *  
1518 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1519 *   {
1520 *   const Tensor<1, dim> phi_j =
1521 *   fe_values[displacement].value(j, q);
1522 *   const Tensor<2, dim> grad_phi_j =
1523 *   fe_values[displacement].gradient(j, q);
1524 *  
1525 * @endcode
1526 *
1527 * calculate the values of the @ref GlossMassMatrix "mass matrix".
1528 *
1529 * @code
1530 *   quadrature_data.mass_coefficient[i][j] =
1531 *   rho_values[q] * xi * phi_i * phi_j;
1532 *  
1533 * @endcode
1534 *
1535 * Loop over the @f$mnkl@f$ indices of the stiffness
1536 * tensor.
1537 *
1538 * @code
1539 *   std::complex<double> stiffness_coefficient = 0;
1540 *   for (unsigned int m = 0; m < dim; ++m)
1541 *   for (unsigned int n = 0; n < dim; ++n)
1542 *   for (unsigned int k = 0; k < dim; ++k)
1543 *   for (unsigned int l = 0; l < dim; ++l)
1544 *   {
1545 * @endcode
1546 *
1547 * Here we calculate the stiffness matrix.
1548 * Note that the stiffness matrix is not
1549 * symmetric because of the PMLs. We use the
1550 * gradient function (see the
1551 * [documentation](https://www.dealii.org/current/doxygen/deal.II/group__vector__valued.html))
1552 * which is a <code>Tensor@<2,dim@></code>.
1553 * The matrix @f$G_{ij}@f$ consists of entries
1554 * @f[
1555 * G_{ij}=
1556 * \frac{\partial\phi_i}{\partial x_j}
1557 * =\partial_j \phi_i
1558 * @f]
1559 * Note the position of the indices @f$i@f$ and
1560 * @f$j@f$ and the notation that we use in this
1561 * tutorial: @f$\partial_j\phi_i@f$. As the
1562 * stiffness tensor is not symmetric, it is
1563 * very easy to make a mistake.
1564 *
1565 * @code
1566 *   stiffness_coefficient +=
1567 *   grad_phi_i[m][n] *
1568 *   (alpha[m][n][k][l] * grad_phi_j[l][k] +
1569 *   beta[m][n][k][l] * grad_phi_j[k][l]);
1570 *   }
1571 *  
1572 * @endcode
1573 *
1574 * We save the value of the stiffness matrix in
1575 * quadrature_data
1576 *
1577 * @code
1578 *   quadrature_data.stiffness_coefficient[i][j] =
1579 *   stiffness_coefficient;
1580 *   }
1581 *  
1582 * @endcode
1583 *
1584 * and the value of the right hand side in
1585 * quadrature_data.
1586 *
1587 * @code
1588 *   quadrature_data.right_hand_side[i] =
1589 *   phi_i * force * fe_values.JxW(q);
1590 *   }
1591 *   }
1592 *  
1593 * @endcode
1594 *
1595 * We loop again over the degrees of freedom of the cells to
1596 * calculate the system matrix. These loops are really quick
1597 * because we have already calculated the stiffness and mass
1598 * matrices, only the value of @f$\omega@f$ changes.
1599 *
1600 * @code
1601 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1602 *   {
1603 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1604 *   {
1605 *   std::complex<double> matrix_sum = 0;
1606 *   matrix_sum += -std::pow(omega, 2) *
1607 *   quadrature_data.mass_coefficient[i][j];
1608 *   matrix_sum += quadrature_data.stiffness_coefficient[i][j];
1609 *   cell_matrix(i, j) += matrix_sum * quadrature_data.JxW;
1610 *   }
1611 *   cell_rhs(i) += quadrature_data.right_hand_side[i];
1612 *   }
1613 *   }
1614 *   cell->get_dof_indices(local_dof_indices);
1615 *   constraints.distribute_local_to_global(cell_matrix,
1616 *   cell_rhs,
1617 *   local_dof_indices,
1618 *   system_matrix,
1619 *   system_rhs);
1620 *   }
1621 *  
1622 *   system_matrix.compress(VectorOperation::add);
1623 *   system_rhs.compress(VectorOperation::add);
1624 *   }
1625 *  
1626 * @endcode
1627 *
1628 *
1629 * <a name="ElasticWavesolve"></a>
1630 * <h4>ElasticWave::solve</h4>
1631 *
1632
1633 *
1634 * This is even more simple than in @ref step_40 "step-40". We use the parallel direct solver
1635 * MUMPS which requires less options than an iterative solver. The drawback is
1636 * that it does not scale very well. It is not straightforward to solve the
1637 * Helmholtz equation with an iterative solver. The shifted Laplacian
1638 * multigrid method is a well known approach to precondition this system, but
1639 * this is beyond the scope of this tutorial.
1640 *
1641 * @code
1642 *   template <int dim>
1643 *   void ElasticWave<dim>::solve()
1644 *   {
1645 *   TimerOutput::Scope t(computing_timer, "solve");
1646 *   LinearAlgebraPETSc::MPI::Vector completely_distributed_solution(
1647 *   locally_owned_dofs, mpi_communicator);
1648 *  
1649 *   SolverControl solver_control;
1650 *   PETScWrappers::SparseDirectMUMPS solver(solver_control, mpi_communicator);
1651 *   solver.solve(system_matrix, completely_distributed_solution, system_rhs);
1652 *  
1653 *   pcout << " Solved in " << solver_control.last_step() << " iterations."
1654 *   << std::endl;
1655 *   constraints.distribute(completely_distributed_solution);
1656 *   locally_relevant_solution = completely_distributed_solution;
1657 *   }
1658 *  
1659 * @endcode
1660 *
1661 *
1662 * <a name="ElasticWaveinitialize_position_vector"></a>
1663 * <h4>ElasticWave::initialize_position_vector</h4>
1664 *
1665
1666 *
1667 * We use this function to calculate the values of the position vector.
1668 *
1669 * @code
1670 *   template <int dim>
1671 *   void ElasticWave<dim>::initialize_probe_positions_vector()
1672 *   {
1673 *   for (unsigned int position_idx = 0;
1674 *   position_idx < parameters.nb_probe_points;
1675 *   ++position_idx)
1676 *   {
1677 * @endcode
1678 *
1679 * Because of the way the operator + and - are overloaded to subtract
1680 * two points, the following has to be done:
1681 * `Point_b<dim> + (-Point_a<dim>)`
1682 *
1683 * @code
1684 *   const Point<dim> p =
1685 *   (position_idx / ((double)(parameters.nb_probe_points - 1))) *
1686 *   (parameters.probe_stop_point + (-parameters.probe_start_point)) +
1687 *   parameters.probe_start_point;
1688 *   probe_positions[position_idx][0] = p[0];
1689 *   probe_positions[position_idx][1] = p[1];
1690 *   if (dim == 3)
1691 *   {
1692 *   probe_positions[position_idx][2] = p[2];
1693 *   }
1694 *   }
1695 *   }
1696 *  
1697 * @endcode
1698 *
1699 *
1700 * <a name="ElasticWavestore_frequency_step_data"></a>
1701 * <h4>ElasticWave::store_frequency_step_data</h4>
1702 *
1703
1704 *
1705 * This function stores in the HDF5 file the measured energy by the probe.
1706 *
1707 * @code
1708 *   template <int dim>
1709 *   void
1710 *   ElasticWave<dim>::store_frequency_step_data(const unsigned int frequency_idx)
1711 *   {
1712 *   TimerOutput::Scope t(computing_timer, "store_frequency_step_data");
1713 *  
1714 * @endcode
1715 *
1716 * We store the displacement in the @f$x@f$ direction; the displacement in the
1717 * @f$y@f$ direction is negligible.
1718 *
1719 * @code
1720 *   const unsigned int probe_displacement_component = 0;
1721 *  
1722 * @endcode
1723 *
1724 * The vector coordinates contains the coordinates in the HDF5 file of the
1725 * points of the probe that are located in locally owned cells. The vector
1726 * displacement_data contains the value of the displacement at these points.
1727 *
1728 * @code
1729 *   std::vector<hsize_t> coordinates;
1730 *   std::vector<std::complex<double>> displacement_data;
1731 *  
1732 *   const auto &mapping = get_default_linear_mapping(triangulation);
1733 *   GridTools::Cache<dim, dim> cache(triangulation, mapping);
1734 *   typename Triangulation<dim, dim>::active_cell_iterator cell_hint{};
1735 *   std::vector<bool> marked_vertices = {};
1736 *   const double tolerance = 1.e-10;
1737 *  
1738 *   for (unsigned int position_idx = 0;
1739 *   position_idx < parameters.nb_probe_points;
1740 *   ++position_idx)
1741 *   {
1742 *   Point<dim> point;
1743 *   for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1744 *   {
1745 *   point[dim_idx] = probe_positions[position_idx][dim_idx];
1746 *   }
1747 *   bool point_in_locally_owned_cell = false;
1748 *   {
1749 *   auto cell_and_ref_point = GridTools::find_active_cell_around_point(
1750 *   cache, point, cell_hint, marked_vertices, tolerance);
1751 *   if (cell_and_ref_point.first.state() == IteratorState::valid)
1752 *   {
1753 *   cell_hint = cell_and_ref_point.first;
1754 *   point_in_locally_owned_cell =
1755 *   cell_and_ref_point.first->is_locally_owned();
1756 *   }
1757 *   }
1758 *   if (point_in_locally_owned_cell)
1759 *   {
1760 * @endcode
1761 *
1762 * Then we can store the values of the displacement in the points of
1763 * the probe in `displacement_data`.
1764 *
1765 * @code
1766 *   Vector<std::complex<double>> tmp_vector(dim);
1767 *   VectorTools::point_value(dof_handler,
1768 *   locally_relevant_solution,
1769 *   point,
1770 *   tmp_vector);
1771 *   coordinates.emplace_back(position_idx);
1772 *   coordinates.emplace_back(frequency_idx);
1773 *   displacement_data.emplace_back(
1774 *   tmp_vector(probe_displacement_component));
1775 *   }
1776 *   }
1777 *  
1778 * @endcode
1779 *
1780 * We write the displacement data in the HDF5 file. The call
1781 * HDF5::DataSet::write_selection() is MPI collective which means that all
1782 * the processes have to participate.
1783 *
1784 * @code
1785 *   if (coordinates.size() > 0)
1786 *   {
1787 *   displacement.write_selection(displacement_data, coordinates);
1788 *   }
1789 * @endcode
1790 *
1791 * Therefore even if the process has no data to write it has to participate
1792 * in the collective call. For this we can use HDF5::DataSet::write_none().
1793 * Note that we have to specify the data type, in this case
1794 * `std::complex<double>`.
1795 *
1796 * @code
1797 *   else
1798 *   {
1799 *   displacement.write_none<std::complex<double>>();
1800 *   }
1801 *  
1802 * @endcode
1803 *
1804 * If the variable `save_vtu_files` in the input file equals `True` then all
1805 * the data will be saved as vtu. The procedure to write `vtu` files has
1806 * been described in @ref step_40 "step-40".
1807 *
1808 * @code
1809 *   if (parameters.save_vtu_files)
1810 *   {
1811 *   std::vector<std::string> solution_names(dim, "displacement");
1812 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
1813 *   interpretation(
1815 *  
1816 *   DataOut<dim> data_out;
1817 *   data_out.add_data_vector(dof_handler,
1818 *   locally_relevant_solution,
1819 *   solution_names,
1820 *   interpretation);
1821 *   Vector<float> subdomain(triangulation.n_active_cells());
1822 *   for (unsigned int i = 0; i < subdomain.size(); ++i)
1823 *   subdomain(i) = triangulation.locally_owned_subdomain();
1824 *   data_out.add_data_vector(subdomain, "subdomain");
1825 *  
1826 *   std::vector<Vector<double>> force(
1827 *   dim, Vector<double>(triangulation.n_active_cells()));
1828 *   std::vector<Vector<double>> pml(
1829 *   dim, Vector<double>(triangulation.n_active_cells()));
1830 *   Vector<double> rho(triangulation.n_active_cells());
1831 *  
1832 *   for (auto &cell : triangulation.active_cell_iterators())
1833 *   {
1834 *   if (cell->is_locally_owned())
1835 *   {
1836 *   for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1837 *   {
1838 *   force[dim_idx](cell->active_cell_index()) =
1839 *   parameters.right_hand_side.value(cell->center(), dim_idx);
1840 *   pml[dim_idx](cell->active_cell_index()) =
1841 *   parameters.pml.value(cell->center(), dim_idx).imag();
1842 *   }
1843 *   rho(cell->active_cell_index()) =
1844 *   parameters.rho.value(cell->center());
1845 *   }
1846 * @endcode
1847 *
1848 * And on the cells that we are not interested in, set the
1849 * respective value to a bogus value in order to make sure that if
1850 * we were somehow wrong about our assumption we would find out by
1851 * looking at the graphical output:
1852 *
1853 * @code
1854 *   else
1855 *   {
1856 *   for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1857 *   {
1858 *   force[dim_idx](cell->active_cell_index()) = -1e+20;
1859 *   pml[dim_idx](cell->active_cell_index()) = -1e+20;
1860 *   }
1861 *   rho(cell->active_cell_index()) = -1e+20;
1862 *   }
1863 *   }
1864 *  
1865 *   for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1866 *   {
1867 *   data_out.add_data_vector(force[dim_idx],
1868 *   "force_" + std::to_string(dim_idx));
1869 *   data_out.add_data_vector(pml[dim_idx],
1870 *   "pml_" + std::to_string(dim_idx));
1871 *   }
1872 *   data_out.add_data_vector(rho, "rho");
1873 *  
1874 *   data_out.build_patches();
1875 *  
1876 *   std::stringstream frequency_idx_stream;
1877 *   const unsigned int nb_number_positions =
1878 *   ((unsigned int)std::log10(parameters.nb_frequency_points)) + 1;
1879 *   frequency_idx_stream << std::setw(nb_number_positions)
1880 *   << std::setfill('0') << frequency_idx;
1881 *   std::string filename = (parameters.simulation_name + "_" +
1882 *   frequency_idx_stream.str() + ".vtu");
1883 *   data_out.write_vtu_in_parallel(filename.c_str(), mpi_communicator);
1884 *   }
1885 *   }
1886 *  
1887 *  
1888 *  
1889 * @endcode
1890 *
1891 *
1892 * <a name="ElasticWaveoutput_results"></a>
1893 * <h4>ElasticWave::output_results</h4>
1894 *
1895
1896 *
1897 * This function writes the datasets that have not already been written.
1898 *
1899 * @code
1900 *   template <int dim>
1901 *   void ElasticWave<dim>::output_results()
1902 *   {
1903 * @endcode
1904 *
1905 * The vectors `frequency` and `position` are the same for all the
1906 * processes. Therefore any of the processes can write the corresponding
1907 * `datasets`. Because the call HDF5::DataSet::write is MPI collective, the
1908 * rest of the processes will have to call HDF5::DataSet::write_none.
1909 *
1910 * @code
1911 *   if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
1912 *   {
1913 *   frequency_dataset.write(frequency);
1914 *   probe_positions_dataset.write(probe_positions);
1915 *   }
1916 *   else
1917 *   {
1918 *   frequency_dataset.write_none<double>();
1919 *   probe_positions_dataset.write_none<double>();
1920 *   }
1921 *   }
1922 *  
1923 *  
1924 *  
1925 * @endcode
1926 *
1927 *
1928 * <a name="ElasticWavesetup_quadrature_cache"></a>
1929 * <h4>ElasticWave::setup_quadrature_cache</h4>
1930 *
1931
1932 *
1933 * We use this function at the beginning of our computations to set up initial
1934 * values of the cache variables. This function has been described in @ref step_18 "step-18".
1935 * There are no differences with the function of @ref step_18 "step-18".
1936 *
1937 * @code
1938 *   template <int dim>
1939 *   void ElasticWave<dim>::setup_quadrature_cache()
1940 *   {
1941 *   triangulation.clear_user_data();
1942 *  
1943 *   {
1944 *   std::vector<QuadratureCache<dim>> tmp;
1945 *   quadrature_cache.swap(tmp);
1946 *   }
1947 *  
1948 *   quadrature_cache.resize(triangulation.n_locally_owned_active_cells() *
1949 *   quadrature_formula.size(),
1950 *   QuadratureCache<dim>(fe.n_dofs_per_cell()));
1951 *   unsigned int cache_index = 0;
1952 *   for (const auto &cell : triangulation.active_cell_iterators())
1953 *   if (cell->is_locally_owned())
1954 *   {
1955 *   cell->set_user_pointer(&quadrature_cache[cache_index]);
1956 *   cache_index += quadrature_formula.size();
1957 *   }
1958 *   Assert(cache_index == quadrature_cache.size(), ExcInternalError());
1959 *   }
1960 *  
1961 *  
1962 *  
1963 * @endcode
1964 *
1965 *
1966 * <a name="ElasticWavefrequency_sweep"></a>
1967 * <h4>ElasticWave::frequency_sweep</h4>
1968 *
1969
1970 *
1971 * For clarity we divide the function `run` of @ref step_40 "step-40" into the functions
1972 * `run` and `frequency_sweep`. In the function `frequency_sweep` we place the
1973 * iteration over the frequency vector.
1974 *
1975 * @code
1976 *   template <int dim>
1977 *   void ElasticWave<dim>::frequency_sweep()
1978 *   {
1979 *   for (unsigned int frequency_idx = 0;
1980 *   frequency_idx < parameters.nb_frequency_points;
1981 *   ++frequency_idx)
1982 *   {
1983 *   pcout << parameters.simulation_name + " frequency idx: "
1984 *   << frequency_idx << '/' << parameters.nb_frequency_points - 1
1985 *   << std::endl;
1986 *  
1987 *  
1988 *  
1989 *   setup_system();
1990 *   if (frequency_idx == 0)
1991 *   {
1992 *   pcout << " Number of active cells : "
1993 *   << triangulation.n_active_cells() << std::endl;
1994 *   pcout << " Number of degrees of freedom : "
1995 *   << dof_handler.n_dofs() << std::endl;
1996 *   }
1997 *  
1998 *   if (frequency_idx == 0)
1999 *   {
2000 * @endcode
2001 *
2002 * Write the simulation parameters only once
2003 *
2004 * @code
2005 *   parameters.data.set_attribute("active_cells",
2006 *   triangulation.n_active_cells());
2007 *   parameters.data.set_attribute("degrees_of_freedom",
2008 *   dof_handler.n_dofs());
2009 *   }
2010 *  
2011 * @endcode
2012 *
2013 * We calculate the frequency and omega values for this particular step.
2014 *
2015 * @code
2016 *   const double current_loop_frequency =
2017 *   (parameters.start_frequency +
2018 *   frequency_idx *
2019 *   (parameters.stop_frequency - parameters.start_frequency) /
2020 *   (parameters.nb_frequency_points - 1));
2021 *   const double current_loop_omega =
2022 *   2 * numbers::PI * current_loop_frequency;
2023 *  
2024 * @endcode
2025 *
2026 * In the first frequency step we calculate the mass and stiffness
2027 * matrices and the right hand side. In the subsequent frequency steps
2028 * we will use those values. This improves considerably the calculation
2029 * time.
2030 *
2031 * @code
2032 *   assemble_system(current_loop_omega,
2033 *   (frequency_idx == 0) ? true : false);
2034 *   solve();
2035 *  
2036 *   frequency[frequency_idx] = current_loop_frequency;
2037 *   store_frequency_step_data(frequency_idx);
2038 *  
2039 *   computing_timer.print_summary();
2040 *   computing_timer.reset();
2041 *   pcout << std::endl;
2042 *   }
2043 *   }
2044 *  
2045 *  
2046 *  
2047 * @endcode
2048 *
2049 *
2050 * <a name="ElasticWaverun"></a>
2051 * <h4>ElasticWave::run</h4>
2052 *
2053
2054 *
2055 * This function is very similar to the one in @ref step_40 "step-40".
2056 *
2057 * @code
2058 *   template <int dim>
2059 *   void ElasticWave<dim>::run()
2060 *   {
2061 *   #ifdef DEBUG
2062 *   pcout << "Debug mode" << std::endl;
2063 *   #else
2064 *   pcout << "Release mode" << std::endl;
2065 *   #endif
2066 *  
2067 *   {
2068 *   Point<dim> p1;
2069 *   p1(0) = -parameters.dimension_x / 2;
2070 *   p1(1) = -parameters.dimension_y / 2;
2071 *   if (dim == 3)
2072 *   {
2073 *   p1(2) = -parameters.dimension_y / 2;
2074 *   }
2075 *   Point<dim> p2;
2076 *   p2(0) = parameters.dimension_x / 2;
2077 *   p2(1) = parameters.dimension_y / 2;
2078 *   if (dim == 3)
2079 *   {
2080 *   p2(2) = parameters.dimension_y / 2;
2081 *   }
2082 *   std::vector<unsigned int> divisions(dim);
2083 *   divisions[0] = int(parameters.dimension_x / parameters.dimension_y);
2084 *   divisions[1] = 1;
2085 *   if (dim == 3)
2086 *   {
2087 *   divisions[2] = 1;
2088 *   }
2090 *   divisions,
2091 *   p1,
2092 *   p2);
2093 *   }
2094 *  
2095 *   triangulation.refine_global(parameters.grid_level);
2096 *  
2097 *   setup_quadrature_cache();
2098 *  
2099 *   initialize_probe_positions_vector();
2100 *  
2101 *   frequency_sweep();
2102 *  
2103 *   output_results();
2104 *   }
2105 *   } // namespace step62
2106 *  
2107 *  
2108 *  
2109 * @endcode
2110 *
2111 *
2112 * <a name="Themainfunction"></a>
2113 * <h4>The main function</h4>
2114 *
2115
2116 *
2117 * The main function is very similar to the one in @ref step_40 "step-40".
2118 *
2119 * @code
2120 *   int main(int argc, char *argv[])
2121 *   {
2122 *   try
2123 *   {
2124 *   using namespace dealii;
2125 *   const unsigned int dim = 2;
2126 *  
2127 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
2128 *  
2129 *   HDF5::File data_file("results.h5",
2131 *   MPI_COMM_WORLD);
2132 *   auto data = data_file.create_group("data");
2133 *  
2134 * @endcode
2135 *
2136 * Each of the simulations (displacement and calibration) is stored in a
2137 * separate HDF5 group:
2138 *
2139 * @code
2140 *   const std::array<std::string, 2> group_names{
2141 *   {"displacement", "calibration"}};
2142 *   for (const std::string &group_name : group_names)
2143 *   {
2144 * @endcode
2145 *
2146 * For each of these two group names, we now create the group and put
2147 * attributes into these groups.
2148 * Specifically, these are:
2149 * - The dimensions of the waveguide (in @f$x@f$ and @f$y@f$ directions)
2150 * - The position of the probe (in @f$x@f$ and @f$y@f$ directions)
2151 * - The number of points in the probe
2152 * - The global refinement level
2153 * - The cavity resonance frequency
2154 * - The number of mirror pairs
2155 * - The material properties
2156 * - The force parameters
2157 * - The PML parameters
2158 * - The frequency parameters
2159 *
2160
2161 *
2162 *
2163 * @code
2164 *   auto group = data.create_group(group_name);
2165 *  
2166 *   group.set_attribute<double>("dimension_x", 2e-5);
2167 *   group.set_attribute<double>("dimension_y", 2e-8);
2168 *   group.set_attribute<double>("probe_pos_x", 8e-6);
2169 *   group.set_attribute<double>("probe_pos_y", 0);
2170 *   group.set_attribute<double>("probe_width_y", 2e-08);
2171 *   group.set_attribute<unsigned int>("nb_probe_points", 5);
2172 *   group.set_attribute<unsigned int>("grid_level", 1);
2173 *   group.set_attribute<double>("cavity_resonance_frequency", 20e9);
2174 *   group.set_attribute<unsigned int>("nb_mirror_pairs", 15);
2175 *  
2176 *   group.set_attribute<double>("poissons_ratio", 0.27);
2177 *   group.set_attribute<double>("youngs_modulus", 270000000000.0);
2178 *   group.set_attribute<double>("material_a_rho", 3200);
2179 *  
2180 *   if (group_name == "displacement")
2181 *   group.set_attribute<double>("material_b_rho", 2000);
2182 *   else
2183 *   group.set_attribute<double>("material_b_rho", 3200);
2184 *  
2185 *   group.set_attribute(
2186 *   "lambda",
2187 *   group.get_attribute<double>("youngs_modulus") *
2188 *   group.get_attribute<double>("poissons_ratio") /
2189 *   ((1 + group.get_attribute<double>("poissons_ratio")) *
2190 *   (1 - 2 * group.get_attribute<double>("poissons_ratio"))));
2191 *   group.set_attribute("mu",
2192 *   group.get_attribute<double>("youngs_modulus") /
2193 *   (2 * (1 + group.get_attribute<double>(
2194 *   "poissons_ratio"))));
2195 *  
2196 *   group.set_attribute<double>("max_force_amplitude", 1e26);
2197 *   group.set_attribute<double>("force_sigma_x", 1e-7);
2198 *   group.set_attribute<double>("force_sigma_y", 1);
2199 *   group.set_attribute<double>("max_force_width_x", 3e-7);
2200 *   group.set_attribute<double>("max_force_width_y", 2e-8);
2201 *   group.set_attribute<double>("force_x_pos", -8e-6);
2202 *   group.set_attribute<double>("force_y_pos", 0);
2203 *  
2204 *   group.set_attribute<bool>("pml_x", true);
2205 *   group.set_attribute<bool>("pml_y", false);
2206 *   group.set_attribute<double>("pml_width_x", 1.8e-6);
2207 *   group.set_attribute<double>("pml_width_y", 5e-7);
2208 *   group.set_attribute<double>("pml_coeff", 1.6);
2209 *   group.set_attribute<unsigned int>("pml_coeff_degree", 2);
2210 *  
2211 *   group.set_attribute<double>("center_frequency", 20e9);
2212 *   group.set_attribute<double>("frequency_range", 0.5e9);
2213 *   group.set_attribute<double>(
2214 *   "start_frequency",
2215 *   group.get_attribute<double>("center_frequency") -
2216 *   group.get_attribute<double>("frequency_range") / 2);
2217 *   group.set_attribute<double>(
2218 *   "stop_frequency",
2219 *   group.get_attribute<double>("center_frequency") +
2220 *   group.get_attribute<double>("frequency_range") / 2);
2221 *   group.set_attribute<unsigned int>("nb_frequency_points", 400);
2222 *  
2223 *   if (group_name == std::string("displacement"))
2224 *   group.set_attribute<std::string>(
2225 *   "simulation_name", std::string("phononic_cavity_displacement"));
2226 *   else
2227 *   group.set_attribute<std::string>(
2228 *   "simulation_name", std::string("phononic_cavity_calibration"));
2229 *  
2230 *   group.set_attribute<bool>("save_vtu_files", false);
2231 *   }
2232 *  
2233 *   {
2234 * @endcode
2235 *
2236 * Displacement simulation. The parameters are read from the
2237 * displacement HDF5 group and the results are saved in the same HDF5
2238 * group.
2239 *
2240 * @code
2241 *   auto displacement = data.open_group("displacement");
2242 *   step62::Parameters<dim> parameters(displacement);
2243 *  
2244 *   step62::ElasticWave<dim> elastic_problem(parameters);
2245 *   elastic_problem.run();
2246 *   }
2247 *  
2248 *   {
2249 * @endcode
2250 *
2251 * Calibration simulation. The parameters are read from the calibration
2252 * HDF5 group and the results are saved in the same HDF5 group.
2253 *
2254 * @code
2255 *   auto calibration = data.open_group("calibration");
2256 *   step62::Parameters<dim> parameters(calibration);
2257 *  
2258 *   step62::ElasticWave<dim> elastic_problem(parameters);
2259 *   elastic_problem.run();
2260 *   }
2261 *   }
2262 *   catch (std::exception &exc)
2263 *   {
2264 *   std::cerr << std::endl
2265 *   << std::endl
2266 *   << "----------------------------------------------------"
2267 *   << std::endl;
2268 *   std::cerr << "Exception on processing: " << std::endl
2269 *   << exc.what() << std::endl
2270 *   << "Aborting!" << std::endl
2271 *   << "----------------------------------------------------"
2272 *   << std::endl;
2273 *  
2274 *   return 1;
2275 *   }
2276 *   catch (...)
2277 *   {
2278 *   std::cerr << std::endl
2279 *   << std::endl
2280 *   << "----------------------------------------------------"
2281 *   << std::endl;
2282 *   std::cerr << "Unknown exception!" << std::endl
2283 *   << "Aborting!" << std::endl
2284 *   << "----------------------------------------------------"
2285 *   << std::endl;
2286 *   return 1;
2287 *   }
2288 *  
2289 *   return 0;
2290 *   }
2291 * @endcode
2292<a name="Results"></a><h1>Results</h1>
2293
2294
2295<a name="Resonancefrequencyandbandgap"></a><h3>Resonance frequency and bandgap</h3>
2296
2297
2298The results are analyzed in the
2299[jupyter notebook](https://github.com/dealii/dealii/blob/master/examples/step-62/step-62.ipynb)
2300with the following code
2301@code{.py}
2302h5_file = h5py.File('results.h5', 'r')
2303data = h5_file['data']
2304
2305# Gaussian function that we use to fit the resonance
2306def resonance_f(freq, freq_m, quality_factor, max_amplitude):
2307 omega = 2 * constants.pi * freq
2308 omega_m = 2 * constants.pi * freq_m
2309 gamma = omega_m / quality_factor
2310 return max_amplitude * omega_m**2 * gamma**2 / (((omega_m**2 - omega**2)**2 + gamma**2 * omega**2))
2311
2312frequency = data['displacement']['frequency'][...]
2313# Average the probe points
2314displacement = np.mean(data['displacement']['displacement'], axis=0)
2315calibration_displacement = np.mean(data['calibration']['displacement'], axis=0)
2316reflection_coefficient = displacement / calibration_displacement
2317reflectivity = (np.abs(np.mean(data['displacement']['displacement'][...]**2, axis=0))/
2318 np.abs(np.mean(data['calibration']['displacement'][...]**2, axis=0)))
2319
2320try:
2321 x_data = frequency
2322 y_data = reflectivity
2323 quality_factor_guess = 1e3
2324 freq_guess = x_data[np.argmax(y_data)]
2325 amplitude_guess = np.max(y_data)
2326 fit_result, covariance = scipy.optimize.curve_fit(resonance_f, x_data, y_data,
2327 [freq_guess, quality_factor_guess, amplitude_guess])
2328 freq_m = fit_result[0]
2329 quality_factor = np.abs(fit_result[1])
2330 max_amplitude = fit_result[2]
2331 y_data_fit = resonance_f(x_data, freq_m, quality_factor, max_amplitude)
2332
2333 fig = plt.figure()
2334 plt.plot(frequency / 1e9, reflectivity, frequency / 1e9, y_data_fit)
2335 plt.xlabel('frequency (GHz)')
2336 plt.ylabel('amplitude^2 (a.u.)')
2337 plt.title('Transmission\n' + 'freq = ' + "%.7g" % (freq_guess / 1e9) + 'GHz Q = ' + "%.6g" % quality_factor)
2338except:
2339 fig = plt.figure()
2340 plt.plot(frequency / 1e9, reflectivity)
2341 plt.xlabel('frequency (GHz)')
2342 plt.ylabel('amplitude^2 (a.u.)')
2343 plt.title('Transmission')
2344
2345fig = plt.figure()
2346plt.plot(frequency / 1e9, np.angle(reflection_coefficient))
2347plt.xlabel('frequency (GHz)')
2348plt.ylabel('phase (rad)')
2349plt.title('Phase (transmission coefficient)\n')
2350
2351plt.show()
2352h5_file.close()
2353@endcode
2354
2355A phononic cavity is characterized by the
2356[resonance frequency](https://en.wikipedia.org/wiki/Resonance) and the
2357[the quality factor](https://en.wikipedia.org/wiki/Q_factor).
2358The quality factor is equal to the ratio between the stored energy in the resonator and the energy
2359dissipated energy per cycle, which is approximately equivalent to the ratio between the
2360resonance frequency and the
2361[full width at half maximum (FWHM)](https://en.wikipedia.org/wiki/Full_width_at_half_maximum).
2362The FWHM is equal to the bandwidth over which the power of vibration is greater than half the
2363power at the resonant frequency.
2364@f[
2365Q = \frac{f_r}{\Delta f} = \frac{\omega_r}{\Delta \omega} =
23662 \pi \times \frac{\text{energy stored}}{\text{energy dissipated per cycle}}
2367@f]
2368
2369The square of the amplitude of the mechanical resonance @f$a^2@f$ as a function of the frequency
2370has a gaussian shape
2371@f[
2372a^2 = a_\textrm{max}^2\frac{\omega^2\Gamma^2}{(\omega_r^2-\omega^2)^2+\Gamma^2\omega^2}
2373@f]
2374where @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.
2375We used the previous equation in the jupyter notebook to fit the mechanical resonance.
2376
2377Given the values we have chosen for the parameters, one could estimate the resonance frequency
2378analytically. Indeed, this is then confirmed by what we get in this program:
2379the phononic superlattice cavity exhibits a mechanical resonance at 20GHz and a quality factor of 5046.
2380The following images show the transmission amplitude and phase as a function of frequency in the
2381vicinity of the resonance frequency:
2382
2383<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.05.png" height="400" />
2384<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.06.png" height="400" />
2385
2386The images above suggest that the periodic structure has its intended effect: It really only lets waves of a very
2387specific frequency pass through, whereas all other waves are reflected. This is of course precisely what one builds
2388these sorts of devices for.
2389But it is not quite this easy. In practice, there is really only a "band gap", i.e., the device blocks waves other than
2390the desired one at 20GHz only within a certain frequency range. Indeed, to find out how large this "gap" is within
2391which waves are blocked, we can extend the frequency range to 16 GHz through the appropriate parameters in the
2392input file. We then obtain the following image:
2393
2394<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.07.png" height="400" />
2395
2396What this image suggests is that in the range of around 18 to around 22 GHz, really only the waves with a frequency
2397of 20 GHz are allowed to pass through, but beyond this range, there are plenty of other frequencies that can pass
2398through the device.
2399
2400<a name="Modeprofile"></a><h3>Mode profile</h3>
2401
2402
2403We can inspect the mode profile with Paraview or VisIt.
2404As we have discussed, at resonance all the mechanical
2405energy is transmitted and the amplitude of motion is amplified inside the cavity.
2406It can be observed that the PMLs are quite effective to truncate the solution.
2407The following image shows the mode profile at resonance:
2408
2409<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.08.png" height="400" />
2410
2411On the other hand, out of resonance all the mechanical energy is
2412reflected. The following image shows the profile at 19.75 GHz.
2413Note the interference between the force pulse and the reflected wave
2414at the position @f$x=-8\mu\textrm{m}@f$.
2415
2416<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.09.png" height="400" />
2417
2418<a name="Experimentalapplications"></a><h3>Experimental applications</h3>
2419
2420
2421Phononic superlattice cavities find application in
2422[quantum optomechanics](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.86.1391).
2423Here we have presented the simulation of a 2D superlattice cavity,
2424but this code can be used as well to simulate "real world" 3D devices such as
2425[micropillar superlattice cavities](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.060101),
2426which are promising candidates to study macroscopic quantum phenomena.
2427The 20GHz mode of a micropillar superlattice cavity is essentially a mechanical harmonic oscillator that is very well isolated
2428from the environment. If the device is cooled down to 20mK in a dilution fridge, the mode would then become a
2429macroscopic quantum harmonic oscillator.
2430
2431
2432<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
2433
2434
2435Instead of setting the parameters in the C++ file we could set the parameters
2436using a python script and save them in the HDF5 file that we will use for
2437the simulations. Then the deal.II program will read the parameters from the
2438HDF5 file.
2439
2440@code{.py}
2441import numpy as np
2442import h5py
2443import matplotlib.pyplot as plt
2444import subprocess
2445import scipy.constants as constants
2446import scipy.optimize
2447
2448# This considerably reduces the size of the svg data
2449plt.rcParams['svg.fonttype'] = 'none'
2450
2451h5_file = h5py.File('results.h5', 'w')
2452data = h5_file.create_group('data')
2453displacement = data.create_group('displacement')
2454calibration = data.create_group('calibration')
2455
2456# Set the parameters
2457for group in [displacement, calibration]:
2458 # Dimensions of the domain
2459 # The waveguide length is equal to dimension_x
2460 group.attrs['dimension_x'] = 2e-5
2461 # The waveguide width is equal to dimension_y
2462 group.attrs['dimension_y'] = 2e-8
2463
2464 # Position of the probe that we use to measure the flux
2465 group.attrs['probe_pos_x'] = 8e-6
2466 group.attrs['probe_pos_y'] = 0
2467 group.attrs['probe_width_y'] = 2e-08
2468
2469 # Number of points in the probe
2470 group.attrs['nb_probe_points'] = 5
2471
2472 # Global refinement
2473 group.attrs['grid_level'] = 1
2474
2475 # Cavity
2476 group.attrs['cavity_resonance_frequency'] = 20e9
2477 group.attrs['nb_mirror_pairs'] = 15
2478
2479 # Material
2480 group.attrs['poissons_ratio'] = 0.27
2481 group.attrs['youngs_modulus'] = 270000000000.0
2482 group.attrs['material_a_rho'] = 3200
2483 if group == displacement:
2484 group.attrs['material_b_rho'] = 2000
2485 else:
2486 group.attrs['material_b_rho'] = 3200
2487 group.attrs['lambda'] = (group.attrs['youngs_modulus'] * group.attrs['poissons_ratio'] /
2488 ((1 + group.attrs['poissons_ratio']) *
2489 (1 - 2 * group.attrs['poissons_ratio'])))
2490 group.attrs['mu']= (group.attrs['youngs_modulus'] / (2 * (1 + group.attrs['poissons_ratio'])))
2491
2492 # Force
2493 group.attrs['max_force_amplitude'] = 1e26
2494 group.attrs['force_sigma_x'] = 1e-7
2495 group.attrs['force_sigma_y'] = 1
2496 group.attrs['max_force_width_x'] = 3e-7
2497 group.attrs['max_force_width_y'] = 2e-8
2498 group.attrs['force_x_pos'] = -8e-6
2499 group.attrs['force_y_pos'] = 0
2500
2501 # PML
2502 group.attrs['pml_x'] = True
2503 group.attrs['pml_y'] = False
2504 group.attrs['pml_width_x'] = 1.8e-6
2505 group.attrs['pml_width_y'] = 5e-7
2506 group.attrs['pml_coeff'] = 1.6
2507 group.attrs['pml_coeff_degree'] = 2
2508
2509 # Frequency sweep
2510 group.attrs['center_frequency'] = 20e9
2511 group.attrs['frequency_range'] = 0.5e9
2512 group.attrs['start_frequency'] = group.attrs['center_frequency'] - group.attrs['frequency_range'] / 2
2513 group.attrs['stop_frequency'] = group.attrs['center_frequency'] + group.attrs['frequency_range'] / 2
2514 group.attrs['nb_frequency_points'] = 400
2515
2516 # Other parameters
2517 if group == displacement:
2518 group.attrs['simulation_name'] = 'phononic_cavity_displacement'
2519 else:
2520 group.attrs['simulation_name'] = 'phononic_cavity_calibration'
2521 group.attrs['save_vtu_files'] = False
2522
2523h5_file.close()
2524@endcode
2525
2526In order to read the HDF5 parameters we have to use the
2527HDF5::File::FileAccessMode::open flag.
2528@code{.py}
2529 HDF5::File data_file("results.h5",
2530 HDF5::File::FileAccessMode::open,
2531 MPI_COMM_WORLD);
2532 auto data = data_file.open_group("data");
2533@endcode
2534 *
2535 *
2536<a name="PlainProg"></a>
2537<h1> The plain program</h1>
2538@include "step-62.cc"
2539*/
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
void write(const Container &data)
Definition hdf5.h:2008
void write_selection(const Container &data, const std::vector< hsize_t > &coordinates)
Definition hdf5.h:2040
void write_none()
Definition hdf5.h:2198
Group create_group(const std::string &name) const
Definition hdf5.cc:381
Definition point.h:112
Point< 2 > first
Definition grid_out.cc:4615
unsigned int level
Definition grid_out.cc:4618
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:439
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition mapping.cc:285
const Event initial
Definition event.cc:65
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={}, const double tolerance=1.e-10)
Definition hdf5.h:346
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:75
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:161
int create_group(const MPI_Comm comm, const MPI_Group &group, const int tag, MPI_Comm *new_comm)
Definition mpi.cc:214
void point_value(const DoFHandler< dim, spacedim > &dof, const VectorType &fe_function, const Point< spacedim, double > &point, Vector< typename VectorType::value_type > &value)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
long double gamma(const unsigned int n)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:71
static constexpr double PI
Definition numbers.h:259
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation