From 7d9d051857d5e8883bd3906779a7a45968282f9c Mon Sep 17 00:00:00 2001 From: Manaswinee Bezbaruah Date: Wed, 2 Feb 2022 20:28:16 -0600 Subject: [PATCH] Step-81: add more text to the results section --- examples/step-81/doc/results.dox | 164 ++++++++++++++++++++++--------- examples/step-81/step-81.cc | 80 ++++++++------- 2 files changed, 162 insertions(+), 82 deletions(-) diff --git a/examples/step-81/doc/results.dox b/examples/step-81/doc/results.dox index db89021b24..dee3123006 100644 --- a/examples/step-81/doc/results.dox +++ b/examples/step-81/doc/results.dox @@ -8,96 +8,168 @@ Number of degrees of freedom: 16640 Program ended with exit code: 0 @endcode -We first run the code without the interface and with the dipole centered at $(0,0)$. To demonstrate the effect of the PML, we first set the strength to 0. Change the following parameters in the .prm file: +

Absorbing boundary conditions and the PML

+ +The following images are the outputs for the imaginary $E_x$ without the interface and with the dipole centered at $(0,0)$. In order to remove the interface, the surface conductivity is set to 0. First, we turn off the absorbing boundary conditions and the PML. Second, we want to see the effect of the PML when absorbing boundary conditions apply. So we set absorbing boundary conditions to true and leave the PML strength to 0. Lastly, we increase the strength of the PML to 4. Change the following in the .prm file: @code -# is there an interface? - set interface boolean = false +# use absorbing boundary conditions? + set absorbing boundary condition boolean = false # position of the dipole set dipole position = 0, 0 # strength of the PML set strength = 0 -@endcode - -The imaginary component of the E_x wave is shown below. We observe that in absence of the PML, there is a lot of distortion and resonance (the real parts will not be generated without a PML). Now, we increase the strength of the PML to 4. -@code -# strength of the PML - set strength = 4 +# surface conductivity between material 1 and material 2 + set sigma = 0, 0; 0, 0| 0, 0; 0, 0 @endcode - -As we see, a much more coherent image is generated with an appropriate PML. + +Following are the output images: - - - + + + +
- Visualization of the solution of step-81 with no interface and PML strength 0 - - Visualization of the solution of step-81 with no interface and PML strength 4 + + Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 0 + + Visualization of the solution of step-81 with no interface, no absorbing boundary conditions and PML strength 0 + + Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4
-Now, let's generate a standing wave by adding an interface at the center. In order to observe this effect, we offset the center of the dipole to $(0, 0.8)$ and set the interface boolean to $true$ +We observe that with absorbing boundary conditions and in absence of the PML, there is a lot of distortion and resonance (the real parts will not be generated without a PML). This is, as we stipulated, due to reflection from infinity. As we see, a much more coherent image is generated with an appropriate PML. -@code -# is there an interface? - set interface boolean = true +

Surface Plasmon Polariton

+Now, let's generate a standing wave by adding an interface at the center. In order to observe this effect, we offset the center of the dipole to $(0, 0.8)$ and set the surface conductivity back to $(0.001, 0.2)$: +@code # position of the dipole set dipole position = 0, 0.8 + +# surface conductivity between material 1 and material 2 + set sigma = 0.001, 0.2; 0, 0| 0, 0; 0.001, 0.2 @endcode -If the PML strength is set to 0, we see the following output (in order from left to right, the images represent the imaginary components of $E_x$ and $E_y$ and the real components of $E_x$ and $E_y$): +Once again, we will visualize the output with absorbing boundary conditions and PML strength 0 and with absorbing boundary conditions and PML strength 4. The following tables are the imaginary part of $E_x$ and the real part of $E_x$. + - - - + + - -
- Visualization of the solution of step-81 with an interface and PML strength 0 - - Visualization of the solution of step-81 with an interface and PML strength 0 + + Visualization of the solution of step-81 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 0 - Visualization of the solution of step-81 with an interface and PML strength 0 - - Visualization of the solution of step-81 with an interface and PML strength 0 + Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4
-Once again, we see significant distortion and resonance but the standing wave is visible. Now, we increase the PML strength to 4 and generate the following images (in order from left to right, the images represent the imaginary components of $E_x$ and $E_y$ and the real components of $E_x$ and $E_y$): + - - - + + - -
- Visualization of the solution of step-81 with an interface and PML strength 4 - - Visualization of the solution of step-81 with an interface and PML strength 4 + + Visualization of the solution of step-81 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 0 - Visualization of the solution of step-81 with an interface and PML strength 4 - - Visualization of the solution of step-81 with an interface and PML strength 4 + Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4
-As expected, the PML removes the distortion and resonance. The standing wave is also dissipating and getting absorbed within the PML, and as we increase the PML strength, the standing wave will dissipate more within the PML ring. +The SPP is confined near the interface that we created, however without absorbing boundary conditions, we don't observe a dissipation effect. On adding the absorbing boundary conditions, we observe distortion and resonance and we still don't notice any dissipation. As expected, the PML removes the distortion and resonance. The standing wave is also dissipating and getting absorbed within the PML, and as we increase the PML strength, the standing wave will dissipate more within the PML ring. +

Notes

+

Real and Complex Matrices

+As is evident from the results, we are splitting our solution matrices into the real and the imaginary components. We started off using the $H^{curl}$ conforming Nédélec Elements, and we made two copies of the Finite Elements in order +to represent the real and the imaginary components of our input (FE_NedelecSZ was used instead of FE_Nedelec to avoid the sign conflicts issues present in traditional Nédélec elements). In the assembly, we create two vectors of dimension $dim$ that assist us in extracting the real and the imaginary components of our finite elements. +

Rotations and Scaling

+As we see in our assembly, our finite element is rotated and scaled as follows: + +@code +const auto phi_i = real_part.value(i, q_point) - 1.0i * imag_part.value(i, q_point); +@endcode + + +This $\phi_i$ variable doesn't need to be scaled in this way, we may choose any arbitrary scaling consents $a$ and $b$. If we choose this scaling, the $\phi_j$ must also be modified with the same scaling, as follows: + +@code +const auto phi_i = a*real_part.value(i, q_point) - + bi * imag_part.value(i, q_point); + +const auto phi_j = a*real_part.value(i, q_point) + + bi * imag_part.value(i, q_point); +@endcode + +Moreover, the cell_rhs need not be the real part of the rhs_value. Say if we modify to take the imaginary part of the computed rhs_value, we must also modify the cell_matrix accordingly to take the imaginary part of temp. However, making these changes to both sides of the equation will not affect our solution, and we will still be able to generate the surface plasmon polariton. + +@code +cell_rhs(i) += rhs_value.imag(); + +cell_matrix(i) += temp.imag(); +@endcode + +

Postprocessing

+We will create a video demonstrating the wave in motion, which is essentially an implementation of $e^{-i\omega t}(Re(E) + i*Im(E))$ as we increment time. This is done by slightly changing the output function to generate a series of .vtk files, which will represent out solution wave as we increment time. Introduce an input variable $t$ in the output_results() class as output_results(unsigned int t). Then change the class itself to the following: + +@code +template +void Maxwell::output_results(unsigned int t) +{ + std::cout << "Running step:" << alpha << std::endl; + DataOut<2> data_out; + data_out.attach_dof_handler(dof_handler); + Vector postprocessed; + postprocessed.reinit(solution); + for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) { + if (i%4 == 0){ + postprocessed[i] = std::cos(2*M_PI* 0.04* t)*solution[i] + - std::sin(2*M_PI * 0.04 * t)*solution[i+1]; + } else if (i%4 == 2) { + postprocessed[i] = std::cos(2*M_PI * 0.04 * t)*solution[i] + - std::sin(2*M_PI * 0.04 * t)*solution[i+1]; + } + } + data_out.add_data_vector(postprocessed, {"E_x","E_y","null0","null1"}); + data_out.build_patches(); + const std::string filename = + "solution-" + Utilities::int_to_string(alpha) + ".vtk"; + std::ofstream output(filename); + data_out.write_vtk(output); + std::cout << "Done running step:" << alpha << std::endl; + +@endcode + +Finally, in the run() function, replace output_results() with +@code +for(int t = 0; t<=100; t++){ + output_results(t); +} +@endcode + +This would generate 100 solution .vtk files, which can be opened in a group on Paraview and then can be saved as an animation. We used FFMPEG to generate gifs. +

Resulting videos

+Following are the resulting videos of our experiments. As we see, TODO +

Possibilities for Extension

+The current program doesn't allow for iterative solvers as the solutions will not converge with an iterative solver. One possible direction for future work is to implement an iterative solver and involve more preconditioners. An advantage of iterative solvers is the more efficient memory usage, and our current memory usage does not allow for a large number of DOFs. +Another possible direction would be to perform Local Mesh Refinement (instead of Global Mesh Refinement). This will also help us visualize more DOFs in a more memory and time efficient way. diff --git a/examples/step-81/step-81.cc b/examples/step-81/step-81.cc index a999e6850c..e318b4d5f1 100644 --- a/examples/step-81/step-81.cc +++ b/examples/step-81/step-81.cc @@ -380,7 +380,7 @@ namespace Step81 unsigned int refinements; unsigned int fe_order; unsigned int quadrature_order; - bool interface_bool; + bool absorbing_boundary; void parse_parameters_callback(); void make_grid(); @@ -408,8 +408,9 @@ namespace Step81 // @sect4{The Constructor} // The Constructor simply consists specifications for the mesh // and the order of the fnite elements. These are editable through - // the .prm file. The interface_bool can be modified to remove the - // interface and the standing wave. + // the .prm file. The absorbing_boundary boolean can be modified to + // remove the absorbing boundary conditions (in which case our boundary + // would be perfectly conducting). template Maxwell::Maxwell() @@ -435,10 +436,10 @@ namespace Step81 quadrature_order, "order of the quadrature"); - interface_bool = true; - add_parameter("interface boolean", - interface_bool, - "is there an interface?"); + absorbing_boundary = true; + add_parameter("absorbing boundary condition", + absorbing_boundary, + "use absorbing boundary conditions?"); } @@ -459,13 +460,18 @@ namespace Step81 GridGenerator::hyper_cube(triangulation, -scaling, scaling); triangulation.refine_global(refinements); - if (interface_bool){ + if (!absorbing_boundary){ + for (auto &face : triangulation.active_face_iterators()) + if (face->at_boundary()) + face->set_boundary_id(1); + }; + for (auto &cell : triangulation.active_cell_iterators()) if (cell->center()[1] > 0.) cell->set_material_id(1); else cell->set_material_id(2); - } + std::cout << "Number of active cells: " << triangulation.n_active_cells() << std::endl; @@ -555,7 +561,7 @@ namespace Step81 update_quadrature_points | update_normal_vectors | update_JxW_values); - + const unsigned int dofs_per_cell = fe->dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -568,20 +574,21 @@ namespace Step81 for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); + FEValuesViews::Vector real_part(fe_values, 0); + FEValuesViews::Vector imag_part(fe_values, dim); + cell_matrix = 0.; cell_rhs = 0.; cell->get_dof_indices(local_dof_indices); + const auto id = cell->material_id(); - FEValuesViews::Vector real_part(fe_values, 0); - FEValuesViews::Vector imag_part(fe_values, dim); const auto &quadrature_points = fe_values.get_quadrature_points(); 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 id = cell->material_id(); const auto inner_radius = perfectly_matched_layer.inner_radius; auto mu_inv = parameters.mu_inv(position, id); @@ -629,21 +636,22 @@ namespace Step81 { if (face->at_boundary()) { - fe_face_values.reinit(cell, face); - FEValuesViews::Vector real_part(fe_face_values, 0); - FEValuesViews::Vector imag_part(fe_face_values, dim); - - for (unsigned int q_point = 0; q_point < n_face_q_points; - ++q_point) - { - const Point position = - quadrature_points[q_point]; - const auto radius = position.norm(); - const auto id = cell->material_id(); - const auto inner_radius = - perfectly_matched_layer.inner_radius; + const auto id = face->boundary_id(); + if (id !=0) + { + fe_face_values.reinit(cell, face); + FEValuesViews::Vector real_part(fe_face_values, 0); + FEValuesViews::Vector imag_part(fe_face_values, dim); + + for (unsigned int q_point = 0; q_point < n_face_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; - if (id == 0){ auto mu_inv = parameters.mu_inv(position, id); auto epsilon = parameters.epsilon(position, id); @@ -682,9 +690,9 @@ namespace Step81 scalar_product((sqrt_prod * phi_j_T), phi_i_T); cell_matrix(i, j) += temp.real(); } /* j */ - } /* i */ - } - } /* q_point */ + } /* i */ + } /* q_point */ + } } else { @@ -760,17 +768,17 @@ namespace Step81 A_direct.vmult(solution, system_rhs); } - // The output is writted into a vtk file with 4 components - template - void Maxwell::output_results() - { +// The output is written into a vtk file with 4 components +template +void Maxwell::output_results() +{ DataOut<2> data_out; data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"}); + data_out.add_data_vector(solution, {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"}); data_out.build_patches(); std::ofstream output("solution.vtk"); data_out.write_vtk(output); - } +} template -- 2.39.5