From 24095e08d2c7b5378956842b0da13614a04c3ec7 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Tue, 24 May 2022 11:13:40 -0500 Subject: [PATCH] Address some more review comments --- doc/doxygen/references.bib | 10 +++ examples/step-81/doc/builds-on | 1 + examples/step-81/doc/intro.dox | 28 ++++---- examples/step-81/doc/results.dox | 60 ++++++++--------- examples/step-81/doc/tooltip | 2 +- examples/step-81/step-81.cc | 111 ++++++++++++++++--------------- 6 files changed, 114 insertions(+), 98 deletions(-) diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index b87479c506..343eeb489f 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -1341,6 +1341,16 @@ SERIES={International Series in Pure and Applied Physics}, } +@ARTICLE{berenger1994, + AUTHOR={Jean-Pierre Bérenger}, + TITLE={A Perfectly Matched Layer for the Absorption of Electromagnetic + Waves}, + JOURNAL={Journal of Computational Physics}, + VOLUME={114}, + PAGES={185-200}, + YEAR={1994}, +} + @ARTICLE{Gopalakrishnan2003, AUTHOR={J. Gopalakrishnan and Pasciak, J. E.}, TITLE={Overlapping Schwarz preconditioners for indefinite time harmonic diff --git a/examples/step-81/doc/builds-on b/examples/step-81/doc/builds-on index 850b582a69..efaf55a221 100644 --- a/examples/step-81/doc/builds-on +++ b/examples/step-81/doc/builds-on @@ -1 +1,2 @@ +step-6 step-8 diff --git a/examples/step-81/doc/intro.dox b/examples/step-81/doc/intro.dox index 3c5ac3d6f0..91247e88c7 100644 --- a/examples/step-81/doc/intro.dox +++ b/examples/step-81/doc/intro.dox @@ -382,19 +382,21 @@ this causes the solution image to be distorted. In order to reduce the resonance and distortion in our solutions, we are implementing a Perfectly Matched Layer (PML) in the scattering configuration. -The concept of a PML was pioneered by Bérenger and it is is an indispensable tool -for truncating unbounded domains for wave equations and often used in the -numerical approximation of scattering problems. It is essentially a thin layer with -modified material parameters placed near the boundary such that all outgoing -electromagnetic waves decay exponentially with no “artificial” reflection due to -truncation of the domain. - -Our PML is essentially a concentric circle with modified material coefficients -($\varepsilon_r, \mu_r, \sigma$). It is located in a small region near the boundary -$\partial\Omega$ and the transformation of the material coordinates is chosen to -be a function of the radial distance $\rho$ from the origin $e_r$. The normal field -$\nu$ of $\Sigma$ is orthogonal to the radial direction $e_r$, which makes -$\mathbf{J}_a \equiv 0$ and $\mathbf{M}_a \equiv 0$ within the PML. +The concept of a PML was pioneered by Bérenger @cite Berenger1994 +and it is is an indispensable tool for truncating unbounded domains for +wave equations and often used in the numerical approximation of scattering +problems. It is essentially a thin layer with modified material parameters +placed near the boundary such that all outgoing electromagnetic waves decay +exponentially with no “artificial” reflection due to truncation of the +domain. + +Our PML is a concentric circle with modified material coefficients +($\varepsilon_r, \mu_r, \sigma$). It is located in a small region near the +boundary $\partial\Omega$ and the transformation of the material +coordinates is chosen to be a function of the radial distance $\rho$ from +the origin $e_r$. The normal field $\nu$ of $\Sigma$ is orthogonal to the +radial direction $e_r$, which makes $\mathbf{J}_a \equiv 0$ and +$\mathbf{M}_a \equiv 0$ within the PML. @htmlonly

diff --git a/examples/step-81/doc/results.dox b/examples/step-81/doc/results.dox index 8542186588..80c3bfd7d8 100644 --- a/examples/step-81/doc/results.dox +++ b/examples/step-81/doc/results.dox @@ -44,12 +44,12 @@ Following are the output images:

Solution with no interface, Dirichlet boundary conditions and PML strength 0.

- + Visualization of the solution of step-81 with no interface, no absorbing boundary conditions and PML strength 0 -

Solution with no interface, absorbing boundary conditions and PML strength 0.

- - - +

Solution with no interface, absorbing boundary conditions and PML strength 0.

+ + + Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4

Solution with no interface, absorbing boundary conditions and PML strength 4.

@@ -82,19 +82,19 @@ of $E_x$. - - - + +
+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0

Solution with an interface, Dirichlet boundary conditions and PML strength 0.

- 		+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 -

Solution with an interface, absorbing boundary conditions and PML strength 0.

+

Solution with an interface, absorbing boundary conditions and PML strength 0.

 Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4 -

Solution with an interface, absorbing boundary conditions and PML strength 4.

+

Solution with an interface, absorbing boundary conditions and PML strength 4.
@@ -102,19 +102,19 @@ of $E_x$. - - - + +
+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 -

Solution with an interface, Dirichlet boundary conditions and PML strength 0.

- 		+

Solution with an interface, Dirichlet boundary conditions and PML strength 0.

+ 		Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 -

Solution with an interface, absorbing boundary conditions and PML strength 0.

+

Solution with an interface, absorbing boundary conditions and PML strength 0.

 Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4 -

Solution with an interface, absorbing boundary conditions and PML strength 4.

+

Solution with an interface, absorbing boundary conditions and PML strength 4.
@@ -130,19 +130,19 @@ the standing wave will dissipate more within the PML ring. Here are some animations to demonstrate the effect of the PML - - - + +
+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0

Solution with an interface, Dirichlet boundary conditions and PML strength 0.

- 		+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 -

Solution with an interface, absorbing boundary conditions and PML strength 0.

+

Solution with an interface, absorbing boundary conditions and PML strength 0.

 Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4 -

Solution with an interface, absorbing boundary conditions and PML strength 4.

+

Solution with an interface, absorbing boundary conditions and PML strength 4.
@@ -150,19 +150,19 @@ Here are some animations to demonstrate the effect of the PML - - - + +
+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0

Solution with an interface, Dirichlet boundary conditions and PML strength 0.

- 		+ Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0 -

Solution with an interface, absorbing boundary conditions and PML strength 0.

+

Solution with an interface, absorbing boundary conditions and PML strength 0.

 Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4 -

Solution with an interface, absorbing boundary conditions and PML strength 4.

+

Solution with an interface, absorbing boundary conditions and PML strength 4.
diff --git a/examples/step-81/doc/tooltip b/examples/step-81/doc/tooltip index ae43afaffd..1ed4f5b8c6 100644 --- a/examples/step-81/doc/tooltip +++ b/examples/step-81/doc/tooltip @@ -1 +1 @@ -Creating a mesh. Refining it. Writing it to a file. +A time-harmonic Maxwell solver with interface jump conditions diff --git a/examples/step-81/step-81.cc b/examples/step-81/step-81.cc index 6b9e964724..0cb2dfeaf2 100644 --- a/examples/step-81/step-81.cc +++ b/examples/step-81/step-81.cc @@ -22,7 +22,8 @@ // The set of include files is quite standard. The most notable include is // the fe/fe_nedelec_sz.h file which allows us to use the FE_NedelecSZ elements. // This is an implementation of the $H^{curl}$ conforming Nédélec Elements -// that resolves the sign conflict issues that arise from parametrization. +// that resolves the sign conflict issues that arise from parametrization +// (for details we refer to the documentation of the FE_NedelecSZ element). #include #include @@ -257,7 +258,8 @@ namespace Step81 class PerfectlyMatchedLayer : public ParameterAcceptor { public: - static_assert(dim == 2, "dim == 2"); /* only works in 2D */ + static_assert(dim == 2, + "The perfectly matched layer is only implemented in 2D."); Parameters parameters; @@ -267,10 +269,6 @@ namespace Step81 PerfectlyMatchedLayer(); - double inner_radius; - double outer_radius; - double strength; - std::complex d(const Point point); std::complex d_bar(const Point point); @@ -285,6 +283,11 @@ namespace Step81 rank2_type b_matrix(const Point point); rank2_type c_matrix(const Point point); + + private: + double inner_radius; + double outer_radius; + double strength; }; @@ -309,11 +312,18 @@ namespace Step81 typename std::complex PerfectlyMatchedLayer::d(const Point point) { - const auto radius = point.norm(); - const double s = - strength * ((radius - inner_radius) * (radius - inner_radius)) / - ((outer_radius - inner_radius) * (outer_radius - inner_radius)); - return 1.0 + 1.0i * s; + const auto radius = point.norm(); + if (radius > inner_radius) + { + const double s = + strength * ((radius - inner_radius) * (radius - inner_radius)) / + ((outer_radius - inner_radius) * (outer_radius - inner_radius)); + return 1.0 + 1.0i * s; + } + else + { + return 1.0; + } } @@ -321,13 +331,21 @@ namespace Step81 typename std::complex PerfectlyMatchedLayer::d_bar(const Point point) { - const auto radius = point.norm(); - const double s_bar = - strength / 3. * - ((radius - inner_radius) * (radius - inner_radius) * - (radius - inner_radius)) / - (radius * (outer_radius - inner_radius) * (outer_radius - inner_radius)); - return 1.0 + 1.0i * s_bar; + const auto radius = point.norm(); + if (radius > inner_radius) + { + const double s_bar = + strength / 3. * + ((radius - inner_radius) * (radius - inner_radius) * + (radius - inner_radius)) / + (radius * (outer_radius - inner_radius) * + (outer_radius - inner_radius)); + return 1.0 + 1.0i * s_bar; + } + else + { + return 1.0; + } } @@ -639,24 +657,19 @@ namespace Step81 for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { const Point &position = quadrature_points[q_point]; - const auto radius = position.norm(); - const auto inner_radius = perfectly_matched_layer.inner_radius; auto mu_inv = parameters.mu_inv(position, id); auto epsilon = parameters.epsilon(position, id); const auto J_a = parameters.J_a(position, id); - if (radius >= inner_radius) - { - auto A = perfectly_matched_layer.a_matrix(position); - auto B = perfectly_matched_layer.b_matrix(position); - auto d = perfectly_matched_layer.d(position); + const auto A = perfectly_matched_layer.a_matrix(position); + const auto B = perfectly_matched_layer.b_matrix(position); + const auto d = perfectly_matched_layer.d(position); - mu_inv = mu_inv / d; - epsilon = invert(A) * epsilon * invert(B); - }; + mu_inv = mu_inv / d; + epsilon = invert(A) * epsilon * invert(B); - for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (const auto i : fe_values.dof_indices()) { const auto phi_i = real_part.value(i, q_point) - 1.0i * imag_part.value(i, q_point); @@ -667,7 +680,7 @@ namespace Step81 (1.0i * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point); cell_rhs(i) += rhs_value.real(); - for (unsigned int j = 0; j < dofs_per_cell; ++j) + for (const auto j : fe_values.dof_indices()) { const auto phi_j = real_part.value(j, q_point) + 1.0i * imag_part.value(j, q_point); @@ -708,34 +721,30 @@ namespace Step81 ++q_point) { const auto &position = quadrature_points[q_point]; - const auto radius = position.norm(); - const auto inner_radius = - perfectly_matched_layer.inner_radius; auto mu_inv = parameters.mu_inv(position, id); auto epsilon = parameters.epsilon(position, id); - if (radius >= inner_radius) - { - auto A = perfectly_matched_layer.a_matrix(position); - auto B = perfectly_matched_layer.b_matrix(position); - auto d = perfectly_matched_layer.d(position); + const auto A = + perfectly_matched_layer.a_matrix(position); + const auto B = + perfectly_matched_layer.b_matrix(position); + const auto d = perfectly_matched_layer.d(position); - mu_inv = mu_inv / d; - epsilon = invert(A) * epsilon * invert(B); - }; + mu_inv = mu_inv / d; + epsilon = invert(A) * epsilon * invert(B); const auto normal = fe_face_values.normal_vector(q_point); - for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (const auto i : fe_face_values.dof_indices()) { const auto phi_i = real_part.value(i, q_point) - 1.0i * imag_part.value(i, q_point); const auto phi_i_T = tangential_part(phi_i, normal); - for (unsigned int j = 0; j < dofs_per_cell; ++j) + for (const auto j : fe_face_values.dof_indices()) { const auto phi_j = real_part.value(j, q_point) + @@ -775,28 +784,22 @@ namespace Step81 ++q_point) { const auto &position = quadrature_points[q_point]; - const auto radius = position.norm(); - const auto inner_radius = - perfectly_matched_layer.inner_radius; auto sigma = parameters.sigma(position, id1, id2); - if (radius >= inner_radius) - { - auto B = perfectly_matched_layer.b_matrix(position); - auto C = perfectly_matched_layer.c_matrix(position); - sigma = invert(C) * sigma * invert(B); - }; + const auto B = perfectly_matched_layer.b_matrix(position); + const auto C = perfectly_matched_layer.c_matrix(position); + sigma = invert(C) * sigma * invert(B); const auto normal = fe_face_values.normal_vector(q_point); - for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (const auto i : fe_face_values.dof_indices()) { const auto phi_i = real_part.value(i, q_point) - 1.0i * imag_part.value(i, q_point); const auto phi_i_T = tangential_part(phi_i, normal); - for (unsigned int j = 0; j < dofs_per_cell; ++j) + for (const auto j : fe_face_values.dof_indices()) { const auto phi_j = real_part.value(j, q_point) + -- 2.39.5