diff --git a/contrib/maxima/maxwellrectcav3d.mac b/contrib/maxima/maxwellrectcav3d.mac
new file mode 100644
index 0000000000000000000000000000000000000000..abda9687d244b25b3cf546bea9f83383e222bc34
--- /dev/null
+++ b/contrib/maxima/maxwellrectcav3d.mac
@@ -0,0 +1,73 @@
+kill(all);
+load("eigen");
+load("itensor");
+load("diag");
+
+/* -------------------------------------------------------------------------- */
+/* Utilities */
+/* -------------------------------------------------------------------------- */
+mu: 1;
+epsilon: 1;
+
+coords:[x,y,z];
+
+curl(x):=crossfunc(lambda([i,j], diff(x[j], coords[i])));
+div(f):=sum(diff(f[j], coords[j]), j, 1, length(coords));
+
+faraday(max_E, max_H):= curl(max_E) - %i * omega * max_H;
+ampere(max_E, max_H):= curl(max_H) + %i * omega * mu * epsilon * max_E;
+/*
+ampere(max_E, max_H):= curl(max_H) + %i * omega * max_E;
+*/
+
+
+crossfunc(f):=makelist(
+ sum(sum(
+ levi_civita([i,j,k])*f(j,k),
+ j,1,3),k,1,3),i,1,3)$
+
+crossprod(a,b):=crossfunc(lambda([j,k], a[j]*b[k]));
+
+/* -------------------------------------------------------------------------- */
+/* Attempts */
+/* -------------------------------------------------------------------------- */
+nabla_t_squared(f):=diff(f, x, 2) + diff(f, y, 2);
+nabla_t(f):= [diff(f, x, 1) , diff(f, y, 1), 0];
+
+/*
+gamma : sqrt((k_y^2 + k_x^2)/ (mu *epsilon)) ;
+*/
+omega : sqrt(k_x^2 + k_y^2 + k_z^2);
+gamma : sqrt((k_y^2 + k_x^2)) ;
+
+psi_cand:E_0*sin(k_x*x)*sin(k_y*y);
+
+
+
+wave_eqn(f, gam_s):=nabla_t_squared(f) + mu*epsilon*gam_s*f;
+gam_s : gamma^2;
+
+/* The _t indicates transverse components (i.e. x and y components only) */
+/*
+E_t(psi):=(-k_z/gam_s)*sin(k_z * z)*Â nabla_t(psi);
+H_t(psi):=crossprod((%i * omega /gam_s)*cos(k_z * z)* [0,0,1], nabla_t(psi));
+*/
+E_t(psi):=(-k_z/gam_s)*sin(k_z * z)*Â nabla_t(psi);
+H_t(psi):=crossprod((%i * omega * epsilon/gam_s)*cos(k_z * z)* [0,0,1], nabla_t(psi));
+
+
+/* These are used as the analytic solution for a rectangular cavity resonator
+ with travsverse magnetic waves */
+
+E : E_t(psi_cand) + [0,0,psi_cand * cos(k_z * z)];
+H :H_t(psi_cand);
+
+Etime : E * %e ^ (- %i * omega * t);
+Htime : H * %e ^ (- %i * omega * t);
+
+
+trigrat(div(E));
+trigrat(div(H));
+trigrat(faraday(E,H));
+trigrat(ampere(E,H));
+
diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index edce0760f260449574bc84d3c5c4c188c009e2f1..ddc7535288cd82c1bbc56cb8decd9b0644cea73a 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -27,7 +27,7 @@ from pyopencl.tools import ( # noqa
import logging
logger = logging.getLogger(__name__)
-from grudge import sym, bind, Discretization
+from grudge import sym, bind, DGDiscretizationWithBoundaries
from pytools.obj_array import join_fields
@@ -44,7 +44,7 @@ def main(write_output=True, order=4):
dt_factor = 4
h = 1/20
- discr = Discretization(cl_ctx, mesh, order=order)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
sym_x = sym.nodes(2)
@@ -68,7 +68,7 @@ def main(write_output=True, order=4):
from grudge.models.advection import VariableCoefficientAdvectionOperator
- discr = Discretization(cl_ctx, mesh, order=order)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
op = VariableCoefficientAdvectionOperator(2, advec_v,
inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
flux_type=flux_type)
diff --git a/examples/advection/weak.py b/examples/advection/weak.py
index 967eb95f594f2b21e216bdf06c17ada799fa9e3c..cdadb49e79e85c08a611f856c54970bf3c417b12 100644
--- a/examples/advection/weak.py
+++ b/examples/advection/weak.py
@@ -27,7 +27,7 @@ from pyopencl.tools import ( # noqa
import logging
logger = logging.getLogger(__name__)
-from grudge import sym, bind, Discretization
+from grudge import sym, bind, DGDiscretizationWithBoundaries
import numpy.linalg as la
@@ -48,7 +48,7 @@ def main(write_output=True, order=4):
dt_factor = 4
h = 1/20
- discr = Discretization(cl_ctx, mesh, order=order)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
c = np.array([0.1,0.1])
norm_c = la.norm(c)
@@ -66,7 +66,7 @@ def main(write_output=True, order=4):
from grudge.models.advection import WeakAdvectionOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
- discr = Discretization(cl_ctx, mesh, order=order)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
op = WeakAdvectionOperator(c,
inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
flux_type=flux_type)
diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py
deleted file mode 100644
index 7ff8884f1d19e46d320e0cc7a8c5301f74d14551..0000000000000000000000000000000000000000
--- a/examples/maxwell/analytic_solutions.py
+++ /dev/null
@@ -1,460 +0,0 @@
-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-
-
-
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-from grudge.tools import \
- cyl_bessel_j, \
- cyl_bessel_j_prime
-from math import sqrt, pi, sin, cos, atan2, acos
-import numpy
-import numpy.linalg as la
-import cmath
-from six.moves import range
-from six.moves import zip
-
-
-
-
-# solution adapters -----------------------------------------------------------
-class RealPartAdapter:
- def __init__(self, adaptee):
- self.adaptee = adaptee
-
- @property
- def shape(self):
- return self.adaptee.shape
-
- def __call__(self, x, el):
- return [xi.real for xi in self.adaptee(x, el)]
-
-class SplitComplexAdapter:
- def __init__(self, adaptee):
- self.adaptee = adaptee
-
- @property
- def shape(self):
- (n,) = self.adaptee.shape
- return (n*2,)
-
- def __call__(self, x, el):
- ad_x = self.adaptee(x, el)
- return [xi.real for xi in ad_x] + [xi.imag for xi in ad_x]
-
-
-
-
-class CylindricalFieldAdapter:
- """Turns a field returned as (r, phi, z) components into
- (x,y,z) components.
- """
- def __init__(self, adaptee):
- self.adaptee = adaptee
-
- @property
- def shape(self):
- return self.adaptee.shape
-
- def __call__(self, x, el):
- xy = x[:2]
- r = la.norm(xy)
- phi = atan2(x[1], x[0])
-
- prev_result = self.adaptee(x, el)
- result = []
- i = 0
- while i < len(prev_result):
- fr, fphi, fz = prev_result[i:i+3]
- result.extend([
- cos(phi)*fr - sin(phi)*fphi, # ex
- sin(phi)*fr + cos(phi)*fphi, # ey
- fz,
- ])
- i += 3
-
- return result
-
-
-
-
-class SphericalFieldAdapter:
- """Turns the adaptee's field C{(Er, Etheta, Ephi)} components into
- C{(Ex,Ey,Ez)} components.
-
- Conventions are those of
- http://mathworld.wolfram.com/SphericalCoordinates.html
- """
- def __init__(self, adaptee):
- self.adaptee = adaptee
-
- @property
- def shape(self):
- return self.adaptee.shape
-
- @staticmethod
- def r_theta_phi_from_xyz(x):
- r = la.norm(x)
- return r, atan2(x[1], x[0]), acos(x[2]/r)
-
- def __call__(self, x, el):
- r, theta, phi = self.r_theta_phi_from_xyz(x)
-
- er = numpy.array([
- cos(theta)*sin(phi),
- sin(theta)*sin(phi),
- cos(phi)])
- etheta = numpy.array([
- -sin(theta),
- cos(theta),
- 0])
- ephi = numpy.array([
- cos(theta)*cos(phi),
- sin(theta)*cos(phi),
- -sin(phi)])
-
- adaptee_result = self.adaptee(x, el)
- result = numpy.empty_like(adaptee_result)
- i = 0
- while i < len(adaptee_result):
- fr, ftheta, fphi = adaptee_result[i:i+3]
- result[i:i+3] = fr*er + ftheta*etheta + fphi*ephi
- i += 3
-
- return result
-
-
-
-
-# actual solutions ------------------------------------------------------------
-class CylindricalCavityMode:
- """A cylindrical TM cavity mode.
-
- Taken from:
- J.D. Jackson, Classical Electrodynamics, Wiley.
- 3rd edition, 2001.
- ch 8.7, p. 368f.
- """
-
- def __init__(self, m, n, p, radius, height, epsilon, mu):
- try:
- from bessel_zeros import bessel_zeros
- except ImportError:
- print("*** You need to generate the bessel root data file.")
- print("*** Execute generate-bessel-zeros.py at the command line.")
- raise
-
- assert m >= 0 and m == int(m)
- assert n >= 1 and n == int(n)
- assert p >= 0 and p == int(p)
- self.m = m
- self.n = n
- self.p = p
- self.phi_sign = 1
-
- R = self.radius = radius
- d = self.height = height
-
- self.epsilon = epsilon
- self.mu = mu
-
- self.t = 0
-
- x_mn = bessel_zeros[m][n-1]
-
- self.omega = 1 / sqrt(mu*epsilon) * sqrt(
- x_mn**2 / R**2
- + p**2 * pi**2 / d**2)
-
- self.gamma_mn = x_mn/R
-
- def set_time(self, t):
- self.t = t
-
- shape = (6,)
-
- def __call__(self, x, el):
- # coordinates -----------------------------------------------------
- xy = x[:2]
- r = sqrt(xy*xy)
- phi = atan2(x[1], x[0])
- z = x[2]
-
- # copy instance variables for easier access -----------------------
- m = self.m
- p = self.p
- gamma = self.gamma_mn
- phi_sign = self.phi_sign
- omega = self.omega
- d = self.height
- epsilon = self.epsilon
-
- # common subexpressions -------------------------------------------
- tdep = cmath.exp(-1j * omega * self.t)
- phi_factor = cmath.exp(phi_sign * 1j * m * phi)
-
- # psi and derivatives ---------------------------------------------
- psi = cyl_bessel_j(m, gamma * r) * phi_factor
- psi_dr = gamma*cyl_bessel_j_prime(m, gamma*r) * phi_factor
- psi_dphi = (cyl_bessel_j(m, gamma * r)
- * 1/r * phi_sign*1j*m * phi_factor)
-
- # field components in polar coordinates ---------------------------
- ez = tdep * cos(p * pi * z / d) * psi
-
- e_transverse_factor = (tdep
- * (-p*pi/(d*gamma**2))
- * sin(p * pi * z / d))
-
- er = e_transverse_factor * psi_dr
- ephi = e_transverse_factor * psi_dphi
-
- hz = 0j
-
- # z x grad psi = z x (psi_x, psi_y) = (-psi_y, psi_x)
- # z x grad psi = z x (psi_r, psi_phi) = (-psi_phi, psi_r)
- h_transverse_factor = (tdep
- * 1j*epsilon*omega/gamma**2
- * cos(p * pi * z / d))
-
- hr = h_transverse_factor * (-psi_dphi)
- hphi = h_transverse_factor * psi_dr
-
- return [er, ephi, ez, hr, hphi, hz]
-
-
-
-
-class RectangularWaveguideMode:
- """A rectangular TM cavity mode."""
-
- def __init__(self, epsilon, mu, mode_indices,
- dimensions=(1,1,1), coefficients=(1,0,0),
- forward_coeff=1, backward_coeff=0):
- for n in mode_indices:
- assert n >= 0 and n == int(n)
- self.mode_indices = mode_indices
- self.dimensions = dimensions
- self.coefficients = coefficients
- self.forward_coeff = forward_coeff
- self.backward_coeff = backward_coeff
-
- self.epsilon = epsilon
- self.mu = mu
-
- self.t = 0
-
- self.factors = [n*pi/a for n, a in zip(self.mode_indices, self.dimensions)]
-
- c = 1/sqrt(mu*epsilon)
- self.k = sqrt(sum(f**2 for f in self.factors))
- self.omega = self.k*c
-
- def set_time(self, t):
- self.t = t
-
- def __call__(self, discr):
- f,g,h = self.factors
- omega = self.omega
- k = self.k
-
- x = discr.nodes[:,0]
- y = discr.nodes[:,1]
- if discr.dimensions > 2:
- z = discr.nodes[:,2]
- else:
- z = discr.volume_zeros()
-
- sx = numpy.sin(f*x)
- sy = numpy.sin(g*y)
- cx = numpy.cos(f*x)
- cy = numpy.cos(g*y)
-
- zdep_add = numpy.exp(1j*h*z)+numpy.exp(-1j*h*z)
- zdep_sub = numpy.exp(1j*h*z)-numpy.exp(-1j*h*z)
-
- tdep = numpy.exp(-1j * omega * self.t)
-
- C = 1j/(f**2+g**2)
-
- result = numpy.zeros((6,len(x)), dtype=zdep_add.dtype)
-
- result[0] = C*f*h*cx*sy*zdep_sub*tdep # ex
- result[1] = C*g*h*sx*cy*zdep_sub*tdep # ey
- result[2] = sx*sy*zdep_add*tdep # ez
- result[3] = -C*g*self.epsilon*omega*sx*cy*zdep_add*tdep # hx
- result[4] = C*f*self.epsilon*omega*cx*sy*zdep_add*tdep
-
- return result
-
-
-
-
-class RectangularCavityMode(RectangularWaveguideMode):
- """A rectangular TM cavity mode."""
-
- def __init__(self, *args, **kwargs):
- if "scale" in kwargs:
- kwargs["forward_coeff"] = scale
- kwargs["backward_coeff"] = scale
- else:
- kwargs["forward_coeff"] = 1
- kwargs["backward_coeff"] = 1
- RectangularWaveguideMode.__init__(self, *args, **kwargs)
-
-
-
-
-
-# dipole solution -------------------------------------------------------------
-class DipoleFarField:
- """Dipole EM field.
-
- See U{http://piki.tiker.net/wiki/Dipole}.
- """
- def __init__(self, q, d, omega, epsilon, mu):
- self.q = q
- self.d = d
-
- self.epsilon = epsilon
- self.mu = mu
- self.c = c = 1/sqrt(mu*epsilon)
-
- self.omega = omega
-
- freq = omega/(2*pi)
- self.wavelength = c/freq
-
- def far_field_mask(self, x, el):
- if la.norm(x) > 0.5*self.wavelength:
- return 1
- else:
- return 0
-
- def set_time(self, t):
- self.t = t
-
- shape = (6,)
-
- def __call__(self, x, el):
- # coordinates ---------------------------------------------------------
- r, theta, phi = SphericalFieldAdapter.r_theta_phi_from_xyz(x)
-
- # copy instance variables for easier access ---------------------------
- q = self.q
- d = self.d
-
- epsilon = self.epsilon
- mu = self.mu
- c = self.c
-
- omega = self.omega
- t = self.t
-
- # field components ----------------------------------------------------
- tdep = omega*t-omega*r/c
- e_r = (2*cos(phi)*q*d) / (4*pi*epsilon) * (
- 1/r**3 * sin(tdep)
- +omega/(c*r**2) * cos(tdep))
- e_phi = (sin(phi)*q*d) / (4*pi*epsilon) * (
- (1/r**3 - omega**2/(c**2*r)*sin(tdep))
- + omega/(c*r**2)*cos(tdep))
- h_theta = (omega*sin(phi)*q*d)/(4*pi) * (
- - omega/(c*r)*sin(tdep)
- + 1/r**2*cos(tdep))
-
- return [e_r, 0, e_phi, 0, h_theta, 0]
-
- def source_modulation(self, t):
- j0 = self.q*self.d*self.omega/self.epsilon
- return j0*cos(self.omega*t)
-
-
-
-
-# analytic solution tools -----------------------------------------------------
-def check_time_harmonic_solution(discr, mode, c_sol):
- from grudge.discretization import bind_nabla, bind_mass_matrix
- from grudge.visualization import SiloVisualizer
- from grudge.silo import SiloFile
- from grudge.tools import dot, cross
- from grudge.silo import DB_VARTYPE_VECTOR
-
- def curl(field):
- return cross(nabla, field)
-
- vis = SiloVisualizer(discr)
-
- nabla = bind_nabla(discr)
- mass = bind_mass_matrix(discr)
-
- def rel_l2_error(err, base):
- def l2_norm(field):
- return sqrt(dot(field, mass*field))
-
- base_l2 = l2_norm(base)
- err_l2 = l2_norm(err)
- if base_l2 == 0:
- if err_l2 == 0:
- return 0.
- else:
- return float("inf")
- else:
- return err_l2/base_l2
-
- dt = 0.1
-
- for step in range(10):
- t = step*dt
- mode.set_time(t)
- fields = discr.interpolate_volume_function(c_sol)
-
- er = fields[0:3]
- hr = fields[3:6]
- ei = fields[6:9]
- hi = fields[9:12]
-
- silo = SiloFile("em-complex-%04d.silo" % step)
- vis.add_to_silo(silo,
- vectors=[
- ("curl_er", curl(er)),
- ("om_hi", -mode.mu*mode.omega*hi),
- ("curl_hr", curl(hr)),
- ("om_ei", mode.epsilon*mode.omega*hi),
- ],
- expressions=[
- ("diff_er", "curl_er-om_hi", DB_VARTYPE_VECTOR),
- ("diff_hr", "curl_hr-om_ei", DB_VARTYPE_VECTOR),
- ],
- write_coarse_mesh=True,
- time=t, step=step
- )
-
- er_res = curl(er) + mode.mu *mode.omega*hi
- ei_res = curl(ei) - mode.mu *mode.omega*hr
- hr_res = curl(hr) - mode.epsilon*mode.omega*ei
- hi_res = curl(hi) + mode.epsilon*mode.omega*er
-
- print("time=%f, rel l2 residual in Re[E]=%g\tIm[E]=%g\tRe[H]=%g\tIm[H]=%g" % (
- t,
- rel_l2_error(er_res, er),
- rel_l2_error(ei_res, ei),
- rel_l2_error(hr_res, hr),
- rel_l2_error(hi_res, hi),
- ))
-
diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py
index 6fb7d54ac4eb69cf87705fce00268f68debfd4ba..85a7c3dec0cef99b4a72fe9c23439c4212c0c32d 100644
--- a/examples/maxwell/cavities.py
+++ b/examples/maxwell/cavities.py
@@ -1,260 +1,131 @@
-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy as np
-
-import logging
-logger = logging.getLogger(__name__)
-
-
-def main(write_output=True, allow_features=None, flux_type_arg=1,
- bdry_flux_type_arg=None, extra_discr_args={}):
- from grudge.mesh.generator import make_cylinder_mesh, make_box_mesh
- from grudge.tools import EOCRecorder, to_obj_array
- from math import sqrt, pi # noqa
- from analytic_solutions import ( # noqa
- RealPartAdapter,
- SplitComplexAdapter,
- CylindricalFieldAdapter,
- CylindricalCavityMode,
- RectangularWaveguideMode,
- RectangularCavityMode)
- from grudge.models.em import MaxwellOperator
-
- logging.basicConfig(level=logging.DEBUG)
-
- from grudge.backends import guess_run_context
- rcon = guess_run_context(allow_features)
-
- epsilon0 = 8.8541878176e-12 # C**2 / (N m**2)
- mu0 = 4*pi*1e-7 # N/A**2.
- epsilon = 1*epsilon0
- mu = 1*mu0
-
- eoc_rec = EOCRecorder()
-
- cylindrical = False
- periodic = False
-
- if cylindrical:
- R = 1
- d = 2
- mode = CylindricalCavityMode(m=1, n=1, p=1,
- radius=R, height=d,
- epsilon=epsilon, mu=mu)
- # r_sol = CylindricalFieldAdapter(RealPartAdapter(mode))
- # c_sol = SplitComplexAdapter(CylindricalFieldAdapter(mode))
-
- if rcon.is_head_rank:
- mesh = make_cylinder_mesh(radius=R, height=d, max_volume=0.01)
- else:
- if periodic:
- mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1))
- periodicity = (False, False, True)
- else:
- periodicity = None
- mode = RectangularCavityMode(epsilon, mu, (1, 2, 2))
-
- if rcon.is_head_rank:
- mesh = make_box_mesh(max_volume=0.001, periodicity=periodicity)
-
- if rcon.is_head_rank:
- mesh_data = rcon.distribute_mesh(mesh)
- else:
- mesh_data = rcon.receive_mesh()
-
- for order in [4, 5, 6]:
- #for order in [1,2,3,4,5,6]:
- extra_discr_args.setdefault("debug", []).extend([
- "cuda_no_plan", "cuda_dump_kernels"])
-
- op = MaxwellOperator(epsilon, mu,
- flux_type=flux_type_arg,
- bdry_flux_type=bdry_flux_type_arg)
-
- discr = rcon.make_discretization(mesh_data, order=order,
- tune_for=op.sym_operator(),
- **extra_discr_args)
-
- from grudge.visualization import VtkVisualizer
- if write_output:
- vis = VtkVisualizer(discr, rcon, "em-%d" % order)
-
- mode.set_time(0)
-
- def get_true_field():
- return discr.convert_volume(
- to_obj_array(mode(discr)
- .real.astype(discr.default_scalar_type).copy()),
- kind=discr.compute_kind)
-
- fields = get_true_field()
+"""Minimal example of a grudge driver."""
- if rcon.is_head_rank:
- print("---------------------------------------------")
- print("order %d" % order)
- print("---------------------------------------------")
- print("#elements=", len(mesh.elements))
+from __future__ import division, print_function
- from grudge.timestep.runge_kutta import LSRK4TimeStepper
- stepper = LSRK4TimeStepper(dtype=discr.default_scalar_type, rcon=rcon)
- #from grudge.timestep.dumka3 import Dumka3TimeStepper
- #stepper = Dumka3TimeStepper(3, dtype=discr.default_scalar_type, rcon=rcon)
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
- # {{{ diagnostics setup
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
- from pytools.log import LogManager, add_general_quantities, \
- add_simulation_quantities, add_run_info
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
- if write_output:
- log_file_name = "maxwell-%d.dat" % order
- else:
- log_file_name = None
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
- logmgr = LogManager(log_file_name, "w", rcon.communicator)
- add_run_info(logmgr)
- add_general_quantities(logmgr)
- add_simulation_quantities(logmgr)
- discr.add_instrumentation(logmgr)
- stepper.add_instrumentation(logmgr)
-
- from pytools.log import IntervalTimer
- vis_timer = IntervalTimer("t_vis", "Time spent visualizing")
- logmgr.add_quantity(vis_timer)
-
- from grudge.log import EMFieldGetter, add_em_quantities
- field_getter = EMFieldGetter(discr, op, lambda: fields)
- add_em_quantities(logmgr, op, field_getter)
-
- logmgr.add_watches(
- ["step.max", "t_sim.max",
- ("W_field", "W_el+W_mag"),
- "t_step.max"]
- )
-
- # }}}
-
- # {{{ timestep loop
-
- rhs = op.bind(discr)
- final_time = 0.5e-9
-
- try:
- from grudge.timestep import times_and_steps
- step_it = times_and_steps(
- final_time=final_time, logmgr=logmgr,
- max_dt_getter=lambda t: op.estimate_timestep(discr,
- stepper=stepper, t=t, fields=fields))
+import numpy as np
+import pyopencl as cl
+from grudge.shortcuts import set_up_rk4
+from grudge import sym, bind, Discretization
- for step, t, dt in step_it:
- if step % 50 == 0 and write_output:
- sub_timer = vis_timer.start_sub_timer()
- e, h = op.split_eh(fields)
- visf = vis.make_file("em-%d-%04d" % (order, step))
- vis.add_data(
- visf,
- [
- ("e",
- discr.convert_volume(e, kind="numpy")),
- ("h",
- discr.convert_volume(h, kind="numpy")),
- ],
- time=t, step=step)
- visf.close()
- sub_timer.stop().submit()
+from grudge.models.em import get_rectangular_cavity_mode
- fields = stepper(fields, t, dt, rhs)
- mode.set_time(final_time)
+def main(dims, write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
- eoc_rec.add_data_point(order,
- discr.norm(fields-get_true_field()))
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(0.0,)*dims,
+ b=(1.0,)*dims,
+ n=(5,)*dims)
- finally:
- if write_output:
- vis.close()
+ discr = Discretization(cl_ctx, mesh, order=order)
- logmgr.close()
- discr.close()
+ if 0:
+ epsilon0 = 8.8541878176e-12 # C**2 / (N m**2)
+ mu0 = 4*np.pi*1e-7 # N/A**2.
+ epsilon = 1*epsilon0
+ mu = 1*mu0
+ else:
+ epsilon = 1
+ mu = 1
- if rcon.is_head_rank:
- print()
- print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
+ from grudge.models.em import MaxwellOperator
+ op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
- # }}}
+ queue = cl.CommandQueue(discr.cl_context)
- assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6
+ if dims == 3:
+ sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2))
+ fields = bind(discr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu)
+ else:
+ sym_mode = get_rectangular_cavity_mode(1, (2, 3))
+ fields = bind(discr, sym_mode)(queue, t=0)
+ # FIXME
+ #dt = op.estimate_rk4_timestep(discr, fields=fields)
-# {{{ entry points for py.test
+ op.check_bc_coverage(mesh)
-from pytools.test import mark_test
+ # print(sym.pretty(op.sym_operator()))
+ bound_op = bind(discr, op.sym_operator())
+ def rhs(t, w):
+ return bound_op(queue, t=t, w=w)
-@mark_test.long
-def test_maxwell_cavities():
- main(write_output=False)
+ if mesh.dim == 2:
+ dt = 0.004
+ elif mesh.dim == 3:
+ dt = 0.002
+ dt_stepper = set_up_rk4("w", dt, fields, rhs)
-@mark_test.long
-def test_maxwell_cavities_lf():
- main(write_output=False, flux_type_arg="lf", bdry_flux_type_arg=1)
+ final_t = dt * 600
+ nsteps = int(final_t/dt)
+ print("dt=%g nsteps=%d" % (dt, nsteps))
-@mark_test.mpi
-@mark_test.long
-def test_maxwell_cavities_mpi():
- from pytools.mpi import run_with_mpi_ranks
- run_with_mpi_ranks(__file__, 2, main,
- write_output=False, allow_features=["mpi"])
+ from grudge.shortcuts import make_visualizer
+ vis = make_visualizer(discr, vis_order=order)
+ step = 0
-def test_cuda():
- try:
- from pycuda.tools import mark_cuda_test
- except ImportError:
- pass
+ norm = bind(discr, sym.norm(2, sym.var("u")))
- yield "SP CUDA Maxwell", mark_cuda_test(
- do_test_maxwell_cavities_cuda), np.float32
- yield "DP CUDA Maxwell", mark_cuda_test(
- do_test_maxwell_cavities_cuda), np.float64
+ from time import time
+ t_last_step = time()
+ e, h = op.split_eh(fields)
-def do_test_maxwell_cavities_cuda(dtype):
- import pytest # noqa
+ if 1:
+ vis.write_vtk_file("fld-%04d.vtu" % step,
+ [
+ ("e", e),
+ ("h", h),
+ ])
- main(write_output=False, allow_features=["cuda"],
- extra_discr_args=dict(
- debug=["cuda_no_plan"],
- default_scalar_type=dtype,
- ))
+ for event in dt_stepper.run(t_end=final_t):
+ if isinstance(event, dt_stepper.StateComputed):
+ assert event.component_id == "w"
-# }}}
+ step += 1
+ print(step, event.t, norm(queue, u=e[0]), norm(queue, u=e[1]),
+ norm(queue, u=h[0]), norm(queue, u=h[1]),
+ time()-t_last_step)
+ if step % 10 == 0:
+ e, h = op.split_eh(event.state_component)
+ vis.write_vtk_file("fld-%04d.vtu" % step,
+ [
+ ("e", e),
+ ("h", h),
+ ])
+ t_last_step = time()
-# entry point -----------------------------------------------------------------
if __name__ == "__main__":
- from pytools.mpi import in_mpi_relaunch
- if in_mpi_relaunch():
- test_maxwell_cavities_mpi()
- else:
- do_test_maxwell_cavities_cuda(np.float32)
+ main(3)
diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py
index 6e1a7ea7368198871e2852e1c7ff43c51ddb75c0..f161d0346c1349e508b5ed5a32737e1cccf61312 100644
--- a/examples/wave/var-propagation-speed.py
+++ b/examples/wave/var-propagation-speed.py
@@ -28,7 +28,7 @@ THE SOFTWARE.
import numpy as np
import pyopencl as cl
from grudge.shortcuts import set_up_rk4
-from grudge import sym, bind, Discretization
+from grudge import sym, bind, DGDiscretizationWithBoundaries
def main(write_output=True, order=4):
@@ -42,7 +42,7 @@ def main(write_output=True, order=4):
b=(0.5,)*dims,
n=(20,)*dims)
- discr = Discretization(cl_ctx, mesh, order=order)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
source_width = 0.05
diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py
index aa119aa506014ff3a1910fb120d64b8e493e4213..d7ebffd33a7732c621d159622804fed064f365b3 100644
--- a/examples/wave/wave-min.py
+++ b/examples/wave/wave-min.py
@@ -28,7 +28,7 @@ THE SOFTWARE.
import numpy as np
import pyopencl as cl
from grudge.shortcuts import set_up_rk4
-from grudge import sym, bind, Discretization
+from grudge import sym, bind, DGDiscretizationWithBoundaries
def main(write_output=True, order=4):
@@ -49,7 +49,7 @@ def main(write_output=True, order=4):
print("%d elements" % mesh.nelements)
- discr = Discretization(cl_ctx, mesh, order=order)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
source_width = 0.05
diff --git a/grudge/__init__.py b/grudge/__init__.py
index c6a5d6b777c5fad20e14cbad4aa2b21a32b8f81b..b854007a48c227d58541a13a972a3a7e25c4e062 100644
--- a/grudge/__init__.py
+++ b/grudge/__init__.py
@@ -24,4 +24,4 @@ THE SOFTWARE.
import grudge.sym # noqa
from grudge.execution import bind # noqa
-from grudge.discretization import Discretization # noqa
+from grudge.discretization import DGDiscretizationWithBoundaries # noqa
diff --git a/grudge/discretization.py b/grudge/discretization.py
index b3596acbe2448cf0367df9ddcf9aee18a41678c2..6f39762c00be545aa387253b9f39ea02f5dbeb8a 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -25,24 +25,18 @@ THE SOFTWARE.
from pytools import memoize_method
from grudge import sym
+import numpy as np
+# FIXME Naming not ideal
class DiscretizationBase(object):
pass
-class Discretization(DiscretizationBase):
+class DGDiscretizationWithBoundaries(DiscretizationBase):
"""
- .. attribute :: volume_discr
-
- .. automethod :: boundary_connection
- .. automethod :: boundary_discr
-
- .. automethod :: interior_faces_connection
- .. automethod :: interior_faces_discr
-
- .. automethod :: all_faces_connection
- .. automethod :: all_faces_discr
+ .. automethod :: discr_from_dd
+ .. automethod :: connection_from_dds
.. autoattribute :: cl_context
.. autoattribute :: dim
@@ -67,8 +61,8 @@ class Discretization(DiscretizationBase):
from meshmode.discretization import Discretization
- self.volume_discr = Discretization(cl_ctx, mesh,
- self.get_group_factory_for_quadrature_tag(sym.QTAG_NONE))
+ self._volume_discr = Discretization(cl_ctx, mesh,
+ self.group_factory_for_quadrature_tag(sym.QTAG_NONE))
# {{{ management of discretization-scoped common subexpressions
@@ -80,40 +74,169 @@ class Discretization(DiscretizationBase):
# }}}
- def get_group_factory_for_quadrature_tag(self, quadrature_tag):
- if quadrature_tag is not sym.QTAG_NONE:
- # FIXME
- raise NotImplementedError("quadrature")
+ @memoize_method
+ def discr_from_dd(self, dd):
+ dd = sym.as_dofdesc(dd)
+
+ qtag = dd.quadrature_tag
+
+ if dd.is_volume():
+ if qtag is not sym.QTAG_NONE:
+ # FIXME
+ raise NotImplementedError("quadrature")
+ return self._volume_discr
+
+ elif dd.domain_tag is sym.FACE_RESTR_ALL:
+ return self._all_faces_discr(qtag)
+ elif dd.domain_tag is sym.FACE_RESTR_INTERIOR:
+ return self._interior_faces_discr(qtag)
+ elif dd.is_boundary():
+ return self._boundary_discr(dd.domain_tag, qtag)
+ else:
+ raise ValueError("DOF desc tag not understood: " + str(dd))
+
+ @memoize_method
+ def connection_from_dds(self, from_dd, to_dd):
+ from_dd = sym.as_dofdesc(from_dd)
+ to_dd = sym.as_dofdesc(to_dd)
+
+ if from_dd.quadrature_tag is not sym.QTAG_NONE:
+ raise ValueError("cannot interpolate *from* a "
+ "(non-interpolatory) quadrature grid")
+
+ to_qtag = to_dd.quadrature_tag
+
+ # {{{ simplify domain + qtag change into chained
+
+ if (from_dd.domain_tag != to_dd.domain_tag
+ and from_dd.quadrature_tag != to_dd.quadrature_tag):
+
+ from meshmode.connection import ChainedDiscretizationConnection
+ intermediate_dd = sym.DOFDesc(to_dd.domain_tag)
+ return ChainedDiscretizationConnection(
+ [
+ # first change domain
+ self.connection_from_dds(
+ from_dd,
+ intermediate_dd),
+
+ # then go to quad grid
+ self.connection_from_dds(
+ intermediate_dd,
+ to_dd
+ )])
+
+ # }}}
+
+ # {{{ generic to-quad
+
+ if (from_dd.domain_tag == to_dd.domain_tag
+ and from_dd.quadrature_tag != to_dd.quadrature_tag):
+ from meshmode.discretization.connection.same_mesh import \
+ make_same_mesh_connection
+
+ return make_same_mesh_connection(
+ self.discr_from_dd(to_dd),
+ self.discr_from_dd(from_dd))
+
+ # }}}
+
+ if from_dd.is_volume():
+ if to_dd.domain_tag is sym.FACE_RESTR_ALL:
+ return self._all_faces_volume_connection(to_qtag)
+ if to_dd.domain_tag is sym.FACE_RESTR_INTERIOR:
+ return self._interior_faces_connection(to_qtag)
+ elif to_dd.is_boundary():
+ assert from_dd.quadrature_tag is sym.QTAG_NONE
+ return self._boundary_connection(to_dd.domain_tag)
+
+ else:
+ raise ValueError("cannot interpolate from volume to: " + str(to_dd))
+
+ elif from_dd.domain_tag is sym.FACE_RESTR_INTERIOR:
+ if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE:
+ return self._all_faces_connection(None)
+ else:
+ raise ValueError(
+ "cannot interpolate from interior faces to: "
+ + str(to_dd))
+
+ elif from_dd.domain_tag is sym.FACE_RESTR_ALL:
+ if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE:
+ return self._all_faces_connection(None)
+
+ elif from_dd.is_boundary():
+ if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE:
+ return self._all_faces_connection(from_dd.domain_tag)
+ else:
+ raise ValueError(
+ "cannot interpolate from interior faces to: "
+ + str(to_dd))
+
+ else:
+ raise ValueError("cannot interpolate from: " + str(from_dd))
+
+ def group_factory_for_quadrature_tag(self, quadrature_tag):
+ """
+ OK to override in user code to control mode/node choice.
+ """
+
+ if quadrature_tag is None:
+ quadrature_tag = sym.QTAG_NONE
from meshmode.discretization.poly_element import \
- PolynomialWarpAndBlendGroupFactory
+ PolynomialWarpAndBlendGroupFactory, \
+ QuadratureSimplexGroupFactory
+
+ if quadrature_tag is not sym.QTAG_NONE:
+ return QuadratureSimplexGroupFactory(order=self.order)
+ else:
+ return PolynomialWarpAndBlendGroupFactory(order=self.order)
+
+ @memoize_method
+ def _quad_volume_discr(self, quadrature_tag):
+ from meshmode.discretization import Discretization
- return PolynomialWarpAndBlendGroupFactory(order=self.order)
+ return Discretization(self._volume_discr.cl_context, self._volume_discr.mesh,
+ self.group_factory_for_quadrature_tag(quadrature_tag))
# {{{ boundary
@memoize_method
- def boundary_connection(self, boundary_tag, quadrature_tag=None):
+ def _boundary_connection(self, boundary_tag):
from meshmode.discretization.connection import make_face_restriction
return make_face_restriction(
- self.volume_discr,
- self.get_group_factory_for_quadrature_tag(quadrature_tag),
+ self._volume_discr,
+ self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
boundary_tag=boundary_tag)
- def boundary_discr(self, boundary_tag, quadrature_tag=None):
- return self.boundary_connection(boundary_tag, quadrature_tag).to_discr
+ @memoize_method
+ def _boundary_discr(self, boundary_tag, quadrature_tag=None):
+ if quadrature_tag is None:
+ quadrature_tag = sym.QTAG_NONE
+
+ if quadrature_tag is sym.QTAG_NONE:
+ return self._boundary_connection(boundary_tag).to_discr
+ else:
+ no_quad_bdry_discr = self.boundary_discr(boundary_tag, sym.QTAG_NONE)
+
+ from meshmode.discretization import Discretization
+ return Discretization(
+ self._volume_discr.cl_context,
+ no_quad_bdry_discr.mesh,
+ self.group_factory_for_quadrature_tag(quadrature_tag))
# }}}
# {{{ interior faces
@memoize_method
- def interior_faces_connection(self, quadrature_tag=None):
+ def _interior_faces_connection(self, quadrature_tag=None):
from meshmode.discretization.connection import (
make_face_restriction, FACE_RESTR_INTERIOR)
return make_face_restriction(
- self.volume_discr,
- self.get_group_factory_for_quadrature_tag(quadrature_tag),
+ self._volume_discr,
+ self.group_factory_for_quadrature_tag(quadrature_tag),
FACE_RESTR_INTERIOR,
# FIXME: This will need to change as soon as we support
@@ -121,8 +244,8 @@ class Discretization(DiscretizationBase):
# types.
per_face_groups=False)
- def interior_faces_discr(self, quadrature_tag=None):
- return self.interior_faces_connection(quadrature_tag).to_discr
+ def _interior_faces_discr(self, quadrature_tag=None):
+ return self._interior_faces_connection(quadrature_tag).to_discr
@memoize_method
def opposite_face_connection(self, quadrature_tag):
@@ -134,19 +257,19 @@ class Discretization(DiscretizationBase):
make_opposite_face_connection
return make_opposite_face_connection(
- self.interior_faces_connection(quadrature_tag))
+ self._interior_faces_connection(quadrature_tag))
# }}}
# {{{ all-faces
@memoize_method
- def all_faces_volume_connection(self, quadrature_tag=None):
+ def _all_faces_volume_connection(self, quadrature_tag=None):
from meshmode.discretization.connection import (
make_face_restriction, FACE_RESTR_ALL)
return make_face_restriction(
- self.volume_discr,
- self.get_group_factory_for_quadrature_tag(quadrature_tag),
+ self._volume_discr,
+ self.group_factory_for_quadrature_tag(quadrature_tag),
FACE_RESTR_ALL,
# FIXME: This will need to change as soon as we support
@@ -154,17 +277,17 @@ class Discretization(DiscretizationBase):
# types.
per_face_groups=False)
- def all_faces_discr(self, quadrature_tag=None):
- return self.all_faces_volume_connection(quadrature_tag).to_discr
+ def _all_faces_discr(self, quadrature_tag=None):
+ return self._all_faces_volume_connection(quadrature_tag).to_discr
@memoize_method
- def all_faces_connection(self, boundary_tag, quadrature_tag=None):
+ def _all_faces_connection(self, boundary_tag):
"""Return a
:class:`meshmode.discretization.connection.DiscretizationConnection`
that goes from either
- :meth:`interior_faces_discr` (if *boundary_tag* is None)
+ :meth:`_interior_faces_discr` (if *boundary_tag* is None)
or
- :meth:`boundary_discr` (if *boundary_tag* is not None)
+ :meth:`_boundary_discr` (if *boundary_tag* is not None)
to a discretization containing all the faces of the volume
discretization.
"""
@@ -172,37 +295,44 @@ class Discretization(DiscretizationBase):
make_face_to_all_faces_embedding
if boundary_tag is None:
- faces_conn = self.interior_faces_connection(quadrature_tag)
+ faces_conn = self._interior_faces_connection()
else:
- faces_conn = self.boundary_connection(boundary_tag, quadrature_tag)
+ faces_conn = self._boundary_connection(boundary_tag)
- return make_face_to_all_faces_embedding(
- faces_conn, self.all_faces_discr(quadrature_tag))
+ return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr())
# }}}
@property
def cl_context(self):
- return self.volume_discr.cl_context
+ return self._volume_discr.cl_context
@property
def dim(self):
- return self.volume_discr.dim
+ return self._volume_discr.dim
@property
def ambient_dim(self):
- return self.volume_discr.ambient_dim
+ return self._volume_discr.ambient_dim
+
+ @property
+ def real_dtype(self):
+ return self._volume_discr.real_dtype
+
+ @property
+ def complex_dtype(self):
+ return self._volume_discr.complex_dtype
@property
def mesh(self):
- return self.volume_discr.mesh
+ return self._volume_discr.mesh
def empty(self, queue=None, dtype=None, extra_dims=None, allocator=None):
- return self.volume_discr.empty(queue, dtype, extra_dims=extra_dims,
+ return self._volume_discr.empty(queue, dtype, extra_dims=extra_dims,
allocator=allocator)
def zeros(self, queue, dtype=None, extra_dims=None, allocator=None):
- return self.volume_discr.zeros(queue, dtype, extra_dims=extra_dims,
+ return self._volume_discr.zeros(queue, dtype, extra_dims=extra_dims,
allocator=allocator)
def is_volume_where(self, where):
@@ -219,6 +349,11 @@ class PointsDiscretization(DiscretizationBase):
def __init__(self, nodes):
self._nodes = nodes
+ self.real_dtype = np.dtype(np.float64)
+ self.complex_dtype = np.dtype({
+ np.float32: np.complex64,
+ np.float64: np.complex128
+ }[self.real_dtype.type])
def ambient_dim(self):
return self._nodes.shape[0]
@@ -234,8 +369,16 @@ class PointsDiscretization(DiscretizationBase):
def nodes(self):
return self._nodes
- @property
- def volume_discr(self):
+ def discr_from_dd(self, dd):
+ dd = sym.as_dofdesc(dd)
+
+ if dd.quadrature_tag is not sym.QTAG_NONE:
+ raise ValueError("quadrature discretization requested from "
+ "PointsDiscretization")
+ if dd.domain_tag is not sym.DTAG_VOLUME_ALL:
+ raise ValueError("non-volume discretization requested from "
+ "PointsDiscretization")
+
return self
diff --git a/grudge/execution.py b/grudge/execution.py
index 5199f8fbbf67550e0b54b0a70ad0a9fc6a9453d3..1f7e60b8804a35e84ab8d3d0c2d8f37b33502130 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -43,45 +43,24 @@ class ExecutionMapper(mappers.Evaluator,
mappers.LocalOpReducerMixin):
def __init__(self, queue, context, bound_op):
super(ExecutionMapper, self).__init__(context)
- self.discr = bound_op.discr
+ self.discrwb = bound_op.discrwb
self.bound_op = bound_op
self.queue = queue
- def get_discr(self, dd):
- qtag = dd.quadrature_tag
- if qtag is None:
- # FIXME: Remove once proper quadrature support arrives
- qtag = sym.QTAG_NONE
-
- if dd.is_volume():
- if qtag is not sym.QTAG_NONE:
- # FIXME
- raise NotImplementedError("quadrature")
- return self.discr.volume_discr
-
- elif dd.domain_tag is sym.FACE_RESTR_ALL:
- return self.discr.all_faces_discr(qtag)
- elif dd.domain_tag is sym.FACE_RESTR_INTERIOR:
- return self.discr.interior_faces_discr(qtag)
- elif dd.is_boundary():
- return self.discr.boundary_discr(dd.domain_tag, qtag)
- else:
- raise ValueError("DOF desc tag not understood: " + str(dd))
-
# {{{ expression mappings -------------------------------------------------
def map_ones(self, expr):
if expr.dd.is_scalar():
return 1
- discr = self.get_discr(expr.dd)
+ discr = self.discrwb.discr_from_dd(expr.dd)
result = discr.empty(self.queue, allocator=self.bound_op.allocator)
result.fill(1)
return result
def map_node_coordinate_component(self, expr):
- discr = self.get_discr(expr.dd)
+ discr = self.discrwb.discr_from_dd(expr.dd)
return discr.nodes()[expr.axis].with_queue(self.queue)
def map_grudge_variable(self, expr):
@@ -89,7 +68,7 @@ class ExecutionMapper(mappers.Evaluator,
value = self.context[expr.name]
if not expr.dd.is_scalar() and isinstance(value, Number):
- discr = self.get_discr(expr.dd)
+ discr = self.discrwb.discr_from_dd(expr.dd)
ary = discr.empty(self.queue)
ary.fill(value)
value = ary
@@ -104,7 +83,7 @@ class ExecutionMapper(mappers.Evaluator,
from numbers import Number
if not dd.is_scalar() and isinstance(value, Number):
- discr = self.get_discr(dd)
+ discr = self.discrwb.discr_from_dd(dd)
ary = discr.empty(self.queue)
ary.fill(value)
value = ary
@@ -116,14 +95,40 @@ class ExecutionMapper(mappers.Evaluator,
# FIXME: Make a way to register functions
- pars = [self.rec(p) for p in expr.parameters]
- if any(isinstance(par, cl.array.Array) for par in pars):
- import pyopencl.clmath as clmath
- func = getattr(clmath, expr.function.name)
- else:
+ args = [self.rec(p) for p in expr.parameters]
+ from numbers import Number
+ representative_arg = args[0]
+ if (
+ isinstance(representative_arg, Number)
+ or (isinstance(representative_arg, np.ndarray)
+ and representative_arg.shape == ())):
func = getattr(np, expr.function.name)
+ return func(*args)
- return func(*pars)
+ cached_name = "map_call_knl_"
+
+ i = Variable("i")
+ func = Variable(expr.function.name)
+ if expr.function.name == "fabs": # FIXME
+ func = Variable("abs")
+ cached_name += "abs"
+ else:
+ cached_name += expr.function.name
+
+ @memoize_in(self.bound_op, cached_name)
+ def knl():
+ knl = lp.make_kernel(
+ "{[i]: 0<=i<n}",
+ [
+ lp.Assignment(Variable("out")[i],
+ func(Variable("a")[i]))
+ ], default_offset=lp.auto)
+ return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
+
+ assert len(args) == 1
+ evt, (out,) = knl()(self.queue, a=args[0])
+
+ return out
def map_nodal_sum(self, op, field_expr):
# FIXME: Could allow array scalars
@@ -200,7 +205,7 @@ class ExecutionMapper(mappers.Evaluator,
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
return lp.tag_inames(knl, dict(k="g.0"))
- discr = self.get_discr(op.dd_in)
+ discr = self.discrwb.discr_from_dd(op.dd_in)
# FIXME: This shouldn't really assume that it's dealing with a volume
# input. What about quadrature? What about boundaries?
@@ -237,46 +242,7 @@ class ExecutionMapper(mappers.Evaluator,
return out
def map_interpolation(self, op, field_expr):
- if op.dd_in.quadrature_tag not in [None, sym.QTAG_NONE]:
- raise ValueError("cannot interpolate *from* a quadrature grid")
-
- dd_in = op.dd_in
- dd_out = op.dd_out
-
- qtag = dd_out.quadrature_tag
- if qtag is None:
- # FIXME: Remove once proper quadrature support arrives
- qtag = sym.QTAG_NONE
-
- if dd_in.is_volume():
- if dd_out.domain_tag is sym.FACE_RESTR_ALL:
- conn = self.discr.all_faces_volume_connection(qtag)
- elif dd_out.domain_tag is sym.FACE_RESTR_INTERIOR:
- conn = self.discr.interior_faces_connection(qtag)
- elif dd_out.is_boundary():
- conn = self.discr.boundary_connection(dd_out.domain_tag, qtag)
- else:
- raise ValueError("cannot interpolate from volume to: " + str(dd_out))
-
- elif dd_in.domain_tag is sym.FACE_RESTR_INTERIOR:
- if dd_out.domain_tag is sym.FACE_RESTR_ALL:
- conn = self.discr.all_faces_connection(None, qtag)
- else:
- raise ValueError(
- "cannot interpolate from interior faces to: "
- + str(dd_out))
-
- elif dd_in.is_boundary():
- if dd_out.domain_tag is sym.FACE_RESTR_ALL:
- conn = self.discr.all_faces_connection(dd_in.domain_tag, qtag)
- else:
- raise ValueError(
- "cannot interpolate from interior faces to: "
- + str(dd_out))
-
- else:
- raise ValueError("cannot interpolate from: " + str(dd_in))
-
+ conn = self.discrwb.connection_from_dds(op.dd_in, op.dd_out)
return conn(self.queue, self.rec(field_expr)).with_queue(self.queue)
def map_opposite_partition_face_swap(self, op, field_expr):
@@ -301,7 +267,8 @@ class ExecutionMapper(mappers.Evaluator,
if qtag is None:
# FIXME: Remove once proper quadrature support arrives
qtag = sym.QTAG_NONE
- return self.discr.opposite_face_connection(qtag)(
+
+ return self.discrwb.opposite_face_connection(qtag)(
self.queue, self.rec(field_expr)).with_queue(self.queue)
# {{{ face mass operator
@@ -327,12 +294,7 @@ class ExecutionMapper(mappers.Evaluator,
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
return lp.tag_inames(knl, dict(k="g.0"))
- qtag = op.dd_in.quadrature_tag
- if qtag is None:
- # FIXME: Remove once proper quadrature support arrives
- qtag = sym.QTAG_NONE
-
- all_faces_conn = self.discr.all_faces_volume_connection(qtag)
+ all_faces_conn = self.discrwb.connection_from_dds("vol", op.dd_in)
all_faces_discr = all_faces_conn.to_discr
vol_discr = all_faces_conn.from_discr
@@ -370,14 +332,13 @@ class ExecutionMapper(mappers.Evaluator,
# }}}
# {{{ code execution functions
-
def map_insn_loopy_kernel(self, insn):
kwargs = {}
kdescr = insn.kernel_descriptor
for name, expr in six.iteritems(kdescr.input_mappings):
kwargs[name] = self.rec(expr)
- discr = self.get_discr(kdescr.governing_dd)
+ discr = self.discrwb.discr_from_dd(kdescr.governing_dd)
for name in kdescr.scalar_args():
v = kwargs[name]
if isinstance(v, (int, float)):
@@ -402,14 +363,14 @@ class ExecutionMapper(mappers.Evaluator,
assignments = []
for name, expr in zip(insn.names, insn.exprs):
value = self.rec(expr)
- self.discr._discr_scoped_subexpr_name_to_value[name] = value
+ self.discrwb._discr_scoped_subexpr_name_to_value[name] = value
assignments.append((name, value))
return assignments, []
def map_insn_assign_from_discr_scoped(self, insn):
return [(insn.name,
- self.discr._discr_scoped_subexpr_name_to_value[insn.name])], []
+ self.discrwb._discr_scoped_subexpr_name_to_value[insn.name])], []
def map_insn_diff_batch_assign(self, insn):
field = self.rec(insn.field)
@@ -419,13 +380,13 @@ class ExecutionMapper(mappers.Evaluator,
# This should be unified with map_elementwise_linear, which should
# be extended to support batching.
- discr = self.get_discr(repr_op.dd_in)
+ discr = self.discrwb.discr_from_dd(repr_op.dd_in)
# FIXME: Enable
# assert repr_op.dd_in == repr_op.dd_out
assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag
- @memoize_in(self.discr, "reference_derivative_knl")
+ @memoize_in(self.discrwb, "reference_derivative_knl")
def knl():
knl = lp.make_kernel(
"""{[imatrix,k,i,j]:
@@ -473,8 +434,8 @@ class ExecutionMapper(mappers.Evaluator,
# {{{ bound operator
class BoundOperator(object):
- def __init__(self, discr, discr_code, eval_code, debug_flags, allocator=None):
- self.discr = discr
+ def __init__(self, discrwb, discr_code, eval_code, debug_flags, allocator=None):
+ self.discrwb = discrwb
self.discr_code = discr_code
self.eval_code = eval_code
self.operator_data_cache = {}
@@ -504,13 +465,16 @@ class BoundOperator(object):
from pytools.obj_array import with_object_array_or_scalar
- # {{{ discr-scope evaluation
+ # {{{ discrwb-scope evaluation
- if any(result_var.name not in self.discr._discr_scoped_subexpr_name_to_value
+ if any(
+ (result_var.name not in
+ self.discrwb._discr_scoped_subexpr_name_to_value)
for result_var in self.discr_code.result):
- # need to do discr-scope evaluation
- discr_eval_context = {}
- self.discr_code.execute(ExecutionMapper(queue, discr_eval_context, self))
+ # need to do discrwb-scope evaluation
+ discrwb_eval_context = {}
+ self.discr_code.execute(
+ ExecutionMapper(queue, discrwb_eval_context, self))
# }}}
@@ -518,34 +482,49 @@ class BoundOperator(object):
for name, var in six.iteritems(context):
new_context[name] = with_object_array_or_scalar(replace_queue, var)
- return self.eval_code.execute(ExecutionMapper(queue, new_context, self))
+ return self.eval_code.execute(
+ ExecutionMapper(queue, new_context, self))
# }}}
# {{{ process_sym_operator function
-def process_sym_operator(sym_operator, post_bind_mapper=None,
- dumper=lambda name, sym_operator: None, mesh=None):
+def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None,
+ dumper=lambda name, sym_operator: None):
import grudge.symbolic.mappers as mappers
dumper("before-bind", sym_operator)
sym_operator = mappers.OperatorBinder()(sym_operator)
- mappers.ErrorChecker(mesh)(sym_operator)
+ mappers.ErrorChecker(discrwb.mesh)(sym_operator)
if post_bind_mapper is not None:
dumper("before-postbind", sym_operator)
sym_operator = post_bind_mapper(sym_operator)
- if mesh is not None:
- dumper("before-empty-flux-killer", sym_operator)
- sym_operator = mappers.EmptyFluxKiller(mesh)(sym_operator)
+ dumper("before-empty-flux-killer", sym_operator)
+ sym_operator = mappers.EmptyFluxKiller(discrwb.mesh)(sym_operator)
dumper("before-cfold", sym_operator)
sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
+ # Work around https://github.com/numpy/numpy/issues/9438
+ #
+ # The idea is that we need 1j as an expression to survive
+ # until code generation time. If it is evaluated and combined
+ # with other constants, we will need to determine its size
+ # (as np.complex64/128) within the expression. But because
+ # of the above numpy bug, sized numbers are not likely to survive
+ # expression building--so that's why we step in here to fix that.
+
+ dumper("before-csize", sym_operator)
+ sym_operator = mappers.ConstantToNumpyConversionMapper(
+ real_type=discrwb.real_dtype.type,
+ complex_type=discrwb.complex_dtype.type,
+ )(sym_operator)
+
# Ordering restriction:
#
# - Must run constant fold before first type inference pass, because zeros,
@@ -566,9 +545,8 @@ def process_sym_operator(sym_operator, post_bind_mapper=None,
# because otherwise it won't differentiate the type of grids (node or quadrature
# grids) that the operators will apply on.
- assert mesh is not None
dumper("before-global-to-reference", sym_operator)
- sym_operator = mappers.GlobalToReferenceMapper(mesh.ambient_dim)(sym_operator)
+ sym_operator = mappers.GlobalToReferenceMapper(discrwb.ambient_dim)(sym_operator)
dumper("before-distributed", sym_operator)
from meshmode.distributed import get_connected_partitions
@@ -614,10 +592,11 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x,
pretty(sym_operator))
stage[0] += 1
- sym_operator = process_sym_operator(sym_operator,
+ sym_operator = process_sym_operator(
+ discr,
+ sym_operator,
post_bind_mapper=post_bind_mapper,
- dumper=dump_optemplate,
- mesh=discr.mesh)
+ dumper=dump_optemplate)
from grudge.symbolic.compiler import OperatorCompiler
discr_code, eval_code = OperatorCompiler(discr)(sym_operator)
diff --git a/grudge/models/em.py b/grudge/models/em.py
index 3a2b70cbac3d00dab20420e060ceb073f6c1447a..f0e44f903d226e28710709651e948293a534d1b8 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -1,11 +1,13 @@
# -*- coding: utf8 -*-
"""grudge operators modelling electromagnetic phenomena."""
-from __future__ import division
-from __future__ import absolute_import
-from six.moves import range
+from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2007, 2010 Andreas Kloeckner, David Powell"
+__copyright__ = """
+Copyright (C) 2007-2017 Andreas Kloeckner
+Copyright (C) 2010 David Powell
+Copyright (C) 2017 Bogdan Enache
+"""
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -27,11 +29,13 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+from six.moves import range
+
from pytools import memoize_method
from grudge.models import HyperbolicOperator
-from grudge.symbolic.primitives import make_common_subexpression as cse
from meshmode.mesh import BTAG_ALL, BTAG_NONE
+from grudge import sym
from pytools.obj_array import join_fields, make_obj_array
# TODO: Check PML
@@ -106,17 +110,8 @@ class MaxwellOperator(HyperbolicOperator):
self.current = current
self.incident_bc_data = incident_bc
- @property
- def c(self):
- from warnings import warn
- warn("MaxwellOperator.c is deprecated", DeprecationWarning)
- if not self.fixed_material:
- raise RuntimeError("Cannot compute speed of light "
- "for non-constant material")
-
- return 1/(self.mu*self.epsilon)**0.5
- def flux(self, flux_type):
+ def flux(self, w):
"""The template for the numerical flux for variable coefficients.
:param flux_type: can be in [0,1] for anything between central and upwind,
@@ -124,88 +119,72 @@ class MaxwellOperator(HyperbolicOperator):
As per Hesthaven and Warburton page 433.
"""
- from grudge.flux import (make_normal, FluxVectorPlaceholder,
- FluxConstantPlaceholder)
- normal = make_normal(self.dimensions)
+ normal = sym.normal(w.dd, self.dimensions)
if self.fixed_material:
- from grudge.tools import count_subset
- w = FluxVectorPlaceholder(count_subset(self.get_eh_subset()))
-
e, h = self.split_eh(w)
- epsilon = FluxConstantPlaceholder(self.epsilon)
- mu = FluxConstantPlaceholder(self.mu)
-
- else:
- from grudge.tools import count_subset
- w = FluxVectorPlaceholder(count_subset(self.get_eh_subset())+2)
+ epsilon = self.epsilon
+ mu = self.mu
- epsilon, mu, e, h = self.split_eps_mu_eh(w)
- Z_int = (mu.int/epsilon.int)**0.5
+ Z_int = (mu/epsilon)**0.5
Y_int = 1/Z_int
- Z_ext = (mu.ext/epsilon.ext)**0.5
+ Z_ext = (mu/epsilon)**0.5
Y_ext = 1/Z_ext
- if flux_type == "lf":
+ if self.flux_type == "lf":
if self.fixed_material:
max_c = (self.epsilon*self.mu)**(-0.5)
- else:
- from grudge.flux import Max
- c_int = (epsilon.int*mu.int)**(-0.5)
- c_ext = (epsilon.ext*mu.ext)**(-0.5)
- max_c = Max(c_int, c_ext) # noqa
return join_fields(
# flux e,
1/2*(
-self.space_cross_h(normal, h.int-h.ext)
# multiplication by epsilon undoes material divisor below
- #-max_c*(epsilon.int*e.int - epsilon.ext*e.ext)
+ #-max_c*(epsilon*e.int - epsilon*e.ext)
),
# flux h
1/2*(
self.space_cross_e(normal, e.int-e.ext)
# multiplication by mu undoes material divisor below
- #-max_c*(mu.int*h.int - mu.ext*h.ext)
+ #-max_c*(mu*h.int - mu*h.ext)
))
- elif isinstance(flux_type, (int, float)):
+ elif isinstance(self.flux_type, (int, float)):
# see doc/maxima/maxwell.mac
return join_fields(
# flux e,
(
-1/(Z_int+Z_ext)*self.space_cross_h(normal,
Z_ext*(h.int-h.ext)
- - flux_type*self.space_cross_e(normal, e.int-e.ext))
+ - self.flux_type*self.space_cross_e(normal, e.int-e.ext))
),
# flux h
(
1/(Y_int + Y_ext)*self.space_cross_e(normal,
Y_ext*(e.int-e.ext)
- + flux_type*self.space_cross_h(normal, h.int-h.ext))
+ + self.flux_type*self.space_cross_h(normal, h.int-h.ext))
),
)
else:
raise ValueError("maxwell: invalid flux_type (%s)"
% self.flux_type)
- def local_derivatives(self, w=None):
+ def local_derivatives(self, w):
"""Template for the spatial derivatives of the relevant components of
:math:`E` and :math:`H`
"""
- e, h = self.split_eh(self.field_placeholder(w))
+ e, h = self.split_eh(w)
+ nabla = sym.nabla(self.dimensions)
def e_curl(field):
return self.space_cross_e(nabla, field)
def h_curl(field):
return self.space_cross_h(nabla, field)
- from grudge.symbolic import make_nabla
- nabla = make_nabla(self.dimensions)
# in conservation form: u_t + A u_x = 0
return join_fields(
@@ -213,62 +192,44 @@ class MaxwellOperator(HyperbolicOperator):
e_curl(e)
)
- def field_placeholder(self, w=None):
- "A placeholder for E and H."
- from grudge.tools import count_subset
- fld_cnt = count_subset(self.get_eh_subset())
- if w is None:
- from grudge.symbolic import make_sym_vector
- w = make_sym_vector("w", fld_cnt)
-
- return w
-
- def pec_bc(self, w=None):
+ def pec_bc(self, w):
"Construct part of the flux operator template for PEC boundary conditions"
- e, h = self.split_eh(self.field_placeholder(w))
+ e, h = self.split_eh(w)
+
+ pec_e = sym.cse(sym.interp("vol", self.pec_tag)(e))
+ pec_h = sym.cse(sym.interp("vol", self.pec_tag)(h))
- from grudge.symbolic import RestrictToBoundary
- pec_e = RestrictToBoundary(self.pec_tag)(e)
- pec_h = RestrictToBoundary(self.pec_tag)(h)
return join_fields(-pec_e, pec_h)
- def pmc_bc(self, w=None):
+ def pmc_bc(self, w):
"Construct part of the flux operator template for PMC boundary conditions"
- e, h = self.split_eh(self.field_placeholder(w))
+ e, h = self.split_eh(w)
- from grudge.symbolic import RestrictToBoundary
- pmc_e = RestrictToBoundary(self.pmc_tag)(e)
- pmc_h = RestrictToBoundary(self.pmc_tag)(h)
+ pmc_e = sym.cse(sym.interp("vol", self.pmc_tag)(e))
+ pmc_h = sym.cse(sym.interp("vol", self.pmc_tag)(h))
return join_fields(pmc_e, -pmc_h)
- def absorbing_bc(self, w=None):
+ def absorbing_bc(self, w):
"""Construct part of the flux operator template for 1st order
absorbing boundary conditions.
"""
- from grudge.symbolic import normal
- absorb_normal = normal(self.absorb_tag, self.dimensions)
+ absorb_normal = sym.normal(self.absorb_tag, self.dimensions)
- from grudge.symbolic import RestrictToBoundary, Field
- e, h = self.split_eh(self.field_placeholder(w))
+ e, h = self.split_eh(w)
if self.fixed_material:
epsilon = self.epsilon
mu = self.mu
- else:
- epsilon = cse(
- RestrictToBoundary(self.absorb_tag)(Field("epsilon")))
- mu = cse(
- RestrictToBoundary(self.absorb_tag)(Field("mu")))
absorb_Z = (mu/epsilon)**0.5
absorb_Y = 1/absorb_Z
- absorb_e = RestrictToBoundary(self.absorb_tag)(e)
- absorb_h = RestrictToBoundary(self.absorb_tag)(h)
+ absorb_e = sym.cse(sym.interp("vol", self.absorb_tag)(e))
+ absorb_h = sym.cse(sym.interp("vol", self.absorb_tag)(h))
bc = join_fields(
absorb_e + 1/2*(self.space_cross_h(absorb_normal, self.space_cross_e(
@@ -281,17 +242,12 @@ class MaxwellOperator(HyperbolicOperator):
return bc
- def incident_bc(self, w=None):
+ def incident_bc(self, w):
"Flux terms for incident boundary conditions"
# NOTE: Untested for inhomogeneous materials, but would usually be
# physically meaningless anyway (are there exceptions to this?)
- e, h = self.split_eh(self.field_placeholder(w))
- if not self.fixed_material:
- from warnings import warn
- if self.incident_tag != BTAG_NONE:
- warn("Incident boundary conditions assume homogeneous"
- " background material, results may be unphysical")
+ e, h = self.split_eh(w)
from grudge.tools import count_subset
fld_cnt = count_subset(self.get_eh_subset())
@@ -301,7 +257,7 @@ class MaxwellOperator(HyperbolicOperator):
if is_zero(incident_bc_data):
return make_obj_array([0]*fld_cnt)
else:
- return cse(-incident_bc_data)
+ return sym.cse(-incident_bc_data)
def sym_operator(self, w=None):
"""The full operator template - the high level description of
@@ -310,23 +266,9 @@ class MaxwellOperator(HyperbolicOperator):
Combines the relevant operator templates for spatial
derivatives, flux, boundary conditions etc.
"""
- w = self.field_placeholder(w)
-
- if self.fixed_material:
- flux_w = w
- else:
- epsilon = self.epsilon
- mu = self.mu
-
- flux_w = join_fields(epsilon, mu, w)
-
- from grudge.symbolic import BoundaryPair, \
- InverseMassOperator, get_flux_operator
-
- flux_op = get_flux_operator(self.flux(self.flux_type))
- bdry_flux_op = get_flux_operator(self.flux(self.bdry_flux_type))
+ from grudge.tools import count_subset
+ w = sym.make_sym_array("w", count_subset(self.get_eh_subset()))
- from grudge.tools.indexing import count_subset
elec_components = count_subset(self.get_eh_subset()[0:3])
mag_components = count_subset(self.get_eh_subset()[3:6])
@@ -334,10 +276,6 @@ class MaxwellOperator(HyperbolicOperator):
# need to check this
material_divisor = (
[self.epsilon]*elec_components+[self.mu]*mag_components)
- else:
- material_divisor = join_fields(
- [epsilon]*elec_components,
- [mu]*mag_components)
tags_and_bcs = [
(self.pec_tag, self.pec_bc(w)),
@@ -346,67 +284,20 @@ class MaxwellOperator(HyperbolicOperator):
(self.incident_tag, self.incident_bc(w)),
]
- def make_flux_bc_vector(tag, bc):
- if self.fixed_material:
- return bc
- else:
- from grudge.symbolic import RestrictToBoundary
- return join_fields(
- cse(RestrictToBoundary(tag)(epsilon)),
- cse(RestrictToBoundary(tag)(mu)),
- bc)
+
+ def flux(pair):
+ return sym.interp(pair.dd, "all_faces")(self.flux(pair))
return (
- self.local_derivatives(w)
- + InverseMassOperator()(
- flux_op(flux_w)
+ - sym.InverseMassOperator()(sym.FaceMassOperator()(
+ flux(sym.int_tpair(w))
+ sum(
- bdry_flux_op(BoundaryPair(
- flux_w, make_flux_bc_vector(tag, bc), tag))
- for tag, bc in tags_and_bcs))
- ) / material_divisor
-
- def bind(self, discr):
- "Convert the abstract operator template into compiled code."
- from grudge.mesh import check_bc_coverage
- check_bc_coverage(discr.mesh, [
- self.pec_tag, self.absorb_tag, self.incident_tag])
-
- compiled_sym_operator = discr.compile(self.op_template())
-
- def rhs(t, w, **extra_context):
- kwargs = {}
- kwargs.update(extra_context)
-
- return compiled_sym_operator(w=w, t=t, **kwargs)
-
- return rhs
-
- def assemble_eh(self, e=None, h=None, discr=None):
- "Combines separate E and H vectors into a single array."
- if discr is None:
- def zero():
- return 0
- else:
- def zero():
- return discr.volume_zeros()
-
- from grudge.tools import count_subset
- e_components = count_subset(self.get_eh_subset()[0:3])
- h_components = count_subset(self.get_eh_subset()[3:6])
-
- def default_fld(fld, comp):
- if fld is None:
- return [zero() for i in range(comp)]
- else:
- return fld
-
- e = default_fld(e, e_components)
- h = default_fld(h, h_components)
+ flux(sym.bv_tpair(tag, w, bc))
+ for tag, bc in tags_and_bcs)
+ ) )) / material_divisor
- return join_fields(e, h)
- assemble_fields = assemble_eh
@memoize_method
def partial_to_eh_subsets(self):
@@ -445,11 +336,8 @@ class MaxwellOperator(HyperbolicOperator):
e_idx, h_idx = self.partial_to_eh_subsets()
e, h = w[e_idx], w[h_idx]
- from grudge.flux import FluxVectorPlaceholder as FVP
- if isinstance(w, FVP):
- return FVP(scalars=e), FVP(scalars=h)
- else:
- return make_obj_array(e), make_obj_array(h)
+ return e,h
+ #return make_obj_array(e), make_obj_array(h)
def get_eh_subset(self):
"""Return a 6-tuple of :class:`bool` objects indicating whether field
@@ -476,6 +364,14 @@ class MaxwellOperator(HyperbolicOperator):
raise ValueError("max_eigenvalue is no longer supported for "
"variable-coefficient problems--use max_eigenvalue_expr")
+ def check_bc_coverage(self, mesh):
+ from meshmode.mesh import check_bc_coverage
+ check_bc_coverage(mesh, [
+ self.pec_tag,
+ self.pmc_tag,
+ self.absorb_tag,
+ self.incident_tag])
+
class TMMaxwellOperator(MaxwellOperator):
"""A 2D TM Maxwell operator with PEC boundaries.
@@ -539,3 +435,61 @@ class SourceFree1DMaxwellOperator(MaxwellOperator):
+
(False, False, True)
)
+
+def get_rectangular_cavity_mode(E_0, mode_indices):
+ """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)
+ if dims != 2 and dims != 3:
+ raise ValueError("Improper mode_indices dimensions")
+ import numpy
+
+ factors = [n*numpy.pi for n in mode_indices]
+
+ kx, ky = factors[0:2]
+ if dims == 3:
+ kz = factors[2]
+
+ 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)
+ if dims == 3:
+ sz = sym.sin(kz*z)
+ cz = sym.cos(kz*z)
+
+ if dims == 2:
+ tfac = sym.ScalarVariable("t") * omega
+
+ result = sym.join_fields(
+ 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,
+ )
+ else:
+ tdep = sym.exp(-1j * omega * sym.ScalarVariable("t"))
+
+ gamma_squared = ky**2 + kx**2
+ result = sym.join_fields(
+ -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,
+ )
+
+ return result
diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py
index 0e6940ff6f120018a8b5c3c9af29de04eb1bddcb..bb12689bf65208c3a48bb299f2a44c6651c5b58b 100644
--- a/grudge/shortcuts.py
+++ b/grudge/shortcuts.py
@@ -44,13 +44,16 @@ def set_up_rk4(field_var_name, dt, fields, rhs, t_start=0):
return dt_stepper
-def make_visualizer(discr, vis_order):
+def make_visualizer(discrwb, vis_order):
from meshmode.discretization.visualization import make_visualizer
- with cl.CommandQueue(discr.cl_context) as queue:
- return make_visualizer(queue, discr.volume_discr, vis_order)
+ with cl.CommandQueue(discrwb.cl_context) as queue:
+ return make_visualizer(queue, discrwb.discr_from_dd("vol"), vis_order)
-def make_boundary_visualizer(discr, vis_order):
+def make_boundary_visualizer(discrwb, vis_order):
from meshmode.discretization.visualization import make_visualizer
- with cl.CommandQueue(discr.cl_context) as queue:
- return make_visualizer(queue, discr.boundary_discr, vis_order)
+ from grudge import sym
+ with cl.CommandQueue(discrwb.cl_context) as queue:
+ return make_visualizer(
+ queue, discrwb.discr_from_dd(sym.BTAG_ALL),
+ vis_order)
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index ac7765ed919f87ed7b8a17099a302263b94ce347..f6e238083266d4ba0c510243b116b2e49b024631 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -852,7 +852,11 @@ class ToLoopyExpressionMapper(mappers.IdentityMapper):
def map_call(self, expr):
if isinstance(expr.function, sym.CFunction):
from pymbolic import var
- return var(expr.function.name)(
+ func_name = expr.function.name
+ if func_name == "fabs":
+ func_name = "abs"
+
+ return var(func_name)(
*[self.rec(par) for par in expr.parameters])
else:
raise NotImplementedError(
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index a0d16d42ad9b1393592193cc1252ef5a096967d3..ddba6c8d6a301f0c482e2e95ae67e887345d49f3 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -34,6 +34,7 @@ import pymbolic.mapper.evaluator
import pymbolic.mapper.dependency
import pymbolic.mapper.substitutor
import pymbolic.mapper.constant_folder
+import pymbolic.mapper.constant_converter
import pymbolic.mapper.flop_counter
from pymbolic.mapper import CSECachingMapperMixin
@@ -893,6 +894,12 @@ class QuadratureUpsamplerChanger(CSECachingMapperMixin, IdentityMapper):
# {{{ simplification / optimization
+class ConstantToNumpyConversionMapper(
+ pymbolic.mapper.constant_converter.ConstantToNumpyConversionMapper,
+ IdentityMapperMixin):
+ pass
+
+
class CommutativeConstantFoldingMapper(
pymbolic.mapper.constant_folder.CommutativeConstantFoldingMapper,
IdentityMapperMixin):
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index c4f6ed659613006337d82bf6e1e477437f2642e0..76bd938c941d5b7e0eea5682f164f188f3e2aff0 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -540,10 +540,10 @@ def norm(p, arg, dd=None):
return sym.CFunction("sqrt")(norm_squared)
elif p == np.Inf:
- result = sym.NodalMax()(sym.CFunction("fabs")(arg))
+ result = sym.NodalMax(dd_in=dd)(sym.CFunction("fabs")(arg))
from pymbolic.primitives import Max
- if isinstance(norm_squared, np.ndarray):
+ if isinstance(result, np.ndarray):
from functools import reduce
result = reduce(Max, result)
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 761fc7a7b6e46b6c50786ceeb5e355c4166a51a6..3281abef525697ee692aacba76a18ba39a574664 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -185,6 +185,7 @@ class DOFDesc(object):
elif isinstance(domain_tag, BTAG_PARTITION):
pass
elif domain_tag in [BTAG_ALL, BTAG_REALLY_ALL, BTAG_NONE]:
+ # FIXME: Should wrap these in DTAG_BOUNDARY
pass
elif isinstance(domain_tag, DTAG_BOUNDARY):
pass
diff --git a/grudge/tools.py b/grudge/tools.py
index cd1cc49f1828ba82f95c0e805262315e5f8e2be6..88d6714364e250d137ee7988558c91138acc1eb7 100644
--- a/grudge/tools.py
+++ b/grudge/tools.py
@@ -25,6 +25,7 @@ THE SOFTWARE.
"""
import numpy as np
+from pytools import levi_civita
def is_zero(x):
@@ -36,3 +37,90 @@ def is_zero(x):
return x == 0
else:
return False
+
+
+class SubsettableCrossProduct:
+ """A cross product that can operate on an arbitrary subsets of its
+ two operands and return an arbitrary subset of its result.
+ """
+
+ full_subset = (True, True, True)
+
+ def __init__(self, op1_subset=full_subset, op2_subset=full_subset,
+ result_subset=full_subset):
+ """Construct a subset-able cross product.
+ :param op1_subset: The subset of indices of operand 1 to be taken into
+ account. Given as a 3-sequence of bools.
+ :param op2_subset: The subset of indices of operand 2 to be taken into
+ account. Given as a 3-sequence of bools.
+ :param result_subset: The subset of indices of the result that are
+ calculated. Given as a 3-sequence of bools.
+ """
+ def subset_indices(subset):
+ return [i for i, use_component in enumerate(subset)
+ if use_component]
+
+ self.op1_subset = op1_subset
+ self.op2_subset = op2_subset
+ self.result_subset = result_subset
+
+ import pymbolic
+ op1 = pymbolic.var("x")
+ op2 = pymbolic.var("y")
+
+ self.functions = []
+ self.component_lcjk = []
+ for i, use_component in enumerate(result_subset):
+ if use_component:
+ this_expr = 0
+ this_component = []
+ for j, j_real in enumerate(subset_indices(op1_subset)):
+ for k, k_real in enumerate(subset_indices(op2_subset)):
+ lc = levi_civita((i, j_real, k_real))
+ if lc != 0:
+ this_expr += lc*op1.index(j)*op2.index(k)
+ this_component.append((lc, j, k))
+ self.functions.append(pymbolic.compile(this_expr,
+ variables=[op1, op2]))
+ self.component_lcjk.append(this_component)
+
+ def __call__(self, x, y, three_mult=None):
+ """Compute the subsetted cross product on the indexables *x* and *y*.
+ :param three_mult: a function of three arguments *sign, xj, yk*
+ used in place of the product *sign*xj*yk*. Defaults to just this
+ product if not given.
+ """
+ from pytools.obj_array import join_fields
+ if three_mult is None:
+ return join_fields(*[f(x, y) for f in self.functions])
+ else:
+ return join_fields(
+ *[sum(three_mult(lc, x[j], y[k]) for lc, j, k in lcjk)
+ for lcjk in self.component_lcjk])
+
+
+cross = SubsettableCrossProduct()
+
+
+def count_subset(subset):
+ from pytools import len_iterable
+ return len_iterable(uc for uc in subset if uc)
+
+
+def partial_to_all_subset_indices(subsets, base=0):
+ """Takes a sequence of bools and generates it into an array of indices
+ to be used to insert the subset into the full set.
+ Example:
+ >>> list(partial_to_all_subset_indices([[False, True, True], [True,False,True]]))
+ [array([0 1]), array([2 3]
+ """
+
+ idx = base
+ for subset in subsets:
+ result = []
+ for is_in in subset:
+ if is_in:
+ result.append(idx)
+ idx += 1
+
+ yield np.array(result, dtype=np.intp)
diff --git a/test/test_grudge.py b/test/test_grudge.py
index bd47f823c36931290c1b7ba9981c5c8b41a7a032..84c46168ede26efe1a42935501d69fba99b7700a 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -37,7 +37,7 @@ from pyopencl.tools import ( # noqa
import logging
logger = logging.getLogger(__name__)
-from grudge import sym, bind, Discretization
+from grudge import sym, bind, DGDiscretizationWithBoundaries
@pytest.mark.parametrize("dim", [2, 3])
@@ -64,7 +64,7 @@ def test_inverse_metric(ctx_factory, dim):
from meshmode.mesh.processing import map_mesh
mesh = map_mesh(mesh, m)
- discr = Discretization(cl_ctx, mesh, order=4)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=4)
sym_op = (
sym.forward_metric_derivative_mat(mesh.dim)
@@ -98,7 +98,7 @@ def test_1d_mass_mat_trig(ctx_factory):
mesh = generate_regular_rect_mesh(a=(-4*np.pi,), b=(9*np.pi,),
n=(17,), order=1)
- discr = Discretization(cl_ctx, mesh, order=8)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=8)
x = sym.nodes(1)
f = bind(discr, sym.cos(x[0])**2)(queue)
@@ -139,7 +139,7 @@ def test_tri_diff_mat(ctx_factory, dim, order=4):
mesh = generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim,
n=(n,)*dim, order=4)
- discr = Discretization(cl_ctx, mesh, order=4)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=4)
nabla = sym.nabla(dim)
for axis in range(dim):
@@ -182,7 +182,7 @@ def test_2d_gauss_theorem(ctx_factory):
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
- discr = Discretization(cl_ctx, mesh, order=2)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=2)
def f(x):
return sym.join_fields(
@@ -272,7 +272,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
from grudge.models.advection import (
StrongAdvectionOperator, WeakAdvectionOperator)
- discr = Discretization(cl_ctx, mesh, order=order)
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
op_class = {
"strong": StrongAdvectionOperator,
"weak": WeakAdvectionOperator,
@@ -333,6 +333,76 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
assert eoc_rec.order_estimate() > order
+@pytest.mark.parametrize("order", [3, 4, 5])
+# test: 'test_convergence_advec(cl._csc, "disk", [0.1, 0.05], "strong", "upwind", 3)'
+def test_convergence_maxwell(ctx_factory, order, visualize=False):
+ """Test whether 3D maxwells actually converges"""
+
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from pytools.convergence import EOCRecorder
+ eoc_rec = EOCRecorder()
+
+ dims = 3
+ ns = [4, 6, 8]
+ for n in ns:
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(0.0,)*dims,
+ b=(1.0,)*dims,
+ n=(n,)*dims)
+
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+
+ epsilon = 1
+ mu = 1
+
+ from grudge.models.em import get_rectangular_cavity_mode
+ sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2))
+
+ analytic_sol = bind(discr, sym_mode)
+ fields = analytic_sol(queue, t=0, epsilon=epsilon, mu=mu)
+
+ from grudge.models.em import MaxwellOperator
+ op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
+ op.check_bc_coverage(mesh)
+ bound_op = bind(discr, op.sym_operator())
+
+ def rhs(t, w):
+ return bound_op(queue, t=t, w=w)
+
+ dt = 0.002
+ final_t = dt * 5
+ nsteps = int(final_t/dt)
+
+ from grudge.shortcuts import set_up_rk4
+ dt_stepper = set_up_rk4("w", dt, fields, rhs)
+
+ print("dt=%g nsteps=%d" % (dt, nsteps))
+
+ norm = bind(discr, sym.norm(2, sym.var("u")))
+
+ step = 0
+ for event in dt_stepper.run(t_end=final_t):
+ if isinstance(event, dt_stepper.StateComputed):
+ assert event.component_id == "w"
+ esc = event.state_component
+
+ step += 1
+ print(step)
+
+ sol = analytic_sol(queue, mu=mu, epsilon=epsilon, t=step * dt)
+ vals = [norm(queue, u=(esc[i] - sol[i])) / norm(queue, u=sol[i]) for i in range(5)] # noqa E501
+ total_error = sum(vals)
+ eoc_rec.add_data_point(1.0/n, total_error)
+
+ print(eoc_rec.pretty_print(abscissa_label="h",
+ error_label="L2 Error"))
+
+ assert eoc_rec.order_estimate() > order
+
+
def test_foreign_points(ctx_factory):
pytest.importorskip("sumpy")
import sumpy.point_calculus as pc