From 3f40665f73903ae3a95fab4904c06e12e8672fe4 Mon Sep 17 00:00:00 2001 From: "[6~" Date: Thu, 21 May 2020 23:32:53 -0500 Subject: [PATCH 1/5] Add eagerly evaluated interface to grudge, plus example --- examples/wave/wave-eager.py | 182 ++++++++++++++++++++++++++++++++++++ grudge/eager.py | 128 +++++++++++++++++++++++++ 2 files changed, 310 insertions(+) create mode 100644 examples/wave/wave-eager.py create mode 100644 grudge/eager.py diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py new file mode 100644 index 00000000..2ec333fc --- /dev/null +++ b/examples/wave/wave-eager.py @@ -0,0 +1,182 @@ +from __future__ import division, print_function + +__copyright__ = "Copyright (C) 2020 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 numpy.linalg as la # noqa +import pyopencl as cl +import pyopencl.array as cla # noqa +import pyopencl.clmath as clmath +from pytools.obj_array import ( + join_fields, + with_object_array_or_scalar) +from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from grudge.eager import EagerDGDiscretization, with_queue +from grudge.symbolic.primitives import TracePair +from grudge.shortcuts import make_visualizer + + +def interior_trace_pair(discr, vec): + i = discr.interp("vol", "int_faces", vec) + e = with_object_array_or_scalar( + lambda el: discr.opposite_face_connection()(el.queue, el), + i) + return TracePair("int_faces", i, e) + + +# {{{ wave equation bits + +def wave_flux(discr, c, w_tpair): + u = w_tpair[0] + v = w_tpair[1:] + + normal = with_queue(u.int.queue, discr.normal(w_tpair.dd)) + + flux_weak = join_fields( + np.dot(v.avg, normal), + normal[0] * u.avg, + normal[1] * u.avg) + + # upwind + v_jump = np.dot(normal, v.int-v.ext) + flux_weak -= join_fields( + 0.5*(u.int-u.ext), + 0.5*normal[0]*v_jump, + 0.5*normal[1]*v_jump, + ) + + flux_strong = join_fields( + np.dot(v.int, normal), + u.int * normal[0], + u.int * normal[1], + ) - flux_weak + + return discr.interp(w_tpair.dd, "all_faces", c*flux_strong) + + +def wave_operator(discr, c, w): + u = w[0] + v = w[1:] + + dir_u = discr.interp("vol", BTAG_ALL, u) + dir_v = discr.interp("vol", BTAG_ALL, v) + dir_bval = join_fields(dir_u, dir_v) + dir_bc = join_fields(-dir_u, dir_v) + + return ( + - join_fields( + -c*discr.div(v), + -c*discr.grad(u) + ) + + # noqa: W504 + discr.inverse_mass( + discr.face_mass( + wave_flux(discr, c=c, w_tpair=interior_trace_pair(discr, w)) + + wave_flux(discr, c=c, w_tpair=TracePair( + BTAG_ALL, dir_bval, dir_bc)) + )) + ) + + +# }}} + + +def rk4_step(y, t, h, f): + k1 = f(t, y) + k2 = f(t+h/2, y + h/2*k1) + k3 = f(t+h/2, y + h/2*k2) + k4 = f(t+h, y + h*k3) + return y + h/6*(k1 + 2*k2 + 2*k3 + k4) + + +def bump(discr, queue, t=0): + source_center = np.array([0.0, 0.05]) + source_width = 0.05 + source_omega = 3 + + nodes = discr.nodes().with_queue(queue) + center_dist = join_fields([ + nodes[0] - source_center[0], + nodes[1] - source_center[1], + ]) + + return ( + np.cos(source_omega*t) + * clmath.exp( + -np.dot(center_dist, center_dist) + / source_width**2)) + + +def main(): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + nel_1d = 16 + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5, -0.5), + b=(0.5, 0.5), + n=(nel_1d, nel_1d)) + + order = 3 + + # no deep meaning here, just a fudge factor + dt = 0.75/(nel_1d*order**2) + + print("%d elements" % mesh.nelements) + + discr = EagerDGDiscretization(cl_ctx, mesh, order=order) + + fields = join_fields( + bump(discr, queue), + [discr.zeros(queue) for i in range(discr.dim)] + ) + + vis = make_visualizer(discr, discr.order+3) + + def rhs(t, w): + return wave_operator(discr, c=1, w=w) + + t = 0 + t_final = 3 + istep = 0 + while t < t_final: + fields = rk4_step(fields, t, dt, rhs) + + if istep % 10 == 0: + print(istep, t, la.norm(fields[0].get())) + vis.write_vtk_file("fld-wave-min-%04d.vtu" % istep, + [ + ("u", fields[0]), + ("v", fields[1:]), + ]) + + t += dt + istep += 1 + + +if __name__ == "__main__": + main() + +# vim: foldmethod=marker diff --git a/grudge/eager.py b/grudge/eager.py new file mode 100644 index 00000000..32e489fb --- /dev/null +++ b/grudge/eager.py @@ -0,0 +1,128 @@ +from __future__ import division, print_function + +__copyright__ = "Copyright (C) 2020 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 +from grudge.discretization import DGDiscretizationWithBoundaries +from pytools import memoize_method +from pytools.obj_array import ( + with_object_array_or_scalar, + is_obj_array) +import pyopencl as cl +import pyopencl.array as cla # noqa +from grudge import sym, bind +from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa + + +def with_queue(queue, ary): + return with_object_array_or_scalar( + lambda x: x.with_queue(queue), ary) + + +def without_queue(ary): + return with_queue(None, ary) + + +class EagerDGDiscretization(DGDiscretizationWithBoundaries): + def interp(self, src, tgt, vec): + if is_obj_array(vec): + return with_object_array_or_scalar( + lambda el: self.interp(src, tgt, el), vec) + + return self.connection_from_dds(src, tgt)(vec.queue, vec) + + def nodes(self): + return self._volume_discr.nodes() + + @memoize_method + def parametrization_derivative(self): + with cl.CommandQueue(self.cl_context) as queue: + fmat = sym.forward_metric_derivative_mat( + self.ambient_dim, self.dim) + result = bind(self, fmat.reshape(-1))(queue) + return result.reshape(*fmat.shape) + + @memoize_method + def vol_jacobian(self): + with cl.CommandQueue(self.cl_context) as queue: + [a, b], [c, d] = with_queue(queue, self.parametrization_derivative()) + return (a*d-b*c).with_queue(None) + + @memoize_method + def inverse_parametrization_derivative(self): + with cl.CommandQueue(self.cl_context) as queue: + [a, b], [c, d] = with_queue(queue, self.parametrization_derivative()) + + result = np.zeros((2, 2), dtype=object) + det = a*d-b*c + result[0, 0] = d/det + result[0, 1] = -b/det + result[1, 0] = -c/det + result[1, 1] = a/det + + return without_queue(result) + + @memoize_method + def _bound_grad(self): + return bind(self, sym.nabla(self.dim) * sym.Variable("u")) + + def grad(self, vec): + return self._bound_grad()(vec.queue, u=vec) + + def div(self, vecs): + return sum( + self.grad(vec_i)[i] for i, vec_i in enumerate(vecs)) + + @memoize_method + def normal(self, dd): + with cl.CommandQueue(self.cl_context) as queue: + surface_discr = self.discr_from_dd(dd) + return without_queue( + bind(self, sym.normal( + dd, surface_discr.ambient_dim, surface_discr.dim))(queue)) + + @memoize_method + def _bound_inverse_mass(self): + return bind(self, sym.InverseMassOperator()(sym.Variable("u"))) + + def inverse_mass(self, vec): + if is_obj_array(vec): + return with_object_array_or_scalar( + lambda el: self.inverse_mass(el), vec) + + return self._bound_inverse_mass()(vec.queue, u=vec) + + @memoize_method + def _bound_face_mass(self): + u = sym.Variable("u", dd=sym.as_dofdesc("all_faces")) + return bind(self, sym.FaceMassOperator()(u)) + + def face_mass(self, vec): + if is_obj_array(vec): + return with_object_array_or_scalar( + lambda el: self.face_mass(el), vec) + + return self._bound_face_mass()(vec.queue, u=vec) + +# vim: foldmethod=marker -- GitLab From 69b1707be8a06d68ada315cacaf2ff90a08be0ea Mon Sep 17 00:00:00 2001 From: "[6~" Date: Thu, 21 May 2020 23:43:44 -0500 Subject: [PATCH 2/5] Delete some more unnecessary code in eager interface --- grudge/eager.py | 30 +----------------------------- 1 file changed, 1 insertion(+), 29 deletions(-) diff --git a/grudge/eager.py b/grudge/eager.py index 32e489fb..5cf380f8 100644 --- a/grudge/eager.py +++ b/grudge/eager.py @@ -23,7 +23,7 @@ THE SOFTWARE. """ -import numpy as np +import numpy as np # noqa from grudge.discretization import DGDiscretizationWithBoundaries from pytools import memoize_method from pytools.obj_array import ( @@ -55,34 +55,6 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries): def nodes(self): return self._volume_discr.nodes() - @memoize_method - def parametrization_derivative(self): - with cl.CommandQueue(self.cl_context) as queue: - fmat = sym.forward_metric_derivative_mat( - self.ambient_dim, self.dim) - result = bind(self, fmat.reshape(-1))(queue) - return result.reshape(*fmat.shape) - - @memoize_method - def vol_jacobian(self): - with cl.CommandQueue(self.cl_context) as queue: - [a, b], [c, d] = with_queue(queue, self.parametrization_derivative()) - return (a*d-b*c).with_queue(None) - - @memoize_method - def inverse_parametrization_derivative(self): - with cl.CommandQueue(self.cl_context) as queue: - [a, b], [c, d] = with_queue(queue, self.parametrization_derivative()) - - result = np.zeros((2, 2), dtype=object) - det = a*d-b*c - result[0, 0] = d/det - result[0, 1] = -b/det - result[1, 0] = -c/det - result[1, 1] = a/det - - return without_queue(result) - @memoize_method def _bound_grad(self): return bind(self, sym.nabla(self.dim) * sym.Variable("u")) -- GitLab From 05308996beb8e8c9dbf6198d92b6ac8ade57bacd Mon Sep 17 00:00:00 2001 From: "[6~" Date: Thu, 21 May 2020 23:55:19 -0500 Subject: [PATCH 3/5] Eager example: towards dim independence --- examples/wave/wave-eager.py | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py index 2ec333fc..7daf2b85 100644 --- a/examples/wave/wave-eager.py +++ b/examples/wave/wave-eager.py @@ -29,7 +29,7 @@ import pyopencl as cl import pyopencl.array as cla # noqa import pyopencl.clmath as clmath from pytools.obj_array import ( - join_fields, + join_fields, make_obj_array, with_object_array_or_scalar) from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.eager import EagerDGDiscretization, with_queue @@ -53,23 +53,25 @@ def wave_flux(discr, c, w_tpair): normal = with_queue(u.int.queue, discr.normal(w_tpair.dd)) + def normal_times(scalar): + # workaround for object array behavior + return make_obj_array([ni*scalar for ni in normal]) + flux_weak = join_fields( np.dot(v.avg, normal), - normal[0] * u.avg, - normal[1] * u.avg) + normal_times(u.avg), + ) # upwind v_jump = np.dot(normal, v.int-v.ext) flux_weak -= join_fields( 0.5*(u.int-u.ext), - 0.5*normal[0]*v_jump, - 0.5*normal[1]*v_jump, + 0.5*normal_times(v_jump), ) flux_strong = join_fields( np.dot(v.int, normal), - u.int * normal[0], - u.int * normal[1], + normal_times(u.int), ) - flux_weak return discr.interp(w_tpair.dd, "all_faces", c*flux_strong) @@ -98,7 +100,6 @@ def wave_operator(discr, c, w): )) ) - # }}} @@ -111,7 +112,7 @@ def rk4_step(y, t, h, f): def bump(discr, queue, t=0): - source_center = np.array([0.0, 0.05]) + source_center = np.array([0.2, 0.35, 0.1])[:discr.dim] source_width = 0.05 source_omega = 3 @@ -132,12 +133,13 @@ def main(): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) + dim = 2 nel_1d = 16 from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh( - a=(-0.5, -0.5), - b=(0.5, 0.5), - n=(nel_1d, nel_1d)) + a=(-0.5,)*dim, + b=(0.5,)*dim, + n=(nel_1d,)*dim) order = 3 @@ -166,7 +168,7 @@ def main(): if istep % 10 == 0: print(istep, t, la.norm(fields[0].get())) - vis.write_vtk_file("fld-wave-min-%04d.vtu" % istep, + vis.write_vtk_file("fld-wave-eager-%04d.vtu" % istep, [ ("u", fields[0]), ("v", fields[1:]), -- GitLab From a43187df2c029dbf329fccb4c71765d29b9f7044 Mon Sep 17 00:00:00 2001 From: "[6~" Date: Fri, 22 May 2020 00:07:19 -0500 Subject: [PATCH 4/5] Get eager wave example to be dim-independent --- examples/wave/wave-eager.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py index 7daf2b85..f80ff2c7 100644 --- a/examples/wave/wave-eager.py +++ b/examples/wave/wave-eager.py @@ -118,8 +118,8 @@ def bump(discr, queue, t=0): nodes = discr.nodes().with_queue(queue) center_dist = join_fields([ - nodes[0] - source_center[0], - nodes[1] - source_center[1], + nodes[i] - source_center[i] + for i in range(discr.dim) ]) return ( @@ -143,8 +143,14 @@ def main(): order = 3 - # no deep meaning here, just a fudge factor - dt = 0.75/(nel_1d*order**2) + if dim == 2: + # no deep meaning here, just a fudge factor + dt = 0.75/(nel_1d*order**2) + elif dim == 3: + # no deep meaning here, just a fudge factor + dt = 0.45/(nel_1d*order**2) + else: + raise ValueError("don't have a stable time step guesstimate") print("%d elements" % mesh.nelements) @@ -155,7 +161,7 @@ def main(): [discr.zeros(queue) for i in range(discr.dim)] ) - vis = make_visualizer(discr, discr.order+3) + vis = make_visualizer(discr, discr.order+3 if dim == 2 else discr.order) def rhs(t, w): return wave_operator(discr, c=1, w=w) -- GitLab From 73c0b9396500f21f64019f3fd232258f15008ef1 Mon Sep 17 00:00:00 2001 From: "[6~" Date: Fri, 22 May 2020 00:14:48 -0500 Subject: [PATCH 5/5] Switch eager example over to weak DG --- examples/wave/wave-eager.py | 17 ++++++----------- grudge/eager.py | 11 +++++++++++ 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py index f80ff2c7..aa9e87ff 100644 --- a/examples/wave/wave-eager.py +++ b/examples/wave/wave-eager.py @@ -69,12 +69,7 @@ def wave_flux(discr, c, w_tpair): 0.5*normal_times(v_jump), ) - flux_strong = join_fields( - np.dot(v.int, normal), - normal_times(u.int), - ) - flux_weak - - return discr.interp(w_tpair.dd, "all_faces", c*flux_strong) + return discr.interp(w_tpair.dd, "all_faces", c*flux_weak) def wave_operator(discr, c, w): @@ -87,12 +82,12 @@ def wave_operator(discr, c, w): dir_bc = join_fields(-dir_u, dir_v) return ( - - join_fields( - -c*discr.div(v), - -c*discr.grad(u) - ) - + # noqa: W504 discr.inverse_mass( + join_fields( + c*discr.weak_div(v), + c*discr.weak_grad(u) + ) + - # noqa: W504 discr.face_mass( wave_flux(discr, c=c, w_tpair=interior_trace_pair(discr, w)) + wave_flux(discr, c=c, w_tpair=TracePair( diff --git a/grudge/eager.py b/grudge/eager.py index 5cf380f8..e64e2fa8 100644 --- a/grudge/eager.py +++ b/grudge/eager.py @@ -66,6 +66,17 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries): return sum( self.grad(vec_i)[i] for i, vec_i in enumerate(vecs)) + @memoize_method + def _bound_weak_grad(self): + return bind(self, sym.stiffness_t(self.dim) * sym.Variable("u")) + + def weak_grad(self, vec): + return self._bound_weak_grad()(vec.queue, u=vec) + + def weak_div(self, vecs): + return sum( + self.weak_grad(vec_i)[i] for i, vec_i in enumerate(vecs)) + @memoize_method def normal(self, dd): with cl.CommandQueue(self.cl_context) as queue: -- GitLab