Line data    Source code

       1              : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
       2              : // SPDX-License-Identifier: Apache-2.0
       3              : 
       4              : #include "farfieldboundaryoperator.hpp"
       5              : 
       6              : #include <set>
       7              : #include "linalg/densematrix.hpp"
       8              : #include "models/materialoperator.hpp"
       9              : #include "utils/communication.hpp"
      10              : #include "utils/geodata.hpp"
      11              : #include "utils/iodata.hpp"
      12              : #include "utils/prettyprint.hpp"
      13              : 
      14              : namespace palace
      15              : {
      16              : 
      17           17 : FarfieldBoundaryOperator::FarfieldBoundaryOperator(const IoData &iodata,
      18              :                                                    const MaterialOperator &mat_op,
      19           17 :                                                    const mfem::ParMesh &mesh)
      20           17 :   : mat_op(mat_op), farfield_attr(SetUpBoundaryProperties(iodata, mesh))
      21              : {
      22              :   // Print out BC info for all farfield attributes.
      23           17 :   if (farfield_attr.Size())
      24              :   {
      25            0 :     Mpi::Print("\nConfiguring Robin absorbing BC (order {:d}) at attributes:\n", order);
      26              :     std::sort(farfield_attr.begin(), farfield_attr.end());
      27            0 :     utils::PrettyPrint(farfield_attr);
      28              :   }
      29           17 : }
      30              : 
      31              : mfem::Array<int>
      32           17 : FarfieldBoundaryOperator::SetUpBoundaryProperties(const IoData &iodata,
      33              :                                                   const mfem::ParMesh &mesh)
      34              : {
      35              :   // Check that impedance boundary attributes have been specified correctly.
      36           17 :   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
      37              :   mfem::Array<int> bdr_attr_marker;
      38           17 :   if (!iodata.boundaries.farfield.empty())
      39              :   {
      40              :     bdr_attr_marker.SetSize(bdr_attr_max);
      41              :     bdr_attr_marker = 0;
      42            0 :     for (auto attr : mesh.bdr_attributes)
      43              :     {
      44            0 :       bdr_attr_marker[attr - 1] = 1;
      45              :     }
      46              :     std::set<int> bdr_warn_list;
      47            0 :     for (auto attr : iodata.boundaries.farfield.attributes)
      48              :     {
      49              :       // MFEM_VERIFY(attr > 0 && attr <= bdr_attr_max,
      50              :       //             "Absorbing boundary attribute tags must be non-negative and correspond
      51              :       //             " "to attributes in the mesh!");
      52              :       // MFEM_VERIFY(bdr_attr_marker[attr - 1],
      53              :       //             "Unknown absorbing boundary attribute " << attr << "!");
      54            0 :       if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
      55              :       {
      56              :         bdr_warn_list.insert(attr);
      57              :       }
      58            0 :       if (!bdr_warn_list.empty())
      59              :       {
      60            0 :         Mpi::Print("\n");
      61            0 :         Mpi::Warning(
      62              :             "Unknown absorbing boundary attributes!\nSolver will just ignore them!");
      63            0 :         utils::PrettyPrint(bdr_warn_list, "Boundary attribute list:");
      64            0 :         Mpi::Print("\n");
      65              :       }
      66              :     }
      67              :   }
      68              : 
      69              :   // Set the order of the farfield boundary condition.
      70           17 :   order = iodata.boundaries.farfield.order;
      71              : 
      72              :   // Mark selected boundary attributes from the mesh as farfield.
      73              :   mfem::Array<int> farfield_bcs;
      74           17 :   farfield_bcs.Reserve(static_cast<int>(iodata.boundaries.farfield.attributes.size()));
      75           17 :   for (auto attr : iodata.boundaries.farfield.attributes)
      76              :   {
      77            0 :     if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
      78              :     {
      79            0 :       continue;  // Can just ignore if wrong
      80              :     }
      81            0 :     farfield_bcs.Append(attr);
      82              :   }
      83           17 :   MFEM_VERIFY(farfield_bcs.Size() == 0 || order < 2 ||
      84              :                   iodata.problem.type == ProblemType::DRIVEN ||
      85              :                   iodata.problem.type == ProblemType::EIGENMODE,
      86              :               "Second-order farfield boundaries are only available for frequency "
      87              :               "domain simulations!");
      88           17 :   return farfield_bcs;
      89              : }
      90              : 
      91            0 : void FarfieldBoundaryOperator::AddDampingBdrCoefficients(double coeff,
      92              :                                                          MaterialPropertyCoefficient &fb)
      93              : {
      94              :   // First-order absorbing boundary condition.
      95            0 :   if (farfield_attr.Size())
      96              :   {
      97            0 :     MaterialPropertyCoefficient invz0_func(mat_op.GetBdrAttributeToMaterial(),
      98            0 :                                            mat_op.GetInvImpedance());
      99            0 :     invz0_func.RestrictCoefficient(mat_op.GetCeedBdrAttributes(farfield_attr));
     100            0 :     fb.AddCoefficient(invz0_func.GetAttributeToMaterial(),
     101              :                       invz0_func.GetMaterialProperties(), coeff);
     102            0 :   }
     103            0 : }
     104              : 
     105            0 : void FarfieldBoundaryOperator::AddExtraSystemBdrCoefficients(
     106              :     double omega, MaterialPropertyCoefficient &dfbr, MaterialPropertyCoefficient &dfbi)
     107              : {
     108              :   // Contribution for second-order absorbing BC. See Jin Section 9.3 for reference. The β
     109              :   // coefficient for the second-order ABC is 1/(2ik+2/r). Taking the radius of curvature
     110              :   // as infinity (plane wave scattering), the r-dependence vanishes and the contribution
     111              :   // is purely imaginary. Multiplying through by μ⁻¹ we get the material coefficient to ω
     112              :   // as 1 / (μ √(με)). Also, this implementation ignores the divergence term ∇⋅Eₜ, as
     113              :   // COMSOL does as well.
     114            0 :   if (farfield_attr.Size() && order > 1)
     115              :   {
     116              :     mfem::DenseTensor muinvc0 =
     117            0 :         linalg::Mult(mat_op.GetInvPermeability(), mat_op.GetLightSpeed());
     118            0 :     MaterialPropertyCoefficient muinvc0_func(mat_op.GetBdrAttributeToMaterial(), muinvc0);
     119            0 :     muinvc0_func.RestrictCoefficient(mat_op.GetCeedBdrAttributes(farfield_attr));
     120              : 
     121              :     // Instead getting the correct normal of farfield boundary elements, just pick the
     122              :     // the first element normal. This is fine as long as the farfield material properties
     123              :     // are not anisotropic.
     124            0 :     mfem::Vector normal(mat_op.SpaceDimension());
     125            0 :     normal = 0.0;
     126            0 :     normal(0) = 1.0;
     127            0 :     muinvc0_func.NormalProjectedCoefficient(normal);
     128              : 
     129            0 :     dfbi.AddCoefficient(muinvc0_func.GetAttributeToMaterial(),
     130              :                         muinvc0_func.GetMaterialProperties(), 0.5 / omega);
     131            0 :   }
     132            0 : }
     133              : 
     134              : }  // namespace palace