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