From ed92043a76ea164b57beb87457aeb6709d03e851 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Fri, 3 Nov 2017 10:01:18 -0500 Subject: [PATCH 1/9] Attempt at solving maxwells analytically --- contrib/maxima/maxwellRectCav3D.mac | 71 +++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 contrib/maxima/maxwellRectCav3D.mac diff --git a/contrib/maxima/maxwellRectCav3D.mac b/contrib/maxima/maxwellRectCav3D.mac new file mode 100644 index 00000000..7a51c157 --- /dev/null +++ b/contrib/maxima/maxwellRectCav3D.mac @@ -0,0 +1,71 @@ +kill(all); +load("eigen"); +load("itensor"); +load("diag"); + +/* -------------------------------------------------------------------------- */ +/* Utilities */ +/* -------------------------------------------------------------------------- */ + +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)); + -- GitLab From 7ea883222231075a2369f50a1d709c22878820a8 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Fri, 3 Nov 2017 11:57:04 -0500 Subject: [PATCH 2/9] Implemented Maxwell's --- examples/maxwell/analytic_solutions.py | 75 ++++-- examples/maxwell/cavities.py | 322 ++++++++----------------- examples/maxwell/testing.py | 161 +++++++++++++ grudge/discretization.py | 6 + grudge/execution.py | 91 +++++-- grudge/models/em.py | 226 +++++------------ grudge/symbolic/compiler.py | 6 +- grudge/symbolic/mappers/__init__.py | 7 + grudge/tools.py | 88 +++++++ 9 files changed, 557 insertions(+), 425 deletions(-) create mode 100644 examples/maxwell/testing.py diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py index 7ff8884f..c7617d91 100644 --- a/examples/maxwell/analytic_solutions.py +++ b/examples/maxwell/analytic_solutions.py @@ -19,15 +19,12 @@ 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 +#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 +#import cmath +from grudge import sym +from six.moves import range, zip @@ -303,22 +300,64 @@ class RectangularWaveguideMode: return result +def get_rectangular_3D_cavity_mode(E_0, mode_indices, dimensions=(1, 1, 1)): + """A rectangular TM cavity mode.""" + kx, ky, kz = factors = [n*numpy.pi/a for n, a in zip(mode_indices, dimensions)] + omega = numpy.sqrt(sum(f**2 for f in factors)) + gamma_squared = ky**2 + kx**2 -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) + nodes = sym.nodes(3) + x = nodes[0] + y = nodes[1] + z = nodes[2] + + 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")) + + 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_2D_cavity_mode(E_0, mode_indices): + """A TM cavity mode. + + Returns an expression depending on *epsilon*, *mu*, and *t*. + """ + kx, ky = factors = [n*numpy.pi for n in mode_indices] + omega = numpy.sqrt(sum(f**2 for f in factors)) + + nodes = sym.nodes(2) + x = nodes[0] + y = nodes[1] + 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, + ) + return result # dipole solution ------------------------------------------------------------- diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index 6fb7d54a..fcea20c6 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -1,260 +1,134 @@ -# 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 - - 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 analytic_solutions import ( + get_rectangular_3D_cavity_mode, + get_rectangular_2D_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_3D_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)) + 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/maxwell/testing.py b/examples/maxwell/testing.py new file mode 100644 index 00000000..12b2f270 --- /dev/null +++ b/examples/maxwell/testing.py @@ -0,0 +1,161 @@ +"""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 +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 + + +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))] + k = np.sqrt(sum(f**2 for f in factors)) + + patch = spy.CalculusPatch((0.4, 0.3, 0.4)) + + 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)) + fields = bind(pdiscr, sym_mode)(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() + + print(frequency_domain_maxwell(patch, np.array([fields[0], fields[1], fields[2]]), np.array([fields[3], fields[4], fields[5]]), k)) + 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__": + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + 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/discretization.py b/grudge/discretization.py index 50e112de..e2f2cfe0 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -25,6 +25,7 @@ THE SOFTWARE. from pytools import memoize_method from grudge import sym +import numpy as np class DiscretizationBase(object): @@ -231,6 +232,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] diff --git a/grudge/execution.py b/grudge/execution.py index ace2dc8b..c4f285cd 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>> 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) -- GitLab From b5a50097845e8a6a215335bccb87f17356045e0f Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Fri, 3 Nov 2017 12:17:24 -0500 Subject: [PATCH 3/9] Satisfied Flake and Removed some dead code --- examples/maxwell/analytic_solutions.py | 419 +------------------------ examples/maxwell/testing.py | 6 +- grudge/discretization.py | 2 +- 3 files changed, 6 insertions(+), 421 deletions(-) diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py index c7617d91..9b8c0661 100644 --- a/examples/maxwell/analytic_solutions.py +++ b/examples/maxwell/analytic_solutions.py @@ -14,290 +14,12 @@ # along with this program. If not, see . - - from __future__ import division from __future__ import absolute_import from __future__ import print_function -#from math import sqrt, pi, sin, cos, atan2, acos import numpy -import numpy.linalg as la -#import cmath from grudge import sym -from six.moves import range, 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 +from six.moves import zip def get_rectangular_3D_cavity_mode(E_0, mode_indices, dimensions=(1, 1, 1)): @@ -358,142 +80,3 @@ def get_rectangular_2D_cavity_mode(E_0, mode_indices): ) return result - - -# 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/testing.py b/examples/maxwell/testing.py index 12b2f270..b9218443 100644 --- a/examples/maxwell/testing.py +++ b/examples/maxwell/testing.py @@ -121,12 +121,14 @@ def sumpy_test_3D(order, cl_ctx, queue): fields = bind(pdiscr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu) for i in range(len(fields)): - if type(fields[i]) == type(0): + if isinstance(fields[i], (int, float)): fields[i] = np.zeros(patch.points.shape[1]) else: fields[i] = fields[i].get() - print(frequency_domain_maxwell(patch, np.array([fields[0], fields[1], fields[2]]), np.array([fields[3], fields[4], fields[5]]), k)) + e = np.array([fields[0], fields[1], fields[2]]) + h = np.array([fields[3], fields[4], fields[5]]) + print(frequency_domain_maxwell(patch, e, h, k)) return diff --git a/grudge/discretization.py b/grudge/discretization.py index e2f2cfe0..817a0bc2 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -236,7 +236,7 @@ class PointsDiscretization(DiscretizationBase): self.complex_dtype = np.dtype({ np.float32: np.complex64, np.float64: np.complex128 -}[self.real_dtype.type]) + }[self.real_dtype.type]) def ambient_dim(self): return self._nodes.shape[0] -- GitLab From 014cd9bce5a30968dd11d79a41451a188fb63114 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Thu, 9 Nov 2017 13:03:27 -0600 Subject: [PATCH 4/9] 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 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 From 2c34a87ce84688c1f9b398cea116a34cbd6abe37 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Thu, 9 Nov 2017 13:14:39 -0600 Subject: [PATCH 5/9] Placate Flake --- test/test_grudge.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 609de951..b14f21d1 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -397,8 +397,6 @@ def test_convergence_maxwell(ctx_factory, order, visualize=False): 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")) -- GitLab From b9139ce1c2af21ae8e6029ba46dcb6d1dd85731d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 21 Nov 2017 00:01:01 -0600 Subject: [PATCH 6/9] Update Maxwell test for DG discr rename --- test/test_grudge.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 68c7b938..84c46168 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -353,7 +353,7 @@ def test_convergence_maxwell(ctx_factory, order, visualize=False): b=(1.0,)*dims, n=(n,)*dims) - discr = Discretization(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) epsilon = 1 mu = 1 -- GitLab From 506ef9dd7d4e14efa96b0d1370e6fd25d3f2a743 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 21 Nov 2017 00:01:15 -0600 Subject: [PATCH 7/9] Delete remnants of analytic_solutions.py --- examples/maxwell/analytic_solutions.py | 82 -------------------------- 1 file changed, 82 deletions(-) delete mode 100644 examples/maxwell/analytic_solutions.py diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py deleted file mode 100644 index 9b8c0661..00000000 --- a/examples/maxwell/analytic_solutions.py +++ /dev/null @@ -1,82 +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 . - - -from __future__ import division -from __future__ import absolute_import -from __future__ import print_function -import numpy -from grudge import sym -from six.moves import zip - - -def get_rectangular_3D_cavity_mode(E_0, mode_indices, dimensions=(1, 1, 1)): - """A rectangular TM cavity mode.""" - kx, ky, kz = factors = [n*numpy.pi/a for n, a in zip(mode_indices, dimensions)] - omega = numpy.sqrt(sum(f**2 for f in factors)) - - gamma_squared = ky**2 + kx**2 - - nodes = sym.nodes(3) - x = nodes[0] - y = nodes[1] - z = nodes[2] - - 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")) - - 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_2D_cavity_mode(E_0, mode_indices): - """A TM cavity mode. - - Returns an expression depending on *epsilon*, *mu*, and *t*. - """ - kx, ky = factors = [n*numpy.pi for n in mode_indices] - omega = numpy.sqrt(sum(f**2 for f in factors)) - - nodes = sym.nodes(2) - x = nodes[0] - y = nodes[1] - - 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, - ) - - return result -- GitLab From a1caee0d830b0bfe714171467354618d8c88c61b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 21 Nov 2017 00:06:50 -0600 Subject: [PATCH 8/9] Delete EM debug leftovers --- examples/maxwell/testing.py | 85 ------------------------------------- 1 file changed, 85 deletions(-) delete mode 100644 examples/maxwell/testing.py diff --git a/examples/maxwell/testing.py b/examples/maxwell/testing.py deleted file mode 100644 index 2ff45a68..00000000 --- a/examples/maxwell/testing.py +++ /dev/null @@ -1,85 +0,0 @@ -"""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 import bind - - -def sumpy_test_3D(order, cl_ctx, queue): - - epsilon = 1 - mu = 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)) - - from grudge.discretization import PointsDiscretization - pdiscr = PointsDiscretization(cl.array.to_device(queue, patch.points)) - - 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)): - if isinstance(fields[i], (int, float)): - fields[i] = np.zeros(patch.points.shape[1]) - else: - fields[i] = fields[i].get() - - e = np.array([fields[0], fields[1], fields[2]]) - h = np.array([fields[3], fields[4], fields[5]]) - print(frequency_domain_maxwell(patch, e, h, k)) - 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__": - cl_ctx = cl.create_some_context() - queue = cl.CommandQueue(cl_ctx) - - print("Doing Sumpy Test") - sumpy_test_3D(4, cl_ctx, queue) -- GitLab From a4e4a990688143c625d6d9f276ecb8506a8f33df Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 21 Nov 2017 00:22:32 -0600 Subject: [PATCH 9/9] Adapt to refactored DGDiscretizationWithBoundaries --- grudge/discretization.py | 8 ++++++++ grudge/execution.py | 12 ++++++------ 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/grudge/discretization.py b/grudge/discretization.py index e75d9661..a1ca6303 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -293,6 +293,14 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): def ambient_dim(self): 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 diff --git a/grudge/execution.py b/grudge/execution.py index a5d81a35..e4aebe5c 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -475,7 +475,7 @@ class BoundOperator(object): # {{{ process_sym_operator function -def process_sym_operator(discr, sym_operator, post_bind_mapper=None, +def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, dumper=lambda name, sym_operator: None): import grudge.symbolic.mappers as mappers @@ -483,14 +483,14 @@ def process_sym_operator(discr, sym_operator, post_bind_mapper=None, dumper("before-bind", sym_operator) sym_operator = mappers.OperatorBinder()(sym_operator) - mappers.ErrorChecker(discr.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) dumper("before-empty-flux-killer", sym_operator) - sym_operator = mappers.EmptyFluxKiller(discr.volume_discr.mesh)(sym_operator) + sym_operator = mappers.EmptyFluxKiller(discrwb.mesh)(sym_operator) dumper("before-cfold", sym_operator) sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator) @@ -506,8 +506,8 @@ def process_sym_operator(discr, sym_operator, post_bind_mapper=None, dumper("before-csize", sym_operator) sym_operator = mappers.ConstantToNumpyConversionMapper( - real_type=discr.volume_discr.real_dtype.type, - complex_type=discr.volume_discr.complex_dtype.type, + real_type=discrwb.real_dtype.type, + complex_type=discrwb.complex_dtype.type, )(sym_operator) # Ordering restriction: @@ -531,7 +531,7 @@ def process_sym_operator(discr, sym_operator, post_bind_mapper=None, # grids) that the operators will apply on. dumper("before-global-to-reference", sym_operator) - sym_operator = mappers.GlobalToReferenceMapper(discr.ambient_dim)(sym_operator) + sym_operator = mappers.GlobalToReferenceMapper(discrwb.ambient_dim)(sym_operator) # Ordering restriction: # -- GitLab