From 0f9ff81856bb3acd4c22b67afd90757e65e71a93 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Sat, 15 May 2021 16:30:43 -0500 Subject: [PATCH] Convert maxwell-related examples and tests --- examples/maxwell/cavities.py | 51 +++++++++++++++++++++--------------- test/test_grudge.py | 32 ++++++++++++---------- 2 files changed, 48 insertions(+), 35 deletions(-) diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index d82a5a5b..03f65362 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -27,17 +27,20 @@ import numpy as np import pyopencl as cl from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, DiscretizationCollection +from grudge import DiscretizationCollection from grudge.models.em import get_rectangular_cavity_mode +import grudge.op as op + STEPS = 60 -def main(dims, write_output=True, order=4): +def main(dims, write_output=False, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) @@ -48,7 +51,7 @@ def main(dims, write_output=True, order=4): b=(1.0,)*dims, nelements_per_axis=(4,)*dims) - discr = DiscretizationCollection(actx, mesh, order=order) + dcoll = DiscretizationCollection(actx, mesh, order=order) if 0: epsilon0 = 8.8541878176e-12 # C**2 / (N m**2) @@ -60,25 +63,30 @@ def main(dims, write_output=True, order=4): mu = 1 from grudge.models.em import MaxwellOperator - op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims) - if dims == 3: - sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2)) - fields = bind(discr, sym_mode)(actx, t=0, epsilon=epsilon, mu=mu) - else: - sym_mode = get_rectangular_cavity_mode(1, (2, 3)) - fields = bind(discr, sym_mode)(actx, t=0) + maxwell_operator = MaxwellOperator( + dcoll, + epsilon, + mu, + flux_type=0.5, + dimensions=dims + ) - # FIXME - #dt = op.estimate_rk4_timestep(discr, fields=fields) + def cavity_mode(x, t=0): + if dims == 3: + return get_rectangular_cavity_mode(actx, x, t, 1, (1, 2, 2)) + else: + return get_rectangular_cavity_mode(actx, x, t, 1, (2, 3)) - op.check_bc_coverage(mesh) + fields = cavity_mode(thaw(actx, op.nodes(dcoll)), t=0) + + # FIXME + # dt = maxwell_operator.estimate_rk4_timestep(dcoll, fields=fields) - # print(sym.pretty(op.sym_operator())) - bound_op = bind(discr, op.sym_operator()) + maxwell_operator.check_bc_coverage(mesh) def rhs(t, w): - return bound_op(t=t, w=w) + return maxwell_operator.operator(t, w) if mesh.dim == 2: dt = 0.004 @@ -93,18 +101,19 @@ def main(dims, write_output=True, order=4): print("dt=%g nsteps=%d" % (dt, nsteps)) from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr) + vis = make_visualizer(dcoll) step = 0 - norm = bind(discr, sym.norm(2, sym.var("u"))) + def norm(u): + return op.norm(dcoll, u, 2) from time import time t_last_step = time() - e, h = op.split_eh(fields) + e, h = maxwell_operator.split_eh(fields) - if 1: + if write_output: vis.write_vtk_file("fld-cavities-%04d.vtu" % step, [ ("e", e), @@ -121,7 +130,7 @@ def main(dims, write_output=True, order=4): norm(u=h[0]), norm(u=h[1]), time()-t_last_step) if step % 10 == 0: - e, h = op.split_eh(event.state_component) + e, h = maxwell_operator.split_eh(event.state_component) vis.write_vtk_file("fld-cavities-%04d.vtu" % step, [ ("e", e), diff --git a/test/test_grudge.py b/test/test_grudge.py index 01359e25..a227e97f 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -814,7 +814,6 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ @pytest.mark.parametrize("order", [3, 4, 5]) def test_convergence_maxwell(actx_factory, order): """Test whether 3D Maxwell's actually converges""" - from grudge import sym, bind actx = actx_factory() @@ -829,24 +828,32 @@ def test_convergence_maxwell(actx_factory, order): b=(1.0,)*dims, nelements_per_axis=(n,)*dims) - discr = DiscretizationCollection(actx, mesh, order=order) + dcoll = DiscretizationCollection(actx, 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(actx, t=0, epsilon=epsilon, mu=mu) + def analytic_sol(x, t=0): + return get_rectangular_cavity_mode(actx, x, t, 1, (1, 2, 2)) + + nodes = thaw(actx, op.nodes(dcoll)) + fields = analytic_sol(nodes, t=0) 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()) + + maxwell_operator = MaxwellOperator( + dcoll, + epsilon, + mu, + flux_type=0.5, + dimensions=dims + ) + maxwell_operator.check_bc_coverage(mesh) def rhs(t, w): - return bound_op(t=t, w=w) + return maxwell_operator.operator(t, w) dt = 0.002 final_t = dt * 5 @@ -857,8 +864,6 @@ def test_convergence_maxwell(actx_factory, order): logger.info("dt %.5e nsteps %5d", 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): @@ -868,9 +873,8 @@ def test_convergence_maxwell(actx_factory, order): step += 1 logger.debug("[%04d] t = %.5e", step, event.t) - sol = analytic_sol(actx, mu=mu, epsilon=epsilon, t=step * dt) - vals = [norm(u=(esc[i] - sol[i])) / norm(u=sol[i]) for i in range(5)] # noqa E501 - total_error = sum(vals) + sol = analytic_sol(nodes, t=step * dt) + total_error = op.norm(dcoll, esc - sol, 2) eoc_rec.add_data_point(1.0/n, total_error) logger.info("\n%s", eoc_rec.pretty_print( -- GitLab