1021 * The speed of sound is defined by
1023 * c = \frac{K_e}{\rho}
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
1030 * K_e = 4\mu\frac{\lambda +\mu}{\lambda+2\mu}
1032 * and
for the three dimensional
case @f$K_e@f$ is
equal to the Young
's modulus.
1034 * K_e = \mu\frac{3\lambda +2\mu}{\lambda+\mu}
1038 * double elastic_constant;
1041 * elastic_constant = 4 * mu * (lambda + mu) / (lambda + 2 * mu);
1043 * else if (dim == 3)
1045 * elastic_constant = mu * (3 * lambda + 2 * mu) / (lambda + mu);
1049 * Assert(false, ExcInternalError());
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;
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"
1066 * where the brown color represents material_a and the green color
1067 * represents material_b.
1070 * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
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))
1079 * const double coefficient =
1081 * (layer_transition_center - average_rho_width / 2)) /
1082 * average_rho_width;
1083 * return (1 - coefficient) * material_a_rho +
1084 * coefficient * material_b_rho;
1090 * Here we define the
1092 * smoothing](https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/)
1093 * which improves the precision of the simulation.
1096 * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
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))
1106 * const double coefficient =
1108 * (layer_transition_center - average_rho_width / 2)) /
1109 * average_rho_width;
1110 * return (1 - coefficient) * material_b_rho +
1111 * coefficient * material_a_rho;
1120 * if (std::abs(p[0]) <= material_a_wavelength / 2)
1122 * return material_a_rho;
1127 * the material_a layers
1130 * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
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))
1140 * return material_a_rho;
1146 * the material_b layers
1149 * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
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))
1159 * return material_b_rho;
1165 * and finally the default is material_a.
1168 * return material_a_rho;
1176 * <a name="TheParametersclassimplementation"></a>
1177 * <h4>The `Parameters` class implementation</h4>
1181 * The constructor reads all the parameters from the HDF5::Group `data` using
1182 * the HDF5::Group::get_attribute() function.
1185 * template <int dim>
1186 * Parameters<dim>::Parameters(HDF5::Group &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)
1215 * <a name="TheQuadratureCacheclassimplementation"></a>
1216 * <h4>The `QuadratureCache` class implementation</h4>
1220 * We need to reserve enough space for the mass and stiffness matrices and the
1221 * right hand side vector.
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)
1237 * <a name="ImplementationoftheElasticWaveclass"></a>
1238 * <h3>Implementation of the `ElasticWave` class</h3>
1243 * <a name="Constructor"></a>
1244 * <h4>Constructor</h4>
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.
1255 * template <int dim>
1256 * ElasticWave<dim>::ElasticWave(const Parameters<dim> ¶meters)
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>(
1270 * std::vector<hsize_t>{parameters.nb_frequency_points}))
1271 * , probe_positions_dataset(parameters.data.template create_dataset<double>(
1273 * std::vector<hsize_t>{parameters.nb_probe_points, dim}))
1275 * parameters.data.template create_dataset<std::complex<double>>(
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,
1283 * TimerOutput::summary,
1284 * TimerOutput::wall_times)
1292 * <a name="ElasticWavesetup_system"></a>
1293 * <h4>ElasticWave::setup_system</h4>
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.
1302 *
template <
int dim>
1303 *
void ElasticWave<dim>::setup_system()
1307 * dof_handler.distribute_dofs(fe);
1309 * locally_owned_dofs = dof_handler.locally_owned_dofs();
1312 * locally_relevant_solution.reinit(locally_owned_dofs,
1313 * locally_relevant_dofs,
1314 * mpi_communicator);
1316 * system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1318 * constraints.clear();
1319 * constraints.reinit(locally_relevant_dofs);
1322 * constraints.close();
1328 * locally_owned_dofs,
1330 * locally_relevant_dofs);
1332 * system_matrix.reinit(locally_owned_dofs,
1333 * locally_owned_dofs,
1335 * mpi_communicator);
1343 * <a name=
"ElasticWaveassemble_system"></a>
1344 * <h4>ElasticWave::assemble_system</h4>
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.
1355 *
template <
int dim>
1356 *
void ElasticWave<dim>::assemble_system(
const double omega,
1357 *
const bool calculate_quadrature_data)
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();
1371 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1375 * Here we store the
value of the right hand side, rho and the PML.
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));
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.
1391 * get_stiffness_tensor<dim>(parameters.lambda, parameters.mu);
1395 * We use the same method of @ref step_20
"step-20" for vector-valued problems.
1400 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1401 *
if (cell->is_locally_owned())
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.
1414 *
if (calculate_quadrature_data)
1416 * fe_values.reinit(cell);
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(),
1422 * parameters.pml.vector_value_list(
1423 * fe_values.get_quadrature_points(), pml_values);
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
1434 * QuadratureCache<dim> *local_quadrature_points_data =
1435 *
reinterpret_cast<QuadratureCache<dim> *
>(cell->user_pointer());
1436 *
Assert(local_quadrature_points_data >= &quadrature_cache.front(),
1438 *
Assert(local_quadrature_points_data <= &quadrature_cache.back(),
1440 *
for (
unsigned int q = 0; q < n_q_points; ++q)
1444 * The quadrature_data variable is used to store the mass and
1445 * stiffness matrices, the right hand side vector and the
value
1449 * QuadratureCache<dim> &quadrature_data =
1450 * local_quadrature_points_data[q];
1454 * Below we declare the force vector and the parameters of the
1455 * PML @f$s@f$ and @f$\xi@f$.
1460 * std::complex<double> xi(1, 0);
1464 * The following block is calculated only in the
first frequency
1468 *
if (calculate_quadrature_data)
1472 * Store the
value of `JxW`.
1475 * quadrature_data.JxW = fe_values.JxW(q);
1477 *
for (
unsigned int component = 0; component < dim; ++component)
1481 * Convert vectors to tensors and calculate xi
1484 * force[component] = rhs_values[q][component];
1485 * s[component] = pml_values[q][component];
1486 * xi *= s[component];
1491 * Here we calculate the @f$\alpha_{mnkl}@f$ and @f$\beta_{mnkl}@f$
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)
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]);
1510 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1513 * fe_values[displacement].value(i, q);
1515 * fe_values[displacement].gradient(i, q);
1517 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1520 * fe_values[displacement].value(j, q);
1522 * fe_values[displacement].gradient(j, q);
1529 * quadrature_data.mass_coefficient[i][j] =
1530 * rho_values[q] * xi * phi_i * phi_j;
1534 * Loop over the @f$mnkl@f$ indices of the stiffness
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)
1546 * Here we calculate the stiffness
matrix.
1547 * Note that the stiffness
matrix is not
1548 *
symmetric because of the PMLs. We use the
1550 * [documentation](https:
1551 * which is a <code>
Tensor@<2,dim@></code>.
1552 * The
matrix @f$G_{ij}@f$ consists of entries
1555 * \frac{\partial\phi_i}{\partial x_j}
1556 * =\partial_j \phi_i
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.
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]);
1577 * quadrature_data.stiffness_coefficient[i][j] =
1578 * stiffness_coefficient;
1583 * and the
value of the right hand side in
1587 * quadrature_data.right_hand_side[i] =
1588 * phi_i * force * fe_values.JxW(q);
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.
1600 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1602 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1604 * std::complex<double> matrix_sum = 0;
1605 * matrix_sum += -
std::pow(omega, 2) *
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;
1610 * cell_rhs(i) += quadrature_data.right_hand_side[i];
1613 * cell->get_dof_indices(local_dof_indices);
1614 * constraints.distribute_local_to_global(
cell_matrix,
1616 * local_dof_indices,
1628 * <a name=
"ElasticWavesolve"></a>
1629 * <h4>ElasticWave::solve</h4>
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.
1641 *
template <
int dim>
1642 *
void ElasticWave<dim>::solve()
1646 * locally_owned_dofs, mpi_communicator);
1650 * solver.solve(system_matrix, completely_distributed_solution, system_rhs);
1652 * pcout <<
" Solved in " << solver_control.
last_step() <<
" iterations."
1654 * constraints.distribute(completely_distributed_solution);
1655 * locally_relevant_solution = completely_distributed_solution;
1661 * <a name=
"ElasticWaveinitialize_position_vector"></a>
1662 * <h4>ElasticWave::initialize_position_vector</h4>
1666 * We use
this function to calculate the
values of the position vector.
1669 *
template <
int dim>
1670 *
void ElasticWave<dim>::initialize_probe_positions_vector()
1672 *
for (
unsigned int position_idx = 0;
1673 * position_idx < parameters.nb_probe_points;
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>)`
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];
1691 * probe_positions[position_idx][2] = p[2];
1699 * <a name=
"ElasticWavestore_frequency_step_data"></a>
1700 * <h4>ElasticWave::store_frequency_step_data</h4>
1704 * This function stores in the
HDF5 file the measured energy by the probe.
1707 *
template <
int dim>
1709 * ElasticWave<dim>::store_frequency_step_data(
const unsigned int frequency_idx)
1715 * We store the displacement in the @f$x@f$ direction; the displacement in the
1716 * @f$y@f$ direction is negligible.
1719 *
const unsigned int probe_displacement_component = 0;
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.
1728 * std::vector<hsize_t> coordinates;
1729 * std::vector<std::complex<double>> displacement_data;
1734 * std::vector<bool> marked_vertices = {};
1735 *
const double tolerance = 1.e-10;
1737 *
for (
unsigned int position_idx = 0;
1738 * position_idx < parameters.nb_probe_points;
1742 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1744 *
point[dim_idx] = probe_positions[position_idx][dim_idx];
1746 *
bool point_in_locally_owned_cell =
false;
1749 * cache,
point, cell_hint, marked_vertices, tolerance);
1752 * cell_hint = cell_and_ref_point.first;
1753 * point_in_locally_owned_cell =
1754 * cell_and_ref_point.first->is_locally_owned();
1757 *
if (point_in_locally_owned_cell)
1761 * Then we can store the
values of the displacement in the points of
1762 * the probe in `displacement_data`.
1767 * locally_relevant_solution,
1770 * coordinates.emplace_back(position_idx);
1771 * coordinates.emplace_back(frequency_idx);
1772 * displacement_data.emplace_back(
1773 * tmp_vector(probe_displacement_component));
1779 * We write the displacement data in the
HDF5 file. The
call
1781 * the processes have to participate.
1784 * if (coordinates.size() > 0)
1786 * displacement.write_selection(displacement_data, coordinates);
1790 * Therefore even
if the process has no data to write it has to participate
1792 * Note that we have to specify the data type, in
this case
1793 * `std::complex<double>`.
1798 * displacement.write_none<std::complex<double>>();
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".
1808 *
if (parameters.save_vtu_files)
1810 * std::vector<std::string> solution_names(dim,
"displacement");
1811 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1817 * locally_relevant_solution,
1821 *
for (
unsigned int i = 0; i < subdomain.size(); ++i)
1825 * std::vector<Vector<double>> force(
1827 * std::vector<Vector<double>> pml(
1833 *
if (cell->is_locally_owned())
1835 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
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();
1842 * rho(cell->active_cell_index()) =
1843 * parameters.rho.value(cell->center());
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:
1855 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1857 * force[dim_idx](cell->active_cell_index()) = -1
e+20;
1858 * pml[dim_idx](cell->active_cell_index()) = -1
e+20;
1860 * rho(cell->active_cell_index()) = -1
e+20;
1864 *
for (
unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
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 * std::string filename = (parameters.simulation_name +
"_" +
1881 * frequency_idx_stream.str() +
".vtu");
1891 * <a name=
"ElasticWaveoutput_results"></a>
1892 * <h4>ElasticWave::output_results</h4>
1896 * This function writes the datasets that have not already been written.
1899 *
template <
int dim>
1900 *
void ElasticWave<dim>::output_results()
1904 * The vectors `frequency` and `position` are the same
for all the
1905 * processes. Therefore any of the processes can write the corresponding
1912 * frequency_dataset.write(frequency);
1913 * probe_positions_dataset.write(probe_positions);
1917 * frequency_dataset.write_none<
double>();
1918 * probe_positions_dataset.write_none<
double>();
1927 * <a name=
"ElasticWavesetup_quadrature_cache"></a>
1928 * <h4>ElasticWave::setup_quadrature_cache</h4>
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".
1937 *
template <
int dim>
1938 *
void ElasticWave<dim>::setup_quadrature_cache()
1943 * std::vector<QuadratureCache<dim>> tmp;
1944 * quadrature_cache.swap(tmp);
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())
1954 * cell->set_user_pointer(&quadrature_cache[cache_index]);
1955 * cache_index += quadrature_formula.size();
1965 * <a name=
"ElasticWavefrequency_sweep"></a>
1966 * <h4>ElasticWave::frequency_sweep</h4>
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.
1975 *
template <
int dim>
1976 *
void ElasticWave<dim>::frequency_sweep()
1978 *
for (
unsigned int frequency_idx = 0;
1979 * frequency_idx < parameters.nb_frequency_points;
1982 * pcout << parameters.simulation_name +
" frequency idx: "
1983 * << frequency_idx <<
'/' << parameters.nb_frequency_points - 1
1989 *
if (frequency_idx == 0)
1991 * pcout <<
" Number of active cells : "
1993 * pcout <<
" Number of degrees of freedom : "
1994 * << dof_handler.n_dofs() << std::endl;
1997 *
if (frequency_idx == 0)
2001 * Write the simulation parameters only once
2004 * parameters.data.set_attribute(
"active_cells",
2006 * parameters.data.set_attribute(
"degrees_of_freedom",
2007 * dof_handler.n_dofs());
2012 * We calculate the frequency and omega
values for this particular step.
2015 *
const double current_loop_frequency =
2016 * (parameters.start_frequency +
2018 * (parameters.stop_frequency - parameters.start_frequency) /
2019 * (parameters.nb_frequency_points - 1));
2020 *
const double current_loop_omega =
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
2031 * assemble_system(current_loop_omega,
2032 * (frequency_idx == 0) ?
true :
false);
2035 * frequency[frequency_idx] = current_loop_frequency;
2036 * store_frequency_step_data(frequency_idx);
2038 * computing_timer.print_summary();
2039 * computing_timer.reset();
2040 * pcout << std::endl;
2049 * <a name=
"ElasticWaverun"></a>
2054 * This function is very similar to the
one in @ref step_40
"step-40".
2057 *
template <
int dim>
2061 * pcout <<
"Debug mode" << std::endl;
2063 * pcout <<
"Release mode" << std::endl;
2068 * p1(0) = -parameters.dimension_x / 2;
2069 * p1(1) = -parameters.dimension_y / 2;
2072 * p1(2) = -parameters.dimension_y / 2;
2075 * p2(0) = parameters.dimension_x / 2;
2076 * p2(1) = parameters.dimension_y / 2;
2079 * p2(2) = parameters.dimension_y / 2;
2081 * std::vector<unsigned int> divisions(dim);
2082 * divisions[0] =
int(parameters.dimension_x / parameters.dimension_y);
2096 * setup_quadrature_cache();
2098 * initialize_probe_positions_vector();
2100 * frequency_sweep();
2111 * <a name=
"Themainfunction"></a>
2112 * <h4>The main function</h4>
2116 * The main function is very similar to the
one in @ref step_40
"step-40".
2119 *
int main(
int argc,
char *argv[])
2123 *
using namespace dealii;
2124 *
const unsigned int dim = 2;
2131 *
auto data = data_file.create_group(
"data");
2135 * Each of the simulations (displacement and calibration) is stored in a
2136 * separate
HDF5 group:
2139 *
const std::vector<std::string> group_names = {
"displacement",
2141 *
for (
auto group_name : group_names)
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
2163 *
auto group = data.create_group(group_name);
2165 * group.set_attribute<
double>(
"dimension_x", 2
e-5);
2166 * group.set_attribute<
double>(
"dimension_y", 2
e-8);
2167 * group.set_attribute<
double>(
"probe_pos_x", 8
e-6);
2168 * group.set_attribute<
double>(
"probe_pos_y", 0);
2169 * group.set_attribute<
double>(
"probe_width_y", 2
e-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);
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);
2179 *
if (group_name == std::string(
"displacement"))
2180 * group.set_attribute<
double>(
"material_b_rho", 2000);
2182 * group.set_attribute<
double>(
"material_b_rho", 3200);
2184 * group.set_attribute(
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"))));
2195 * group.set_attribute<
double>(
"max_force_amplitude", 1e26);
2196 * group.set_attribute<
double>(
"force_sigma_x", 1
e-7);
2197 * group.set_attribute<
double>(
"force_sigma_y", 1);
2198 * group.set_attribute<
double>(
"max_force_width_x", 3
e-7);
2199 * group.set_attribute<
double>(
"max_force_width_y", 2
e-8);
2200 * group.set_attribute<
double>(
"force_x_pos", -8
e-6);
2201 * group.set_attribute<
double>(
"force_y_pos", 0);
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", 5
e-7);
2207 * group.set_attribute<
double>(
"pml_coeff", 1.6);
2208 * group.set_attribute<
unsigned int>(
"pml_coeff_degree", 2);
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>(
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);
2222 *
if (group_name == std::string(
"displacement"))
2223 * group.set_attribute<std::string>(
2224 *
"simulation_name", std::string(
"phononic_cavity_displacement"));
2226 * group.set_attribute<std::string>(
2227 *
"simulation_name", std::string(
"phononic_cavity_calibration"));
2229 * group.set_attribute<
bool>(
"save_vtu_files",
false);
2235 * Displacement simulation. The parameters are read from the
2236 * displacement
HDF5 group and the results are saved in the same
HDF5
2240 *
auto displacement = data.open_group(
"displacement");
2241 * step62::Parameters<dim> parameters(displacement);
2243 * step62::ElasticWave<dim> elastic_problem(parameters);
2244 * elastic_problem.run();
2250 * Calibration simulation. The parameters are read from the calibration
2251 *
HDF5 group and the results are saved in the same
HDF5 group.
2254 *
auto calibration = data.open_group(
"calibration");
2255 * step62::Parameters<dim> parameters(calibration);
2257 * step62::ElasticWave<dim> elastic_problem(parameters);
2258 * elastic_problem.run();
2261 *
catch (std::exception &exc)
2263 * std::cerr << std::endl
2265 * <<
"----------------------------------------------------"
2267 * std::cerr <<
"Exception on processing: " << std::endl
2268 * << exc.what() << std::endl
2269 * <<
"Aborting!" << std::endl
2270 * <<
"----------------------------------------------------"
2277 * std::cerr << std::endl
2279 * <<
"----------------------------------------------------"
2281 * std::cerr <<
"Unknown exception!" << std::endl
2282 * <<
"Aborting!" << std::endl
2283 * <<
"----------------------------------------------------"
2291<a name=
"Results"></a><h1>Results</h1>
2294<a name=
"Resonancefrequencyandbandgap"></a><h3>Resonance frequency and bandgap</h3>
2297The results are analyzed in the
2298[jupyter notebook](https:
2299with the following code
2301h5_file = h5py.File(
'results.h5',
'r')
2302data = h5_file[
'data']
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))
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)))
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)
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)
2339 plt.plot(frequency / 1e9, reflectivity)
2340 plt.xlabel(
'frequency (GHz)')
2341 plt.ylabel(
'amplitude^2 (a.u.)')
2342 plt.title(
'Transmission')
2345plt.plot(frequency / 1e9, np.
angle(reflection_coefficient))
2346plt.xlabel(
'frequency (GHz)')
2347plt.ylabel(
'phase (rad)')
2348plt.title(
'Phase (transmission coefficient)\n')
2354A phononic cavity is characterized by the
2355[resonance frequency](https:
2356[the quality factor](https:
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:
2361The FWHM is
equal to the bandwidth over which the power of vibration is greater than half the
2362power at the resonant frequency.
2364Q = \frac{f_r}{\Delta f} = \frac{\omega_r}{\Delta \omega} =
23652 \pi \times \frac{\text{energy stored}}{\text{energy dissipated per cycle}}
2368The square of the amplitude of the mechanical resonance @f$a^2@f$ as a function of the frequency
2371a^2 = a_\textrm{
max}^2\frac{\omega^2\Gamma^2}{(\omega_r^2-\omega^2)^2+\Gamma^2\omega^2}
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.
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:
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" />
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:
2393<img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.07.png" height=
"400" />
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
2399<a name=
"Modeprofile"></a><h3>Mode profile</h3>
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:
2408<img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.08.png" height=
"400" />
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$.
2415<img alt=
"Phononic superlattice cavity" src=
"https://www.dealii.org/images/steps/developer/step-62.09.png" height=
"400" />
2417<a name=
"Experimentalapplications"></a><h3>Experimental applications</h3>
2420Phononic superlattice cavities find application in
2421[quantum optomechanics](https:
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:
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.
2431<a name=
"Possibilitiesforextensions"></a><h3>Possibilities
for extensions</h3>
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
2442import matplotlib.pyplot as plt
2444import scipy.constants as constants
2445import scipy.optimize
2447# This considerably reduces the size of the svg data
2448plt.rcParams[
'svg.fonttype'] =
'none'
2450h5_file = h5py.
File(
'results.h5',
'w')
2456for group in [displacement, calibration]:
2457 # Dimensions of the domain
2458 # The waveguide length is equal to dimension_x
2459 group.attrs[
'dimension_x'] = 2
e-5
2460 # The waveguide width is equal to dimension_y
2461 group.attrs[
'dimension_y'] = 2
e-8
2463 # Position of the probe that we use to measure the flux
2464 group.attrs[
'probe_pos_x'] = 8
e-6
2465 group.attrs[
'probe_pos_y'] = 0
2466 group.attrs[
'probe_width_y'] = 2
e-08
2468 # Number of points in the probe
2469 group.attrs[
'nb_probe_points'] = 5
2472 group.attrs[
'grid_level'] = 1
2475 group.attrs[
'cavity_resonance_frequency'] = 20e9
2476 group.attrs[
'nb_mirror_pairs'] = 15
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
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'])))
2492 group.attrs[
'max_force_amplitude'] = 1e26
2493 group.attrs[
'force_sigma_x'] = 1
e-7
2494 group.attrs[
'force_sigma_y'] = 1
2495 group.attrs[
'max_force_width_x'] = 3
e-7
2496 group.attrs[
'max_force_width_y'] = 2
e-8
2497 group.attrs[
'force_x_pos'] = -8
e-6
2498 group.attrs[
'force_y_pos'] = 0
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'] = 5
e-7
2505 group.attrs[
'pml_coeff'] = 1.6
2506 group.attrs[
'pml_coeff_degree'] = 2
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
2516 if group == displacement:
2517 group.attrs[
'simulation_name'] =
'phononic_cavity_displacement'
2519 group.attrs[
'simulation_name'] =
'phononic_cavity_calibration'
2520 group.attrs[
'save_vtu_files'] = False
2525In order to read the
HDF5 parameters we have to use the
2531 auto data = data_file.open_group(
"data");
2535<a name=
"PlainProg"></a>
2536<h1> The plain program</h1>
2537@include
"step-62.cc"
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=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
virtual void build_patches(const unsigned int n_subdivisions=0)
void write(const Container &data)
void write_selection(const Container &data, const std::vector< hsize_t > &coordinates)
Group create_group(const std::string &name) const
unsigned int last_step() const
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
__global__ void set(Number *val, const Number s, const size_type N)
void write_vtu_in_parallel(const std::string &filename, const MPI_Comm &comm) const
#define Assert(cond, exc)
std::string to_string(const T &t)
static ::ExceptionBase & ExcInternalError()
void loop(ITERATOR begin, typename identity< ITERATOR >::type 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())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
@ component_is_part_of_vector
Expression log10(const Expression &x)
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)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
static const types::blas_int one
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.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
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)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
void run(const Iterator &begin, const typename identity< Iterator >::type &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)
static constexpr double PI
::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