From: Manaswinee Bezbaruah Date: Thu, 27 Jan 2022 07:24:59 +0000 (-0600) Subject: Step-81: result section updated, minor typos fixed X-Git-Tag: v9.4.0-rc1~136^2~19 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c75701fd7d597b17614c5acba3afacdff03a215f;p=dealii.git Step-81: result section updated, minor typos fixed --- diff --git a/examples/step-81/doc/results.dox b/examples/step-81/doc/results.dox index a5b174b029..db89021b24 100644 --- a/examples/step-81/doc/results.dox +++ b/examples/step-81/doc/results.dox @@ -1,16 +1,103 @@

Results

-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.
-TODO: add the images
-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.
-TODO: add the images
-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.
-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. + + + + + + + +
+ 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 +
+ +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$): + + + + + + + + + + +
+ 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 an interface and PML strength 0 + + Visualization of the solution of step-81 with an interface and PML strength 0 +
+ +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 an interface and PML strength 4 + + Visualization of the solution of step-81 with an interface 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. + + + + + + diff --git a/examples/step-81/step-81.cc b/examples/step-81/step-81.cc index 9d6b7fe62c..a999e6850c 100644 --- a/examples/step-81/step-81.cc +++ b/examples/step-81/step-81.cc @@ -179,7 +179,7 @@ namespace Step81 dipole_position = Point(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 Maxwell::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>(FE_NedelecSZ(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 void Maxwell::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(2 * dim), + 0, /* boundary id */ + constraints); + VectorTools::project_boundary_values_curl_conforming_l2( + dof_handler, + dim, /* imaginary part */ + dealii::ZeroFunction(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 void Maxwell::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 void Maxwell::output_results() {