hedge.models.em

1 # -*- coding: utf8 -*- 2 """Hedge operators modelling electromagnetic phenomena.""" 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 from pytools import memoize_method 26 27 import hedge.mesh 28 from hedge.models import TimeDependentOperator

29 30 31 32 33 -class MaxwellOperator(TimeDependentOperator):

34 """A 3D Maxwell operator with PEC boundaries. 35 36 Field order is [Ex Ey Ez Hx Hy Hz]. 37 """ 38 39 _default_dimensions = 3 40

41 - def __init__(self, epsilon, mu, 42 flux_type, 43 bdry_flux_type=None, 44 pec_tag=hedge.mesh.TAG_ALL, 45 absorb_tag=hedge.mesh.TAG_NONE, 46 incident_tag=hedge.mesh.TAG_NONE, 47 incident_bc=None, current=None, dimensions=None):

48 """ 49 @arg flux_type: can be in [0,1] for anything between central and upwind, 50 or "lf" for Lax-Friedrichs. 51 """ 52 53 self.dimensions = dimensions or self._default_dimensions 54 55 space_subset = [True]*self.dimensions + [False]*(3-self.dimensions) 56 57 e_subset = self.get_eh_subset()[0:3] 58 h_subset = self.get_eh_subset()[3:6] 59 60 from hedge.tools import SubsettableCrossProduct 61 self.space_cross_e = SubsettableCrossProduct( 62 op1_subset=space_subset, 63 op2_subset=e_subset, 64 result_subset=h_subset) 65 self.space_cross_h = SubsettableCrossProduct( 66 op1_subset=space_subset, 67 op2_subset=h_subset, 68 result_subset=e_subset) 69 70 from math import sqrt 71 72 self.epsilon = epsilon 73 self.mu = mu 74 self.c = 1/sqrt(mu*epsilon) 75 76 self.Z = sqrt(mu/epsilon) 77 self.Y = 1/self.Z 78 79 self.flux_type = flux_type 80 if bdry_flux_type is None: 81 self.bdry_flux_type = flux_type 82 else: 83 self.bdry_flux_type = bdry_flux_type 84 85 self.pec_tag = pec_tag 86 self.absorb_tag = absorb_tag 87 self.incident_tag = incident_tag 88 89 self.current = current 90 self.incident_bc_data = incident_bc

91

92 - def flux(self, flux_type):

93 from hedge.flux import make_normal, FluxVectorPlaceholder 94 from hedge.tools import join_fields 95 96 normal = make_normal(self.dimensions) 97 98 from hedge.tools import count_subset 99 w = FluxVectorPlaceholder(count_subset(self.get_eh_subset())) 100 e, h = self.split_eh(w) 101 102 if flux_type == "lf": 103 return join_fields( 104 # flux e, 105 1/2*( 106 -1/self.epsilon*self.space_cross_h(normal, h.int-h.ext) 107 -self.c/2*(e.int-e.ext) 108 ), 109 # flux h 110 1/2*( 111 1/self.mu*self.space_cross_e(normal, e.int-e.ext) 112 -self.c/2*(h.int-h.ext)) 113 ) 114 elif isinstance(flux_type, (int, float)): 115 # see doc/maxima/maxwell.mac 116 return join_fields( 117 # flux e, 118 1/self.epsilon*( 119 -1/2*self.space_cross_h(normal, 120 h.int-h.ext 121 -flux_type/self.Z*self.space_cross_e( 122 normal, e.int-e.ext)) 123 ), 124 # flux h 125 1/self.mu*( 126 1/2*self.space_cross_e(normal, 127 e.int-e.ext 128 +flux_type/(self.Y)*self.space_cross_h( 129 normal, h.int-h.ext)) 130 ), 131 ) 132 else: 133 raise ValueError, "maxwell: invalid flux_type (%s)" % self.flux_type

134

135 - def local_derivatives(self, w=None):

136 e, h = self.split_eh(self.field_placeholder(w)) 137 138 def e_curl(field): 139 return self.space_cross_e(nabla, field)

140 141 def h_curl(field): 142 return self.space_cross_h(nabla, field)

143 144 from hedge.optemplate import make_nabla 145 from hedge.tools import join_fields, count_subset 146 147 nabla = make_nabla(self.dimensions) 148 149 if self.current is not None: 150 from hedge.optemplate import make_vector_field 151 j = make_vector_field("j", 152 count_subset(self.get_eh_subset()[:3])) 153 else: 154 j = 0 155 156 # in conservation form: u_t + A u_x = 0 157 return join_fields( 158 1/self.epsilon * (j - h_curl(h)), 159 1/self.mu * e_curl(e), 160 ) 161 162 from hedge.optemplate import pair_with_boundary, \ 163 InverseMassOperator, get_flux_operator, \ 164 BoundarizeOperator 165

166 - def field_placeholder(self, w=None):

167 from hedge.tools import count_subset 168 fld_cnt = count_subset(self.get_eh_subset()) 169 if w is None: 170 from hedge.optemplate import make_vector_field 171 w = make_vector_field("w", fld_cnt) 172 173 return w

174

175 - def pec_bc(self, w=None):

176 e, h = self.split_eh(self.field_placeholder(w)) 177 from hedge.tools import join_fields 178 179 from hedge.optemplate import BoundarizeOperator 180 pec_e = BoundarizeOperator(self.pec_tag) * e 181 pec_h = BoundarizeOperator(self.pec_tag) * h 182 return join_fields(-pec_e, pec_h)

183

184 - def absorbing_bc(self, w=None):

185 e, h = self.split_eh(self.field_placeholder(w)) 186 187 from hedge.optemplate import make_normal 188 absorb_normal = make_normal(self.absorb_tag, self.dimensions) 189 190 from hedge.optemplate import BoundarizeOperator 191 absorb_e = BoundarizeOperator(self.absorb_tag) * e 192 absorb_h = BoundarizeOperator(self.absorb_tag) * h 193 absorb_w = BoundarizeOperator(self.absorb_tag) * w 194 195 from hedge.tools import join_fields 196 return absorb_w + 1/2*join_fields( 197 self.space_cross_h(absorb_normal, self.space_cross_e( 198 absorb_normal, absorb_e)) 199 - self.Z*self.space_cross_h(absorb_normal, absorb_h), 200 self.space_cross_e(absorb_normal, self.space_cross_h( 201 absorb_normal, absorb_h)) 202 + self.Y*self.space_cross_e(absorb_normal, absorb_e) 203 )

204

205 - def incident_bc(self, w=None):

206 e, h = self.split_eh(self.field_placeholder(w)) 207 208 from hedge.tools import count_subset 209 fld_cnt = count_subset(self.get_eh_subset()) 210 211 if self.incident_bc_data is not None: 212 from hedge.optemplate import make_common_subexpression 213 return make_common_subexpression( 214 -make_vector_field("incident_bc", fld_cnt)) 215 else: 216 from hedge.tools import make_obj_array 217 return make_obj_array([0]*fld_cnt)

218

219 - def op_template(self, w=None):

220 w = self.field_placeholder(w) 221 222 from hedge.optemplate import pair_with_boundary, \ 223 InverseMassOperator, get_flux_operator 224 225 flux_op = get_flux_operator(self.flux(self.flux_type)) 226 bdry_flux_op = get_flux_operator(self.flux(self.bdry_flux_type)) 227 228 return - self.local_derivatives(w) \ 229 + InverseMassOperator()*( 230 flux_op * w 231 +bdry_flux_op * pair_with_boundary( 232 w, self.pec_bc(w), self.pec_tag) 233 +bdry_flux_op * pair_with_boundary( 234 w, self.absorbing_bc(w), self.absorb_tag) 235 +bdry_flux_op * pair_with_boundary( 236 w, self.incident_bc(w), self.incident_tag) 237 )

238

239 - def bind(self, discr, **extra_context):

240 from hedge.mesh import check_bc_coverage 241 check_bc_coverage(discr.mesh, [ 242 self.pec_tag, self.absorb_tag, self.incident_tag]) 243 244 compiled_op_template = discr.compile(self.op_template()) 245 246 from hedge.tools import full_to_subset_indices 247 e_indices = full_to_subset_indices(self.get_eh_subset()[0:3]) 248 all_indices = full_to_subset_indices(self.get_eh_subset()) 249 250 def rhs(t, w): 251 if self.current is not None: 252 j = self.current.volume_interpolant(t, discr)[e_indices] 253 else: 254 j = 0 255 256 if self.incident_bc_data is not None: 257 incident_bc_data = self.incident_bc_data.boundary_interpolant( 258 t, discr, self.incident_tag)[all_indices] 259 else: 260 incident_bc_data = 0 261 262 return compiled_op_template( 263 w=w, j=j, incident_bc=incident_bc_data, **extra_context)

264 265 return rhs 266

267 - def assemble_eh(self, e=None, h=None, discr=None):

268 if discr is None: 269 def zero(): return 0 270 else: 271 def zero(): return discr.volume_zeros() 272 273 from hedge.tools import count_subset 274 e_components = count_subset(self.get_eh_subset()[0:3]) 275 h_components = count_subset(self.get_eh_subset()[3:6]) 276 277 def default_fld(fld, comp): 278 if fld is None: 279 return [zero() for i in xrange(comp)] 280 else: 281 return fld

282 283 e = default_fld(e, e_components) 284 h = default_fld(h, h_components) 285 286 from hedge.tools import join_fields 287 return join_fields(e, h) 288 289 @memoize_method

290 - def partial_to_eh_subsets(self):

291 e_subset = self.get_eh_subset()[0:3] 292 h_subset = self.get_eh_subset()[3:6] 293 294 from hedge.tools import partial_to_all_subset_indices 295 return tuple(partial_to_all_subset_indices( 296 [e_subset, h_subset]))

297

298 - def split_eh(self, w):

299 e_idx, h_idx = self.partial_to_eh_subsets() 300 e, h = w[e_idx], w[h_idx] 301 302 from hedge.flux import FluxVectorPlaceholder as FVP 303 if isinstance(w, FVP): 304 return FVP(scalars=e), FVP(scalars=h) 305 else: 306 from hedge.tools import make_obj_array as moa 307 return moa(e), moa(h)

308

309 - def get_eh_subset(self):

310 """Return a 6-tuple of C{bool}s indicating whether field components 311 are to be computed. The fields are numbered in the order specified 312 in the class documentation. 313 """ 314 return 6*(True,)

315

316 - def max_eigenvalue(self):

317 """Return the largest eigenvalue of Maxwell's equations as a hyperbolic system.""" 318 from math import sqrt 319 return 1/sqrt(self.mu*self.epsilon)

320

321 322 323 324 -class TMMaxwellOperator(MaxwellOperator):

325 """A 2D TM Maxwell operator with PEC boundaries. 326 327 Field order is [Ez Hx Hy]. 328 """ 329 330 _default_dimensions = 2 331

332 - def get_eh_subset(self):

333 return ( 334 (False,False,True) # only ez 335 + 336 (True,True,False) # hx and hy 337 )

338

339 340 341 342 -class TEMaxwellOperator(MaxwellOperator):

343 """A 2D TE Maxwell operator. 344 345 Field order is [Ex Ey Hz]. 346 """ 347 348 _default_dimensions = 2 349

350 - def get_eh_subset(self):

351 return ( 352 (True,True,False) # ex and ey 353 + 354 (False,False,True) # only hz 355 )

356

357 -class TE1DMaxwellOperator(MaxwellOperator):

358 """A 1D TE Maxwell operator. 359 360 Field order is [Ex Ey Hz]. 361 """ 362 363 _default_dimensions = 1 364

365 - def get_eh_subset(self):

366 return ( 367 (True,True,False) 368 + 369 (False,False,True) 370 )

371

372 373 -class SourceFree1DMaxwellOperator(MaxwellOperator):

374 """A 1D TE Maxwell operator. 375 376 Field order is [Ey Hz]. 377 """ 378 379 _default_dimensions = 1 380

381 - def get_eh_subset(self):

382 return ( 383 (False,True,False) 384 + 385 (False,False,True) 386 )

387

Source Code for Module hedge.models.em