From 2b6b435048e63ac9faa92ca60d4f901230f788aa Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 6 Mar 2017 18:20:24 -0600 Subject: [PATCH 01/34] Wave is broken? --- examples/wave/var2.py | 129 +++++++++++++++++++++++++++++++++ examples/wave/weak.py | 125 ++++++++++++++++++++++++++++++++ grudge/models/wave.py | 162 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 416 insertions(+) create mode 100644 examples/wave/var2.py create mode 100644 examples/wave/weak.py diff --git a/examples/wave/var2.py b/examples/wave/var2.py new file mode 100644 index 00000000..772b9a8e --- /dev/null +++ b/examples/wave/var2.py @@ -0,0 +1,129 @@ +"""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 StrongWaveOperator + from meshmode.mesh import BTAG_ALL, BTAG_NONE + op = StrongWaveOperator( + sym.If(sym.Comparison( + np.dot(sym_x, sym_x), "<", 0.4**2), + 1, 0.5), + 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="central") + + 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 00000000..9df822ab --- /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/models/wave.py b/grudge/models/wave.py index 98f0b440..36655e38 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -180,6 +180,168 @@ class StrongWaveOperator(HyperbolicOperator): 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) + return flux_weak + + ##Below is trash + + 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 = 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) + + # }}} -- GitLab From 0a91ae4c0eb319d0f8dbaf60b81179960ab09e14 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 14 Mar 2017 02:47:58 -0500 Subject: [PATCH 02/34] Added upwinding, maybe --- grudge/models/wave.py | 253 +++++++++++++++++++----------------------- 1 file changed, 114 insertions(+), 139 deletions(-) diff --git a/grudge/models/wave.py b/grudge/models/wave.py index 36655e38..5901cd23 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -231,25 +231,27 @@ class WeakWaveOperator(HyperbolicOperator): flux_weak = join_fields( np.dot(v.avg, normal), u.avg * normal) - return flux_weak + #return -self.c*flux_weak - ##Below is trash + u_brac = (u.int - u.ext) * normal + v_brac = np.dot(v.int, normal) - np.dot(v.ext, normal) if self.flux_type == "central": - pass + return -self.c*flux_weak 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) + return -self.c * (flux_weak + 0.5 * join_fields(u_brac, v_brac)) + #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 + #flux_strong = join_fields( + #np.dot(v.int, normal), + #u.int * normal) - flux_weak - return self.c*flux_strong + #return self.c*flux_strong def sym_operator(self): d = self.ambient_dim @@ -320,12 +322,12 @@ class WeakWaveOperator(HyperbolicOperator): ) - -sym.FaceMassOperator()( flux(sym.int_tpair(w)) + - 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 @@ -347,182 +349,155 @@ class WeakWaveOperator(HyperbolicOperator): # {{{ variable-velocity -class VariableVelocityStrongWaveOperator(HyperbolicOperator): - r"""This operator discretizes the wave equation - :math:`\partial_t^2 u = c^2 \Delta u`. +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 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. - \partial_t v - c \nabla u = 0 + :math:`c` is assumed to be constant across all space. """ - def __init__( - self, c, ambient_dim, source=0, + 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, - time_sign=1, - diffusion_coeff=None, - diffusion_scheme=CentralSecondDerivative() - ): - """*c* and *source* are optemplate expressions. - """ + radiation_tag=BTAG_NONE): assert isinstance(ambient_dim, int) self.c = c - self.time_sign = time_sign self.ambient_dim = ambient_dim - self.source = source + 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.flux_type = flux_type + self.dirichlet_bc_f = dirichlet_bc_f - self.diffusion_coeff = diffusion_coeff - self.diffusion_scheme = diffusion_scheme + self.flux_type = flux_type - # {{{ flux ---------------------------------------------------------------- - def flux(self, w, c_vol): + def flux(self, w): u = w[0] v = w[1:] - normal = sym.normal(w.tag, self.ambient_dim) + normal = sym.normal(w.dd, self.ambient_dim) - c = sym.RestrictToBoundary(w.tag)(c_vol) + flux_weak = join_fields( + np.dot(v.avg, normal), + u.avg * normal) + return -self.c*flux_weak - flux = self.time_sign*1/2*join_fields( - c.ext * np.dot(v.ext, normal) - - c.int * np.dot(v.int, normal), - normal*(c.ext*u.ext - c.int*u.int)) if self.flux_type == "central": - pass - elif self.flux_type == "upwind": - flux += join_fields( - c.ext*u.ext - c.int*u.int, - c.ext*normal*np.dot(normal, v.ext) - - c.int*normal*np.dot(normal, v.int) - ) - else: - raise ValueError("invalid flux type '%s'" % self.flux_type) - - return flux - - # }}} - - def bind_characteristic_velocity(self, discr): - from grudge.symbolic.operators import ElementwiseMaxOperator - - compiled = discr.compile(ElementwiseMaxOperator()(self.c)) + return -self.c*flux_weak + #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) - def do(t, w, **extra_context): - return compiled(t=t, w=w, **extra_context) + #flux_strong = join_fields( + #np.dot(v.int, normal), + #u.int * normal) - flux_weak - return do + #return self.c*flux_strong - def sym_operator(self, with_sensor=False): + 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 - dir_c = sym.RestrictToBoundary(self.dirichlet_tag) * self.c - dir_u = sym.RestrictToBoundary(self.dirichlet_tag) * u - dir_v = sym.RestrictToBoundary(self.dirichlet_tag) * v - - dir_bc = join_fields(dir_c, -dir_u, dir_v) - - # Neumann - neu_c = sym.RestrictToBoundary(self.neumann_tag) * self.c - neu_u = sym.RestrictToBoundary(self.neumann_tag) * u - neu_v = sym.RestrictToBoundary(self.neumann_tag) * v - - neu_bc = join_fields(neu_c, neu_u, -neu_v) + # boundary conditions ------------------------------------------------- - # Radiation - rad_normal = sym.make_normal(self.radiation_tag, d) + # 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.") - rad_c = sym.RestrictToBoundary(self.radiation_tag) * self.c - rad_u = sym.RestrictToBoundary(self.radiation_tag) * u - rad_v = sym.RestrictToBoundary(self.radiation_tag) * v + 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) - rad_bc = join_fields( - rad_c, - 0.5*(rad_u - self.time_sign*np.dot(rad_normal, rad_v)), - 0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.time_sign*rad_u) - ) + 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") - # {{{ diffusion ------------------------------------------------------- - from pytools.obj_array import with_object_array_or_scalar + # radiation BCs ------------------------------------------------------- + rad_normal = sym.normal(self.radiation_tag, d) - def make_diffusion(arg): - if with_sensor or ( - self.diffusion_coeff is not None and self.diffusion_coeff != 0): - if self.diffusion_coeff is None: - diffusion_coeff = 0 - else: - diffusion_coeff = self.diffusion_coeff + rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u)) + rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v)) - if with_sensor: - diffusion_coeff += sym.Field("sensor") + 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") - from grudge.second_order import SecondDerivativeTarget + # entire operator ----------------------------------------------------- + nabla = sym.nabla(d) - # strong_form here allows the reuse the value of grad u. - grad_tgt = SecondDerivativeTarget( - self.ambient_dim, strong_form=True, - operand=arg) + def flux(pair): + return sym.interp(pair.dd, "all_faces")(self.flux(pair)) - self.diffusion_scheme.grad(grad_tgt, bc_getter=None, - dirichlet_tags=[], neumann_tags=[]) - div_tgt = SecondDerivativeTarget( - self.ambient_dim, strong_form=False, - operand=diffusion_coeff*grad_tgt.minv_all) + #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)) + #))) - self.diffusion_scheme.div(div_tgt, - bc_getter=None, - dirichlet_tags=[], neumann_tags=[]) + 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) + ) - return div_tgt.minv_all - else: - return 0 - # }}} + - 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)) - # entire operator ----------------------------------------------------- - nabla = sym.nabla(d) - flux_op = sym.get_flux_operator(self.flux()) + ) ) - return ( - - join_fields( - - self.time_sign*self.c*np.dot(nabla, v) - make_diffusion(u) - + self.source, + result[0] += self.source_f - -self.time_sign*self.c*(nabla*u) - with_object_array_or_scalar( - make_diffusion, v) - ) - + - sym.InverseMassOperator() * ( - flux_op(flux_w) - + flux_op(sym.BoundaryPair(flux_w, dir_bc, self.dirichlet_tag)) - + flux_op(sym.BoundaryPair(flux_w, neu_bc, self.neumann_tag)) - + flux_op(sym.BoundaryPair(flux_w, rad_bc, self.radiation_tag)) - )) + return result def check_bc_coverage(self, mesh): from meshmode.mesh import check_bc_coverage @@ -531,9 +506,9 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator): self.neumann_tag, self.radiation_tag]) - def max_eigenvalue_expr(self): - import grudge.symbolic as sym - return sym.NodalMax()(sym.CFunction("fabs")(self.c)) + def max_eigenvalue(self, t, fields=None, discr=None): + return abs(self.c) + # }}} -- GitLab From 93b65969ed5fbb978bcb4e824f15009555337476 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 20 Mar 2017 22:11:53 -0500 Subject: [PATCH 03/34] Made variable coefficient wave equation work --- examples/wave/var2.py | 11 +++++----- grudge/execution.py | 11 ++++++++++ grudge/models/wave.py | 49 ++++++++++++++++++++++--------------------- 3 files changed, 42 insertions(+), 29 deletions(-) diff --git a/examples/wave/var2.py b/examples/wave/var2.py index 772b9a8e..24719059 100644 --- a/examples/wave/var2.py +++ b/examples/wave/var2.py @@ -51,13 +51,14 @@ def main(write_output=True, order=4): 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 StrongWaveOperator + from grudge.models.wave import VariableCoefficientWeakWaveOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE - op = StrongWaveOperator( - sym.If(sym.Comparison( - np.dot(sym_x, sym_x), "<", 0.4**2), - 1, 0.5), + op = VariableCoefficientWeakWaveOperator(c, discr.dim, source_f=( sym.sin(source_omega*sym_t) diff --git a/grudge/execution.py b/grudge/execution.py index 0241c82b..b88b0fe2 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -119,6 +119,17 @@ class ExecutionMapper(mappers.Evaluator, then = self.rec(expr.then) else_ = self.rec(expr.else_) + if (type(then) is not pyopencl.array.Array): + fill_val = then + then = cl.array.empty(self.queue,bool_crit.shape, dtype=np.float32) #what type should this be lol? + then.fill(fill_val) + + if (type(else_) is not pyopencl.array.Array): + fill_val = else_ + else_ = cl.array.empty(self.queue,bool_crit.shape, dtype=np.float32) + else_.fill(fill_val) + + result = cl.array.empty_like(then, queue=self.queue, allocator=self.bound_op.allocator) cl.array.if_positive(bool_crit, then, else_, out=result, diff --git a/grudge/models/wave.py b/grudge/models/wave.py index 5901cd23..6ea1d923 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -231,16 +231,16 @@ class WeakWaveOperator(HyperbolicOperator): flux_weak = join_fields( np.dot(v.avg, normal), u.avg * normal) - #return -self.c*flux_weak - u_brac = (u.int - u.ext) * normal - v_brac = np.dot(v.int, normal) - np.dot(v.ext, 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)))) # see doc/notes/grudge-notes.tm - return -self.c * (flux_weak + 0.5 * join_fields(u_brac, v_brac)) + # THis is terrible #flux_weak -= self.sign*join_fields( #0.5*(u.int-u.ext), #0.5*(normal * np.dot(normal, v.int-v.ext))) @@ -379,10 +379,9 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): self.ambient_dim = ambient_dim self.source_f = source_f - #if self.c > 0: - #self.sign = 1 - #else: - self.sign = -1 + self.sign = sym.If(sym.Comparison( + self.c, ">", 0), + 1, -1) self.dirichlet_tag = dirichlet_tag self.neumann_tag = neumann_tag @@ -397,14 +396,16 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): v = w[1:] normal = sym.normal(w.dd, self.ambient_dim) + surf_c = sym.interp("vol", u.dd)(self.c) + flux_weak = join_fields( np.dot(v.avg, normal), u.avg * normal) - return -self.c*flux_weak + return -surf_c*flux_weak if self.flux_type == "central": - return -self.c*flux_weak + return -surf_c*flux_weak #elif self.flux_type == "upwind": # see doc/notes/grudge-notes.tm #flux_weak -= self.sign*join_fields( @@ -456,8 +457,8 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): 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) + 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 ----------------------------------------------------- @@ -467,19 +468,19 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): return sym.interp(pair.dd, "all_faces")(self.flux(pair)) - #result = ( - #- join_fields( - #-self.c*np.dot(nabla, v), - #-self.c*(nabla*u) + #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.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)) - #))) + + + #- 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 = sym.InverseMassOperator()( join_fields( -- GitLab From edba3f813a9b48d2452e47e3f51cb324ed557041 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 22 Mar 2017 17:27:12 -0500 Subject: [PATCH 04/34] Doing EM stuff --- examples/maxwell/cavities.py | 53 ++----- examples/maxwell/cavities2.py | 261 ++++++++++++++++++++++++++++++++++ examples/maxwell/cavities3.py | 117 +++++++++++++++ examples/wave/var2.py | 2 +- grudge/execution.py | 6 +- grudge/models/em.py | 66 ++------- grudge/models/wave.py | 30 +--- 7 files changed, 410 insertions(+), 125 deletions(-) create mode 100644 examples/maxwell/cavities2.py create mode 100644 examples/maxwell/cavities3.py diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index a387bb26..dc842914 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -20,14 +20,17 @@ from __future__ import absolute_import from __future__ import print_function import numpy as np -import logging -logger = logging.getLogger(__name__) +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, 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 + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + from math import sqrt, pi # noqa from analytic_solutions import ( # noqa RealPartAdapter, @@ -36,21 +39,16 @@ def main(write_output=True, allow_features=None, flux_type_arg=1, 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 + cylindrical = False #True is unsupported because meshmode periodic = False if cylindrical: @@ -62,36 +60,29 @@ def main(write_output=True, allow_features=None, flux_type_arg=1, # 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) + #mesh = make_cylinder_mesh(radius=R, height=d, max_volume=0.01) else: if periodic: - mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1)) + mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1)) #What is the point of this if it is overwritten? 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) + 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]: + discr = Discretization(mesh_data, order=order, + tune_for=op.sym_operator(), + **extra_discr_args) #for order in [1,2,3,4,5,6]: - extra_discr_args.setdefault("debug", []).extend([ - "cuda_no_plan", "cuda_dump_kernels"]) + from grudge.models.em import MaxwellOperator 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: @@ -188,20 +179,6 @@ def main(write_output=True, allow_features=None, flux_type_arg=1, eoc_rec.add_data_point(order, discr.norm(fields-get_true_field())) - finally: - if write_output: - vis.close() - - logmgr.close() - discr.close() - - if rcon.is_head_rank: - print() - print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) - - # }}} - - assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6 # {{{ entry points for py.test diff --git a/examples/maxwell/cavities2.py b/examples/maxwell/cavities2.py new file mode 100644 index 00000000..a387bb26 --- /dev/null +++ b/examples/maxwell/cavities2.py @@ -0,0 +1,261 @@ +# grudge - the Hybrid'n'Easy DG Environment +# 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() + + 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) + + 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() + + if rcon.is_head_rank: + print() + print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) + + # }}} + + assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6 + + +# {{{ entry points for py.test + +from pytools.test import mark_test + + +@mark_test.long +def test_maxwell_cavities(): + main(write_output=False) + + +@mark_test.long +def test_maxwell_cavities_lf(): + main(write_output=False, flux_type_arg="lf", bdry_flux_type_arg=1) + + +@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"]) + + +def test_cuda(): + try: + from pycuda.tools import mark_cuda_test + except ImportError: + pass + + 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 + + +def do_test_maxwell_cavities_cuda(dtype): + import pytest # noqa + + main(write_output=False, allow_features=["cuda"], + extra_discr_args=dict( + debug=["cuda_no_plan"], + default_scalar_type=dtype, + )) + +# }}} + + +# 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) diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py new file mode 100644 index 00000000..a78e5b6c --- /dev/null +++ b/examples/maxwell/cavities3.py @@ -0,0 +1,117 @@ +"""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) + + + epsilon0 = 8.8541878176e-12 # C**2 / (N m**2) + mu0 = 4*pi*1e-7 # N/A**2. + epsilon = 1*epsilon0 + mu = 1*mu0 + + + from grudge.models.wave import StrongWaveOperator + from meshmode.mesh import BTAG_ALL, BTAG_NONE + op = MaxwellOperator(epsilon,mu, + flux_type=flux_type_arg) + + + 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, + [ + ("e", event.state_component[0]), + ("h", event.state_component[1:]), + ]) + t_last_step = time() + + +if __name__ == "__main__": + main() diff --git a/examples/wave/var2.py b/examples/wave/var2.py index 24719059..e07acd44 100644 --- a/examples/wave/var2.py +++ b/examples/wave/var2.py @@ -68,7 +68,7 @@ def main(write_output=True, order=4): dirichlet_tag=BTAG_NONE, neumann_tag=BTAG_NONE, radiation_tag=BTAG_ALL, - flux_type="central") + flux_type="upwind") queue = cl.CommandQueue(discr.cl_context) from pytools.obj_array import join_fields diff --git a/grudge/execution.py b/grudge/execution.py index b88b0fe2..6390d99b 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -119,14 +119,14 @@ class ExecutionMapper(mappers.Evaluator, then = self.rec(expr.then) else_ = self.rec(expr.else_) - if (type(then) is not pyopencl.array.Array): + if (type(then) is not pyopencl.array.Array): # is instance, poarenthesis fill_val = then - then = cl.array.empty(self.queue,bool_crit.shape, dtype=np.float32) #what type should this be lol? + then = cl.array.empty(self.queue,bool_crit.shape, dtype=np.float64) #what type should this be? then.fill(fill_val) if (type(else_) is not pyopencl.array.Array): fill_val = else_ - else_ = cl.array.empty(self.queue,bool_crit.shape, dtype=np.float32) + else_ = cl.array.empty(self.queue,bool_crit.shape, dtype=np.float64) else_.fill(fill_val) diff --git a/grudge/models/em.py b/grudge/models/em.py index 3a2b70cb..a1e80976 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -106,17 +106,8 @@ class MaxwellOperator(HyperbolicOperator): self.current = current self.incident_bc_data = incident_bc - @property - def c(self): - from warnings import warn - warn("MaxwellOperator.c is deprecated", DeprecationWarning) - if not self.fixed_material: - raise RuntimeError("Cannot compute speed of light " - "for non-constant material") - - return 1/(self.mu*self.epsilon)**0.5 - def flux(self, flux_type): + def flux(self, flux_type, w): """The template for the numerical flux for variable coefficients. :param flux_type: can be in [0,1] for anything between central and upwind, @@ -124,38 +115,22 @@ class MaxwellOperator(HyperbolicOperator): As per Hesthaven and Warburton page 433. """ - from grudge.flux import (make_normal, FluxVectorPlaceholder, - FluxConstantPlaceholder) - normal = make_normal(self.dimensions) + normal = sym.normal(w.dd, self.dimensions) if self.fixed_material: - from grudge.tools import count_subset - w = FluxVectorPlaceholder(count_subset(self.get_eh_subset())) - e, h = self.split_eh(w) - epsilon = FluxConstantPlaceholder(self.epsilon) - mu = FluxConstantPlaceholder(self.mu) - - else: - from grudge.tools import count_subset - w = FluxVectorPlaceholder(count_subset(self.get_eh_subset())+2) - - epsilon, mu, e, h = self.split_eps_mu_eh(w) + epsilon = self.elpsion + mu = self.mu - Z_int = (mu.int/epsilon.int)**0.5 + Z_int = (mu/epsilon)**0.5 Y_int = 1/Z_int - Z_ext = (mu.ext/epsilon.ext)**0.5 + Z_ext = (mu/epsilon)**0.5 Y_ext = 1/Z_ext if flux_type == "lf": if self.fixed_material: max_c = (self.epsilon*self.mu)**(-0.5) - else: - from grudge.flux import Max - c_int = (epsilon.int*mu.int)**(-0.5) - c_ext = (epsilon.ext*mu.ext)**(-0.5) - max_c = Max(c_int, c_ext) # noqa return join_fields( # flux e, @@ -195,7 +170,7 @@ class MaxwellOperator(HyperbolicOperator): :math:`E` and :math:`H` """ - e, h = self.split_eh(self.field_placeholder(w)) + e, h = self.split_eh(w) def e_curl(field): return self.space_cross_e(nabla, field) @@ -205,7 +180,7 @@ class MaxwellOperator(HyperbolicOperator): from grudge.symbolic import make_nabla - nabla = make_nabla(self.dimensions) + nabla = sym.nabla(self.dimensions) # in conservation form: u_t + A u_x = 0 return join_fields( @@ -213,16 +188,6 @@ class MaxwellOperator(HyperbolicOperator): e_curl(e) ) - def field_placeholder(self, w=None): - "A placeholder for E and H." - from grudge.tools import count_subset - fld_cnt = count_subset(self.get_eh_subset()) - if w is None: - from grudge.symbolic import make_sym_vector - w = make_sym_vector("w", fld_cnt) - - return w - def pec_bc(self, w=None): "Construct part of the flux operator template for PEC boundary conditions" e, h = self.split_eh(self.field_placeholder(w)) @@ -310,21 +275,12 @@ class MaxwellOperator(HyperbolicOperator): Combines the relevant operator templates for spatial derivatives, flux, boundary conditions etc. """ - w = self.field_placeholder(w) + w = sym.make_sum_array("w", self.dimensions) # probs wrong # dimensions if self.fixed_material: flux_w = w - else: - epsilon = self.epsilon - mu = self.mu - flux_w = join_fields(epsilon, mu, w) - from grudge.symbolic import BoundaryPair, \ - InverseMassOperator, get_flux_operator - - flux_op = get_flux_operator(self.flux(self.flux_type)) - bdry_flux_op = get_flux_operator(self.flux(self.bdry_flux_type)) from grudge.tools.indexing import count_subset elec_components = count_subset(self.get_eh_subset()[0:3]) @@ -334,10 +290,6 @@ class MaxwellOperator(HyperbolicOperator): # need to check this material_divisor = ( [self.epsilon]*elec_components+[self.mu]*mag_components) - else: - material_divisor = join_fields( - [epsilon]*elec_components, - [mu]*mag_components) tags_and_bcs = [ (self.pec_tag, self.pec_bc(w)), diff --git a/grudge/models/wave.py b/grudge/models/wave.py index 6ea1d923..df17df4c 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -239,19 +239,6 @@ class WeakWaveOperator(HyperbolicOperator): 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)))) - # see doc/notes/grudge-notes.tm - # THis is terrible - #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 @@ -406,19 +393,10 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): if self.flux_type == "central": return -surf_c*flux_weak - #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 + if 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)))) def sym_operator(self): d = self.ambient_dim -- GitLab From 07c3e0144d52f6a88dce550ade5c4164ab637a95 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 28 Mar 2017 08:02:11 -0500 Subject: [PATCH 05/34] Copy pasted cross product in --- grudge/tools.py | 78 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) diff --git a/grudge/tools.py b/grudge/tools.py index cd1cc49f..fa329cad 100644 --- a/grudge/tools.py +++ b/grudge/tools.py @@ -36,3 +36,81 @@ def is_zero(x): return x == 0 else: return False + +def levi_civita(tuple): + """Compute an entry of the Levi-Civita tensor for the indices *tuple*. + Only three-tuples are supported for now. + """ + if len(tuple) == 3: + if tuple in [(0,1,2), (2,0,1), (1,2,0)]: + return 1 + elif tuple in [(2,1,0), (0,2,1), (1,0,2)]: + return -1 + else: + return 0 + else: + raise NotImplementedError + + +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() -- GitLab From 531f27e7d4ecc5ec668daf49b2dedd84b08b6d80 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 4 Apr 2017 10:37:09 -0500 Subject: [PATCH 06/34] Made em compile --- examples/maxwell/cavities3.py | 10 +-- grudge/models/em.py | 123 ++++++++++++++-------------------- grudge/tools.py | 24 +++++++ 3 files changed, 81 insertions(+), 76 deletions(-) diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py index a78e5b6c..0913be6e 100644 --- a/examples/maxwell/cavities3.py +++ b/examples/maxwell/cavities3.py @@ -35,7 +35,7 @@ def main(write_output=True, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - dims = 2 + dims = 3 from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh( a=(-0.5,)*dims, @@ -46,21 +46,21 @@ def main(write_output=True, order=4): epsilon0 = 8.8541878176e-12 # C**2 / (N m**2) - mu0 = 4*pi*1e-7 # N/A**2. + mu0 = 4*np.pi*1e-7 # N/A**2. epsilon = 1*epsilon0 mu = 1*mu0 - from grudge.models.wave import StrongWaveOperator + from grudge.models.em import MaxwellOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE op = MaxwellOperator(epsilon,mu, - flux_type=flux_type_arg) + flux_type=1) 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)]) + [discr.zeros(queue) for i in range(discr.dim*2 - 1)]) # FIXME #dt = op.estimate_rk4_timestep(discr, fields=fields) diff --git a/grudge/models/em.py b/grudge/models/em.py index a1e80976..ad826564 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -30,8 +30,8 @@ THE SOFTWARE. from pytools import memoize_method from grudge.models import HyperbolicOperator -from grudge.symbolic.primitives import make_common_subexpression as cse from meshmode.mesh import BTAG_ALL, BTAG_NONE +from grudge import sym from pytools.obj_array import join_fields, make_obj_array # TODO: Check PML @@ -107,7 +107,7 @@ class MaxwellOperator(HyperbolicOperator): self.incident_bc_data = incident_bc - def flux(self, flux_type, w): + def flux(self, w): """The template for the numerical flux for variable coefficients. :param flux_type: can be in [0,1] for anything between central and upwind, @@ -120,15 +120,16 @@ class MaxwellOperator(HyperbolicOperator): if self.fixed_material: e, h = self.split_eh(w) - epsilon = self.elpsion + epsilon = self.epsilon mu = self.mu + Z_int = (mu/epsilon)**0.5 Y_int = 1/Z_int Z_ext = (mu/epsilon)**0.5 Y_ext = 1/Z_ext - if flux_type == "lf": + if self.flux_type == "lf": if self.fixed_material: max_c = (self.epsilon*self.mu)**(-0.5) @@ -145,42 +146,41 @@ class MaxwellOperator(HyperbolicOperator): # multiplication by mu undoes material divisor below #-max_c*(mu.int*h.int - mu.ext*h.ext) )) - elif isinstance(flux_type, (int, float)): + elif isinstance(self.flux_type, (int, float)): # see doc/maxima/maxwell.mac return join_fields( # flux e, ( -1/(Z_int+Z_ext)*self.space_cross_h(normal, Z_ext*(h.int-h.ext) - - flux_type*self.space_cross_e(normal, e.int-e.ext)) + - self.flux_type*self.space_cross_e(normal, e.int-e.ext)) ), # flux h ( 1/(Y_int + Y_ext)*self.space_cross_e(normal, Y_ext*(e.int-e.ext) - + flux_type*self.space_cross_h(normal, h.int-h.ext)) + + self.flux_type*self.space_cross_h(normal, h.int-h.ext)) ), ) else: raise ValueError("maxwell: invalid flux_type (%s)" % self.flux_type) - def local_derivatives(self, w=None): + def local_derivatives(self, w): """Template for the spatial derivatives of the relevant components of :math:`E` and :math:`H` """ e, h = self.split_eh(w) + nabla = sym.nabla(self.dimensions) def e_curl(field): return self.space_cross_e(nabla, field) def h_curl(field): return self.space_cross_h(nabla, field) - from grudge.symbolic import make_nabla - nabla = sym.nabla(self.dimensions) # in conservation form: u_t + A u_x = 0 return join_fields( @@ -188,52 +188,44 @@ class MaxwellOperator(HyperbolicOperator): e_curl(e) ) - def pec_bc(self, w=None): + def pec_bc(self, w): "Construct part of the flux operator template for PEC boundary conditions" - e, h = self.split_eh(self.field_placeholder(w)) + e, h = self.split_eh(w) + + pec_e = sym.cse(sym.interp("vol", self.pec_tag)(e)) + pec_h = sym.cse(sym.interp("vol", self.pec_tag)(h)) - from grudge.symbolic import RestrictToBoundary - pec_e = RestrictToBoundary(self.pec_tag)(e) - pec_h = RestrictToBoundary(self.pec_tag)(h) return join_fields(-pec_e, pec_h) - def pmc_bc(self, w=None): + def pmc_bc(self, w): "Construct part of the flux operator template for PMC boundary conditions" - e, h = self.split_eh(self.field_placeholder(w)) + e, h = self.split_eh(w) - from grudge.symbolic import RestrictToBoundary - pmc_e = RestrictToBoundary(self.pmc_tag)(e) - pmc_h = RestrictToBoundary(self.pmc_tag)(h) + pmc_e = sym.cse(sym.interp("vol", self.pmc_tag)(e)) + pmc_h = sym.cse(sym.interp("vol", self.pmc_tag)(h)) return join_fields(pmc_e, -pmc_h) - def absorbing_bc(self, w=None): + def absorbing_bc(self, w): """Construct part of the flux operator template for 1st order absorbing boundary conditions. """ - from grudge.symbolic import normal - absorb_normal = normal(self.absorb_tag, self.dimensions) + absorb_normal = sym.normal(self.absorb_tag, self.dimensions) - from grudge.symbolic import RestrictToBoundary, Field - e, h = self.split_eh(self.field_placeholder(w)) + e, h = self.split_eh(w) if self.fixed_material: epsilon = self.epsilon mu = self.mu - else: - epsilon = cse( - RestrictToBoundary(self.absorb_tag)(Field("epsilon"))) - mu = cse( - RestrictToBoundary(self.absorb_tag)(Field("mu"))) absorb_Z = (mu/epsilon)**0.5 absorb_Y = 1/absorb_Z - absorb_e = RestrictToBoundary(self.absorb_tag)(e) - absorb_h = RestrictToBoundary(self.absorb_tag)(h) + absorb_e = sym.cse(sym.interp("vol", self.absorb_tag)(e)) + absorb_h = sym.cse(sym.interp("vol", self.absorb_tag)(h)) bc = join_fields( absorb_e + 1/2*(self.space_cross_h(absorb_normal, self.space_cross_e( @@ -246,17 +238,12 @@ class MaxwellOperator(HyperbolicOperator): return bc - def incident_bc(self, w=None): + def incident_bc(self, w): "Flux terms for incident boundary conditions" # NOTE: Untested for inhomogeneous materials, but would usually be # physically meaningless anyway (are there exceptions to this?) - e, h = self.split_eh(self.field_placeholder(w)) - if not self.fixed_material: - from warnings import warn - if self.incident_tag != BTAG_NONE: - warn("Incident boundary conditions assume homogeneous" - " background material, results may be unphysical") + e, h = self.split_eh(w) from grudge.tools import count_subset fld_cnt = count_subset(self.get_eh_subset()) @@ -266,7 +253,7 @@ class MaxwellOperator(HyperbolicOperator): if is_zero(incident_bc_data): return make_obj_array([0]*fld_cnt) else: - return cse(-incident_bc_data) + return sym.cse(-incident_bc_data) def sym_operator(self, w=None): """The full operator template - the high level description of @@ -275,14 +262,9 @@ class MaxwellOperator(HyperbolicOperator): Combines the relevant operator templates for spatial derivatives, flux, boundary conditions etc. """ - w = sym.make_sum_array("w", self.dimensions) # probs wrong # dimensions - - if self.fixed_material: - flux_w = w - - + from grudge.tools import count_subset + w = sym.make_sym_array("w", count_subset(self.get_eh_subset())) - from grudge.tools.indexing import count_subset elec_components = count_subset(self.get_eh_subset()[0:3]) mag_components = count_subset(self.get_eh_subset()[3:6]) @@ -293,30 +275,24 @@ class MaxwellOperator(HyperbolicOperator): tags_and_bcs = [ (self.pec_tag, self.pec_bc(w)), - (self.pmc_tag, self.pmc_bc(w)), - (self.absorb_tag, self.absorbing_bc(w)), - (self.incident_tag, self.incident_bc(w)), + #(self.pmc_tag, self.pmc_bc(w)), + #(self.absorb_tag, self.absorbing_bc(w)), + #(self.incident_tag, self.incident_bc(w)), ] - def make_flux_bc_vector(tag, bc): - if self.fixed_material: - return bc - else: - from grudge.symbolic import RestrictToBoundary - return join_fields( - cse(RestrictToBoundary(tag)(epsilon)), - cse(RestrictToBoundary(tag)(mu)), - bc) + + def flux(pair): + return sym.interp(pair.dd, "all_faces")(self.flux(pair)) return ( - self.local_derivatives(w) - + InverseMassOperator()( - flux_op(flux_w) - + sum( - bdry_flux_op(BoundaryPair( - flux_w, make_flux_bc_vector(tag, bc), tag)) - for tag, bc in tags_and_bcs)) - ) / material_divisor + + sym.InverseMassOperator()( + flux(sym.int_tpair(w)) + #+ sum( + #flux(sym.bv_tpair(tag, w, bc)) + #for tag, bc in tags_and_bcs)) + )) / material_divisor + def bind(self, discr): "Convert the abstract operator template into compiled code." @@ -397,11 +373,8 @@ class MaxwellOperator(HyperbolicOperator): e_idx, h_idx = self.partial_to_eh_subsets() e, h = w[e_idx], w[h_idx] - from grudge.flux import FluxVectorPlaceholder as FVP - if isinstance(w, FVP): - return FVP(scalars=e), FVP(scalars=h) - else: - return make_obj_array(e), make_obj_array(h) + return e,h + #return make_obj_array(e), make_obj_array(h) def get_eh_subset(self): """Return a 6-tuple of :class:`bool` objects indicating whether field @@ -428,6 +401,14 @@ class MaxwellOperator(HyperbolicOperator): raise ValueError("max_eigenvalue is no longer supported for " "variable-coefficient problems--use max_eigenvalue_expr") + def check_bc_coverage(self, mesh): + from meshmode.mesh import check_bc_coverage + check_bc_coverage(mesh, [ + self.pec_tag, + self.pmc_tag, + self.absorb_tag, + self.incident_tag]) + class TMMaxwellOperator(MaxwellOperator): """A 2D TM Maxwell operator with PEC boundaries. diff --git a/grudge/tools.py b/grudge/tools.py index fa329cad..cf5aea0b 100644 --- a/grudge/tools.py +++ b/grudge/tools.py @@ -114,3 +114,27 @@ class SubsettableCrossProduct: 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) + + -- GitLab From 3c63a45f719fbfb2b71091745676a0286c62a94c Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 11 Apr 2017 10:47:33 -0500 Subject: [PATCH 07/34] Changed to simpler flux --- examples/maxwell/cavities3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py index 0913be6e..cae35f76 100644 --- a/examples/maxwell/cavities3.py +++ b/examples/maxwell/cavities3.py @@ -54,7 +54,7 @@ def main(write_output=True, order=4): from grudge.models.em import MaxwellOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE op = MaxwellOperator(epsilon,mu, - flux_type=1) + flux_type="lf") queue = cl.CommandQueue(discr.cl_context) -- GitLab From 93e7bdbcbfbef4c2552a588ad6e6b2ab95ad60dc Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 17 Apr 2017 19:21:49 -0500 Subject: [PATCH 08/34] Made EM compile --- examples/maxwell/analytic_solutions.py | 6 +++--- examples/maxwell/cavities3.py | 8 ++++++-- grudge/models/em.py | 12 ++++++------ 3 files changed, 15 insertions(+), 11 deletions(-) diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py index 152130f9..eaae8cdc 100644 --- a/examples/maxwell/analytic_solutions.py +++ b/examples/maxwell/analytic_solutions.py @@ -20,9 +20,9 @@ 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 grudge.tools import \ + #cyl_bessel_j, \ + #cyl_bessel_j_prime from math import sqrt, pi, sin, cos, atan2, acos import numpy import numpy.linalg as la diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py index cae35f76..b06632dc 100644 --- a/examples/maxwell/cavities3.py +++ b/examples/maxwell/cavities3.py @@ -30,6 +30,9 @@ import pyopencl as cl from grudge.shortcuts import set_up_rk4 from grudge import sym, bind, Discretization +from analytic_solutions import ( # noqa + RectangularWaveguideMode, + RectangularCavityMode) def main(write_output=True, order=4): cl_ctx = cl.create_some_context() @@ -59,8 +62,9 @@ def main(write_output=True, order=4): 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*2 - 1)]) + fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr) + #fields = join_fields(discr.zeros(queue), + #[discr.zeros(queue) for i in range(discr.dim*2 - 1)]) # FIXME #dt = op.estimate_rk4_timestep(discr, fields=fields) diff --git a/grudge/models/em.py b/grudge/models/em.py index ad826564..e728dd78 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -275,9 +275,9 @@ class MaxwellOperator(HyperbolicOperator): tags_and_bcs = [ (self.pec_tag, self.pec_bc(w)), - #(self.pmc_tag, self.pmc_bc(w)), - #(self.absorb_tag, self.absorbing_bc(w)), - #(self.incident_tag, self.incident_bc(w)), + (self.pmc_tag, self.pmc_bc(w)), + (self.absorb_tag, self.absorbing_bc(w)), + (self.incident_tag, self.incident_bc(w)), ] @@ -286,12 +286,12 @@ class MaxwellOperator(HyperbolicOperator): return ( - self.local_derivatives(w) - + sym.InverseMassOperator()( + + sym.InverseMassOperator()(sym.FaceMassOperator()( flux(sym.int_tpair(w)) #+ sum( #flux(sym.bv_tpair(tag, w, bc)) - #for tag, bc in tags_and_bcs)) - )) / material_divisor + #for tag, bc in tags_and_bcs) + ))) / material_divisor def bind(self, discr): -- GitLab From d1e275c015177eae0a86fa3cfd237acaa3ac4db5 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 1 May 2017 21:06:56 -0500 Subject: [PATCH 09/34] progress with EM --- examples/maxwell/analytic_solutions.py | 88 ++++++++++++++++++++++++-- examples/maxwell/cavities3.py | 11 +++- 2 files changed, 92 insertions(+), 7 deletions(-) diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py index eaae8cdc..6e57f7ff 100644 --- a/examples/maxwell/analytic_solutions.py +++ b/examples/maxwell/analytic_solutions.py @@ -273,11 +273,11 @@ class RectangularWaveguideMode: 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] + nodes = discr.nodes() + x = nodes[:,0] + y = nodes[:,1] + if discr.dim > 2: + z = nodes[:,2] else: z = discr.volume_zeros() @@ -318,6 +318,84 @@ class RectangularCavityMode(RectangularWaveguideMode): kwargs["backward_coeff"] = 1 RectangularWaveguideMode.__init__(self, *args, **kwargs) +class SymRectangularWaveguideMode: + """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, nodes): + f,g,h = self.factors + omega = self.omega + k = self.k + + x = nodes[0] + y = nodes[1] + if nodes.shape[0] > 2: + z = nodes[2] + else: + z = 0 + from grudge import sym + + sx = sym.sin(f*x) + sy = sym.sin(g*y) + cx = sym.cos(f*x) + cy = sym.cos(g*y) + + zdep_add = sym.exp(1j*h*z)+sym.exp(-1j*h*z) + zdep_sub = sym.exp(1j*h*z)-sym.exp(-1j*h*z) + + tdep = numpy.exp(-1j * omega * self.t) + + C = 1j/(f**2+g**2) + + from pytools.obj_array import join_fields + result = join_fields( C*f*h*cx*sy*zdep_sub*tdep # ex + , C*g*h*sx*cy*zdep_sub*tdep # ey + , sx*sy*zdep_add*tdep # ez + , -C*g*self.epsilon*omega*sx*cy*zdep_add*tdep # hx + , C*f*self.epsilon*omega*cx*sy*zdep_add*tdep + , C*f*self.epsilon*omega*cx*sy*zdep_add*tdep + ) # THIS IS BAD, last line should be all zeros + + return result + + + + +class SymRectangularCavityMode(SymRectangularWaveguideMode): + """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) diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py index b06632dc..226f74ba 100644 --- a/examples/maxwell/cavities3.py +++ b/examples/maxwell/cavities3.py @@ -32,7 +32,9 @@ from grudge import sym, bind, Discretization from analytic_solutions import ( # noqa RectangularWaveguideMode, - RectangularCavityMode) + RectangularCavityMode, + SymRectangularWaveguideMode, + SymRectangularCavityMode) def main(write_output=True, order=4): cl_ctx = cl.create_some_context() @@ -62,7 +64,12 @@ def main(write_output=True, order=4): queue = cl.CommandQueue(discr.cl_context) from pytools.obj_array import join_fields - fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr) + + mode = SymRectangularCavityMode(epsilon, mu, (1, 2, 2)) + fields = bind(discr, mode(sym.nodes(dims)))(queue, t=0) + #fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr.volume_discr) + + #fields = (discr.volume_discr) #fields = join_fields(discr.zeros(queue), #[discr.zeros(queue) for i in range(discr.dim*2 - 1)]) -- GitLab From 23ea07b2e7238cc2d37eb788f8dbb0060f30499a Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Fri, 12 May 2017 11:54:16 -0500 Subject: [PATCH 10/34] Broke everything --- examples/maxwell/testing.py | 82 +++++++++++++++++++++++++++++++++++++ grudge/execution.py | 39 +++++++++++++++--- 2 files changed, 115 insertions(+), 6 deletions(-) create mode 100644 examples/maxwell/testing.py diff --git a/examples/maxwell/testing.py b/examples/maxwell/testing.py new file mode 100644 index 00000000..96de8f24 --- /dev/null +++ b/examples/maxwell/testing.py @@ -0,0 +1,82 @@ +"""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 + +from analytic_solutions import ( # noqa + RectangularWaveguideMode, + RectangularCavityMode, + SymRectangularWaveguideMode, + SymRectangularCavityMode) + +def main(write_output=True, order=4): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + + from grudge.models.em import MaxwellOperator + from meshmode.mesh import BTAG_ALL, BTAG_NONE + + + queue = cl.CommandQueue(cl_ctx) + from pytools.obj_array import join_fields + + from pymbolic.primitives import Variable as var + import pymbolic.primitives as p + import loopy as lp + + + nom = "abs" + i = var("i") + func = var(nom) + if("fabs" == nom): + func = var("abs") + + def knl(): + knl = lp.make_kernel( + "{[i]: 0<=i Date: Fri, 12 May 2017 12:40:00 -0500 Subject: [PATCH 11/34] nothing --- grudge/execution.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index 82549719..e35ed30a 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -113,8 +113,8 @@ class ExecutionMapper(mappers.Evaluator, i = var("i") func = var(expr.function.name) - #if(expr.function.name == "fabs"): - #func = var("abs") + if(expr.function.name == "fabs"): + func = var("abs") @memoize_in(self.bound_op, "map_call_knl") def knl(): -- GitLab From 9e315eef8596371a984daf00ee2221a6c21fe5cb Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 15 May 2017 16:07:46 -0500 Subject: [PATCH 12/34] Got openCL functions to use loopy --- examples/maxwell/cavities3.py | 10 ++++++++-- grudge/execution.py | 19 +++++++++---------- grudge/models/em.py | 6 +++--- 3 files changed, 20 insertions(+), 15 deletions(-) diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py index 226f74ba..c69e05f7 100644 --- a/examples/maxwell/cavities3.py +++ b/examples/maxwell/cavities3.py @@ -45,7 +45,7 @@ def main(write_output=True, order=4): mesh = generate_regular_rect_mesh( a=(-0.5,)*dims, b=(0.5,)*dims, - n=(8,)*dims) + n=(4,)*dims) discr = Discretization(cl_ctx, mesh, order=order) @@ -107,6 +107,12 @@ def main(write_output=True, order=4): from time import time t_last_step = time() + vis.write_vtk_file("fld-%04d.vtu" % step, + [ + ("e", fields[0:2]), + ("h", fields[2:]), + ]) + for event in dt_stepper.run(t_end=final_t): if isinstance(event, dt_stepper.StateComputed): assert event.component_id == "w" @@ -115,7 +121,7 @@ def main(write_output=True, order=4): print(step, event.t, norm(queue, u=event.state_component[0]), time()-t_last_step) - if step % 10 == 0: + if step % 1 == 0: vis.write_vtk_file("fld-%04d.vtu" % step, [ ("e", event.state_component[0]), diff --git a/grudge/execution.py b/grudge/execution.py index e35ed30a..160c5b1e 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -98,25 +98,24 @@ 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: - #func = getattr(np, expr.function.name) - - #HELD TOGETHER BY GLUE AND DUCT TAPE - if(pars[0].shape == ()): + if pars[0].shape == (): import pyopencl.clmath as clmath func = getattr(clmath, expr.function.name) return func(*pars) + cached_name = "map_call_knl_" + i = var("i") func = var(expr.function.name) - if(expr.function.name == "fabs"): + if expr.function.name == "fabs": #FIXME func = var("abs") + cached_name += "abs" + else: + cached_name += expr.function.name + - @memoize_in(self.bound_op, "map_call_knl") + @memoize_in(self.bound_op, cached_name) def knl(): knl = lp.make_kernel( "{[i]: 0<=i Date: Tue, 16 May 2017 14:44:09 -0500 Subject: [PATCH 13/34] Cleaned some stuff up --- examples/maxwell/cavities3.py | 14 +++++---- grudge/models/em.py | 53 ++++------------------------------- 2 files changed, 14 insertions(+), 53 deletions(-) diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py index c69e05f7..f79563ff 100644 --- a/examples/maxwell/cavities3.py +++ b/examples/maxwell/cavities3.py @@ -87,9 +87,9 @@ def main(write_output=True, order=4): return bound_op(queue, t=t, w=w) if mesh.dim == 2: - dt = 0.04 + dt = 0.004 elif mesh.dim == 3: - dt = 0.02 + dt = 0.002 dt_stepper = set_up_rk4("w", dt, fields, rhs) @@ -107,10 +107,11 @@ def main(write_output=True, order=4): from time import time t_last_step = time() + e,h = op.split_eh(fields) vis.write_vtk_file("fld-%04d.vtu" % step, [ - ("e", fields[0:2]), - ("h", fields[2:]), + ("e", e), + ("h", h), ]) for event in dt_stepper.run(t_end=final_t): @@ -122,10 +123,11 @@ def main(write_output=True, order=4): print(step, event.t, norm(queue, u=event.state_component[0]), time()-t_last_step) if step % 1 == 0: + e,h = op.split_eh(event.state_component) vis.write_vtk_file("fld-%04d.vtu" % step, [ - ("e", event.state_component[0]), - ("h", event.state_component[1:]), + ("e", e), + ("h", h), ]) t_last_step = time() diff --git a/grudge/models/em.py b/grudge/models/em.py index 434f7c00..ff491bef 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -138,13 +138,13 @@ class MaxwellOperator(HyperbolicOperator): 1/2*( -self.space_cross_h(normal, h.int-h.ext) # multiplication by epsilon undoes material divisor below - #-max_c*(epsilon.int*e.int - epsilon.ext*e.ext) + #-max_c*(epsilon*e.int - epsilon*e.ext) ), # flux h 1/2*( self.space_cross_e(normal, e.int-e.ext) # multiplication by mu undoes material divisor below - #-max_c*(mu.int*h.int - mu.ext*h.ext) + #-max_c*(mu*h.int - mu*h.ext) )) elif isinstance(self.flux_type, (int, float)): # see doc/maxima/maxwell.mac @@ -288,53 +288,12 @@ class MaxwellOperator(HyperbolicOperator): - self.local_derivatives(w) + sym.InverseMassOperator()(sym.FaceMassOperator()( flux(sym.int_tpair(w)) - + sum( - flux(sym.bv_tpair(tag, w, bc)) - for tag, bc in tags_and_bcs) - ))) / material_divisor + #+ sum( + #flux(sym.bv_tpair(tag, w, bc)) + #for tag, bc in tags_and_bcs) + ) )) / material_divisor - def bind(self, discr): - "Convert the abstract operator template into compiled code." - from grudge.mesh import check_bc_coverage - check_bc_coverage(discr.mesh, [ - self.pec_tag, self.absorb_tag, self.incident_tag]) - - compiled_sym_operator = discr.compile(self.op_template()) - - def rhs(t, w, **extra_context): - kwargs = {} - kwargs.update(extra_context) - - return compiled_sym_operator(w=w, t=t, **kwargs) - - return rhs - - def assemble_eh(self, e=None, h=None, discr=None): - "Combines separate E and H vectors into a single array." - if discr is None: - def zero(): - return 0 - else: - def zero(): - return discr.volume_zeros() - - from grudge.tools import count_subset - e_components = count_subset(self.get_eh_subset()[0:3]) - h_components = count_subset(self.get_eh_subset()[3:6]) - - def default_fld(fld, comp): - if fld is None: - return [zero() for i in range(comp)] - else: - return fld - - e = default_fld(e, e_components) - h = default_fld(h, h_components) - - return join_fields(e, h) - - assemble_fields = assemble_eh @memoize_method def partial_to_eh_subsets(self): -- GitLab From f75ccd4ce2cfab3ed682b14c33b45ff6cf26ba50 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Thu, 1 Jun 2017 20:58:49 -0700 Subject: [PATCH 14/34] Made time steps smaller --- examples/maxwell/cavities3.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py index f79563ff..0040755c 100644 --- a/examples/maxwell/cavities3.py +++ b/examples/maxwell/cavities3.py @@ -87,9 +87,9 @@ def main(write_output=True, order=4): return bound_op(queue, t=t, w=w) if mesh.dim == 2: - dt = 0.004 + dt = 0.00004 elif mesh.dim == 3: - dt = 0.002 + dt = 0.00002 dt_stepper = set_up_rk4("w", dt, fields, rhs) @@ -108,11 +108,11 @@ def main(write_output=True, order=4): t_last_step = time() e,h = op.split_eh(fields) - vis.write_vtk_file("fld-%04d.vtu" % step, - [ - ("e", e), - ("h", h), - ]) + #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): @@ -124,11 +124,11 @@ def main(write_output=True, order=4): time()-t_last_step) if step % 1 == 0: e,h = op.split_eh(event.state_component) - vis.write_vtk_file("fld-%04d.vtu" % step, - [ - ("e", e), - ("h", h), - ]) + #vis.write_vtk_file("fld-%04d.vtu" % step, + #[ + #("e", e), + #("h", h), + #]) t_last_step = time() -- GitLab From c5355dcfbfb6ea73afd53f6ddfab6934608aae1b Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sun, 11 Jun 2017 16:14:57 -0700 Subject: [PATCH 15/34] Added some debug logs --- examples/maxwell/cavities3.py | 6 ++++-- grudge/execution.py | 35 ++++++++++++++++++++++++++++++++++- grudge/shortcuts.py | 15 +++++++++++++++ 3 files changed, 53 insertions(+), 3 deletions(-) diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py index 0040755c..2a756fa6 100644 --- a/examples/maxwell/cavities3.py +++ b/examples/maxwell/cavities3.py @@ -28,6 +28,7 @@ THE SOFTWARE. import numpy as np import pyopencl as cl from grudge.shortcuts import set_up_rk4 +from grudge.shortcuts import set_up_forward_euler from grudge import sym, bind, Discretization from analytic_solutions import ( # noqa @@ -91,10 +92,11 @@ def main(write_output=True, order=4): elif mesh.dim == 3: dt = 0.00002 - dt_stepper = set_up_rk4("w", dt, fields, rhs) + dt_stepper = set_up_forward_euler("w", dt, fields, rhs) - final_t = 10 + final_t = dt * 0.9 nsteps = int(final_t/dt) + print("dt=%g nsteps=%d" % (dt, nsteps)) from grudge.shortcuts import make_visualizer diff --git a/grudge/execution.py b/grudge/execution.py index 160c5b1e..16f04042 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -226,6 +226,7 @@ class ExecutionMapper(mappers.Evaluator, knl()(self.queue, mat=matrix, result=grp.view(result), vec=grp.view(field)) + return result def map_elementwise_max(self, op, field_expr): @@ -358,11 +359,43 @@ class ExecutionMapper(mappers.Evaluator, # }}} # {{{ code execution functions + # DEBUG stuff + def get_max(self, arr): + from pymbolic.primitives import Variable as var + i = var("i") + + def knl2(): + knl = lp.make_kernel( + "{[i]: 0<=i"+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 -- GitLab From 680ffd29b2045c874b09f2e60f99c7e134e8b7df Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 17 Jul 2017 00:58:58 -0700 Subject: [PATCH 16/34] Removed scratchwork cavities, and re-enabled all the fluxes on EM --- examples/maxwell/cavities.py | 295 +++++++++++----------------------- examples/maxwell/cavities2.py | 261 ------------------------------ examples/maxwell/cavities3.py | 138 ---------------- grudge/models/em.py | 6 +- 4 files changed, 101 insertions(+), 599 deletions(-) delete mode 100644 examples/maxwell/cavities2.py delete mode 100644 examples/maxwell/cavities3.py diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index 0c5a33f5..2a756fa6 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -1,237 +1,138 @@ -# 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 +"""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.shortcuts import set_up_forward_euler from grudge import sym, bind, Discretization +from analytic_solutions import ( # noqa + RectangularWaveguideMode, + RectangularCavityMode, + SymRectangularWaveguideMode, + SymRectangularCavityMode) -def main(write_output=True, allow_features=None, flux_type_arg=1, - bdry_flux_type_arg=None, extra_discr_args={}): +def main(write_output=True, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - from math import sqrt, pi # noqa - from analytic_solutions import ( # noqa - RealPartAdapter, - SplitComplexAdapter, - CylindricalFieldAdapter, - CylindricalCavityMode, - RectangularWaveguideMode, - RectangularCavityMode) + dims = 3 + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*dims, + b=(0.5,)*dims, + n=(4,)*dims) + discr = Discretization(cl_ctx, mesh, order=order) epsilon0 = 8.8541878176e-12 # C**2 / (N m**2) - mu0 = 4*pi*1e-7 # N/A**2. + mu0 = 4*np.pi*1e-7 # N/A**2. epsilon = 1*epsilon0 mu = 1*mu0 - cylindrical = False #True is unsupported because meshmode - 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)) - - #mesh = make_cylinder_mesh(radius=R, height=d, max_volume=0.01) - else: - if periodic: - mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1)) #What is the point of this if it is overwritten? - periodicity = (False, False, True) - else: - periodicity = None - mode = RectangularCavityMode(epsilon, mu, (1, 2, 2)) - - mesh = make_box_mesh(max_volume=0.001, periodicity=periodicity) - - - for order in [4, 5, 6]: - discr = Discretization(mesh_data, order=order, - tune_for=op.sym_operator(), - **extra_discr_args) - #for order in [1,2,3,4,5,6]: - from grudge.models.em import MaxwellOperator - - op = MaxwellOperator(epsilon, mu, - flux_type=flux_type_arg, - bdry_flux_type=bdry_flux_type_arg) - - - 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) - - 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) + from grudge.models.em import MaxwellOperator + from meshmode.mesh import BTAG_ALL, BTAG_NONE + op = MaxwellOperator(epsilon,mu, + flux_type="lf") + - eoc_rec.add_data_point(order, - discr.norm(fields-get_true_field())) + queue = cl.CommandQueue(discr.cl_context) + from pytools.obj_array import join_fields + mode = SymRectangularCavityMode(epsilon, mu, (1, 2, 2)) + fields = bind(discr, mode(sym.nodes(dims)))(queue, t=0) + #fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr.volume_discr) + #fields = (discr.volume_discr) + #fields = join_fields(discr.zeros(queue), + #[discr.zeros(queue) for i in range(discr.dim*2 - 1)]) -# {{{ entry points for py.test + # FIXME + #dt = op.estimate_rk4_timestep(discr, fields=fields) -from pytools.test import mark_test + op.check_bc_coverage(mesh) + # print(sym.pretty(op.sym_operator())) + bound_op = bind(discr, op.sym_operator()) + print(bound_op) + # 1/0 -@mark_test.long -def test_maxwell_cavities(): - main(write_output=False) + def rhs(t, w): + return bound_op(queue, t=t, w=w) + if mesh.dim == 2: + dt = 0.00004 + elif mesh.dim == 3: + dt = 0.00002 -@mark_test.long -def test_maxwell_cavities_lf(): - main(write_output=False, flux_type_arg="lf", bdry_flux_type_arg=1) + dt_stepper = set_up_forward_euler("w", dt, fields, rhs) + final_t = dt * 0.9 + nsteps = int(final_t/dt) -@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"]) + print("dt=%g nsteps=%d" % (dt, nsteps)) + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) -def test_cuda(): - try: - from pycuda.tools import mark_cuda_test - except ImportError: - pass + step = 0 - 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 + norm = bind(discr, sym.norm(2, sym.var("u"))) + from time import time + t_last_step = time() -def do_test_maxwell_cavities_cuda(dtype): - import pytest # noqa + e,h = op.split_eh(fields) + #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=event.state_component[0]), + time()-t_last_step) + if step % 1 == 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/cavities2.py b/examples/maxwell/cavities2.py deleted file mode 100644 index a387bb26..00000000 --- a/examples/maxwell/cavities2.py +++ /dev/null @@ -1,261 +0,0 @@ -# grudge - the Hybrid'n'Easy DG Environment -# 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() - - 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) - - 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() - - if rcon.is_head_rank: - print() - print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) - - # }}} - - assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6 - - -# {{{ entry points for py.test - -from pytools.test import mark_test - - -@mark_test.long -def test_maxwell_cavities(): - main(write_output=False) - - -@mark_test.long -def test_maxwell_cavities_lf(): - main(write_output=False, flux_type_arg="lf", bdry_flux_type_arg=1) - - -@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"]) - - -def test_cuda(): - try: - from pycuda.tools import mark_cuda_test - except ImportError: - pass - - 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 - - -def do_test_maxwell_cavities_cuda(dtype): - import pytest # noqa - - main(write_output=False, allow_features=["cuda"], - extra_discr_args=dict( - debug=["cuda_no_plan"], - default_scalar_type=dtype, - )) - -# }}} - - -# 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) diff --git a/examples/maxwell/cavities3.py b/examples/maxwell/cavities3.py deleted file mode 100644 index 2a756fa6..00000000 --- a/examples/maxwell/cavities3.py +++ /dev/null @@ -1,138 +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 -from grudge.shortcuts import set_up_rk4 -from grudge.shortcuts import set_up_forward_euler -from grudge import sym, bind, Discretization - -from analytic_solutions import ( # noqa - RectangularWaveguideMode, - RectangularCavityMode, - SymRectangularWaveguideMode, - SymRectangularCavityMode) - -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.5,)*dims, - b=(0.5,)*dims, - n=(4,)*dims) - - discr = Discretization(cl_ctx, mesh, order=order) - - - epsilon0 = 8.8541878176e-12 # C**2 / (N m**2) - mu0 = 4*np.pi*1e-7 # N/A**2. - epsilon = 1*epsilon0 - mu = 1*mu0 - - - from grudge.models.em import MaxwellOperator - from meshmode.mesh import BTAG_ALL, BTAG_NONE - op = MaxwellOperator(epsilon,mu, - flux_type="lf") - - - queue = cl.CommandQueue(discr.cl_context) - from pytools.obj_array import join_fields - - mode = SymRectangularCavityMode(epsilon, mu, (1, 2, 2)) - fields = bind(discr, mode(sym.nodes(dims)))(queue, t=0) - #fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr.volume_discr) - - #fields = (discr.volume_discr) - #fields = join_fields(discr.zeros(queue), - #[discr.zeros(queue) for i in range(discr.dim*2 - 1)]) - - # 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.00004 - elif mesh.dim == 3: - dt = 0.00002 - - dt_stepper = set_up_forward_euler("w", dt, fields, rhs) - - final_t = dt * 0.9 - 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) - #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=event.state_component[0]), - time()-t_last_step) - if step % 1 == 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() - - -if __name__ == "__main__": - main() diff --git a/grudge/models/em.py b/grudge/models/em.py index ff491bef..8241c3a2 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -288,9 +288,9 @@ class MaxwellOperator(HyperbolicOperator): - self.local_derivatives(w) + sym.InverseMassOperator()(sym.FaceMassOperator()( flux(sym.int_tpair(w)) - #+ sum( - #flux(sym.bv_tpair(tag, w, bc)) - #for tag, bc in tags_and_bcs) + + sum( + flux(sym.bv_tpair(tag, w, bc)) + for tag, bc in tags_and_bcs) ) )) / material_divisor -- GitLab From fc42169d7c8d37b3ff3e2d86eb0da2756216c975 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 18 Jul 2017 14:18:40 -0500 Subject: [PATCH 17/34] Add constant sizing mapper --- grudge/execution.py | 36 ++++++++++++++++++++--------- grudge/symbolic/mappers/__init__.py | 7 ++++++ 2 files changed, 32 insertions(+), 11 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index b198c076..e9baa5d6 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -530,27 +530,41 @@ class BoundOperator(object): # {{{ process_sym_operator function -def process_sym_operator(sym_operator, post_bind_mapper=None, - dumper=lambda name, sym_operator: None, mesh=None): +def process_sym_operator(discr, sym_operator, post_bind_mapper=None, + dumper=lambda name, sym_operator: None): import grudge.symbolic.mappers as mappers dumper("before-bind", sym_operator) sym_operator = mappers.OperatorBinder()(sym_operator) - mappers.ErrorChecker(mesh)(sym_operator) + mappers.ErrorChecker(discr.mesh)(sym_operator) if post_bind_mapper is not None: dumper("before-postbind", sym_operator) sym_operator = post_bind_mapper(sym_operator) - if mesh is not None: - dumper("before-empty-flux-killer", sym_operator) - sym_operator = mappers.EmptyFluxKiller(mesh)(sym_operator) + dumper("before-empty-flux-killer", sym_operator) + sym_operator = mappers.EmptyFluxKiller(discr.volume_discr.mesh)(sym_operator) dumper("before-cfold", sym_operator) sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator) + # Work around https://github.com/numpy/numpy/issues/9438 + # + # The idea is that we need 1j as an expression to survive + # until code generation time. If it is evaluated and combined + # with other constants, we will need to determine its size + # (as np.complex64/128) within the expression. But because + # of the above numpy bug, sized numbers are not likely to survive + # expression building--so that's why we step in here to fix that. + + dumper("before-csize", sym_operator) + sym_operator = mappers.ConstantToNumpyConversionMapper( + real_type=discr.volume_discr.real_dtype.type, + complex_type=discr.volume_discr.complex_dtype.type, + )(sym_operator) + # Ordering restriction: # # - Must run constant fold before first type inference pass, because zeros, @@ -571,9 +585,8 @@ def process_sym_operator(sym_operator, post_bind_mapper=None, # because otherwise it won't differentiate the type of grids (node or quadrature # grids) that the operators will apply on. - assert mesh is not None dumper("before-global-to-reference", sym_operator) - sym_operator = mappers.GlobalToReferenceMapper(mesh.ambient_dim)(sym_operator) + sym_operator = mappers.GlobalToReferenceMapper(discr.ambient_dim)(sym_operator) # Ordering restriction: # @@ -614,10 +627,11 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x, pretty(sym_operator)) stage[0] += 1 - sym_operator = process_sym_operator(sym_operator, + sym_operator = process_sym_operator( + discr, + sym_operator, post_bind_mapper=post_bind_mapper, - dumper=dump_optemplate, - mesh=discr.mesh) + dumper=dump_optemplate) from grudge.symbolic.compiler import OperatorCompiler discr_code, eval_code = OperatorCompiler(discr)(sym_operator) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 60b489ce..9eb8e2b3 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): -- GitLab From a4b16e3c53f753ed8d0544ad5a0f6b231069acec Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 18 Jul 2017 14:18:59 -0500 Subject: [PATCH 18/34] Clean up EM header, add Bogdan to Copyright --- grudge/models/em.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/grudge/models/em.py b/grudge/models/em.py index 8241c3a2..7b0def8b 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -1,11 +1,13 @@ # -*- coding: utf8 -*- """grudge operators modelling electromagnetic phenomena.""" -from __future__ import division -from __future__ import absolute_import -from six.moves import range +from __future__ import division, absolute_import -__copyright__ = "Copyright (C) 2007, 2010 Andreas Kloeckner, David Powell" +__copyright__ = """ +Copyright (C) 2007-2017 Andreas Kloeckner +Copyright (C) 2010 David Powell +Copyright (C) 2017 Bogdan Enache +""" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -27,6 +29,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from six.moves import range + from pytools import memoize_method from grudge.models import HyperbolicOperator -- GitLab From e255e0e8ba138720e76934d7baf4bea63de13aae Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 18 Jul 2017 14:20:01 -0500 Subject: [PATCH 19/34] Fix scalar expr eval, plus increased Flake8 happiness --- grudge/execution.py | 50 ++++++++++++++++++++++----------------------- 1 file changed, 24 insertions(+), 26 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index e9baa5d6..87188134 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -85,46 +85,46 @@ class ExecutionMapper(mappers.Evaluator, return self.context[expr.name] def map_call(self, expr): - from pymbolic.primitives import Variable as var - import pymbolic.primitives as p - import pymbolic.functions as f - assert isinstance(expr.function, var) + from pymbolic.primitives import Variable + assert isinstance(expr.function, Variable) # FIXME: Make a way to register functions - pars = [self.rec(p) for p in expr.parameters] - if pars[0].shape == (): - import pyopencl.clmath as clmath - func = getattr(clmath, expr.function.name) - return func(*pars) + 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 = var("i") - func = var(expr.function.name) - if expr.function.name == "fabs": #FIXME - func = var("abs") + 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 Date: Tue, 18 Jul 2017 14:20:36 -0500 Subject: [PATCH 20/34] Make to-loopy conversion understand 'fabs' --- grudge/symbolic/compiler.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index ac7765ed..f6e23808 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( -- GitLab From d8ab43f85d78320916663e16310379e5b7b5193a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 18 Jul 2017 14:21:04 -0500 Subject: [PATCH 21/34] Clean up Maxwell cavities solution/example --- examples/maxwell/analytic_solutions.py | 306 ++++++++----------------- examples/maxwell/cavities.py | 30 ++- 2 files changed, 105 insertions(+), 231 deletions(-) diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py index 73c3e901..22fe4394 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,169 +216,67 @@ 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 - - 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 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. - def __call__(self, discr): - f,g,h = self.factors - omega = self.omega - k = self.k - nodes = discr.nodes() - x = nodes[:,0] - y = nodes[:,1] - if discr.dim > 2: - z = nodes[:,2] - else: - z = discr.volume_zeros() - - sx = numpy.sin(f*x) - sy = numpy.sin(g*y) - cx = numpy.cos(f*x) - cy = numpy.cos(g*y) - - zdep_add = numpy.exp(1j*h*z)+numpy.exp(-1j*h*z) - zdep_sub = numpy.exp(1j*h*z)-numpy.exp(-1j*h*z) - - tdep = numpy.exp(-1j * omega * self.t) - - C = 1j/(f**2+g**2) - - result = numpy.zeros((6,len(x)), dtype=zdep_add.dtype) - - result[0] = C*f*h*cx*sy*zdep_sub*tdep # ex - result[1] = C*g*h*sx*cy*zdep_sub*tdep # ey - result[2] = sx*sy*zdep_add*tdep # ez - result[3] = -C*g*self.epsilon*omega*sx*cy*zdep_add*tdep # hx - result[4] = C*f*self.epsilon*omega*cx*sy*zdep_add*tdep - - return result - - - - -class RectangularCavityMode(RectangularWaveguideMode): - """A rectangular TM cavity mode.""" - - def __init__(self, *args, **kwargs): - if "scale" in kwargs: - kwargs["forward_coeff"] = scale - kwargs["backward_coeff"] = scale - else: - kwargs["forward_coeff"] = 1 - kwargs["backward_coeff"] = 1 - RectangularWaveguideMode.__init__(self, *args, **kwargs) - -class SymRectangularWaveguideMode: - """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, nodes): - f,g,h = self.factors - omega = self.omega - k = self.k + Returns an expression depending on *epsilon*, *mu*, and *t*. + """ - x = nodes[0] - y = nodes[1] - if nodes.shape[0] > 2: - z = nodes[2] - else: - z = 0 - from grudge import sym + f, g, h = factors = [n*np.pi/a for n, a in zip(mode_indices, dimensions)] + k = np.sqrt(sum(f**2 for f in factors)) - sx = sym.sin(f*x) - sy = sym.sin(g*y) - cx = sym.cos(f*x) - cy = sym.cos(g*y) + epsilon = sym.ScalarVariable("epsilon") + mu = sym.ScalarVariable("mu") - zdep_add = sym.exp(1j*h*z)+sym.exp(-1j*h*z) - zdep_sub = sym.exp(1j*h*z)-sym.exp(-1j*h*z) + c = sym.cse(1/sym.sqrt(mu*epsilon), "c") + omega = k*c - tdep = numpy.exp(-1j * omega * self.t) + nodes = sym.nodes(ambient_dim) + x = nodes[0] + y = nodes[1] + if nodes.shape[0] > 2: + z = nodes[2] + else: + z = 0 - C = 1j/(f**2+g**2) + sx = sym.sin(f*x) + sy = sym.sin(g*y) + cx = sym.cos(f*x) + cy = sym.cos(g*y) - from pytools.obj_array import join_fields - result = join_fields( C*f*h*cx*sy*zdep_sub*tdep # ex - , C*g*h*sx*cy*zdep_sub*tdep # ey - , sx*sy*zdep_add*tdep # ez - , -C*g*self.epsilon*omega*sx*cy*zdep_add*tdep # hx - , C*f*self.epsilon*omega*cx*sy*zdep_add*tdep - , C*f*self.epsilon*omega*cx*sy*zdep_add*tdep - ) # THIS IS BAD, last line should be all 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) - return result + tdep = sym.exp(-1j * omega * sym.ScalarVariable("t")) + const = 1j/(f**2+g**2) + result = sym.join_fields( + const*f*h*cx*sy*zdep_sub*tdep, # ex + const*g*h*sx*cy*zdep_sub*tdep, # ey + sx*sy*zdep_add*tdep, # ez + -const*g*epsilon*omega*sx*cy*zdep_add*tdep, # hx + const*f*epsilon*omega*cx*sy*zdep_add*tdep, + const*f*epsilon*omega*cx*sy*zdep_add*tdep, + ) # THIS IS BAD, last line should be all zeros -class SymRectangularCavityMode(SymRectangularWaveguideMode): - """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) + return result +# }}} +# {{{ dipole solution -# dipole solution ------------------------------------------------------------- class DipoleFarField: """Dipole EM field. @@ -462,16 +339,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) @@ -510,15 +385,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 ) @@ -536,3 +411,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 2a756fa6..3f9dec59 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -27,15 +27,16 @@ THE SOFTWARE. import numpy as np import pyopencl as cl -from grudge.shortcuts import set_up_rk4 -from grudge.shortcuts import set_up_forward_euler +from grudge.shortcuts import set_up_rk4, set_up_forward_euler # noqa: F401 from grudge import sym, bind, Discretization from analytic_solutions import ( # noqa - RectangularWaveguideMode, - RectangularCavityMode, - SymRectangularWaveguideMode, - SymRectangularCavityMode) + get_rectangular_waveguide_mode, + #RectangularCavityMode, + #RectangularWaveguideMode, + #RectangularCavityMode + ) + def main(write_output=True, order=4): cl_ctx = cl.create_some_context() @@ -50,24 +51,19 @@ def main(write_output=True, order=4): discr = Discretization(cl_ctx, mesh, order=order) - epsilon0 = 8.8541878176e-12 # C**2 / (N m**2) mu0 = 4*np.pi*1e-7 # N/A**2. epsilon = 1*epsilon0 mu = 1*mu0 - + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.models.em import MaxwellOperator - from meshmode.mesh import BTAG_ALL, BTAG_NONE - op = MaxwellOperator(epsilon,mu, - flux_type="lf") - + op = MaxwellOperator(epsilon, mu, flux_type="lf") queue = cl.CommandQueue(discr.cl_context) - from pytools.obj_array import join_fields - mode = SymRectangularCavityMode(epsilon, mu, (1, 2, 2)) - fields = bind(discr, mode(sym.nodes(dims)))(queue, t=0) + sym_mode = get_rectangular_waveguide_mode(mesh.ambient_dim, (1, 2, 2)) + fields = bind(discr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu) #fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr.volume_discr) #fields = (discr.volume_discr) @@ -109,7 +105,7 @@ def main(write_output=True, order=4): from time import time t_last_step = time() - e,h = op.split_eh(fields) + e, h = op.split_eh(fields) #vis.write_vtk_file("fld-%04d.vtu" % step, #[ #("e", e), @@ -125,7 +121,7 @@ def main(write_output=True, order=4): print(step, event.t, norm(queue, u=event.state_component[0]), time()-t_last_step) if step % 1 == 0: - e,h = op.split_eh(event.state_component) + e, h = op.split_eh(event.state_component) #vis.write_vtk_file("fld-%04d.vtu" % step, #[ #("e", e), -- GitLab From 8823bfb306d2d2825123202eaae226eaeb678e19 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 26 Jul 2017 16:35:41 -0500 Subject: [PATCH 22/34] Voodoo EM flux sign flip --- grudge/models/em.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/models/em.py b/grudge/models/em.py index 7b0def8b..f50f393a 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -290,7 +290,7 @@ class MaxwellOperator(HyperbolicOperator): return ( - self.local_derivatives(w) - + sym.InverseMassOperator()(sym.FaceMassOperator()( + - sym.InverseMassOperator()(sym.FaceMassOperator()( flux(sym.int_tpair(w)) + sum( flux(sym.bv_tpair(tag, w, bc)) -- GitLab From c21968570d5877f1f85a989b5f1b5139a8c0e0ef Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 26 Jul 2017 16:36:15 -0500 Subject: [PATCH 23/34] Cavities test debugging hackery --- examples/maxwell/cavities.py | 48 ++++++++++++++++++++---------------- 1 file changed, 27 insertions(+), 21 deletions(-) diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index 3f9dec59..6f9f9819 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -42,23 +42,27 @@ def main(write_output=True, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - dims = 3 + 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=(4,)*dims) + n=(20,)*dims) discr = Discretization(cl_ctx, mesh, order=order) - epsilon0 = 8.8541878176e-12 # C**2 / (N m**2) - mu0 = 4*np.pi*1e-7 # N/A**2. - epsilon = 1*epsilon0 - mu = 1*mu0 + 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 MaxwellOperator - op = MaxwellOperator(epsilon, mu, flux_type="lf") + op = MaxwellOperator(epsilon, mu, flux_type="lf", dimensions=dims) queue = cl.CommandQueue(discr.cl_context) @@ -84,13 +88,13 @@ def main(write_output=True, order=4): return bound_op(queue, t=t, w=w) if mesh.dim == 2: - dt = 0.00004 + dt = 0.0004 elif mesh.dim == 3: - dt = 0.00002 + dt = 0.0002 dt_stepper = set_up_forward_euler("w", dt, fields, rhs) - final_t = dt * 0.9 + final_t = dt * 300 nsteps = int(final_t/dt) print("dt=%g nsteps=%d" % (dt, nsteps)) @@ -106,11 +110,13 @@ def main(write_output=True, order=4): t_last_step = time() e, h = op.split_eh(fields) - #vis.write_vtk_file("fld-%04d.vtu" % step, - #[ - #("e", e), - #("h", h), - #]) + + if 1: + 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): @@ -120,13 +126,13 @@ def main(write_output=True, order=4): print(step, event.t, norm(queue, u=event.state_component[0]), time()-t_last_step) - if step % 1 == 0: + if step % 10 == 0: e, h = op.split_eh(event.state_component) - #vis.write_vtk_file("fld-%04d.vtu" % step, - #[ - #("e", e), - #("h", h), - #]) + vis.write_vtk_file("fld-%04d.vtu" % step, + [ + ("e", e), + ("h", h), + ]) t_last_step = time() -- GitLab From f120cc13042c50f1ae6f91e00b5499b47f915f8f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 15 Aug 2017 16:54:57 -0500 Subject: [PATCH 24/34] Allow elwise functions to take arrays with offsets --- grudge/execution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/execution.py b/grudge/execution.py index 95200d87..c4f285cd 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -143,7 +143,7 @@ class ExecutionMapper(mappers.Evaluator, [ lp.Assignment(Variable("out")[i], func(Variable("a")[i])) - ]) + ], default_offset=lp.auto) return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0") assert len(args) == 1 -- GitLab From eb42e7fe68b8b0246fe1deb5e6c7b5b697fc2f8b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 15 Aug 2017 17:30:24 -0500 Subject: [PATCH 25/34] Flake8 fixes to pytools --- grudge/tools.py | 36 +++++++++++------------------------- 1 file changed, 11 insertions(+), 25 deletions(-) diff --git a/grudge/tools.py b/grudge/tools.py index cf5aea0b..88d67143 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): @@ -37,20 +38,6 @@ def is_zero(x): else: return False -def levi_civita(tuple): - """Compute an entry of the Levi-Civita tensor for the indices *tuple*. - Only three-tuples are supported for now. - """ - if len(tuple) == 3: - if tuple in [(0,1,2), (2,0,1), (1,2,0)]: - return 1 - elif tuple in [(2,1,0), (0,2,1), (1,0,2)]: - return -1 - else: - return 0 - else: - raise NotImplementedError - class SubsettableCrossProduct: """A cross product that can operate on an arbitrary subsets of its @@ -59,14 +46,15 @@ class SubsettableCrossProduct: full_subset = (True, True, True) - def __init__(self, op1_subset=full_subset, op2_subset=full_subset, result_subset=full_subset): + 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. + :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) @@ -111,14 +99,14 @@ class SubsettableCrossProduct: 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. @@ -136,5 +124,3 @@ def partial_to_all_subset_indices(subsets, base=0): idx += 1 yield np.array(result, dtype=np.intp) - - -- GitLab From 8b4c78bc3cf8d6b90b14a00efa0786c3194e3784 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 15 Aug 2017 17:52:51 -0500 Subject: [PATCH 26/34] Flake8 fixes to shortcuts --- grudge/shortcuts.py | 1 + 1 file changed, 1 insertion(+) diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py index 34e8f39d..84fd8b3f 100644 --- a/grudge/shortcuts.py +++ b/grudge/shortcuts.py @@ -43,6 +43,7 @@ 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 -- GitLab From 2e77ccc3ab6f7651919d18cb5bd07094676f690c Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 28 Aug 2017 03:54:43 -0700 Subject: [PATCH 27/34] Attempted to add rectangular cavity mode --- examples/maxwell/analytic_solutions.py | 72 +++++++++++++++++++++----- examples/maxwell/cavities.py | 7 +-- 2 files changed, 64 insertions(+), 15 deletions(-) diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py index 22fe4394..5458feb0 100644 --- a/examples/maxwell/analytic_solutions.py +++ b/examples/maxwell/analytic_solutions.py @@ -231,7 +231,7 @@ def get_rectangular_waveguide_mode( Returns an expression depending on *epsilon*, *mu*, and *t*. """ - f, g, h = factors = [n*np.pi/a for n, a in zip(mode_indices, dimensions)] + 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)) epsilon = sym.ScalarVariable("epsilon") @@ -248,35 +248,83 @@ def get_rectangular_waveguide_mode( else: z = 0 - sx = sym.sin(f*x) - sy = sym.sin(g*y) - cx = sym.cos(f*x) - cy = sym.cos(g*y) + sx = sym.sin(kx*x) + sy = sym.sin(ky*y) + cx = sym.cos(kx*x) + cy = sym.cos(ky*y) zdep_add = sym.exp(1j*h*z)+sym.exp(-1j*h*z) zdep_sub = sym.exp(1j*h*z)-sym.exp(-1j*h*z) tdep = sym.exp(-1j * omega * sym.ScalarVariable("t")) - const = 1j/(f**2+g**2) + const = 1j/(kx**2+ky**2) result = sym.join_fields( - const*f*h*cx*sy*zdep_sub*tdep, # ex - const*g*h*sx*cy*zdep_sub*tdep, # ey + const*kx*h*cx*sy*zdep_sub*tdep, # ex + const*ky*h*sx*cy*zdep_sub*tdep, # ey sx*sy*zdep_add*tdep, # ez - -const*g*epsilon*omega*sx*cy*zdep_add*tdep, # hx - const*f*epsilon*omega*cx*sy*zdep_add*tdep, - const*f*epsilon*omega*cx*sy*zdep_add*tdep, - ) # THIS IS BAD, last line should be all zeros + -const*ky*epsilon*omega*sx*cy*zdep_add*tdep, # hx + const*kx*epsilon*omega*cx*sy*zdep_add*tdep, + 0, + ) return result # }}} +def get_rectangular_cavity_mode(TildeHatA, + 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 + + 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")) + + result = sym.join_fields( + TildeHatA*cx*sy*sz*tdep, # ex + TildeHatA*sx*cy*cz*tdep * ky / kx, # ey + -TildeHatA * (kx/kz + ky * ky / (kx*kz)) * sx*sy*cz*tdep, # ez + + 1j*TildeHatA*((kx + ky*ky/kx)*(ky/kz) + ky*kz/kx)*sx*cy*cz*tdep/omega, # hx + -1j*TildeHatA*(kz + (kx + ky*ky/kx))*cx*sy*cz*tdep/omega, + 0, + ) + + return result # {{{ dipole solution + class DipoleFarField: """Dipole EM field. diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index 6f9f9819..0bcfbe40 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -31,7 +31,7 @@ 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_waveguide_mode, + get_rectangular_cavity_mode, #RectangularCavityMode, #RectangularWaveguideMode, #RectangularCavityMode @@ -66,7 +66,7 @@ def main(write_output=True, order=4): queue = cl.CommandQueue(discr.cl_context) - sym_mode = get_rectangular_waveguide_mode(mesh.ambient_dim, (1, 2, 2)) + 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) #fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr.volume_discr) @@ -124,7 +124,8 @@ def main(write_output=True, order=4): step += 1 - print(step, event.t, norm(queue, u=event.state_component[0]), + 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) -- GitLab From 6f66b7ff46bfe49f33ce5cc6109c88866f0f4b03 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 22 Sep 2017 17:01:18 -0500 Subject: [PATCH 28/34] Implement PointsDiscretization abomination --- grudge/discretization.py | 33 ++++++++++++++++++++++++++++++++- test/test_grudge.py | 16 ++++++++++++++++ 2 files changed, 48 insertions(+), 1 deletion(-) diff --git a/grudge/discretization.py b/grudge/discretization.py index 8cbacd89..50e112de 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -27,7 +27,11 @@ from pytools import memoize_method from grudge import sym -class Discretization(object): +class DiscretizationBase(object): + pass + + +class Discretization(DiscretizationBase): """ .. attribute :: volume_discr @@ -220,4 +224,31 @@ 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 + + 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/test/test_grudge.py b/test/test_grudge.py index 0ceaf8fe..bd47f823 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()' -- GitLab From ec76b11fd8cbfc86884a52aecb5859dacb602d8b Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Thu, 28 Sep 2017 23:55:49 -0500 Subject: [PATCH 29/34] Progress with sumpy analytic solution check --- examples/maxwell/analytic_solutions.py | 58 ++++++ examples/maxwell/cavities.py | 25 ++- examples/maxwell/testing.py | 249 ++++++++++++++++++++++--- grudge/discretization.py | 6 + 4 files changed, 294 insertions(+), 44 deletions(-) diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py index 5458feb0..14a962f0 100644 --- a/examples/maxwell/analytic_solutions.py +++ b/examples/maxwell/analytic_solutions.py @@ -322,6 +322,64 @@ def get_rectangular_cavity_mode(TildeHatA, 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] + + 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 diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index 0bcfbe40..628ddb9f 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -32,6 +32,7 @@ from grudge import sym, bind, Discretization from analytic_solutions import ( # noqa get_rectangular_cavity_mode, + get_rectangular_2D_cavity_mode, #RectangularCavityMode, #RectangularWaveguideMode, #RectangularCavityMode @@ -45,8 +46,8 @@ def main(write_output=True, order=4): dims = 2 from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh( - a=(-0.5,)*dims, - b=(0.5,)*dims, + a=(0.0,)*dims, + b=(1.0,)*dims, n=(20,)*dims) discr = Discretization(cl_ctx, mesh, order=order) @@ -62,17 +63,13 @@ def main(write_output=True, order=4): from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.models.em import MaxwellOperator - op = MaxwellOperator(epsilon, mu, flux_type="lf", dimensions=dims) + 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)) - fields = bind(discr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu) - #fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr.volume_discr) - - #fields = (discr.volume_discr) - #fields = join_fields(discr.zeros(queue), - #[discr.zeros(queue) for i in range(discr.dim*2 - 1)]) + #sym_mode = get_rectangular_cavity_mode(1,mesh.ambient_dim, (1, 2, 2)) + sym_mode = get_rectangular_2D_cavity_mode(2,3) + fields = bind(discr, sym_mode)(queue, t=0) # FIXME #dt = op.estimate_rk4_timestep(discr, fields=fields) @@ -88,13 +85,13 @@ def main(write_output=True, order=4): return bound_op(queue, t=t, w=w) if mesh.dim == 2: - dt = 0.0004 + dt = 0.004 elif mesh.dim == 3: - dt = 0.0002 + dt = 0.002 - dt_stepper = set_up_forward_euler("w", dt, fields, rhs) + dt_stepper = set_up_rk4("w", dt, fields, rhs) - final_t = dt * 300 + final_t = dt * 600 nsteps = int(final_t/dt) print("dt=%g nsteps=%d" % (dt, nsteps)) diff --git a/examples/maxwell/testing.py b/examples/maxwell/testing.py index 96de8f24..b69b3410 100644 --- a/examples/maxwell/testing.py +++ b/examples/maxwell/testing.py @@ -27,56 +27,245 @@ THE SOFTWARE. import numpy as np import pyopencl as cl -from grudge.shortcuts import set_up_rk4 +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 - RectangularWaveguideMode, - RectangularCavityMode, - SymRectangularWaveguideMode, - SymRectangularCavityMode) + get_rectangular_cavity_mode, + get_rectangular_2D_cavity_mode, + get_rectangular_2D_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(write_output=True, order=4): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dims = 2 + + 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 + + for ns in [8,16, 32]: + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(0.0,)*dims, + b=(1.0,)*dims, + n=(ns,)*dims) + + discr = Discretization(cl_ctx, mesh, order=order) + + 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) + + 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) + fields = bind(discr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu) + + sym_mode_d = get_rectangular_2D_cavity_mode_deriv(1,1) + fields_d = bind(discr, sym_mode_d)(queue, t=0, epsilon=epsilon, mu=mu) + + + from grudge.tools import count_subset + w = sym.make_sym_array("w", count_subset(op.get_eh_subset())) + bound_op = bind(discr, op.local_derivatives(w)) + + dt = 0.0002 + #print(fields[2]) + + + + 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(write_output=True, order=4): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dims = 2 + + 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.5,0.5)) + + + #discr = Discretization(cl_ctx, mesh, order=order) + + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.models.em import MaxwellOperator - from meshmode.mesh import BTAG_ALL, BTAG_NONE - + op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims) queue = cl.CommandQueue(cl_ctx) - from pytools.obj_array import join_fields - from pymbolic.primitives import Variable as var - import pymbolic.primitives as p - import loopy as lp + from grudge.discretization import PointsDiscretization + pdiscr = PointsDiscretization(cl.array.to_device(queue, patch.points)) + #sym_mode = get_rectangular_cavity_mode(1,mesh.ambient_dim, (1, 2, 2)) + sym_mode = get_rectangular_2D_cavity_mode(1,1) + fields = bind(pdiscr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu) - nom = "abs" - i = var("i") - func = var(nom) - if("fabs" == nom): - func = var("abs") + 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) - def knl(): - knl = lp.make_kernel( - "{[i]: 0<=i Date: Tue, 3 Oct 2017 15:27:29 -0500 Subject: [PATCH 30/34] Add missing .get() in analytic solution test --- examples/maxwell/testing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/maxwell/testing.py b/examples/maxwell/testing.py index b69b3410..9a1a6d28 100644 --- a/examples/maxwell/testing.py +++ b/examples/maxwell/testing.py @@ -248,7 +248,7 @@ def sumpy_test(write_output=True, order=4): #print(fields[2]) - print( patch.dy(fields[2])) + print( patch.dy(fields[2].get())) print( fields_d[2]) -- GitLab From 42cebd7f7502e47e69b12b79f13e5135837b556c Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 11 Oct 2017 16:32:38 -0500 Subject: [PATCH 31/34] Attempt to change 3D analytics --- examples/maxwell/analytic_solutions.py | 59 ++++++++++++++++++++++++++ examples/maxwell/cavities.py | 11 +++-- examples/maxwell/testing.py | 41 +++++++++++++++--- 3 files changed, 101 insertions(+), 10 deletions(-) diff --git a/examples/maxwell/analytic_solutions.py b/examples/maxwell/analytic_solutions.py index 14a962f0..637f0ac7 100644 --- a/examples/maxwell/analytic_solutions.py +++ b/examples/maxwell/analytic_solutions.py @@ -309,6 +309,7 @@ def get_rectangular_cavity_mode(TildeHatA, cz = sym.cos(kz*z) tdep = sym.exp(-1j * omega * sym.ScalarVariable("t")) + #tdep = 1 result = sym.join_fields( TildeHatA*cx*sy*sz*tdep, # ex @@ -320,6 +321,64 @@ def get_rectangular_cavity_mode(TildeHatA, 0, ) + #result = sym.join_fields( + #-kx * kz * TildeHatA*cx*sy*sz*tdep / (k**2), # ex + #-ky * kz * TildeHatA*sx*cy*sz*tdep / (k**2), # ey + #TildeHatA * sx*sy*cz*tdep, # ez + + #1j * omega * ky*TildeHatA*sx*cy*cz*tdep / (k**2), # hx + #-1j * omega * kx*TildeHatA*cx*sy*cz*tdep / (k**2), + #0, + #) + + return result + +def get_rectangular_3D_cavity_mode_deriv(TildeHatA, + 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 + + 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")) + + result = sym.join_fields( + -1j * omega * TildeHatA*cx*sy*sz*tdep, # ex + -1j * omega * TildeHatA*sx*cy*cz*tdep * ky / kx, # ey + 1j * omega * TildeHatA * (kx/kz + ky * ky / (kx*kz)) * sx*sy*cz*tdep, # ez + + TildeHatA*((kx + ky*ky/kx)*(ky/kz) + ky*kz/kx)*sx*cy*cz*tdep, # hx + -TildeHatA*(kz + (kx + ky*ky/kx))*cx*sy*cz*tdep, + 0, + ) + return result def get_rectangular_2D_cavity_mode(m,n ): diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index 628ddb9f..ad2a4478 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -43,12 +43,12 @@ def main(write_output=True, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - dims = 2 + 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=(20,)*dims) + n=(5,)*dims) discr = Discretization(cl_ctx, mesh, order=order) @@ -67,8 +67,11 @@ def main(write_output=True, order=4): 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(2,3) + if dims == 3: + sym_mode = get_rectangular_cavity_mode(1,mesh.ambient_dim, (1, 2, 2)) + else: + sym_mode = get_rectangular_2D_cavity_mode(2,3) + fields = bind(discr, sym_mode)(queue, t=0) # FIXME diff --git a/examples/maxwell/testing.py b/examples/maxwell/testing.py index 9a1a6d28..dfa04a1a 100644 --- a/examples/maxwell/testing.py +++ b/examples/maxwell/testing.py @@ -35,6 +35,7 @@ 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 @@ -216,7 +217,7 @@ def sumpy_test(write_output=True, order=4): c = 1/np.sqrt(mu*epsilon) omega = k*c - patch = spy.CalculusPatch((0.5,0.5)) + patch = spy.CalculusPatch((0.4,0.3, 0.4)) #discr = Discretization(cl_ctx, mesh, order=order) @@ -230,11 +231,12 @@ def sumpy_test(write_output=True, order=4): from grudge.discretization import PointsDiscretization pdiscr = PointsDiscretization(cl.array.to_device(queue, patch.points)) - #sym_mode = get_rectangular_cavity_mode(1,mesh.ambient_dim, (1, 2, 2)) - sym_mode = get_rectangular_2D_cavity_mode(1,1) + sym_mode = get_rectangular_cavity_mode(1,3, (1, 2, 2)) + #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) + #sym_mode_d = get_rectangular_2D_cavity_mode_deriv(1,1) + 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) #fields[0] = np.zeros(fields[3].shape) @@ -247,10 +249,21 @@ def sumpy_test(write_output=True, order=4): #bound_op = bind(discr, op.local_derivatives(w)) #print(fields[2]) + 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( patch.dy(fields[2].get())) - print( fields_d[2]) + # 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())) + # 3D + #print( fields_d[3].get() + patch.dy(fields[2].get()) - patch.dz(fields[1].get())) + 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) @@ -263,6 +276,22 @@ def sumpy_test(write_output=True, order=4): #print(norm(queue, u=(-1j * omega* fields[5] - bound_op(queue, w=fields)[5]))/norm(queue, u =(-1j * omega * fields[5]))) 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(False,4) sumpy_test(False,4) -- GitLab From 58d01782775b82f76c70402030679ae450a50acb Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 25 Oct 2017 03:12:30 -0500 Subject: [PATCH 32/34] Attempt to come up with an analytic solution to maxwell's cube resonator --- contrib/maxima/em-modes.mac | 157 ++++++++++++++++++++++++++++++++++++ contrib/maxima/maxwell2.mac | 84 +++++++++++++++++++ 2 files changed, 241 insertions(+) create mode 100644 contrib/maxima/em-modes.mac create mode 100644 contrib/maxima/maxwell2.mac diff --git a/contrib/maxima/em-modes.mac b/contrib/maxima/em-modes.mac new file mode 100644 index 00000000..6482ec17 --- /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 00000000..5c4dcba4 --- /dev/null +++ b/contrib/maxima/maxwell2.mac @@ -0,0 +1,84 @@ +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); +*/ +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, k_x = m * %pi, k_y = n*%pi, k_z = p*%pi; +H : H, k_x = m * %pi, k_y = n*%pi, k_z = p*%pi; +*/ + +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])); + -- GitLab From 88095cedf31143e7d40a4d30b5dce874fe8d2e83 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sat, 28 Oct 2017 12:05:46 -0500 Subject: [PATCH 33/34] Got maxwells to converge --- contrib/maxima/maxwell2.mac | 10 +- examples/maxwell/analytic_solutions.py | 75 ++++++++---- examples/maxwell/cavities.py | 3 +- examples/maxwell/testing.py | 160 +++++++++++++++---------- 4 files changed, 164 insertions(+), 84 deletions(-) diff --git a/contrib/maxima/maxwell2.mac b/contrib/maxima/maxwell2.mac index 5c4dcba4..13b6f9dd 100644 --- a/contrib/maxima/maxwell2.mac +++ b/contrib/maxima/maxwell2.mac @@ -42,6 +42,9 @@ crossprod(a,b):=crossfunc(lambda([j,k], a[j]*b[k])); /* 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); @@ -61,14 +64,16 @@ 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, k_x = m * %pi, k_y = n*%pi, k_z = p*%pi; -H : H, k_x = m * %pi, k_y = n*%pi, k_z = p*%pi; +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])); @@ -81,4 +86,5 @@ 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 637f0ac7..2bf0ee46 100644 --- a/examples/maxwell/analytic_solutions.py +++ b/examples/maxwell/analytic_solutions.py @@ -274,7 +274,7 @@ def get_rectangular_waveguide_mode( # }}} -def get_rectangular_cavity_mode(TildeHatA, +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, @@ -292,6 +292,7 @@ def get_rectangular_cavity_mode(TildeHatA, 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] @@ -311,29 +312,39 @@ def get_rectangular_cavity_mode(TildeHatA, tdep = sym.exp(-1j * omega * sym.ScalarVariable("t")) #tdep = 1 - result = sym.join_fields( - TildeHatA*cx*sy*sz*tdep, # ex - TildeHatA*sx*cy*cz*tdep * ky / kx, # ey - -TildeHatA * (kx/kz + ky * ky / (kx*kz)) * sx*sy*cz*tdep, # ez + #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*TildeHatA*((kx + ky*ky/kx)*(ky/kz) + ky*kz/kx)*sx*cy*cz*tdep/omega, # hx - -1j*TildeHatA*(kz + (kx + ky*ky/kx))*cx*sy*cz*tdep/omega, - 0, - ) + #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 * TildeHatA*cx*sy*sz*tdep / (k**2), # ex - #-ky * kz * TildeHatA*sx*cy*sz*tdep / (k**2), # ey - #TildeHatA * sx*sy*cz*tdep, # ez + #-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*TildeHatA*sx*cy*cz*tdep / (k**2), # hx - #-1j * omega * kx*TildeHatA*cx*sy*cz*tdep / (k**2), + #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(TildeHatA, +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, @@ -342,7 +353,6 @@ def get_rectangular_3D_cavity_mode_deriv(TildeHatA, 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)) @@ -351,6 +361,7 @@ def get_rectangular_3D_cavity_mode_deriv(TildeHatA, 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] @@ -368,19 +379,41 @@ def get_rectangular_3D_cavity_mode_deriv(TildeHatA, 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 * TildeHatA*cx*sy*sz*tdep, # ex - -1j * omega * TildeHatA*sx*cy*cz*tdep * ky / kx, # ey - 1j * omega * TildeHatA * (kx/kz + ky * ky / (kx*kz)) * sx*sy*cz*tdep, # ez + 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 - TildeHatA*((kx + ky*ky/kx)*(ky/kz) + ky*kz/kx)*sx*cy*cz*tdep, # hx - -TildeHatA*(kz + (kx + ky*ky/kx))*cx*sy*cz*tdep, + -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. diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index ad2a4478..8877fab8 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -69,10 +69,11 @@ def main(write_output=True, order=4): 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) - fields = bind(discr, sym_mode)(queue, t=0) # FIXME #dt = op.estimate_rk4_timestep(discr, fields=fields) diff --git a/examples/maxwell/testing.py b/examples/maxwell/testing.py index dfa04a1a..09f96058 100644 --- a/examples/maxwell/testing.py +++ b/examples/maxwell/testing.py @@ -145,68 +145,93 @@ def main(write_output=True, order=4): print(norm(queue, u=( event.state_component[5] - sol[5]))/ norm(queue, u=( sol[5]))) return -def analytic_test(write_output=True, order=4): +def analytic_test(dims = 3, n = 8, 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=(n,)*dims) - epsilon = 1 - mu = 1 + discr = Discretization(cl_ctx, mesh, order=order) - 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)) + 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 - c = 1/np.sqrt(mu*epsilon) - omega = k*c + 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) - for ns in [8,16, 32]: - from meshmode.mesh.generation import generate_regular_rect_mesh - mesh = generate_regular_rect_mesh( - a=(0.0,)*dims, - b=(1.0,)*dims, - n=(ns,)*dims) + queue = cl.CommandQueue(discr.cl_context) - discr = Discretization(cl_ctx, mesh, order=order) + 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)) - 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) + analytic_sol = bind(discr, sym_mode) + fields = analytic_sol(queue, t=0, epsilon=epsilon, mu=mu) - 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) - fields = bind(discr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu) + op.check_bc_coverage(mesh) - sym_mode_d = get_rectangular_2D_cavity_mode_deriv(1,1) - fields_d = bind(discr, sym_mode_d)(queue, t=0, epsilon=epsilon, mu=mu) + 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) - from grudge.tools import count_subset - w = sym.make_sym_array("w", count_subset(op.get_eh_subset())) - bound_op = bind(discr, op.local_derivatives(w)) + if mesh.dim == 2: + dt = 0.002 + elif mesh.dim == 3: + dt = 0.002 - dt = 0.0002 - #print(fields[2]) + 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"))) - 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(write_output=True, order=4): + 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 = 2 + dims = 3 epsilon = 1 mu = 1 @@ -220,35 +245,17 @@ def sumpy_test(write_output=True, order=4): patch = spy.CalculusPatch((0.4,0.3, 0.4)) - #discr = Discretization(cl_ctx, mesh, order=order) - - 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) - 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)) - #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) 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) - #fields[0] = np.zeros(fields[3].shape) - #fields[1] = np.zeros(fields[3].shape) - #fields[5] = np.zeros(fields[3].shape) - #fields = np.vstack(fields) - - #from grudge.tools import count_subset - #w = sym.make_sym_array("w", count_subset(op.get_eh_subset())) - #bound_op = bind(discr, op.local_derivatives(w)) - - #print(fields[2]) for i in range(len(fields)): if type(fields[i]) == type(0): fields[i] = np.zeros(patch.points.shape[1]) @@ -259,9 +266,11 @@ def sumpy_test(write_output=True, order=4): #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)) @@ -276,6 +285,36 @@ def sumpy_test(write_output=True, order=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 @@ -293,8 +332,9 @@ def frequency_domain_maxwell(cpatch, e, h, k): return ( resid_faraday, resid_ampere, resid_div_e, resid_div_h) if __name__ == "__main__": - #analytic_test(False,4) - sumpy_test(False,4) + analytic_test(3, 4, False,4) + analytic_test(3, 8, False,4) + #sumpy_test_3D(False,4) #main() -- GitLab From d322c001d7103ee57dbfe505413c94c73dad1be6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 29 Oct 2017 13:46:53 +0100 Subject: [PATCH 34/34] Placate new Flake8 --- grudge/discretization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/discretization.py b/grudge/discretization.py index e2f2cfe0..d578e636 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