diff --git a/grudge/models/em.py b/grudge/models/em.py
index 3f4bc13d42a017fcee36ba7b05237c512ca11451..a3fc5b323cde4f46764758ea3f704eb30b2f908c 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -29,10 +29,15 @@ THE SOFTWARE.
from pytools import memoize_method
from grudge.models import HyperbolicOperator
+from grudge.symbolic.primitives import TracePair
+
+from meshmode.dof_array import thaw
from meshmode.mesh import BTAG_ALL, BTAG_NONE
-from grudge import sym
+
from pytools.obj_array import flat_obj_array, make_obj_array
+import grudge.op as op
+
class MaxwellOperator(HyperbolicOperator):
"""A strong-form 3D Maxwell operator which supports fixed or variable
@@ -43,7 +48,7 @@ class MaxwellOperator(HyperbolicOperator):
_default_dimensions = 3
- def __init__(self, epsilon, mu,
+ def __init__(self, dcoll, epsilon, mu,
flux_type,
bdry_flux_type=None,
pec_tag=BTAG_ALL,
@@ -66,6 +71,7 @@ class MaxwellOperator(HyperbolicOperator):
boundary condition
"""
+ self.dcoll = dcoll
self.dimensions = dimensions or self._default_dimensions
space_subset = [True]*self.dimensions + [False]*(3-self.dimensions)
@@ -103,7 +109,7 @@ class MaxwellOperator(HyperbolicOperator):
self.current = current
self.incident_bc_data = incident_bc
- def flux(self, w):
+ def flux(self, wtpair):
"""The numerical flux for variable coefficients.
:param flux_type: can be in [0,1] for anything between central and upwind,
@@ -112,10 +118,10 @@ class MaxwellOperator(HyperbolicOperator):
As per Hesthaven and Warburton page 433.
"""
- normal = sym.normal(w.dd, self.dimensions)
+ normal = thaw(self.dcoll._setup_actx, op.normal(self.dcoll, wtpair.dd))
if self.fixed_material:
- e, h = self.split_eh(w)
+ e, h = self.split_eh(wtpair)
epsilon = self.epsilon
mu = self.mu
@@ -168,13 +174,18 @@ class MaxwellOperator(HyperbolicOperator):
e, h = self.split_eh(w)
- nabla = sym.nabla(self.dimensions)
+ # Object array of derivative operators
+ nabla = flat_obj_array(
+ [_Dx(self.dcoll, i) for i in range(self.dimensions)]
+ )
def e_curl(field):
- return self.space_cross_e(nabla, field)
+ return self.space_cross_e(nabla, field,
+ three_mult=lambda lc, x, y: lc * (x * y))
def h_curl(field):
- return self.space_cross_h(nabla, field)
+ return self.space_cross_h(nabla, field,
+ three_mult=lambda lc, x, y: lc * (x * y))
# in conservation form: u_t + A u_x = 0
return flat_obj_array(
@@ -187,8 +198,8 @@ class MaxwellOperator(HyperbolicOperator):
"""
e, h = self.split_eh(w)
- pec_e = sym.cse(sym.project("vol", self.pec_tag)(e))
- pec_h = sym.cse(sym.project("vol", self.pec_tag)(h))
+ pec_e = op.project(self.dcoll, "vol", self.pec_tag, e)
+ pec_h = op.project(self.dcoll, "vol", self.pec_tag, h)
return flat_obj_array(-pec_e, pec_h)
@@ -197,8 +208,8 @@ class MaxwellOperator(HyperbolicOperator):
"""
e, h = self.split_eh(w)
- pmc_e = sym.cse(sym.project("vol", self.pmc_tag)(e))
- pmc_h = sym.cse(sym.project("vol", self.pmc_tag)(h))
+ pmc_e = op.project(self.dcoll, "vol", self.pmc_tag, e)
+ pmc_h = op.project(self.dcoll, "vol", self.pmc_tag, h)
return flat_obj_array(pmc_e, -pmc_h)
@@ -207,7 +218,8 @@ class MaxwellOperator(HyperbolicOperator):
absorbing boundary conditions.
"""
- absorb_normal = sym.normal(self.absorb_tag, self.dimensions)
+ absorb_normal = thaw(self.dcoll._setup_actx,
+ op.normal(self.dcoll, dd=self.absorb_tag))
e, h = self.split_eh(w)
@@ -218,8 +230,8 @@ class MaxwellOperator(HyperbolicOperator):
absorb_Z = (mu/epsilon)**0.5 # noqa: N806
absorb_Y = 1/absorb_Z # noqa: N806
- absorb_e = sym.cse(sym.project("vol", self.absorb_tag)(e))
- absorb_h = sym.cse(sym.project("vol", self.absorb_tag)(h))
+ absorb_e = op.project(self.dcoll, "vol", self.absorb_tag, e)
+ absorb_h = op.project(self.dcoll, "vol", self.absorb_tag, h)
bc = flat_obj_array(
absorb_e + 1/2*(self.space_cross_h(absorb_normal, self.space_cross_e(
@@ -247,9 +259,9 @@ class MaxwellOperator(HyperbolicOperator):
if is_zero(incident_bc_data):
return make_obj_array([0]*fld_cnt)
else:
- return sym.cse(-incident_bc_data)
+ return -incident_bc_data
- def sym_operator(self, w=None):
+ def operator(self, t, w):
"""The full operator template - the high level description of
the Maxwell operator.
@@ -257,7 +269,6 @@ class MaxwellOperator(HyperbolicOperator):
derivatives, flux, boundary conditions etc.
"""
from grudge.tools import count_subset
- w = sym.make_sym_array("w", count_subset(self.get_eh_subset()))
elec_components = count_subset(self.get_eh_subset()[0:3])
mag_components = count_subset(self.get_eh_subset()[3:6])
@@ -274,17 +285,27 @@ class MaxwellOperator(HyperbolicOperator):
(self.incident_tag, self.incident_bc(w)),
]
+ dcoll = self.dcoll
+
def flux(pair):
- return sym.project(pair.dd, "all_faces")(self.flux(pair))
+ return op.project(dcoll, pair.dd, "all_faces", self.flux(pair))
+
+ def bv_tpair(tag, w, bc):
+ return TracePair(tag, interior=op.project(dcoll, "vol", tag, w),
+ exterior=bc)
return (
- - self.local_derivatives(w)
- - sym.InverseMassOperator()(sym.FaceMassOperator()(
- flux(sym.int_tpair(w))
- + sum(
- flux(sym.bv_tpair(tag, w, bc))
- for tag, bc in tags_and_bcs)
- ))) / material_divisor
+ - self.local_derivatives(w)
+ - op.inverse_mass(
+ dcoll,
+ op.face_mass(
+ dcoll,
+ flux(op.interior_trace_pair(dcoll, w))
+ + sum(flux(bv_tpair(tag, w, bc))
+ for tag, bc in tags_and_bcs)
+ )
+ )
+ ) / material_divisor
@memoize_method
def partial_to_eh_subsets(self):
@@ -322,10 +343,8 @@ class MaxwellOperator(HyperbolicOperator):
if self.fixed_material:
return 1/sqrt(self.epsilon*self.mu) # a number
else:
- import grudge.symbolic as sym
- return sym.NodalMax("vol")(
- 1 / sym.sqrt(self.epsilon * self.mu)
- )
+ actx = self.dcoll._setup_actx
+ return op.nodal_maximum(1 / actx.np.sqrt(self.epsilon * self.mu))
def max_eigenvalue(self, t, fields=None, discr=None, context=None):
if context is None:
@@ -345,6 +364,17 @@ class MaxwellOperator(HyperbolicOperator):
self.incident_tag])
+# NOTE: Hack for getting the derivative operators to play nice
+# with grudge.tools.SubsettableCrossProduct
+class _Dx:
+ def __init__(self, dcoll, i):
+ self.dcoll = dcoll
+ self.i = i
+
+ def __mul__(self, other):
+ return op.local_d_dx(self.dcoll, self.i, other)
+
+
class TMMaxwellOperator(MaxwellOperator):
"""A 2D TM Maxwell operator with PEC boundaries.
@@ -405,7 +435,7 @@ class SourceFree1DMaxwellOperator(MaxwellOperator):
)
-def get_rectangular_cavity_mode(E_0, mode_indices): # noqa: N803
+def get_rectangular_cavity_mode(actx, nodes, t, E_0, mode_indices): # noqa: N803
"""A rectangular TM cavity mode for a rectangle / cube
with one corner at the origin and the other at (1,1[,1])."""
dims = len(mode_indices)
@@ -421,43 +451,44 @@ def get_rectangular_cavity_mode(E_0, mode_indices): # noqa: N803
omega = numpy.sqrt(sum(f**2 for f in factors))
- nodes = sym.nodes(dims)
x = nodes[0]
y = nodes[1]
if dims == 3:
z = nodes[2]
- sx = sym.sin(kx*x)
- cx = sym.cos(kx*x)
- sy = sym.sin(ky*y)
- cy = sym.cos(ky*y)
+ zeros = 0*x
+ sx = actx.np.sin(kx*x)
+ cx = actx.np.cos(kx*x)
+ sy = actx.np.sin(ky*y)
+ cy = actx.np.cos(ky*y)
if dims == 3:
- sz = sym.sin(kz*z)
- cz = sym.cos(kz*z)
+ sz = actx.np.sin(kz*z)
+ cz = actx.np.cos(kz*z)
if dims == 2:
- tfac = sym.ScalarVariable("t") * omega
+ tfac = t * omega
result = flat_obj_array(
- 0,
- 0,
- sym.sin(kx * x) * sym.sin(ky * y) * sym.cos(tfac), # ez
- -ky * sym.sin(kx * x) * sym.cos(ky * y) * sym.sin(tfac) / omega, # hx
- kx * sym.cos(kx * x) * sym.sin(ky * y) * sym.sin(tfac) / omega, # hy
- 0,
- )
+ zeros,
+ zeros,
+ actx.np.sin(kx * x) * actx.np.sin(ky * y) * actx.np.cos(tfac), # ez
+ (-ky * actx.np.sin(kx * x) * actx.np.cos(ky * y)
+ * actx.np.sin(tfac) / omega), # hx
+ (kx * actx.np.cos(kx * x) * actx.np.sin(ky * y)
+ * actx.np.sin(tfac) / omega), # hy
+ zeros,
+ )
else:
- tdep = sym.exp(-1j * omega * sym.ScalarVariable("t"))
+ tdep = numpy.exp(-1j * omega * t)
gamma_squared = ky**2 + kx**2
result = flat_obj_array(
-kx * kz * E_0*cx*sy*sz*tdep / gamma_squared, # ex
-ky * kz * E_0*sx*cy*sz*tdep / gamma_squared, # ey
E_0 * sx*sy*cz*tdep, # ez
-
-1j * omega * ky*E_0*sx*cy*cz*tdep / gamma_squared, # hx
1j * omega * kx*E_0*cx*sy*cz*tdep / gamma_squared,
- 0,
- )
+ zeros,
+ )
return result