diff --git a/contrib/maxima/em-modes.mac b/contrib/maxima/em-modes.mac
new file mode 100644
index 0000000000000000000000000000000000000000..6482ec1774c2f2f0ba63db7b9f69406f601e53f1
--- /dev/null
+++ b/contrib/maxima/em-modes.mac
@@ -0,0 +1,157 @@
+kill(all);
+load("eigen");
+load("itensor");
+load("diag");
+
+/* -------------------------------------------------------------------------- */
+/* Utilities */
+/* -------------------------------------------------------------------------- */
+
+curl(x):=crossfunc(lambda([i,j], diff(x[j], coords[i])));
+div(f):=sum(diff(f[j], coords[j]), j, 1, length(coords));
+
+assert(condition):=if not condition then error("Assertion violated") else true$
+
+norm_2_squared(v):=v.v;
+
+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]));
+
+/* -------------------------------------------------------------------------- */
+/* Variable declarations */
+/* -------------------------------------------------------------------------- */
+
+coords:[x,y,z];
+allvars:append([t],coords);
+
+mu: 1;
+epsilon: 1;
+c: 1/sqrt(mu*epsilon);
+
+epsilon_0: 1;
+mu_0: 1;
+c_0: 1/sqrt(mu_0*epsilon_0);
+
+max_B(max_H):= mu*max_H;
+max_D(max_E):= epsilon*max_E;
+
+/* SI conventions, assumed time dep: exp(%i*omega*t) */
+faraday(max_E, max_H):= curl(max_E) + %i * omega * max_B(max_H);
+ampere(max_E, max_H):= curl(max_B(max_H)) - %i * omega / c_0**2 * max_E;
+
+div_e(max_E, max_H):= div(max_E);
+div_h(max_E, max_H):= div(max_H);
+
+maxwell_pde(max_E, max_H):=append(
+ faraday(max_E, max_H),
+ ampere(max_E, max_H),
+ [div_e(max_E, max_H), div_h(max_E, max_H)]);
+
+/*
+spatial_deps:[
+ exp(%i*m*x)*exp(%i*n*y),
+ exp(%i*m*x)*exp(-%i*n*y),
+ exp(-%i*m*x)*exp(%i*n*y),
+ exp(-%i*m*x)*exp(-%i*n*y)
+ ];
+*/
+
+spatial_deps:[
+ exp(+%i*m*x)*exp(+%i*n*y)*exp(+%i*l*z),
+ exp(+%i*m*x)*exp(+%i*n*y)*exp(-%i*l*z),
+ exp(+%i*m*x)*exp(-%i*n*y)*exp(+%i*l*z),
+ exp(+%i*m*x)*exp(-%i*n*y)*exp(-%i*l*z)
+
+ /*
+ exp(-%i*m*x)*exp(+%i*n*y)*exp(+%i*l*z),
+ exp(-%i*m*x)*exp(+%i*n*y)*exp(-%i*l*z),
+ exp(-%i*m*x)*exp(-%i*n*y)*exp(+%i*l*z),
+ exp(-%i*m*x)*exp(-%i*n*y)*exp(-%i*l*z)
+ */
+ ];
+
+declare(m, integer, n, integer, l, integer);
+
+get_maxwell_solution(spatial_dep):=block([
+ max_B, max_D, coeffs, max_E, max_H, faraday, ampere, div_e, div_h, maxwell_pde, soln
+ ],
+ max_B: mu*max_H,
+ max_D: epsilon*max_E,
+
+ coeffs: [
+ Exr, Eyr, Ezr, Hxr, Hyr, Hzr,
+ Exi, Eyi, Ezi, Hxi, Hyi, Hzi
+ ],
+
+ max_E: [
+ (Exr+%i*Exi)*spatial_dep,
+ (Eyr+%i*Eyi)*spatial_dep,
+ (Ezr+%i*Ezi)*spatial_dep
+ ],
+ max_H: [
+ (Hxr+%i*Hxi)*spatial_dep,
+ (Hyr+%i*Hyi)*spatial_dep,
+ (Hzr+%i*Hzi)*spatial_dep
+ ],
+
+ soln:solve(
+ maxwell_pde(max_E, max_H),
+ append(coeffs, [omega])),
+
+ family1: soln[1],
+ omega_eq: family1[length(family1)],
+ assert(lhs(omega_eq) = omega),
+
+ [subst(family1, [max_E, max_H]), rhs(omega_eq)]
+ );
+
+maxwell_solutions:makelist(
+ get_maxwell_solution(spatial_deps[i]),
+ i, 1, length(spatial_deps));
+
+omegas:makelist(
+ maxwell_solutions[i][2],
+ i, 1, length(maxwell_solutions));
+display(omegas);
+
+max_E:ratsimp(sum(
+ maxwell_solutions[i][1][1],
+ i, 1, length(maxwell_solutions)));
+max_H:ratsimp(sum(
+ maxwell_solutions[i][1][2],
+ i, 1, length(maxwell_solutions)));
+
+print("Check Maxwell:");
+print(ratsimp(subst([omega=omegas[1]],maxwell_pde(max_E,max_H))));
+
+pec_bcs:append(
+ realpart(crossprod([-1,0,0], subst(x=0, max_E))),
+ realpart(crossprod([1,0,0], subst(x=%pi, max_E))),
+ realpart(crossprod([0,-1,0], subst(y=0, max_E))),
+ realpart(crossprod([0,1,0], subst(y=%pi, max_E))),
+ [
+ realpart([-1,0,0].subst(x=0, max_H)),
+ realpart([1,0,0].subst(x=%pi, max_H)),
+ realpart([0,-1,0].subst(y=0, max_H)),
+ realpart([0,1,0].subst(y=%pi, max_H))
+ ]);
+
+freevars: sublist(
+ listofvars([max_E, max_H]),
+ lambda([rvar], substring(string(rvar),1,2) = "%"));
+
+ev_pec_bcs:append(
+ subst([x=0, y=0, z=0], pec_bcs),
+ subst([x=0, y=0, z=%pi], pec_bcs),
+ subst([x=0, y=%pi, z=%pi], pec_bcs)
+ );
+
+/*
+Fails:
+
+pec_soln:linsolve(pec_bcs, freevars);
+*/
diff --git a/contrib/maxima/maxwell2.mac b/contrib/maxima/maxwell2.mac
new file mode 100644
index 0000000000000000000000000000000000000000..13b6f9ddf723288aae163e5e0312934130059b8d
--- /dev/null
+++ b/contrib/maxima/maxwell2.mac
@@ -0,0 +1,90 @@
+kill(all);
+load("eigen");
+load("itensor");
+load("diag");
+
+/* -------------------------------------------------------------------------- */
+/* Utilities */
+/* -------------------------------------------------------------------------- */
+
+coords:[x,y,z];
+mu: 1;
+epsilon: 1;
+c: 1/sqrt(mu*epsilon);
+
+epsilon_0: 1;
+mu_0: 1;
+c_0: 1/sqrt(mu_0*epsilon_0);
+max_B(max_H):= mu*max_H;
+max_D(max_E):= epsilon*max_E;
+
+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_B(max_H);
+ampere(max_E, max_H):= curl(max_B(max_H)) + %i * omega / c_0**2 * max_E;
+
+assert(condition):=if not condition then error("Assertion violated") else true$
+
+norm_2_squared(v):=v.v;
+
+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 */
+/* -------------------------------------------------------------------------- */
+
+/*
+psi_cand:E_0*sin(k_x*x)*sin(k_y*y);
+*/
+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);
+
+nabla_t_squared(f):=diff(f, x, 2) + diff(f, y, 2);
+nabla_t(f):= [diff(f, x, 1) , diff(f, y, 1), 0];
+
+
+wave_eqn(f, gam_s):=nabla_t_squared(f) + gam_s*f;
+gam_s : gamma^2;
+/*
+gam_s : omega^2 - k_z^2;
+gam_s : k_y^2 + k_x^2;
+*/
+omega2 : sqrt(k_y^2 + k_x^2 + k_z^2);
+
+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 : E_t(psi_cand) + [0,0,psi_cand * cos(k_z * z)];
+H :H_t(psi_cand);
+
+/*
+E : E, omega = sqrt(k_x^2 + k_y^2 + k_z^2), gamma = sqrt(k_y^2 + k_x^2);
+H : H, omega = sqrt(k_x^2 + k_y^2 + k_z^2), gamma = sqrt(k_y^2 + k_x^2);
+*/
+
+Etime : E * %e ^ (- %i * omega * t);
+Htime : H * %e ^ (- %i * omega * t);
+
+/*
+print(solve([div(E),div(H),faraday(E,H),ampere(E, H)], [gamma,omega]));
+print(solve(div(E), [gamma,omega]));
+print(solve(div(H), [gamma,omega]));
+print(solve(faraday(E,H), [gamma,omega]));
+print(solve(ampere(E,H), [gamma,omega]));
+
+
+print(solve([div(Etime),div(Htime),faraday(Etime,Htime),ampere(Etime, Htime)], [gamma,omega]));
+print(solve(div(Etime), [gamma,omega]));
+print(solve(div(Htime), [gamma,omega]));
+print(solve(faraday(Etime,Htime), [gamma,omega]));
+print(solve(ampere(Etime,Htime), [gamma,omega]));
+*/
+
diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py
index 7ff8884f1d19e46d320e0cc7a8c5301f74d14551..2bf0ee4668f10a850e785f0b47e811011016dd73 100644
--- a/examples/maxwell/analytic_solutions.py
+++ b/examples/maxwell/analytic_solutions.py
@@ -1,64 +1,40 @@
-# 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 .
-
-
-
-
-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
+from __future__ import division, absolute_import, print_function
+
+__copyright__ = """
+Copyright (C) 2007-17 Andreas Kloeckner
+Copyright (C) 2017 Bogdan Enache
+"""
+
+__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:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+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.
+"""
+
+
+from six.moves import range, zip
+import numpy as np
import numpy.linalg as la
-import cmath
-from six.moves import range
-from six.moves import zip
+from grudge import sym
-# 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]
-
-
-
+# {{{ cylindrical field adapter
class CylindricalFieldAdapter:
"""Turns a field returned as (r, phi, z) components into
@@ -90,8 +66,10 @@ class CylindricalFieldAdapter:
return result
+# }}}
+# {{{ spherical field adapter
class SphericalFieldAdapter:
"""Turns the adaptee's field C{(Er, Etheta, Ephi)} components into
@@ -138,10 +116,11 @@ class SphericalFieldAdapter:
return result
+# }}}
+# {{{ cylindrical cavity mode
-# actual solutions ------------------------------------------------------------
class CylindricalCavityMode:
"""A cylindrical TM cavity mode.
@@ -150,7 +129,7 @@ class CylindricalCavityMode:
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
@@ -211,7 +190,7 @@ class CylindricalCavityMode:
# 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)
+ psi_dphi = (cyl_bessel_j(m, gamma * r)
* 1/r * phi_sign*1j*m * phi_factor)
# field components in polar coordinates ---------------------------
@@ -237,91 +216,265 @@ class CylindricalCavityMode:
return [er, ephi, ez, hr, hphi, hz]
+# }}}
+# {{{ symmetric rectangular waveguide mode
-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
+def get_rectangular_waveguide_mode(
+ ambient_dim, mode_indices,
+ dimensions=(1, 1, 1), coefficients=(1, 0, 0),
+ forward_coeff=1, backward_coeff=0,
+ ):
+ """A TM cavity mode.
- self.epsilon = epsilon
- self.mu = mu
+ Returns an expression depending on *epsilon*, *mu*, and *t*.
+ """
- self.t = 0
+ kx, ky, h = factors = [n*np.pi/a for n, a in zip(mode_indices, dimensions)]
+ k = np.sqrt(sum(f**2 for f in factors))
- self.factors = [n*pi/a for n, a in zip(self.mode_indices, self.dimensions)]
+ epsilon = sym.ScalarVariable("epsilon")
+ mu = sym.ScalarVariable("mu")
- c = 1/sqrt(mu*epsilon)
- self.k = sqrt(sum(f**2 for f in self.factors))
- self.omega = self.k*c
+ c = sym.cse(1/sym.sqrt(mu*epsilon), "c")
+ omega = k*c
- def set_time(self, t):
- self.t = t
+ nodes = sym.nodes(ambient_dim)
+ x = nodes[0]
+ y = nodes[1]
+ if nodes.shape[0] > 2:
+ z = nodes[2]
+ else:
+ z = 0
- def __call__(self, discr):
- f,g,h = self.factors
- omega = self.omega
- k = self.k
+ sx = sym.sin(kx*x)
+ sy = sym.sin(ky*y)
+ cx = sym.cos(kx*x)
+ cy = sym.cos(ky*y)
- x = discr.nodes[:,0]
- y = discr.nodes[:,1]
- if discr.dimensions > 2:
- z = discr.nodes[:,2]
- else:
- z = discr.volume_zeros()
+ zdep_add = sym.exp(1j*h*z)+sym.exp(-1j*h*z)
+ zdep_sub = sym.exp(1j*h*z)-sym.exp(-1j*h*z)
- sx = numpy.sin(f*x)
- sy = numpy.sin(g*y)
- cx = numpy.cos(f*x)
- cy = numpy.cos(g*y)
+ tdep = sym.exp(-1j * omega * sym.ScalarVariable("t"))
- zdep_add = numpy.exp(1j*h*z)+numpy.exp(-1j*h*z)
- zdep_sub = numpy.exp(1j*h*z)-numpy.exp(-1j*h*z)
+ const = 1j/(kx**2+ky**2)
- tdep = numpy.exp(-1j * omega * self.t)
+ result = sym.join_fields(
+ const*kx*h*cx*sy*zdep_sub*tdep, # ex
+ const*ky*h*sx*cy*zdep_sub*tdep, # ey
+ sx*sy*zdep_add*tdep, # ez
- C = 1j/(f**2+g**2)
+ -const*ky*epsilon*omega*sx*cy*zdep_add*tdep, # hx
+ const*kx*epsilon*omega*cx*sy*zdep_add*tdep,
+ 0,
+ )
- result = numpy.zeros((6,len(x)), dtype=zdep_add.dtype)
+ return result
- 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
+def get_rectangular_cavity_mode(E_0,
+ ambient_dim, mode_indices,
+ dimensions=(1, 1, 1), coefficients=(1, 0, 0),
+ forward_coeff=1, backward_coeff=0,
+ ):
+ """A TM cavity mode.
+
+ Returns an expression depending on *epsilon*, *mu*, and *t*.
+ """
+
+ kx, ky, kz = factors = [n*np.pi/a for n, a in zip(mode_indices, dimensions)]
+ k = np.sqrt(sum(f**2 for f in factors))
+
+ epsilon = sym.ScalarVariable("epsilon")
+ mu = sym.ScalarVariable("mu")
+
+ c = sym.cse(1/sym.sqrt(mu*epsilon), "c")
+ omega = k*c
+ gamma_squared = ky**2 + kx**2
+
+ nodes = sym.nodes(ambient_dim)
+ x = nodes[0]
+ y = nodes[1]
+ if nodes.shape[0] > 2:
+ z = nodes[2]
+ else:
+ z = 0
+
+ sx = sym.sin(kx*x)
+ sy = sym.sin(ky*y)
+ sz = sym.sin(kz*z)
+ cx = sym.cos(kx*x)
+ cy = sym.cos(ky*y)
+ cz = sym.cos(kz*z)
+
+ tdep = sym.exp(-1j * omega * sym.ScalarVariable("t"))
+ #tdep = 1
+
+ #result = sym.join_fields(
+ #E_0*cx*sy*sz*tdep, # ex
+ #E_0*sx*cy*cz*tdep * ky / kx, # ey
+ #-E_0 * (kx/kz + ky * ky / (kx*kz)) * sx*sy*cz*tdep, # ez
+
+ #1j*E_0*((kx + ky*ky/kx)*(ky/kz) + ky*kz/kx)*sx*cy*cz*tdep/omega, # hx
+ #-1j*E_0*(kz + (kx + ky*ky/kx))*cx*sy*cz*tdep/omega,
+ #0,
+ #)
+
+ #result = sym.join_fields(
+ #-kx * kz * E_0*cx*sy*sz*tdep / (k**2), # ex
+ #-ky * kz * E_0*sx*cy*sz*tdep / (k**2), # ey
+ #E_0 * sx*sy*cz*tdep, # ez
+
+ #1j * omega * ky*E_0*sx*cy*cz*tdep / (k**2), # hx
+ #-1j * omega * kx*E_0*cx*sy*cz*tdep / (k**2),
+ #0,
+ #)
+
+ 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
+
+def get_rectangular_3D_cavity_mode_deriv(E_0,
+ ambient_dim, mode_indices,
+ dimensions=(1, 1, 1), coefficients=(1, 0, 0),
+ forward_coeff=1, backward_coeff=0,
+ ):
+ """A TM cavity mode.
+
+ Returns an expression depending on *epsilon*, *mu*, and *t*.
+ """
+ kx, ky, kz = factors = [n*np.pi/a for n, a in zip(mode_indices, dimensions)]
+ k = np.sqrt(sum(f**2 for f in factors))
+
+ epsilon = sym.ScalarVariable("epsilon")
+ mu = sym.ScalarVariable("mu")
+
+ c = sym.cse(1/sym.sqrt(mu*epsilon), "c")
+ omega = k*c
+ gamma_squared = ky**2 + kx**2
+
+ nodes = sym.nodes(ambient_dim)
+ x = nodes[0]
+ y = nodes[1]
+ if nodes.shape[0] > 2:
+ z = nodes[2]
+ else:
+ z = 0
+
+ sx = sym.sin(kx*x)
+ sy = sym.sin(ky*y)
+ sz = sym.sin(kz*z)
+ cx = sym.cos(kx*x)
+ cy = sym.cos(ky*y)
+ cz = sym.cos(kz*z)
+
+ tdep = sym.exp(-1j * omega * sym.ScalarVariable("t"))
+ #tdep = 1
+
+ #result = sym.join_fields(
+ #E_0*cx*sy*sz*tdep, # ex
+ #E_0*sx*cy*cz*tdep * ky / kx, # ey
+ #-E_0 * (kx/kz + ky * ky / (kx*kz)) * sx*sy*cz*tdep, # ez
+
+ #1j*E_0*((kx + ky*ky/kx)*(ky/kz) + ky*kz/kx)*sx*cy*cz*tdep/omega, # hx
+ #-1j*E_0*(kz + (kx + ky*ky/kx))*cx*sy*cz*tdep/omega,
+ #0,
+ #)
+
+ #result = sym.join_fields(
+ #-kx * kz * E_0*cx*sy*sz*tdep / (k**2), # ex
+ #-ky * kz * E_0*sx*cy*sz*tdep / (k**2), # ey
+ #E_0 * sx*sy*cz*tdep, # ez
+
+ #1j * omega * ky*E_0*sx*cy*cz*tdep / (k**2), # hx
+ #-1j * omega * kx*E_0*cx*sy*cz*tdep / (k**2),
+ #0,
+ #)
+
+ result = sym.join_fields(
+ 1j*omega*kx * kz * E_0*cx*sy*sz*tdep / gamma_squared, # ex
+ 1j * omega * ky * kz * E_0*sx*cy*sz*tdep / gamma_squared, # ey
+ -1j * omega * E_0 * sx*sy*cz*tdep, # ez
+
+ -omega * omega * ky*E_0*sx*cy*cz*tdep / gamma_squared, # hx
+ omega * omega * kx*E_0*cx*sy*cz*tdep / gamma_squared,
+ 0,
+ )
+
+ return result
+
+
+def get_rectangular_2D_cavity_mode(m,n ):
+ """A TM cavity mode.
+
+ Returns an expression depending on *epsilon*, *mu*, and *t*.
+ """
+ omega = np.pi * np.sqrt(m*m+n*n)
+ nodes = sym.nodes(2)
+ x = nodes[0]
+ y = nodes[1]
-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)
+ xfac = m * np.pi * x
+ yfac = n * np.pi * y
+
+ tfac = sym.ScalarVariable("t") * omega
+ result = sym.join_fields(
+ 0,
+ 0,
+ sym.sin(xfac) * sym.sin(yfac) * sym.cos(tfac), # ez
+ -np.pi * n * sym.sin(xfac) * sym.cos(yfac) * sym.sin(tfac) / omega, # hx
+ np.pi * m * sym.cos(xfac) * sym.sin(yfac) * sym.sin(tfac) / omega, # hy
+ 0,
+ )
+ return result
+def get_rectangular_2D_cavity_mode_deriv(m,n ):
+ """A TM cavity mode.
+
+ Returns an expression depending on *epsilon*, *mu*, and *t*.
+ """
+
+ omega = np.pi * np.sqrt(m*m+n*n)
+
+
+ nodes = sym.nodes(2)
+ x = nodes[0]
+ y = nodes[1]
+
+ xfac = m * np.pi * x
+ yfac = n * np.pi * y
+
+ tfac = sym.ScalarVariable("t") * omega
+
+ result = sym.join_fields(
+ 0,
+ 0,
+ -omega * sym.sin(xfac) * sym.sin(yfac) * sym.sin(tfac), # ez
+ -np.pi * n * sym.sin(xfac) * sym.cos(yfac) * sym.cos(tfac) , # hx
+ np.pi * m * sym.cos(xfac) * sym.sin(yfac) * sym.cos(tfac) , # hy
+ 0,
+ )
+
+ return result
+
+# {{{ dipole solution
-# dipole solution -------------------------------------------------------------
class DipoleFarField:
"""Dipole EM field.
@@ -384,16 +537,14 @@ class DipoleFarField:
j0 = self.q*self.d*self.omega/self.epsilon
return j0*cos(self.omega*t)
+# }}}
+# {{{ check_time_harmonic_solution
-# 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)
@@ -432,15 +583,15 @@ def check_time_harmonic_solution(discr, mode, c_sol):
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),
+ ("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),
- ],
+ ("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
)
@@ -458,3 +609,6 @@ def check_time_harmonic_solution(discr, mode, c_sol):
rel_l2_error(hi_res, hi),
))
+# }}}
+
+# vim: foldmethod=marker
diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py
index 6fb7d54ac4eb69cf87705fce00268f68debfd4ba..8877fab8b72598111809e475d643d33b6b789f27 100644
--- a/examples/maxwell/cavities.py
+++ b/examples/maxwell/cavities.py
@@ -1,260 +1,142 @@
-# 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 .
-
-
-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
+"""Minimal example of a grudge driver."""
- logging.basicConfig(level=logging.DEBUG)
+from __future__ import division, print_function
- from grudge.backends import guess_run_context
- rcon = guess_run_context(allow_features)
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
- epsilon0 = 8.8541878176e-12 # C**2 / (N m**2)
- mu0 = 4*pi*1e-7 # N/A**2.
- epsilon = 1*epsilon0
- mu = 1*mu0
+__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:
- eoc_rec = EOCRecorder()
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
- cylindrical = False
- periodic = False
+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.
+"""
- 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)
+import numpy as np
+import pyopencl as cl
+from grudge.shortcuts import set_up_rk4, set_up_forward_euler # noqa: F401
+from grudge import sym, bind, Discretization
+
+from analytic_solutions import ( # noqa
+ get_rectangular_cavity_mode,
+ get_rectangular_2D_cavity_mode,
+ #RectangularCavityMode,
+ #RectangularWaveguideMode,
+ #RectangularCavityMode
+ )
+
+
+def main(write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dims = 3
+ 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)
+
+ discr = Discretization(cl_ctx, mesh, order=order)
+
+ 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:
- 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()
-
- if rcon.is_head_rank:
- print("---------------------------------------------")
- print("order %d" % order)
- print("---------------------------------------------")
- print("#elements=", len(mesh.elements))
-
- 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)
-
- # {{{ diagnostics setup
-
- from pytools.log import LogManager, add_general_quantities, \
- add_simulation_quantities, add_run_info
-
- if write_output:
- log_file_name = "maxwell-%d.dat" % order
- else:
- log_file_name = None
-
- 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)
+ epsilon = 1
+ mu = 1
- 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))
-
- 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()
-
- fields = stepper(fields, t, dt, rhs)
-
- mode.set_time(final_time)
-
- eoc_rec.add_data_point(order,
- discr.norm(fields-get_true_field()))
-
- finally:
- if write_output:
- vis.close()
-
- logmgr.close()
- discr.close()
+ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa
+ from grudge.models.em import MaxwellOperator
+ op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
- if rcon.is_head_rank:
- print()
- print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
+ queue = cl.CommandQueue(discr.cl_context)
- # }}}
+ if dims == 3:
+ sym_mode = get_rectangular_cavity_mode(1,mesh.ambient_dim, (1, 2, 2))
+ fields = bind(discr, sym_mode)(queue, t=0, epsilon = epsilon, mu=mu)
+ else:
+ sym_mode = get_rectangular_2D_cavity_mode(2,3)
+ fields = bind(discr, sym_mode)(queue, t=0)
- assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6
+ # 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())
+ print(bound_op)
+ # 1/0
+ 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()
diff --git a/examples/maxwell/testing.py b/examples/maxwell/testing.py
new file mode 100644
index 0000000000000000000000000000000000000000..09f960588d80fa75c8bae4de2cd8705482275ba0
--- /dev/null
+++ b/examples/maxwell/testing.py
@@ -0,0 +1,340 @@
+"""Minimal example of a grudge driver."""
+
+from __future__ import division, print_function
+
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+
+__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:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+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.
+"""
+
+
+import numpy as np
+import pyopencl as cl
+import sumpy.point_calculus as spy
+from grudge.shortcuts import set_up_rk4, set_up_forward_euler # noqa: F401
+from grudge import sym, bind, Discretization
+
+from analytic_solutions import ( # noqa
+ get_rectangular_cavity_mode,
+ get_rectangular_2D_cavity_mode,
+ get_rectangular_2D_cavity_mode_deriv,
+ get_rectangular_3D_cavity_mode_deriv,
+ #RectangularCavityMode,
+ #RectangularWaveguideMode,
+ #RectangularCavityMode
+ )
+
+
+def main(write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dims = 2
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(0.0,)*dims,
+ b=(1.0,)*dims,
+ n=(32,)*dims)
+
+ discr = Discretization(cl_ctx, mesh, order=order)
+
+ 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
+
+ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa
+ from grudge.models.em import TMMaxwellOperator, MaxwellOperator
+ op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
+
+ queue = cl.CommandQueue(discr.cl_context)
+
+ #sym_mode = get_rectangular_cavity_mode(1,mesh.ambient_dim, (1, 2, 2))
+ sym_mode = get_rectangular_2D_cavity_mode(1,1)
+ analytic_sol = bind(discr, sym_mode)
+ fields = analytic_sol(queue, t=0, epsilon=epsilon, mu=mu)
+
+ # FIXME
+ #dt = op.estimate_rk4_timestep(discr, fields=fields)
+
+ op.check_bc_coverage(mesh)
+
+ from grudge.tools import count_subset
+ bound_op = bind(discr, op.sym_operator())
+
+ def rhs(t, w):
+ return bound_op(queue, t=t, w=w)
+
+ if mesh.dim == 2:
+ dt = 0.002
+ elif mesh.dim == 3:
+ dt = 0.002
+
+ dt_stepper = set_up_rk4("w", dt, fields, rhs)
+
+ final_t = dt * 300
+ nsteps = int(final_t/dt)
+
+ print("dt=%g nsteps=%d" % (dt, nsteps))
+
+ from grudge.shortcuts import make_visualizer
+ vis = make_visualizer(discr, vis_order=order)
+
+ step = 0
+
+ norm = bind(discr, sym.norm(2, sym.var("u")))
+
+ from time import time
+ t_last_step = time()
+
+ e, h = op.split_eh(fields)
+
+ if 0:
+ vis.write_vtk_file("fld-%04d.vtu" % step,
+ [
+ ("e", e),
+ ("h", h),
+ ])
+
+ 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:
+ print(step)
+ #e, h = op.split_eh(event.state_component)
+ #vis.write_vtk_file("fld-%04d.vtu" % step,
+ #[
+ #("e", e),
+ #("h", h),
+ #])
+ t_last_step = time()
+ if step == 20:
+ sol = analytic_sol(queue, mu=mu, epsilon=epsilon, t= step * dt)
+ print(norm(queue, u=( event.state_component[0] - sol[0]))/ norm(queue, u=( sol[0])))
+ print(norm(queue, u=( event.state_component[1] - sol[1]))/ norm(queue, u=( sol[1])))
+ print(norm(queue, u=( event.state_component[2] - sol[2]))/ norm(queue, u=( sol[2])))
+ print(norm(queue, u=( event.state_component[3] - sol[3]))/ norm(queue, u=( sol[3])))
+ print(norm(queue, u=( event.state_component[4] - sol[4]))/ norm(queue, u=( sol[4])))
+ print(norm(queue, u=( event.state_component[5] - sol[5]))/ norm(queue, u=( sol[5])))
+ return
+
+def analytic_test(dims = 3, n = 8, write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ 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 = Discretization(cl_ctx, mesh, order=order)
+
+ 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
+
+ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa
+ from grudge.models.em import TMMaxwellOperator, MaxwellOperator
+ op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
+
+ queue = cl.CommandQueue(discr.cl_context)
+
+ if dims == 2:
+ sym_mode = get_rectangular_2D_cavity_mode(1,1)
+ else:
+ sym_mode = get_rectangular_cavity_mode(1,mesh.ambient_dim, (1, 2, 2))
+
+ analytic_sol = bind(discr, sym_mode)
+ fields = analytic_sol(queue, t=0, epsilon=epsilon, mu=mu)
+
+
+ op.check_bc_coverage(mesh)
+
+ from grudge.tools import count_subset
+ bound_op = bind(discr, op.sym_operator())
+
+ def rhs(t, w):
+ return bound_op(queue, t=t, w=w)
+
+ if mesh.dim == 2:
+ dt = 0.002
+ elif mesh.dim == 3:
+ dt = 0.002
+
+ dt_stepper = set_up_rk4("w", dt, fields, rhs)
+
+ final_t = dt * 300
+ nsteps = int(final_t/dt)
+
+ print("dt=%g nsteps=%d" % (dt, nsteps))
+
+ step = 0
+
+ norm = bind(discr, sym.norm(2, sym.var("u")))
+
+
+ e, h = op.split_eh(fields)
+
+ for event in dt_stepper.run(t_end=final_t):
+ if isinstance(event, dt_stepper.StateComputed):
+ assert event.component_id == "w"
+
+ step += 1
+
+ if step % 10 == 0:
+ print(step)
+
+ if step == 20:
+ sol = analytic_sol(queue, mu=mu, epsilon=epsilon, t= step * dt)
+ print(norm(queue, u=( event.state_component[0] - sol[0]))/ norm(queue, u=( sol[0])))
+ print(norm(queue, u=( event.state_component[1] - sol[1]))/ norm(queue, u=( sol[1])))
+ print(norm(queue, u=( event.state_component[2] - sol[2]))/ norm(queue, u=( sol[2])))
+ print(norm(queue, u=( event.state_component[3] - sol[3]))/ norm(queue, u=( sol[3])))
+ print(norm(queue, u=( event.state_component[4] - sol[4]))/ norm(queue, u=( sol[4])))
+ print(norm(queue, u=( event.state_component[5] - sol[5]))/ norm(queue, u=( sol[5])))
+ return
+
+def sumpy_test_3D(write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dims = 3
+
+ epsilon = 1
+ mu = 1
+
+ kx, ky, kz = factors = [n*np.pi/a for n, a in zip((1,2,2), (1,1,1))]
+ k = np.sqrt(sum(f**2 for f in factors))
+
+ c = 1/np.sqrt(mu*epsilon)
+ omega = k*c
+
+ patch = spy.CalculusPatch((0.4,0.3, 0.4))
+
+
+ queue = cl.CommandQueue(cl_ctx)
+
+ from grudge.discretization import PointsDiscretization
+ pdiscr = PointsDiscretization(cl.array.to_device(queue, patch.points))
+
+ sym_mode = get_rectangular_cavity_mode(1,3, (1, 2, 2))
+ fields = bind(pdiscr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu)
+
+ sym_mode_d = get_rectangular_3D_cavity_mode_deriv(1,3, (1, 2, 2))
+ fields_d = bind(pdiscr, sym_mode_d)(queue, t=0, epsilon=epsilon, mu=mu)
+
+ for i in range(len(fields)):
+ if type(fields[i]) == type(0):
+ fields[i] = np.zeros(patch.points.shape[1])
+ else:
+ fields[i] = fields[i].get()
+
+ # 2D
+ #print( patch.dy(fields[2].get()) + fields_d[3].get())
+ #print( patch.dx(fields[2].get()) - fields_d[4].get())
+ #print( fields_d[2].get() + patch.dy(fields[3].get()) - patch.dx(fields[4].get()))
+ #print( fields_d[2] + patch.dy(fields[3]) - patch.dx(fields[4]))
+
+ # 3D
+ #print( fields_d[3].get() + patch.dy(fields[2].get()) - patch.dz(fields[1].get()))
+ print( fields_d[3].get() + patch.dy(fields[2]) - patch.dz(fields[1]))
+ print( frequency_domain_maxwell(patch, np.array([fields[0], fields[1], fields[2]]), np.array([fields[3], fields[4], fields[5]]), k))
+ #print( frequency_domain_maxwell(patch, [fields[0].get(), fields[1].get(), fields[2].get()], [fields[3].get(), fields[4].get(), np.zeros(fields[4].get().size)], k))
+
+
+ #ans = bound_op(queue, w = fields)
+ #norm = bind(discr, sym.norm(2, sym.var("u")))
+ #print(norm(queue, u=(-1j * omega* fields[0] - bound_op(queue, w=fields)[0]))/norm(queue, u =(-1j * omega * fields[0])))
+ #print(norm(queue, u=(-1j * omega* fields[1] - bound_op(queue, w=fields)[1]))/norm(queue, u =(-1j * omega * fields[1])))
+ #print(norm(queue, u= (fields_d[2] - ans[2])))
+ #print(norm(queue, u=(fields_d[3] - ans[3])))
+ #print(norm(queue, u=(fields_d[4] - ans[4])))
+ #print(norm(queue, u=(-1j * omega* fields[5] - bound_op(queue, w=fields)[5]))/norm(queue, u =(-1j * omega * fields[5])))
+ return
+
+def sumpy_test_2D(write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dims = 2
+ epsilon = 1
+ mu = 1
+
+
+ patch = spy.CalculusPatch((0.4,0.3))
+
+
+ queue = cl.CommandQueue(cl_ctx)
+
+ from grudge.discretization import PointsDiscretization
+ pdiscr = PointsDiscretization(cl.array.to_device(queue, patch.points))
+
+ sym_mode = get_rectangular_2D_cavity_mode(1,1)
+ fields = bind(pdiscr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu)
+
+ sym_mode_d = get_rectangular_2D_cavity_mode_deriv(1,1)
+ fields_d = bind(pdiscr, sym_mode_d)(queue, t=0, epsilon=epsilon, mu=mu)
+
+
+ print( patch.dy(fields[2].get()) + fields_d[3].get())
+ print( patch.dx(fields[2].get()) - fields_d[4].get())
+ print( fields_d[2].get() + patch.dy(fields[3].get()) - patch.dx(fields[4].get()))
+
+ return
+
+def frequency_domain_maxwell(cpatch, e, h, k):
+ mu = 1
+ epsilon = 1
+ c = 1/np.sqrt(mu*epsilon)
+ omega = k*c
+ b = mu*h
+ d = epsilon*e
+ # https://en.wikipedia.org/w/index.php?title=Maxwell%27s_equations&oldid=798940325#Macroscopic_formulation
+ # assumed time dependence exp(-1j*omega*t)
+ # Agrees with Jackson, Third Ed., (8.16)
+ resid_faraday = np.vstack(cpatch.curl(e)) - 1j * omega * b
+ resid_ampere = np.vstack( cpatch.curl(h)) + 1j * omega * d
+ resid_div_e = cpatch.div(e)
+ resid_div_h = cpatch.div(h)
+ return ( resid_faraday, resid_ampere, resid_div_e, resid_div_h)
+
+if __name__ == "__main__":
+ analytic_test(3, 4, False,4)
+ analytic_test(3, 8, False,4)
+ #sumpy_test_3D(False,4)
+
+ #main()
+
diff --git a/examples/wave/var2.py b/examples/wave/var2.py
new file mode 100644
index 0000000000000000000000000000000000000000..e07acd44db99b45a131452c64a4f330e55088292
--- /dev/null
+++ b/examples/wave/var2.py
@@ -0,0 +1,130 @@
+"""Minimal example of a grudge driver."""
+
+from __future__ import division, print_function
+
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+
+__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:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+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.
+"""
+
+
+import numpy as np
+import pyopencl as cl
+from grudge.shortcuts import set_up_rk4
+from grudge import sym, bind, Discretization
+
+
+def main(write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dims = 2
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(-0.5,)*dims,
+ b=(0.5,)*dims,
+ n=(8,)*dims)
+
+ discr = Discretization(cl_ctx, mesh, order=order)
+
+ source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
+ source_width = 0.05
+ source_omega = 3
+
+ sym_x = sym.nodes(mesh.dim)
+ sym_source_center_dist = sym_x - source_center
+ sym_t = sym.ScalarVariable("t")
+ c = sym.If(sym.Comparison(
+ np.dot(sym_x, sym_x), "<", 0.2),
+ -0.1, -0.2)
+ #c = np.dot(sym_x,sym_x)
+
+ from grudge.models.wave import VariableCoefficientWeakWaveOperator
+ from meshmode.mesh import BTAG_ALL, BTAG_NONE
+ op = VariableCoefficientWeakWaveOperator(c,
+ discr.dim,
+ source_f=(
+ sym.sin(source_omega*sym_t)
+ * sym.exp(
+ -np.dot(sym_source_center_dist, sym_source_center_dist)
+ / source_width**2)),
+ dirichlet_tag=BTAG_NONE,
+ neumann_tag=BTAG_NONE,
+ radiation_tag=BTAG_ALL,
+ flux_type="upwind")
+
+ queue = cl.CommandQueue(discr.cl_context)
+ from pytools.obj_array import join_fields
+ fields = join_fields(discr.zeros(queue),
+ [discr.zeros(queue) for i in range(discr.dim)])
+
+ # FIXME
+ #dt = op.estimate_rk4_timestep(discr, fields=fields)
+
+ op.check_bc_coverage(mesh)
+
+ # print(sym.pretty(op.sym_operator()))
+ bound_op = bind(discr, op.sym_operator())
+ print(bound_op)
+ # 1/0
+
+ def rhs(t, w):
+ return bound_op(queue, t=t, w=w)
+
+ if mesh.dim == 2:
+ dt = 0.04
+ elif mesh.dim == 3:
+ dt = 0.02
+
+ dt_stepper = set_up_rk4("w", dt, fields, rhs)
+
+ final_t = 10
+ nsteps = int(final_t/dt)
+ print("dt=%g nsteps=%d" % (dt, nsteps))
+
+ from grudge.shortcuts import make_visualizer
+ vis = make_visualizer(discr, vis_order=order)
+
+ step = 0
+
+ norm = bind(discr, sym.norm(2, sym.var("u")))
+
+ from time import time
+ t_last_step = time()
+
+ 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=event.state_component[0]),
+ time()-t_last_step)
+ if step % 10 == 0:
+ vis.write_vtk_file("fld-%04d.vtu" % step,
+ [
+ ("u", event.state_component[0]),
+ ("v", event.state_component[1:]),
+ ])
+ t_last_step = time()
+
+
+if __name__ == "__main__":
+ main()
diff --git a/examples/wave/weak.py b/examples/wave/weak.py
new file mode 100644
index 0000000000000000000000000000000000000000..9df822abdf16425e5f3eb25d481b3c0117cd9087
--- /dev/null
+++ b/examples/wave/weak.py
@@ -0,0 +1,125 @@
+"""Minimal example of a grudge driver."""
+
+from __future__ import division, print_function
+
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+
+__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:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+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.
+"""
+
+
+import numpy as np
+import pyopencl as cl
+from grudge.shortcuts import set_up_rk4
+from grudge import sym, bind, Discretization
+
+
+def main(write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dims = 2
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(-0.5,)*dims,
+ b=(0.5,)*dims,
+ n=(8,)*dims)
+
+ discr = Discretization(cl_ctx, mesh, order=order)
+
+ source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
+ source_width = 0.05
+ source_omega = 3
+
+ sym_x = sym.nodes(mesh.dim)
+ sym_source_center_dist = sym_x - source_center
+ sym_t = sym.ScalarVariable("t")
+
+ from grudge.models.wave import WeakWaveOperator
+ from meshmode.mesh import BTAG_ALL, BTAG_NONE
+ op = WeakWaveOperator(-0.1, discr.dim,
+ source_f=(
+ sym.sin(source_omega*sym_t)
+ * sym.exp(
+ -np.dot(sym_source_center_dist, sym_source_center_dist)
+ / source_width**2)),
+ dirichlet_tag=BTAG_NONE,
+ neumann_tag=BTAG_NONE,
+ radiation_tag=BTAG_ALL,
+ flux_type="upwind")
+
+ queue = cl.CommandQueue(discr.cl_context)
+ from pytools.obj_array import join_fields
+ fields = join_fields(discr.zeros(queue),
+ [discr.zeros(queue) for i in range(discr.dim)])
+
+ # FIXME
+ #dt = op.estimate_rk4_timestep(discr, fields=fields)
+
+ op.check_bc_coverage(mesh)
+
+ # print(sym.pretty(op.sym_operator()))
+ bound_op = bind(discr, op.sym_operator())
+ print(bound_op)
+ # 1/0
+
+ def rhs(t, w):
+ return bound_op(queue, t=t, w=w)
+
+ if mesh.dim == 2:
+ dt = 0.04
+ elif mesh.dim == 3:
+ dt = 0.02
+
+ dt_stepper = set_up_rk4("w", dt, fields, rhs)
+
+ final_t = 10
+ nsteps = int(final_t/dt)
+ print("dt=%g nsteps=%d" % (dt, nsteps))
+
+ from grudge.shortcuts import make_visualizer
+ vis = make_visualizer(discr, vis_order=order)
+
+ step = 0
+
+ norm = bind(discr, sym.norm(2, sym.var("u")))
+
+ from time import time
+ t_last_step = time()
+
+ 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=event.state_component[0]),
+ time()-t_last_step)
+ if step % 10 == 0:
+ vis.write_vtk_file("fld-%04d.vtu" % step,
+ [
+ ("u", event.state_component[0]),
+ ("v", event.state_component[1:]),
+ ])
+ t_last_step = time()
+
+
+if __name__ == "__main__":
+ main()
diff --git a/grudge/discretization.py b/grudge/discretization.py
index 8cbacd8981f4c353b35b3f325b465b0030b5be97..d578e636df68b12205561bb0be424ef228b5372d 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -25,9 +25,14 @@ THE SOFTWARE.
from pytools import memoize_method
from grudge import sym
+import numpy as np
-class Discretization(object):
+class DiscretizationBase(object):
+ pass
+
+
+class Discretization(DiscretizationBase):
"""
.. attribute :: volume_discr
@@ -220,4 +225,36 @@ class Discretization(object):
or where == sym.VTAG_ALL)
+class PointsDiscretization(DiscretizationBase):
+ """Implements just enough of the discretization interface to be
+ able to smuggle some points into :func:`bind`.
+ """
+
+ 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]
+
+ @property
+ def mesh(self):
+ return self
+
+ @property
+ def nnodes(self):
+ return self._nodes.shape[-1]
+
+ def nodes(self):
+ return self._nodes
+
+ @property
+ def volume_discr(self):
+ return self
+
+
# vim: foldmethod=marker
diff --git a/grudge/execution.py b/grudge/execution.py
index ace2dc8bd1403035e217b5936cf3eca302f358c6..c4f285cd6c49a7e6a9fb9eacabdb9a030b453607 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -116,14 +116,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)
+
+ 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 0:
- self.sign = 1
- else:
- self.sign = -1
-
- self.dirichlet_tag = dirichlet_tag
- self.neumann_tag = neumann_tag
- self.radiation_tag = radiation_tag
-
- self.dirichlet_bc_f = dirichlet_bc_f
-
- self.flux_type = flux_type
-
- def flux(self, w):
- u = w[0]
- v = w[1:]
- normal = sym.normal(w.dd, self.ambient_dim)
-
- flux_weak = join_fields(
- np.dot(v.avg, normal),
- u.avg * normal)
-
- if self.flux_type == "central":
- pass
- elif self.flux_type == "upwind":
- # see doc/notes/grudge-notes.tm
- flux_weak -= self.sign*join_fields(
- 0.5*(u.int-u.ext),
- 0.5*(normal * np.dot(normal, v.int-v.ext)))
- else:
- raise ValueError("invalid flux type '%s'" % self.flux_type)
-
- flux_strong = join_fields(
- np.dot(v.int, normal),
- u.int * normal) - flux_weak
-
- return self.c*flux_strong
-
- def sym_operator(self):
- d = self.ambient_dim
-
- w = sym.make_sym_array("w", d+1)
- u = w[0]
- v = w[1:]
-
- # boundary conditions -------------------------------------------------
-
- # dirichlet BCs -------------------------------------------------------
- dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
- dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
- if self.dirichlet_bc_f:
- # FIXME
- from warnings import warn
- warn("Inhomogeneous Dirichlet conditions on the wave equation "
- "are still having issues.")
-
- dir_g = sym.Field("dir_bc_u")
- dir_bc = join_fields(2*dir_g - dir_u, dir_v)
- else:
- dir_bc = join_fields(-dir_u, dir_v)
-
- dir_bc = sym.cse(dir_bc, "dir_bc")
-
- # neumann BCs ---------------------------------------------------------
- neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
- neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
- neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc")
-
- # radiation BCs -------------------------------------------------------
- rad_normal = sym.normal(self.radiation_tag, d)
-
- rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
- rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
-
- rad_bc = sym.cse(join_fields(
- 0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
- 0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.sign*rad_u)
- ), "rad_bc")
-
- # entire operator -----------------------------------------------------
- nabla = sym.nabla(d)
-
- def flux(pair):
- return sym.interp(pair.dd, "all_faces")(
- self.flux(pair))
-
- result = (
- - join_fields(
- -self.c*np.dot(nabla, v),
- -self.c*(nabla*u)
- )
- +
- sym.InverseMassOperator()(
- sym.FaceMassOperator()(
- flux(sym.int_tpair(w))
- + flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
- + flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
- + flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
- )))
-
- result[0] += self.source_f
-
- return result
-
- def check_bc_coverage(self, mesh):
- from meshmode.mesh import check_bc_coverage
- check_bc_coverage(mesh, [
- self.dirichlet_tag,
- self.neumann_tag,
- self.radiation_tag])
-
- def max_eigenvalue(self, t, fields=None, discr=None):
- return abs(self.c)
-
-
-class WeakWaveOperator(HyperbolicOperator):
- """This operator discretizes the wave equation
- :math:`\\partial_t^2 u = c^2 \\Delta u`.
-
- To be precise, we discretize the hyperbolic system
-
- .. math::
-
- \partial_t u - c \\nabla \\cdot v = 0
-
- \partial_t v - c \\nabla u = 0
-
- The sign of :math:`v` determines whether we discretize the forward or the
- backward wave equation.
-
- :math:`c` is assumed to be constant across all space.
- """
-
- def __init__(self, c, ambient_dim, source_f=0,
- flux_type="upwind",
- dirichlet_tag=BTAG_ALL,
- dirichlet_bc_f=0,
- neumann_tag=BTAG_NONE,
- radiation_tag=BTAG_NONE):
- assert isinstance(ambient_dim, int)
-
- self.c = c
- self.ambient_dim = ambient_dim
- self.source_f = source_f
-
- if self.c > 0:
- self.sign = 1
- else:
- self.sign = -1
-
- self.dirichlet_tag = dirichlet_tag
- self.neumann_tag = neumann_tag
- self.radiation_tag = radiation_tag
-
- self.dirichlet_bc_f = dirichlet_bc_f
-
- self.flux_type = flux_type
-
- def flux(self, w):
- u = w[0]
- v = w[1:]
- normal = sym.normal(w.dd, self.ambient_dim)
-
- flux_weak = join_fields(
- np.dot(v.avg, normal),
- u.avg * normal)
-
- if self.flux_type == "central":
- return -self.c*flux_weak
- elif self.flux_type == "upwind":
- return -self.c*(flux_weak + self.sign*join_fields(
- 0.5*(u.int-u.ext),
- 0.5*(normal * np.dot(normal, v.int-v.ext))))
- else:
- raise ValueError("invalid flux type '%s'" % self.flux_type)
-
- def sym_operator(self):
- d = self.ambient_dim
-
- w = sym.make_sym_array("w", d+1)
- u = w[0]
- v = w[1:]
-
- # boundary conditions -------------------------------------------------
-
- # dirichlet BCs -------------------------------------------------------
- dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
- dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
- if self.dirichlet_bc_f:
- # FIXME
- from warnings import warn
- warn("Inhomogeneous Dirichlet conditions on the wave equation "
- "are still having issues.")
-
- dir_g = sym.Field("dir_bc_u")
- dir_bc = join_fields(2*dir_g - dir_u, dir_v)
- else:
- dir_bc = join_fields(-dir_u, dir_v)
-
- dir_bc = sym.cse(dir_bc, "dir_bc")
-
- # neumann BCs ---------------------------------------------------------
- neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
- neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
- neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc")
-
- # radiation BCs -------------------------------------------------------
- rad_normal = sym.normal(self.radiation_tag, d)
-
- rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
- rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
-
- rad_bc = sym.cse(join_fields(
- 0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
- 0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.sign*rad_u)
- ), "rad_bc")
-
- # entire operator -----------------------------------------------------
- def flux(pair):
- return sym.interp(pair.dd, "all_faces")(self.flux(pair))
-
- result = sym.InverseMassOperator()(
- join_fields(
- -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v),
- -self.c*(sym.stiffness_t(self.ambient_dim)*u)
- )
-
-
- - sym.FaceMassOperator()(flux(sym.int_tpair(w))
- + flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
- + flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
- + flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
-
- ))
-
- result[0] += self.source_f
-
- return result
-
- def check_bc_coverage(self, mesh):
- from meshmode.mesh import check_bc_coverage
- check_bc_coverage(mesh, [
- self.dirichlet_tag,
- self.neumann_tag,
- self.radiation_tag])
-
- def max_eigenvalue(self, t, fields=None, discr=None):
- return abs(self.c)
-
-
-# }}}
-
-
-# {{{ variable-velocity
-
-class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
- """This operator discretizes the wave equation
- :math:`\\partial_t^2 u = c^2 \\Delta u`.
-
- To be precise, we discretize the hyperbolic system
-
- .. math::
-
- \partial_t u - c \\nabla \\cdot v = 0
-
- \partial_t v - c \\nabla u = 0
-
- The sign of :math:`v` determines whether we discretize the forward or the
- backward wave equation.
-
- :math:`c` is assumed to be constant across all space.
- """
-
- def __init__(self, c, ambient_dim, source_f=0,
- flux_type="upwind",
- dirichlet_tag=BTAG_ALL,
- dirichlet_bc_f=0,
- neumann_tag=BTAG_NONE,
- radiation_tag=BTAG_NONE):
- assert isinstance(ambient_dim, int)
-
- self.c = c
- self.ambient_dim = ambient_dim
- self.source_f = source_f
-
- self.sign = sym.If(sym.Comparison(
- self.c, ">", 0),
- np.int32(1), np.int32(-1))
-
- self.dirichlet_tag = dirichlet_tag
- self.neumann_tag = neumann_tag
- self.radiation_tag = radiation_tag
-
- self.dirichlet_bc_f = dirichlet_bc_f
-
- self.flux_type = flux_type
-
- def flux(self, w):
- c = w[0]
- u = w[1]
- v = w[2:]
- normal = sym.normal(w.dd, self.ambient_dim)
-
- if self.flux_type == "central":
- return -0.5 * join_fields(
- np.dot(v.int*c.int + v.ext*c.ext, normal),
- (u.int * c.int + u.ext*c.ext) * normal)
-
- elif self.flux_type == "upwind":
- return -0.5 * join_fields(np.dot(normal, c.ext * v.ext + c.int * v.int)
- + c.ext*u.ext - c.int * u.int,
- normal * (np.dot(normal, c.ext * v.ext - c.int * v.int)
- + c.ext*u.ext + c.int * u.int))
-
- else:
- raise ValueError("invalid flux type '%s'" % self.flux_type)
-
- def sym_operator(self):
- d = self.ambient_dim
-
- w = sym.make_sym_array("w", d+1)
- u = w[0]
- v = w[1:]
- flux_w = join_fields(self.c, w)
-
- # boundary conditions -------------------------------------------------
-
- # dirichlet BCs -------------------------------------------------------
- dir_c = sym.cse(sym.interp("vol", self.dirichlet_tag)(self.c))
- dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
- dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
- if self.dirichlet_bc_f:
- # FIXME
- from warnings import warn
- warn("Inhomogeneous Dirichlet conditions on the wave equation "
- "are still having issues.")
-
- dir_g = sym.Field("dir_bc_u")
- dir_bc = join_fields(dir_c, 2*dir_g - dir_u, dir_v)
- else:
- dir_bc = join_fields(dir_c, -dir_u, dir_v)
-
- dir_bc = sym.cse(dir_bc, "dir_bc")
-
- # neumann BCs ---------------------------------------------------------
- neu_c = sym.cse(sym.interp("vol", self.neumann_tag)(self.c))
- neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
- neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
- neu_bc = sym.cse(join_fields(neu_c, neu_u, -neu_v), "neu_bc")
-
- # radiation BCs -------------------------------------------------------
- rad_normal = sym.normal(self.radiation_tag, d)
-
- rad_c = sym.cse(sym.interp("vol", self.radiation_tag)(self.c))
- rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
- rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
-
- rad_bc = sym.cse(join_fields(rad_c,
- 0.5*(rad_u - sym.interp("vol", self.radiation_tag)(self.sign)
- * np.dot(rad_normal, rad_v)),
- 0.5*rad_normal*(np.dot(rad_normal, rad_v)
- - sym.interp("vol", self.radiation_tag)(self.sign)*rad_u)
- ), "rad_bc")
-
- # entire operator -----------------------------------------------------
- def flux(pair):
- return sym.interp(pair.dd, "all_faces")(self.flux(pair))
-
- result = sym.InverseMassOperator()(
- join_fields(
- -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v),
- -self.c*(sym.stiffness_t(self.ambient_dim)*u)
- )
-
-
- - sym.FaceMassOperator()(flux(sym.int_tpair(flux_w))
- + flux(sym.bv_tpair(self.dirichlet_tag, flux_w, dir_bc))
- + flux(sym.bv_tpair(self.neumann_tag, flux_w, neu_bc))
- + flux(sym.bv_tpair(self.radiation_tag, flux_w, rad_bc))
-
- ))
-
- result[0] += self.source_f
-
- return result
-
- def check_bc_coverage(self, mesh):
- from meshmode.mesh import check_bc_coverage
- check_bc_coverage(mesh, [
- self.dirichlet_tag,
- self.neumann_tag,
- self.radiation_tag])
-
- def max_eigenvalue(self, t, fields=None, discr=None):
- return sym.NodalMax()(sym.CFunction("fabs")(self.c))
-
-# }}}
-
-
-# vim: foldmethod=marker
diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py
index 0e6940ff6f120018a8b5c3c9af29de04eb1bddcb..84fd8b3fd7a01217e4e30f1e1df525a912581c05 100644
--- a/grudge/shortcuts.py
+++ b/grudge/shortcuts.py
@@ -44,6 +44,22 @@ def set_up_rk4(field_var_name, dt, fields, rhs, t_start=0):
return dt_stepper
+def set_up_forward_euler(field_var_name, dt, fields, rhs, t_start=0):
+ from leap.rk import ForwardEulerMethod
+ from dagrt.codegen import PythonCodeGenerator
+
+ dt_method = ForwardEulerMethod(component_id=field_var_name)
+ dt_code = dt_method.generate()
+ dt_stepper_class = PythonCodeGenerator("TimeStep").get_class(dt_code)
+ dt_stepper = dt_stepper_class({""+dt_method.component_id: rhs})
+
+ dt_stepper.set_up(
+ t_start=t_start, dt_start=dt,
+ context={dt_method.component_id: fields})
+
+ return dt_stepper
+
+
def make_visualizer(discr, vis_order):
from meshmode.discretization.visualization import make_visualizer
with cl.CommandQueue(discr.cl_context) as queue:
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 60b489cecb50d6f69c562ecd7a14e4a0e786ad07..9eb8e2b3bcb18c822335a8a6a884daf903ae3efc 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
@@ -812,6 +813,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/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 0ceaf8fec8339734e93dbdbfd12bdd1e603f171f..bd47f823c36931290c1b7ba9981c5c8b41a7a032 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -333,6 +333,22 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
assert eoc_rec.order_estimate() > order
+def test_foreign_points(ctx_factory):
+ pytest.importorskip("sumpy")
+ import sumpy.point_calculus as pc
+
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ dim = 2
+ cp = pc.CalculusPatch(np.zeros(dim))
+
+ from grudge.discretization import PointsDiscretization
+ pdiscr = PointsDiscretization(cl.array.to_device(queue, cp.points))
+
+ bind(pdiscr, sym.nodes(dim)**2)(queue)
+
+
# You can test individual routines by typing
# $ python test_grudge.py 'test_routine()'