1
2 """Models describing absorbing boundary layers."""
3
4 from __future__ import division
5
6 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
7
8 __license__ = """
9 This program is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with this program. If not, see U{http://www.gnu.org/licenses/}.
21 """
22
23
24
25
26 import numpy
27 import numpy.linalg as la
28
29 from pytools import memoize_method, Record
30 from hedge.models.em import \
31 MaxwellOperator, \
32 TMMaxwellOperator, \
33 TEMaxwellOperator
39 """Implements a PML as in
40
41 [1] S. Abarbanel and D. Gottlieb, "On the construction and analysis of absorbing
42 layers in CEM," Applied Numerical Mathematics, vol. 27, 1998, S. 331-340.
43 (eq 3.7-3.11)
44
45 [2] E. Turkel and A. Yefet, "Absorbing PML
46 boundary layers for wave-like equations,"
47 Applied Numerical Mathematics, vol. 27,
48 1998, S. 533-557.
49 (eq. 4.10)
50
51 [3] Abarbanel, D. Gottlieb, and J.S. Hesthaven, "Long Time Behavior of the
52 Perfectly Matched Layer Equations in Computational Electromagnetics,"
53 Journal of Scientific Computing, vol. 17, Dez. 2002, S. 405-422.
54
55 Generalized to 3D in doc/maxima/abarbanel-pml.mac.
56 """
57
59 __slots__ = ["sigma", "sigma_prime", "tau"]
60
61
63 return self.__class__(
64 **dict((name, f(getattr(self, name)))
65 for name in self.fields))
66
70
72 sub_e, sub_h, sub_p, sub_q = self.split_ehpq(w)
73
74 e_subset = self.get_eh_subset()[0:3]
75 h_subset = self.get_eh_subset()[3:6]
76 dim_subset = (True,) * self.dimensions + (False,) * (3-self.dimensions)
77
78 def pad_vec(v, subset):
79 result = numpy.zeros((3,), dtype=object)
80 result[numpy.array(subset, dtype=bool)] = v
81 return result
82
83 from hedge.optemplate import make_vector_field
84 sig = pad_vec(
85 make_vector_field("sigma", self.dimensions),
86 dim_subset)
87 sig_prime = pad_vec(
88 make_vector_field("sigma_prime", self.dimensions),
89 dim_subset)
90 if self.add_decay:
91 tau = pad_vec(
92 make_vector_field("tau", self.dimensions),
93 dim_subset)
94 else:
95 tau = numpy.zeros((3,))
96
97 e = pad_vec(sub_e, e_subset)
98 h = pad_vec(sub_h, h_subset)
99 p = pad_vec(sub_p, dim_subset)
100 q = pad_vec(sub_q, dim_subset)
101
102 rhs = numpy.zeros(12, dtype=object)
103
104 for mx in range(3):
105 my = (mx+1) % 3
106 mz = (mx+2) % 3
107
108 from hedge.tools import levi_civita
109 assert levi_civita((mx,my,mz)) == 1
110
111 rhs[mx] += -sig[my]/self.epsilon*(2*e[mx]+p[mx]) - 2*tau[my]/self.epsilon*e[mx]
112 rhs[my] += -sig[mx]/self.epsilon*(2*e[my]+p[my]) - 2*tau[mx]/self.epsilon*e[my]
113 rhs[3+mz] += 1/(self.epsilon*self.mu) * (
114 sig_prime[mx] * q[mx] - sig_prime[my] * q[my])
115
116 rhs[6+mx] += sig[my]/self.epsilon*e[mx]
117 rhs[6+my] += sig[mx]/self.epsilon*e[my]
118 rhs[9+mx] += -sig[mx]/self.epsilon*q[mx] - (e[my] + e[mz])
119
120 from hedge.tools import full_to_subset_indices
121 sub_idx = full_to_subset_indices(e_subset+h_subset+dim_subset+dim_subset)
122
123 return rhs[sub_idx]
124
137
138 - def bind(self, discr, coefficients):
143
144 - def assemble_ehpq(self, e=None, h=None, p=None, q=None, discr=None):
145 if discr is None:
146 def zero(): return 0
147 else:
148 def zero(): return discr.volume_zeros()
149
150 from hedge.tools import count_subset
151 e_components = count_subset(self.get_eh_subset()[0:3])
152 h_components = count_subset(self.get_eh_subset()[3:6])
153
154 def default_fld(fld, comp):
155 if fld is None:
156 return [zero() for i in xrange(comp)]
157 else:
158 return fld
159
160 e = default_fld(e, e_components)
161 h = default_fld(h, h_components)
162 p = default_fld(p, self.dimensions)
163 q = default_fld(q, self.dimensions)
164
165 from hedge.tools import join_fields
166 return join_fields(e, h, p, q)
167
168 @memoize_method
170 e_subset = self.get_eh_subset()[0:3]
171 h_subset = self.get_eh_subset()[3:6]
172
173 dim_subset = [True] * self.dimensions + [False] * (3-self.dimensions)
174
175 from hedge.tools import partial_to_all_subset_indices
176 return tuple(partial_to_all_subset_indices(
177 [e_subset, h_subset, dim_subset, dim_subset]))
178
180 e_idx, h_idx, p_idx, q_idx = self.partial_to_ehpq_subsets()
181 e, h, p, q = w[e_idx], w[h_idx], w[p_idx], w[q_idx]
182
183 from hedge.flux import FluxVectorPlaceholder as FVP
184 if isinstance(w, FVP):
185 return FVP(scalars=e), FVP(scalars=h)
186 else:
187 from hedge.tools import make_obj_array as moa
188 return moa(e), moa(h), moa(p), moa(q)
189
190
193 assert o_min < i_min <= i_max < o_max
194
195 if o_min != i_min:
196 l_dist = (i_min - node_coord) / (i_min-o_min)
197 l_dist_prime = discr.volume_zeros(kind="numpy", dtype=node_coord.dtype)
198 l_dist_prime[l_dist >= 0] = -1 / (i_min-o_min)
199 l_dist[l_dist < 0] = 0
200 else:
201 l_dist = l_dist_prime = numpy.zeros_like(node_coord)
202
203 if i_max != o_max:
204 r_dist = (node_coord - i_max) / (o_max-i_max)
205 r_dist_prime = discr.volume_zeros(kind="numpy", dtype=node_coord.dtype)
206 r_dist_prime[r_dist >= 0] = 1 / (o_max-i_max)
207 r_dist[r_dist < 0] = 0
208 else:
209 r_dist = r_dist_prime = numpy.zeros_like(node_coord)
210
211 l_plus_r = l_dist+r_dist
212 return l_plus_r**exponent, \
213 (l_dist_prime+r_dist_prime)*exponent*l_plus_r**(exponent-1), \
214 l_plus_r
215
216 - def coefficients_from_boxes(self, discr,
217 inner_bbox, outer_bbox=None,
218 magnitude=None, tau_magnitude=None,
219 exponent=None, dtype=None):
220 if outer_bbox is None:
221 outer_bbox = discr.mesh.bounding_box()
222
223 if exponent is None:
224 exponent = 2
225
226 if magnitude is None:
227 magnitude = 20
228
229 if tau_magnitude is None:
230 tau_magnitude = 0.4
231
232
233 from math import sqrt
234 magnitude = magnitude*sqrt(self.epsilon/self.mu)
235 tau_magnitude = tau_magnitude*sqrt(self.epsilon/self.mu)
236
237 i_min, i_max = inner_bbox
238 o_min, o_max = outer_bbox
239
240 from hedge.tools import make_obj_array
241
242 nodes = discr.nodes
243 if dtype is not None:
244 nodes = nodes.astype(dtype)
245
246 sigma, sigma_prime, tau = zip(*[self._construct_scalar_coefficients(
247 discr, nodes[:,i],
248 i_min[i], i_max[i], o_min[i], o_max[i],
249 exponent)
250 for i in range(discr.dimensions)])
251
252 return self.PMLCoefficients(
253 sigma=magnitude*make_obj_array(sigma),
254 sigma_prime=magnitude*make_obj_array(sigma_prime),
255 tau=tau_magnitude*make_obj_array(tau))
256
257 - def coefficients_from_width(self, discr, width,
258 magnitude=None, tau_magnitude=None, exponent=None,
259 dtype=None):
265
273
278