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 "surfaceconductivityoperator.hpp"
5 :
6 : #include <set>
7 : #include "models/materialoperator.hpp"
8 : #include "utils/communication.hpp"
9 : #include "utils/geodata.hpp"
10 : #include "utils/iodata.hpp"
11 : #include "utils/prettyprint.hpp"
12 :
13 : namespace palace
14 : {
15 :
16 : using namespace std::complex_literals;
17 :
18 17 : SurfaceConductivityOperator::SurfaceConductivityOperator(const IoData &iodata,
19 : const MaterialOperator &mat_op,
20 17 : const mfem::ParMesh &mesh)
21 17 : : mat_op(mat_op)
22 : {
23 : // Print out BC info for all finite conductivity boundary attributes.
24 17 : SetUpBoundaryProperties(iodata, mesh);
25 17 : PrintBoundaryInfo(iodata, mesh);
26 17 : }
27 :
28 17 : void SurfaceConductivityOperator::SetUpBoundaryProperties(const IoData &iodata,
29 : const mfem::ParMesh &mesh)
30 : {
31 : // Check that conductivity boundary attributes have been specified correctly.
32 17 : int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
33 : mfem::Array<int> bdr_attr_marker;
34 17 : if (!iodata.boundaries.conductivity.empty())
35 : {
36 : mfem::Array<int> conductivity_marker(bdr_attr_max);
37 : bdr_attr_marker.SetSize(bdr_attr_max);
38 : bdr_attr_marker = 0;
39 : conductivity_marker = 0;
40 0 : for (auto attr : mesh.bdr_attributes)
41 : {
42 0 : bdr_attr_marker[attr - 1] = 1;
43 : }
44 : std::set<int> bdr_warn_list;
45 0 : for (const auto &data : iodata.boundaries.conductivity)
46 : {
47 0 : for (auto attr : data.attributes)
48 : {
49 0 : MFEM_VERIFY(!conductivity_marker[attr - 1],
50 : "Multiple definitions of conductivity boundary properties for boundary "
51 : "attribute "
52 : << attr << "!");
53 0 : conductivity_marker[attr - 1] = 1;
54 : // MFEM_VERIFY(attr > 0 && attr <= bdr_attr_max,
55 : // "Conductivity boundary attribute tags must be non-negative and "
56 : // "correspond to attributes in the mesh!");
57 : // MFEM_VERIFY(bdr_attr_marker[attr - 1],
58 : // "Unknown conductivity boundary attribute " << attr << "!");
59 0 : if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
60 : {
61 : bdr_warn_list.insert(attr);
62 : }
63 : }
64 : }
65 0 : if (!bdr_warn_list.empty())
66 : {
67 0 : Mpi::Print("\n");
68 0 : Mpi::Warning(
69 : "Unknown conductivity boundary attributes!\nSolver will just ignore them!");
70 0 : utils::PrettyPrint(bdr_warn_list, "Boundary attribute list:");
71 0 : Mpi::Print("\n");
72 : }
73 : }
74 :
75 : // Finite conductivity boundaries are defined using the user provided surface conductivity
76 : // and optionally conductor thickness.
77 17 : boundaries.reserve(iodata.boundaries.conductivity.size());
78 17 : for (const auto &data : iodata.boundaries.conductivity)
79 : {
80 0 : MFEM_VERIFY(data.sigma > 0.0 && data.mu_r > 0.0,
81 : "Conductivity boundary has no conductivity or no "
82 : "permeability defined!");
83 0 : MFEM_VERIFY(data.h >= 0.0, "Conductivity boundary should have non-negative thickness!");
84 0 : auto &bdr = boundaries.emplace_back();
85 0 : bdr.sigma = data.sigma;
86 0 : bdr.mu = data.mu_r;
87 0 : bdr.h = data.h;
88 0 : if (data.external)
89 : {
90 : // External surfaces have twice the effective thickness since the BC is applied at one
91 : // side.
92 0 : bdr.h *= 2.0;
93 : }
94 0 : bdr.attr_list.Reserve(static_cast<int>(data.attributes.size()));
95 0 : for (auto attr : data.attributes)
96 : {
97 0 : if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
98 : {
99 0 : continue; // Can just ignore if wrong
100 : }
101 0 : bdr.attr_list.Append(attr);
102 : }
103 : }
104 17 : MFEM_VERIFY(boundaries.empty() || iodata.problem.type == ProblemType::DRIVEN ||
105 : iodata.problem.type == ProblemType::EIGENMODE,
106 : "Finite conductivity boundaries are only available for frequency "
107 : "domain simulations!");
108 17 : }
109 :
110 17 : void SurfaceConductivityOperator::PrintBoundaryInfo(const IoData &iodata,
111 : const mfem::ParMesh &mesh)
112 : {
113 17 : if (boundaries.empty())
114 : {
115 : return;
116 : }
117 0 : Mpi::Print("\nConfiguring Robin finite conductivity BC at attributes:\n");
118 0 : for (const auto &bdr : boundaries)
119 : {
120 0 : for (auto attr : bdr.attr_list)
121 : {
122 0 : Mpi::Print(" {:d}: σ = {:.3e} S/m", attr,
123 0 : iodata.units.Dimensionalize<Units::ValueType::CONDUCTIVITY>(bdr.sigma));
124 0 : if (bdr.h > 0.0)
125 : {
126 0 : Mpi::Print(", h = {:.3e} m",
127 0 : iodata.units.Dimensionalize<Units::ValueType::LENGTH>(bdr.h));
128 : }
129 0 : Mpi::Print(", n = ({:+.1f})\n", fmt::join(mesh::GetSurfaceNormal(mesh, attr), ","));
130 : }
131 : }
132 : }
133 :
134 17 : mfem::Array<int> SurfaceConductivityOperator::GetAttrList() const
135 : {
136 : mfem::Array<int> attr_list;
137 17 : for (const auto &bdr : boundaries)
138 : {
139 : attr_list.Append(bdr.attr_list);
140 : }
141 17 : return attr_list;
142 : }
143 :
144 0 : void SurfaceConductivityOperator::AddExtraSystemBdrCoefficients(
145 : double omega, MaterialPropertyCoefficient &fbr, MaterialPropertyCoefficient &fbi)
146 : {
147 : // If the provided conductor thickness is empty (zero), prescribe a surface impedance
148 : // (1+i)/σδ, where δ is the skin depth. If it is nonzero, use a finite thickness
149 : // modification which correctly produces the DC limit when h << δ. See the Ansys HFSS
150 : // user manual section titled "Surface Impedance Boundary Condition for Metal Traces of
151 : // Finite Thickness."
152 0 : for (const auto &bdr : boundaries)
153 : {
154 0 : if (std::abs(bdr.sigma) > 0.0)
155 : {
156 0 : double delta = std::sqrt(2.0 / (bdr.mu * bdr.sigma * omega));
157 0 : std::complex<double> Z = 1.0 / (bdr.sigma * delta);
158 : Z.imag(Z.real());
159 0 : if (bdr.h > 0.0)
160 : {
161 0 : double nu = bdr.h / delta;
162 0 : double den = std::cosh(nu) - std::cos(nu);
163 0 : Z.real(Z.real() * (std::sinh(nu) + std::sin(nu)) / den);
164 0 : Z.imag(Z.imag() * (std::sinh(nu) - std::sin(nu)) / den);
165 : }
166 : // The BC term has coefficient iω/Z (like for standard lumped surface impedance).
167 : std::complex<double> s(1i * omega / Z);
168 0 : fbr.AddMaterialProperty(mat_op.GetCeedBdrAttributes(bdr.attr_list), s.real());
169 0 : fbi.AddMaterialProperty(mat_op.GetCeedBdrAttributes(bdr.attr_list), s.imag());
170 : }
171 : }
172 0 : }
173 :
174 : } // namespace palace
|
