From 014cd9bce5a30968dd11d79a41451a188fb63114 Mon Sep 17 00:00:00 2001
From: Bogdan Enache <enache2@illinois.edu>
Date: Thu, 9 Nov 2017 13:03:27 -0600
Subject: [PATCH] Cleaned Maxwell-related code, and added a convergence pytest
---
...wellRectCav3D.mac => maxwellrectcav3d.mac} | 2 +
examples/maxwell/cavities.py | 9 +-
examples/maxwell/testing.py | 86 +------------------
grudge/execution.py | 18 ----
grudge/models/em.py | 58 +++++++++++++
test/test_grudge.py | 72 ++++++++++++++++
6 files changed, 139 insertions(+), 106 deletions(-)
rename contrib/maxima/{maxwellRectCav3D.mac => maxwellrectcav3d.mac} (99%)
diff --git a/contrib/maxima/maxwellRectCav3D.mac b/contrib/maxima/maxwellrectcav3d.mac
similarity index 99%
rename from contrib/maxima/maxwellRectCav3D.mac
rename to contrib/maxima/maxwellrectcav3d.mac
index 7a51c157..abda9687 100644
--- a/contrib/maxima/maxwellRectCav3D.mac
+++ b/contrib/maxima/maxwellrectcav3d.mac
@@ -6,6 +6,8 @@ load("diag");
/* -------------------------------------------------------------------------- */
/* Utilities */
/* -------------------------------------------------------------------------- */
+mu: 1;
+epsilon: 1;
coords:[x,y,z];
diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py
index fcea20c6..85a7c3de 100644
--- a/examples/maxwell/cavities.py
+++ b/examples/maxwell/cavities.py
@@ -30,10 +30,7 @@ import pyopencl as cl
from grudge.shortcuts import set_up_rk4
from grudge import sym, bind, Discretization
-from analytic_solutions import (
- get_rectangular_3D_cavity_mode,
- get_rectangular_2D_cavity_mode,
- )
+from grudge.models.em import get_rectangular_cavity_mode
def main(dims, write_output=True, order=4):
@@ -63,10 +60,10 @@ def main(dims, write_output=True, order=4):
queue = cl.CommandQueue(discr.cl_context)
if dims == 3:
- sym_mode = get_rectangular_3D_cavity_mode(1, (1, 2, 2))
+ 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_2D_cavity_mode(1, (2, 3))
+ sym_mode = get_rectangular_cavity_mode(1, (2, 3))
fields = bind(discr, sym_mode)(queue, t=0)
# FIXME
diff --git a/examples/maxwell/testing.py b/examples/maxwell/testing.py
index b9218443..2ff45a68 100644
--- a/examples/maxwell/testing.py
+++ b/examples/maxwell/testing.py
@@ -28,80 +28,7 @@ THE SOFTWARE.
import numpy as np
import pyopencl as cl
import sumpy.point_calculus as spy
-from grudge.shortcuts import set_up_rk4
-from grudge import sym, bind, Discretization
-
-from analytic_solutions import (
- get_rectangular_3D_cavity_mode,
- get_rectangular_2D_cavity_mode,
- )
-
-
-def analytic_test(dims, n, order, cl_ctx, queue):
-
- 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
-
- if dims == 2:
- sym_mode = get_rectangular_2D_cavity_mode(1, (1, 2))
- else:
- sym_mode = get_rectangular_3D_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 * 300
- nsteps = int(final_t/dt)
-
- 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
-
- if step % 10 == 0:
- print(step)
-
- if step == 20:
- sol = analytic_sol(queue, mu=mu, epsilon=epsilon, t=step * dt)
- print(norm(queue, u=(esc[0] - sol[0])) / norm(queue, u=sol[0]))
- print(norm(queue, u=(esc[1] - sol[1])) / norm(queue, u=sol[1]))
- print(norm(queue, u=(esc[2] - sol[2])) / norm(queue, u=sol[2]))
- print(norm(queue, u=(esc[3] - sol[3])) / norm(queue, u=sol[3]))
- print(norm(queue, u=(esc[4] - sol[4])) / norm(queue, u=sol[4]))
- print(norm(queue, u=(esc[5] - sol[5])) / norm(queue, u=sol[5]))
- return
+from grudge import bind
def sumpy_test_3D(order, cl_ctx, queue):
@@ -109,7 +36,7 @@ def sumpy_test_3D(order, cl_ctx, queue):
epsilon = 1
mu = 1
- kx, ky, kz = factors = [n*np.pi/a for n, a in zip((1, 2, 2), (1, 1, 1))]
+ kx, ky, kz = factors = [n*np.pi for n in (1, 2, 2)]
k = np.sqrt(sum(f**2 for f in factors))
patch = spy.CalculusPatch((0.4, 0.3, 0.4))
@@ -117,7 +44,8 @@ def sumpy_test_3D(order, cl_ctx, queue):
from grudge.discretization import PointsDiscretization
pdiscr = PointsDiscretization(cl.array.to_device(queue, patch.points))
- sym_mode = get_rectangular_3D_cavity_mode(1, (1, 2, 2))
+ from grudge.models.em import get_rectangular_cavity_mode
+ sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2))
fields = bind(pdiscr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu)
for i in range(len(fields)):
@@ -155,9 +83,3 @@ if __name__ == "__main__":
print("Doing Sumpy Test")
sumpy_test_3D(4, cl_ctx, queue)
- print("Doing 2D")
- analytic_test(2, 8, 4, cl_ctx, queue)
- analytic_test(2, 16, 4, cl_ctx, queue)
- print("Doing 3D")
- analytic_test(3, 4, 4, cl_ctx, queue)
- analytic_test(3, 8, 4, cl_ctx, queue)
diff --git a/grudge/execution.py b/grudge/execution.py
index c4f285cd..28735582 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -382,24 +382,6 @@ class ExecutionMapper(mappers.Evaluator,
# }}}
# {{{ code execution functions
- # DEBUG stuff
- def get_max(self, arr):
- from pymbolic.primitives import Variable
- i = Variable("i")
-
- def knl2():
- knl = lp.make_kernel(
- "{[i]: 0<=i<n}",
- [
- lp.Assignment(Variable("out")[i],
- Variable("abs")(Variable("a")[i]))
- ])
- return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
-
- evt, (out,) = knl2()(self.queue, a=arr)
-
- return pyopencl.array.max(out)
-
def map_insn_loopy_kernel(self, insn):
kwargs = {}
kdescr = insn.kernel_descriptor
diff --git a/grudge/models/em.py b/grudge/models/em.py
index f50f393a..f0e44f90 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -435,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/test/test_grudge.py b/test/test_grudge.py
index bd47f823..609de951 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -333,6 +333,78 @@ 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 = Discretization(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
--
GitLab