]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Step-81: result section updated, minor typos fixed
authorManaswinee Bezbaruah <bezba004@tamu.edu>
Thu, 27 Jan 2022 07:24:59 +0000 (01:24 -0600)
committerMatthias Maier <tamiko@43-1.org>
Thu, 26 May 2022 05:16:02 +0000 (00:16 -0500)
examples/step-81/doc/results.dox
examples/step-81/step-81.cc

index a5b174b02916a86758d8449a4b7330591c75d445..db89021b243bf42a6bdb5782ad2980bf5991bfe4 100644 (file)
 <h1>Results</h1>
 
-Using the above code (without the PML) and a forcing term of a Hertzian dipole
-at the center, we have generated the following solution wave. The two complex
-plane solutions are followed by the two solutions in the real plane. Furthermore,
-a significant resonance is observed, causing distorted images and necessitating
-a PML boundary condition. <br />
-TODO: add the images <br />
-The solution of the same problem with a PML of strength 2 of radii 8 and 10
-is shown below (in the same order or complex and real solutions). Clearly,
-the PML significantly reduces the distortion in our solution.<br />
-TODO: add the images<br />
-Additionally, an interface is added at y = 0, and by adjusting the surface
-conductivity value and the position of our dipole, we arrive at a standing
-wave.<br />
-TODO: add the images
+The solution is written to a .vtk file with four components. These are the real and imaginary parts of the $E_x$ and $E_y$ solution waves. With the current setup, the output should read
+
+@code
+Number of active cells: 4096
+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:
+
+@code
+# is there an interface?
+  set interface 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
+@endcode
+
+As we see, a much more coherent image is generated with an appropriate PML.
+
+<table width="80%" align="center">
+  <tr>
+    <td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_nointerface.png" alt="Visualization of the solution of step-81 with no interface and PML strength 0" height="500">
+    </td>
+    <td></td>
+<td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_nointerface.png" alt="Visualization of the solution of step-81 with no interface and PML strength 4" height="500">
+    </td>
+  </tr>
+</table>
+
+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$
+
+@code
+# is there an interface?
+  set interface boolean = true
+
+# position of the dipole
+  set dipole position    = 0, 0.8
+@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$):
+<table width="80%" align="center">
+  <tr>
+    <td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_imagEx.png" alt="Visualization of the solution of step-81 with an interface and PML strength 0" height="300">
+    </td>
+    <td></td>
+<td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_imagEy.png" alt="Visualization of the solution of step-81 with  an interface and PML strength 0" height="300">
+    </td>
+    <td></td>
+    <td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_realEx.png" alt="Visualization of the solution of step-81 with  an interface and PML strength 0" height="300">
+    </td>
+    <td></td>
+<td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML0_realEy.png" alt="Visualization of the solution of step-81 with an interface and PML strength 0" height="300">
+    </td>
+  </tr>
+</table>
+
+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$):
+<table width="80%" align="center">
+  <tr>
+    <td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_imagEx.png" alt="Visualization of the solution of step-81 with an interface and PML strength 4" height="300">
+    </td>
+    <td></td>
+<td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_imagEy.png" alt="Visualization of the solution of step-81 with  an interface and PML strength 4" height="300">
+    </td>
+    <td></td>
+    <td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_realEx.png" alt="Visualization of the solution of step-81 with  an interface and PML strength 4" height="300">
+    </td>
+    <td></td>
+<td align="center">
+      <img src="/Users/manaswineebezbaruah/workspace/dealii/examples/step-81/images/PML4_realEy.png" alt="Visualization of the solution of step-81 with an interface and PML strength 4" height="300">
+    </td>
+  </tr>
+</table>
+
+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.
+
+
+
+
+
+
index 9d6b7fe62ceed70137b16b7ba180cf4d22fbea75..a999e6850c2d15d753b0954443fbfefe92e3826c 100644 (file)
@@ -179,7 +179,7 @@ namespace Step81
     dipole_position = Point<dim>(0., 0.8);
     add_parameter("dipole position",
                   dipole_position,
-                  "posititon of the dipole");
+                  "position of the dipole");
 
     dipole_orientation = Tensor<1, dim, double>{{0., 1.}};
     add_parameter("dipole orientation",
@@ -380,7 +380,7 @@ namespace Step81
     unsigned int refinements;
     unsigned int fe_order;
     unsigned int quadrature_order;
-    unsigned int n_outputs;
+    bool interface_bool;
 
     void parse_parameters_callback();
     void make_grid();
@@ -408,7 +408,8 @@ 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 .prm file. The interface_bool can be modified to remove the
+  // interface and the standing wave.
 
   template <int dim>
   Maxwell<dim>::Maxwell()
@@ -421,7 +422,7 @@ namespace Step81
     scaling = 20;
     add_parameter("scaling", scaling, "scale of the hypercube geometry");
 
-    refinements = 10;
+    refinements = 8;
     add_parameter("refinements",
                   refinements,
                   "number of refinements of the geometry");
@@ -433,6 +434,11 @@ namespace Step81
     add_parameter("quadrature order",
                   quadrature_order,
                   "order of the quadrature");
+        
+    interface_bool = true;
+    add_parameter("interface boolean",
+                    interface_bool,
+                    "is there an interface?");
   }
 
 
@@ -442,18 +448,24 @@ namespace Step81
     fe = std::make_unique<FESystem<dim>>(FE_NedelecSZ<dim>(fe_order), 2);
   }
 
-  // Make the mesh for the domain and generate the triangulation required
+  // Make the mesh for the domain and generate the triangulation required.
+  // Additionally, there is an interface added here to visualize
+  // a standing wave. To generate a solution without any interface,
+  // comment out lines 455-459.
+
   template <int dim>
   void Maxwell<dim>::make_grid()
   {
     GridGenerator::hyper_cube(triangulation, -scaling, scaling);
     triangulation.refine_global(refinements);
-
+    
+    if (interface_bool){
     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;
@@ -472,9 +484,26 @@ namespace Step81
 
     solution.reinit(dof_handler.n_dofs());
     system_rhs.reinit(dof_handler.n_dofs());
+
     constraints.clear();
+
     DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+    VectorTools::project_boundary_values_curl_conforming_l2(
+      dof_handler,
+      0, /* real part */
+      dealii::ZeroFunction<dim>(2 * dim),
+      0, /* boundary id */
+      constraints);
+    VectorTools::project_boundary_values_curl_conforming_l2(
+      dof_handler,
+      dim, /* imaginary part */
+      dealii::ZeroFunction<dim>(2 * dim),
+      0, /* boundary id */
+      constraints);
+
     constraints.close();
+
     DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
     DoFTools::make_sparsity_pattern(dof_handler,
                                     dsp,
@@ -614,47 +643,48 @@ namespace Step81
                     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);
 
-                    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_tensor(position);
-
-                        mu_inv  = mu_inv / d;
-                        epsilon = invert(A) * epsilon * invert(B);
-                      };
+                        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_tensor(position);
 
-                    const auto normal = fe_face_values.normal_vector(q_point);
+                            mu_inv  = mu_inv / d;
+                            epsilon = invert(A) * epsilon * invert(B);
+                          };
 
-                    for (unsigned int i = 0; i < dofs_per_cell; ++i)
-                      {
-                        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);
+                        const auto normal = fe_face_values.normal_vector(q_point);
 
-                        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                        for (unsigned int i = 0; i < dofs_per_cell; ++i)
                           {
-                            const auto phi_j =
-                              real_part.value(j, q_point) +
-                              1.0i * imag_part.value(j, q_point);
-                            const auto phi_j_T =
-                              tangential_part(phi_j, normal) *
-                              fe_face_values.JxW(q_point);
-
-                            const auto prod      = mu_inv * epsilon;
-                            const auto sqrt_prod = prod;
-
-                            const auto temp =
-                              -1.0i *
-                              scalar_product((sqrt_prod * phi_j_T), phi_i_T);
-                            cell_matrix(i, j) += temp.real();
-                          } /* j */
-                      }     /* i */
-                  }         /* q_point */
+                            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)
+                              {
+                                const auto phi_j =
+                                  real_part.value(j, q_point) +
+                                  1.0i * imag_part.value(j, q_point);
+                                const auto phi_j_T =
+                                  tangential_part(phi_j, normal) *
+                                  fe_face_values.JxW(q_point);
+
+                                const auto prod      = mu_inv * epsilon;
+                                const auto sqrt_prod = prod;
+
+                                const auto temp =
+                                  -1.0i *
+                                  scalar_product((sqrt_prod * phi_j_T), phi_i_T);
+                                cell_matrix(i, j) += temp.real();
+                              } /* j */
+                          }     /* i */
+                        }
+                    }         /* q_point */
               }
             else
               {
@@ -721,7 +751,7 @@ namespace Step81
       }
   };
 
-
+  // We use a direct solver from the SparseDirectUMFPACK to solve the system
   template <int dim>
   void Maxwell<dim>::solve()
   {
@@ -730,7 +760,7 @@ namespace Step81
     A_direct.vmult(solution, system_rhs);
   }
 
-
+  // The output is writted into a vtk file with 4 components
   template <int dim>
   void Maxwell<dim>::output_results()
   {

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.