From 7982284d48bbc8702fd6438b858e4e0c5ff3121a Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 4 Dec 2017 10:26:37 -0600 Subject: [PATCH 01/30] In the middle of fluxy quadratures, 2 --- examples/advection/var-velocity.py | 7 ++- examples/quad/flux_quad.py | 91 +++++++++++++++++++++++++++ examples/quad/mass_quad.py | 90 +++++++++++++++++++++++++++ examples/quad/quad.py | 80 ++++++++++++++++++++++++ examples/quad/stiffness_quad.py | 99 ++++++++++++++++++++++++++++++ grudge/discretization.py | 33 ++++++++-- grudge/execution.py | 59 +++++++++++------- grudge/models/advection.py | 2 +- grudge/symbolic/operators.py | 85 ++++++++++++++++++++----- grudge/symbolic/primitives.py | 11 +++- 10 files changed, 505 insertions(+), 52 deletions(-) create mode 100644 examples/quad/flux_quad.py create mode 100644 examples/quad/mass_quad.py create mode 100644 examples/quad/quad.py create mode 100644 examples/quad/stiffness_quad.py diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index ddc75352..0a69967e 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -37,12 +37,13 @@ def main(write_output=True, order=4): dim = 2 + n_divs = 15 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=(15, 15), order=order) + n=(n_divs, n_divs), order=order) - dt_factor = 4 - h = 1/20 + dt_factor = 5 + h = 1/n_divs discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) diff --git a/examples/quad/flux_quad.py b/examples/quad/flux_quad.py new file mode 100644 index 00000000..3d7f9a32 --- /dev/null +++ b/examples/quad/flux_quad.py @@ -0,0 +1,91 @@ +# 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 . + + +import numpy as np +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + + + +def main(order=1): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + + from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh + #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + n=(3, 3), order=7) + + + + + pquad_dd = sym.DOFDesc("vol", "product") + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) + + + magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL, "product") + #magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL) + def fluxy_stuff(u): + return sym.interp(u.dd, magic_thing)(u.int) + + + ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0) + nds = sym.nodes(dim) + op = sym.FaceMassOperator(magic_thing, "vol")(fluxy_stuff(sym.int_tpair(nds, "product"))) + #op = (sym.int_tpair(nds, "product").int) + #op = sym.FaceMassOperator()(fluxy_stuff(sym.int_tpair(nds))) + print(op) + + print(sym.pretty(op)) + u = bind(discr, op)(queue, t=0) + + + print(u) + + + + + #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ]) + #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ]) + + + + + + + +if __name__ == "__main__": + main(4) + + diff --git a/examples/quad/mass_quad.py b/examples/quad/mass_quad.py new file mode 100644 index 00000000..05c306bd --- /dev/null +++ b/examples/quad/mass_quad.py @@ -0,0 +1,90 @@ +# 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 . + + +import numpy as np +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + + + +def main(order=1): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + + from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh + #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + n=(3, 3), order=7) + + + + + pquad_dd = sym.DOFDesc("vol", "product") + to_pquad = sym.interp("vol", pquad_dd) + + def f(x): + x = x[0] / x[0] + #return sym.MassOperator()(x) + return sym.MassOperator(pquad_dd, "vol")(to_pquad(x)) + #return np.dot(from_pquad_stiffness_t, to_pquad(x)) + #return sym.interp("vol", pquad_dd)(sym.sin(3*x[0])) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) + + + + + + + + ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0) + u = bind(discr, f(sym.nodes(dim)))(queue, t=0) + + print(np.dot(ones.T, u)) + + + + + #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ]) + #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ]) + + + + + + + +if __name__ == "__main__": + main(4) + + diff --git a/examples/quad/quad.py b/examples/quad/quad.py new file mode 100644 index 00000000..67c41fa2 --- /dev/null +++ b/examples/quad/quad.py @@ -0,0 +1,80 @@ +# 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 . + + +import numpy as np +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + + + +def main(order=1): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + + from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh + #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + n=(3, 3), order=7) + + + + + pquad_dd = sym.DOFDesc("vol", "product") + + def f(x): + #return sym.sin(3*x[0]) + return sym.interp("vol", pquad_dd)(sym.sin(3*x[0])) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) + + u = bind(discr, f(sym.nodes(dim)))(queue, t=0) + + + + from meshmode.discretization.visualization import make_visualizer + #vis = make_visualizer(queue, discr.quad_volume_discr['product'], 2*order) + #vis = make_visualizer(queue, discr.volume_discr, order) + + #vis.write_vtk_file("fld-000.vtu", [ ("u", u) ]) + print(u) + + + + + + + +if __name__ == "__main__": + main(4) + + diff --git a/examples/quad/stiffness_quad.py b/examples/quad/stiffness_quad.py new file mode 100644 index 00000000..ff334c4e --- /dev/null +++ b/examples/quad/stiffness_quad.py @@ -0,0 +1,99 @@ +# 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 . + + +import numpy as np +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + + + +def main(order=1): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + + from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh + #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + n=(3, 3), order=7) + + + + + pquad_dd = sym.DOFDesc("vol", "product") + to_pquad = sym.interp("vol", pquad_dd) + from_pquad_stiffness_t = sym.stiffness_t(dim, pquad_dd, "vol") + from_pquad_mass = sym.MassOperator() + + def lin(x): + return x[0] + + def non_lin(x): + return x[1] * x[0] * x[0] + 3* x[0] * x[0] - 4 + + def f(x, fun): + #return sym.MassOperator()(x) + #return np.dot(sym.stiffness_t(2), to_pquad(fun(x))) + return np.dot(from_pquad_stiffness_t, fun(to_pquad(x))) + + def reg(x, fun): + #return sym.MassOperator()(x) + #return np.dot(sym.stiffness_t(2), to_pquad(fun(x))) + return np.dot(sym.stiffness_t(2), fun(x)) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) + + + #ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0) + nds = sym.nodes(dim) + #u = bind(discr, to_pquad(nds))(queue, t=0) + #u = bind(discr, reg(nds, lin))(queue, t=0) + u = bind(discr, f(nds, lin))(queue, t=0) + + print(u) + + + + + #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ]) + #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ]) + + + + + + + +if __name__ == "__main__": + main(4) + + diff --git a/grudge/discretization.py b/grudge/discretization.py index a1ca6303..05360ecd 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -82,8 +82,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): if dd.is_volume(): if qtag is not sym.QTAG_NONE: - # FIXME - raise NotImplementedError("quadrature") + return self._quad_volume_discr(qtag) return self._volume_discr elif dd.domain_tag is sym.FACE_RESTR_ALL: @@ -120,6 +119,12 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): self._boundary_to_quad_connection( to_dd.domain_tag.tag, to_dd.quadrature_tag)]) + elif to_dd.is_volume(): + from meshmode.discretization.connection import make_same_mesh_connection + to_discr = self._quad_volume_discr(to_dd.quadrature_tag) + from_discr = self._volume_discr + return make_same_mesh_connection(to_discr, from_discr) + else: raise ValueError("cannot interpolate from volume to: " + str(to_dd)) @@ -130,6 +135,24 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): raise ValueError( "cannot interpolate from interior faces to: " + str(to_dd)) + elif from_dd.domain_tag is sym.FACE_RESTR_ALL: + if to_dd.domain_tag is sym.FACE_RESTR_ALL and qtag is not sym.QTAG_NONE: + from meshmode.discretization.connection import make_face_restriction, FACE_RESTR_ALL + return make_face_restriction( + self._all_faces_discr(), + self.group_factory_for_quadrature_tag(qtag), + FACE_RESTR_ALL, + + # FIXME: This will need to change as soon as we support + # pyramids or other elements with non-identical face + # types. + per_face_groups=False) + + #return self._all_faces_connection(None, qtag) + else: + raise ValueError( + "cannot interpolate from interior faces to: " + + str(to_dd)) elif from_dd.is_boundary(): if to_dd.domain_tag is sym.FACE_RESTR_ALL and qtag is sym.QTAG_NONE: @@ -155,7 +178,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): QuadratureSimplexGroupFactory if quadrature_tag is not sym.QTAG_NONE: - return QuadratureSimplexGroupFactory(order=self.order) + return QuadratureSimplexGroupFactory(order=self.quad_min_degrees[quadrature_tag]) else: return PolynomialWarpAndBlendGroupFactory(order=self.order) @@ -259,7 +282,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): return self._all_faces_volume_connection(quadrature_tag).to_discr @memoize_method - def _all_faces_connection(self, boundary_tag): + def _all_faces_connection(self, boundary_tag, to_qtag=None): """Return a :class:`meshmode.discretization.connection.DiscretizationConnection` that goes from either @@ -277,7 +300,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): else: faces_conn = self._boundary_connection(boundary_tag) - return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr()) + return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr(to_qtag)) # }}} diff --git a/grudge/execution.py b/grudge/execution.py index e4aebe5c..40dcb19a 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -198,36 +198,42 @@ class ExecutionMapper(mappers.Evaluator, knl = lp.make_kernel( """{[k,i,j]: 0<=k Date: Mon, 4 Dec 2017 10:26:37 -0600 Subject: [PATCH 02/30] In the middle of fluxy quadratures, 2 --- examples/advection/var-velocity.py | 7 ++- examples/quad/flux_quad.py | 91 +++++++++++++++++++++++++++ examples/quad/mass_quad.py | 90 +++++++++++++++++++++++++++ examples/quad/quad.py | 80 ++++++++++++++++++++++++ examples/quad/stiffness_quad.py | 99 ++++++++++++++++++++++++++++++ grudge/discretization.py | 33 ++++++++-- grudge/execution.py | 59 +++++++++++------- grudge/models/advection.py | 2 +- grudge/symbolic/operators.py | 85 ++++++++++++++++++++----- grudge/symbolic/primitives.py | 11 +++- 10 files changed, 505 insertions(+), 52 deletions(-) create mode 100644 examples/quad/flux_quad.py create mode 100644 examples/quad/mass_quad.py create mode 100644 examples/quad/quad.py create mode 100644 examples/quad/stiffness_quad.py diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index ddc75352..0a69967e 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -37,12 +37,13 @@ def main(write_output=True, order=4): dim = 2 + n_divs = 15 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=(15, 15), order=order) + n=(n_divs, n_divs), order=order) - dt_factor = 4 - h = 1/20 + dt_factor = 5 + h = 1/n_divs discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) diff --git a/examples/quad/flux_quad.py b/examples/quad/flux_quad.py new file mode 100644 index 00000000..3d7f9a32 --- /dev/null +++ b/examples/quad/flux_quad.py @@ -0,0 +1,91 @@ +# 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 . + + +import numpy as np +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + + + +def main(order=1): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + + from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh + #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + n=(3, 3), order=7) + + + + + pquad_dd = sym.DOFDesc("vol", "product") + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) + + + magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL, "product") + #magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL) + def fluxy_stuff(u): + return sym.interp(u.dd, magic_thing)(u.int) + + + ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0) + nds = sym.nodes(dim) + op = sym.FaceMassOperator(magic_thing, "vol")(fluxy_stuff(sym.int_tpair(nds, "product"))) + #op = (sym.int_tpair(nds, "product").int) + #op = sym.FaceMassOperator()(fluxy_stuff(sym.int_tpair(nds))) + print(op) + + print(sym.pretty(op)) + u = bind(discr, op)(queue, t=0) + + + print(u) + + + + + #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ]) + #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ]) + + + + + + + +if __name__ == "__main__": + main(4) + + diff --git a/examples/quad/mass_quad.py b/examples/quad/mass_quad.py new file mode 100644 index 00000000..05c306bd --- /dev/null +++ b/examples/quad/mass_quad.py @@ -0,0 +1,90 @@ +# 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 . + + +import numpy as np +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + + + +def main(order=1): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + + from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh + #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + n=(3, 3), order=7) + + + + + pquad_dd = sym.DOFDesc("vol", "product") + to_pquad = sym.interp("vol", pquad_dd) + + def f(x): + x = x[0] / x[0] + #return sym.MassOperator()(x) + return sym.MassOperator(pquad_dd, "vol")(to_pquad(x)) + #return np.dot(from_pquad_stiffness_t, to_pquad(x)) + #return sym.interp("vol", pquad_dd)(sym.sin(3*x[0])) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) + + + + + + + + ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0) + u = bind(discr, f(sym.nodes(dim)))(queue, t=0) + + print(np.dot(ones.T, u)) + + + + + #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ]) + #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ]) + + + + + + + +if __name__ == "__main__": + main(4) + + diff --git a/examples/quad/quad.py b/examples/quad/quad.py new file mode 100644 index 00000000..67c41fa2 --- /dev/null +++ b/examples/quad/quad.py @@ -0,0 +1,80 @@ +# 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 . + + +import numpy as np +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + + + +def main(order=1): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + + from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh + #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + n=(3, 3), order=7) + + + + + pquad_dd = sym.DOFDesc("vol", "product") + + def f(x): + #return sym.sin(3*x[0]) + return sym.interp("vol", pquad_dd)(sym.sin(3*x[0])) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) + + u = bind(discr, f(sym.nodes(dim)))(queue, t=0) + + + + from meshmode.discretization.visualization import make_visualizer + #vis = make_visualizer(queue, discr.quad_volume_discr['product'], 2*order) + #vis = make_visualizer(queue, discr.volume_discr, order) + + #vis.write_vtk_file("fld-000.vtu", [ ("u", u) ]) + print(u) + + + + + + + +if __name__ == "__main__": + main(4) + + diff --git a/examples/quad/stiffness_quad.py b/examples/quad/stiffness_quad.py new file mode 100644 index 00000000..ff334c4e --- /dev/null +++ b/examples/quad/stiffness_quad.py @@ -0,0 +1,99 @@ +# 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 . + + +import numpy as np +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + + + +def main(order=1): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + + from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh + #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + n=(3, 3), order=7) + + + + + pquad_dd = sym.DOFDesc("vol", "product") + to_pquad = sym.interp("vol", pquad_dd) + from_pquad_stiffness_t = sym.stiffness_t(dim, pquad_dd, "vol") + from_pquad_mass = sym.MassOperator() + + def lin(x): + return x[0] + + def non_lin(x): + return x[1] * x[0] * x[0] + 3* x[0] * x[0] - 4 + + def f(x, fun): + #return sym.MassOperator()(x) + #return np.dot(sym.stiffness_t(2), to_pquad(fun(x))) + return np.dot(from_pquad_stiffness_t, fun(to_pquad(x))) + + def reg(x, fun): + #return sym.MassOperator()(x) + #return np.dot(sym.stiffness_t(2), to_pquad(fun(x))) + return np.dot(sym.stiffness_t(2), fun(x)) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) + + + #ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0) + nds = sym.nodes(dim) + #u = bind(discr, to_pquad(nds))(queue, t=0) + #u = bind(discr, reg(nds, lin))(queue, t=0) + u = bind(discr, f(nds, lin))(queue, t=0) + + print(u) + + + + + #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ]) + #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ]) + + + + + + + +if __name__ == "__main__": + main(4) + + diff --git a/grudge/discretization.py b/grudge/discretization.py index 2d15c4a7..5a9cea3e 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -82,8 +82,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): if dd.is_volume(): if qtag is not sym.QTAG_NONE: - # FIXME - raise NotImplementedError("quadrature") + return self._quad_volume_discr(qtag) return self._volume_discr elif dd.domain_tag is sym.FACE_RESTR_ALL: @@ -150,6 +149,12 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): assert from_dd.quadrature_tag is sym.QTAG_NONE return self._boundary_connection(to_dd.domain_tag) + elif to_dd.is_volume(): + from meshmode.discretization.connection import make_same_mesh_connection + to_discr = self._quad_volume_discr(to_dd.quadrature_tag) + from_discr = self._volume_discr + return make_same_mesh_connection(to_discr, from_discr) + else: raise ValueError("cannot interpolate from volume to: " + str(to_dd)) @@ -160,6 +165,24 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): raise ValueError( "cannot interpolate from interior faces to: " + str(to_dd)) + elif from_dd.domain_tag is sym.FACE_RESTR_ALL: + if to_dd.domain_tag is sym.FACE_RESTR_ALL and qtag is not sym.QTAG_NONE: + from meshmode.discretization.connection import make_face_restriction, FACE_RESTR_ALL + return make_face_restriction( + self._all_faces_discr(), + self.group_factory_for_quadrature_tag(qtag), + FACE_RESTR_ALL, + + # FIXME: This will need to change as soon as we support + # pyramids or other elements with non-identical face + # types. + per_face_groups=False) + + #return self._all_faces_connection(None, qtag) + else: + raise ValueError( + "cannot interpolate from interior faces to: " + + str(to_dd)) elif from_dd.domain_tag is sym.FACE_RESTR_ALL: if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE: @@ -189,7 +212,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): QuadratureSimplexGroupFactory if quadrature_tag is not sym.QTAG_NONE: - return QuadratureSimplexGroupFactory(order=self.order) + return QuadratureSimplexGroupFactory(order=self.quad_min_degrees[quadrature_tag]) else: return PolynomialWarpAndBlendGroupFactory(order=self.order) @@ -281,7 +304,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): return self._all_faces_volume_connection(quadrature_tag).to_discr @memoize_method - def _all_faces_connection(self, boundary_tag): + def _all_faces_connection(self, boundary_tag, to_qtag=None): """Return a :class:`meshmode.discretization.connection.DiscretizationConnection` that goes from either @@ -299,7 +322,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): else: faces_conn = self._boundary_connection(boundary_tag) - return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr()) + return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr(to_qtag)) # }}} diff --git a/grudge/execution.py b/grudge/execution.py index e4aebe5c..40dcb19a 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -198,36 +198,42 @@ class ExecutionMapper(mappers.Evaluator, knl = lp.make_kernel( """{[k,i,j]: 0<=k Date: Mon, 4 Dec 2017 13:13:05 -0600 Subject: [PATCH 03/30] cleaned up discretization.py a tad --- grudge/discretization.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/grudge/discretization.py b/grudge/discretization.py index c554e614..e559220e 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -148,19 +148,11 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): elif to_dd.is_boundary(): assert from_dd.quadrature_tag is sym.QTAG_NONE return self._boundary_connection(to_dd.domain_tag) - - elif to_dd.is_volume(): - from meshmode.discretization.connection import make_same_mesh_connection - to_discr = self._quad_volume_discr(to_dd.quadrature_tag) - from_discr = self._volume_discr - return make_same_mesh_connection(to_discr, from_discr) - elif to_dd.is_volume(): from meshmode.discretization.connection import make_same_mesh_connection to_discr = self._quad_volume_discr(to_dd.quadrature_tag) from_discr = self._volume_discr return make_same_mesh_connection(to_discr, from_discr) - else: raise ValueError("cannot interpolate from volume to: " + str(to_dd)) -- GitLab From a7bdf749a9e646e425a2b151238b6114d6fae20a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 4 Dec 2017 21:52:07 -0600 Subject: [PATCH 04/30] Further clean up DGDiscretizationWithBoundaries.{discr_from_dd,connection_from_dds} --- grudge/discretization.py | 152 ++++++++++++++++----------------------- 1 file changed, 60 insertions(+), 92 deletions(-) diff --git a/grudge/discretization.py b/grudge/discretization.py index e559220e..2babc3a4 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -85,12 +85,23 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): return self._quad_volume_discr(qtag) return self._volume_discr - elif dd.domain_tag is sym.FACE_RESTR_ALL: - return self._all_faces_discr(qtag) + if qtag is not sym.QTAG_NONE: + no_quad_discr = self.discr_from_dd(sym.DOFDesc(dd.domain_tag)) + + from meshmode.discretization import Discretization + return Discretization( + self._volume_discr.cl_context, + no_quad_discr.mesh, + self.group_factory_for_quadrature_tag(qtag)) + + assert qtag is sym.QTAG_NONE + + if dd.domain_tag is sym.FACE_RESTR_ALL: + return self._all_faces_volume_connection().to_discr elif dd.domain_tag is sym.FACE_RESTR_INTERIOR: - return self._interior_faces_discr(qtag) + return self._interior_faces_connection().to_discr elif dd.is_boundary(): - return self._boundary_discr(dd.domain_tag, qtag) + return self._boundary_connection(dd.domain_tag).to_discr else: raise ValueError("DOF desc tag not understood: " + str(dd)) @@ -99,19 +110,33 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): from_dd = sym.as_dofdesc(from_dd) to_dd = sym.as_dofdesc(to_dd) - if from_dd.quadrature_tag is not sym.QTAG_NONE: - raise ValueError("cannot interpolate *from* a " - "(non-interpolatory) quadrature grid") - to_qtag = to_dd.quadrature_tag + if ( + not from_dd.is_volume() + and from_dd.quadrature_tag == to_dd.quadrature_tag + and to_dd.domain_tag is sym.FACE_RESTR_ALL): + faces_conn = self.connection_from_dds( + sym.DOFDesc("vol"), + sym.DOFDesc(from_dd.domain_tag)) + + from meshmode.discretization.connection import \ + make_face_to_all_faces_embedding + + return make_face_to_all_faces_embedding( + faces_conn, self.discr_from_dd(to_dd), + self.discr_from_dd(from_dd)) + # {{{ simplify domain + qtag change into chained - if (from_dd.domain_tag != to_dd.domain_tag - and from_dd.quadrature_tag != to_dd.quadrature_tag): + if ( + from_dd.domain_tag != to_dd.domain_tag + and from_dd.quadrature_tag is sym.QTAG_NONE + and to_dd.quadrature_tag is not sym.QTAG_NONE): - from meshmode.connection import ChainedDiscretizationConnection - intermediate_dd = sym.DOFDesc(from_dd.domain_tag) + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection + intermediate_dd = sym.DOFDesc(to_dd.domain_tag) return ChainedDiscretizationConnection( [ # first change domain @@ -129,8 +154,10 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): # {{{ generic to-quad - if (from_dd.domain_tag == to_dd.domain_tag - and from_dd.quadrature_tag != to_dd.quadrature_tag): + if ( + from_dd.domain_tag == to_dd.domain_tag + and from_dd.quadrature_tag is sym.QTAG_NONE + and to_dd.quadrature_tag is not sym.QTAG_NONE): from meshmode.discretization.connection.same_mesh import \ make_same_mesh_connection @@ -140,41 +167,29 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): # }}} + if from_dd.quadrature_tag is not sym.QTAG_NONE: + raise ValueError("cannot interpolate *from* a " + "(non-interpolatory) quadrature grid") + + assert to_qtag is sym.QTAG_NONE + if from_dd.is_volume(): if to_dd.domain_tag is sym.FACE_RESTR_ALL: - return self._all_faces_volume_connection(to_qtag) + return self._all_faces_volume_connection() if to_dd.domain_tag is sym.FACE_RESTR_INTERIOR: - return self._interior_faces_connection(to_qtag) + return self._interior_faces_connection() elif to_dd.is_boundary(): assert from_dd.quadrature_tag is sym.QTAG_NONE return self._boundary_connection(to_dd.domain_tag) elif to_dd.is_volume(): - from meshmode.discretization.connection import make_same_mesh_connection + from meshmode.discretization.connection import \ + make_same_mesh_connection to_discr = self._quad_volume_discr(to_dd.quadrature_tag) from_discr = self._volume_discr return make_same_mesh_connection(to_discr, from_discr) - else: - raise ValueError("cannot interpolate from volume to: " + str(to_dd)) - - elif from_dd.domain_tag is sym.FACE_RESTR_INTERIOR: - if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE: - return self._all_faces_connection(None) - else: - raise ValueError( - "cannot interpolate from interior faces to: " - + str(to_dd)) - - elif from_dd.domain_tag is sym.FACE_RESTR_ALL: - if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE: - return self._all_faces_connection(None) - elif from_dd.is_boundary(): - if to_dd.domain_tag is sym.FACE_RESTR_ALL and to_qtag is sym.QTAG_NONE: - return self._all_faces_connection(from_dd.domain_tag) else: - raise ValueError( - "cannot interpolate from interior faces to: " - + str(to_dd)) + raise ValueError("cannot interpolate from volume to: " + str(to_dd)) else: raise ValueError("cannot interpolate from: " + str(from_dd)) @@ -192,7 +207,8 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): QuadratureSimplexGroupFactory if quadrature_tag is not sym.QTAG_NONE: - return QuadratureSimplexGroupFactory(order=self.quad_min_degrees[quadrature_tag]) + return QuadratureSimplexGroupFactory( + order=self.quad_min_degrees[quadrature_tag]) else: return PolynomialWarpAndBlendGroupFactory(order=self.order) @@ -213,33 +229,17 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): self.group_factory_for_quadrature_tag(sym.QTAG_NONE), boundary_tag=boundary_tag) - @memoize_method - def _boundary_discr(self, boundary_tag, quadrature_tag=None): - if quadrature_tag is None: - quadrature_tag = sym.QTAG_NONE - - if quadrature_tag is sym.QTAG_NONE: - return self._boundary_connection(boundary_tag).to_discr - else: - no_quad_bdry_discr = self.boundary_discr(boundary_tag, sym.QTAG_NONE) - - from meshmode.discretization import Discretization - return Discretization( - self._volume_discr.cl_context, - no_quad_bdry_discr.mesh, - self.group_factory_for_quadrature_tag(quadrature_tag)) - # }}} # {{{ interior faces @memoize_method - def _interior_faces_connection(self, quadrature_tag=None): + def _interior_faces_connection(self): from meshmode.discretization.connection import ( make_face_restriction, FACE_RESTR_INTERIOR) return make_face_restriction( self._volume_discr, - self.group_factory_for_quadrature_tag(quadrature_tag), + self.group_factory_for_quadrature_tag(sym.QTAG_NONE), FACE_RESTR_INTERIOR, # FIXME: This will need to change as soon as we support @@ -247,32 +247,24 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): # types. per_face_groups=False) - def _interior_faces_discr(self, quadrature_tag=None): - return self._interior_faces_connection(quadrature_tag).to_discr - @memoize_method - def opposite_face_connection(self, quadrature_tag): - if quadrature_tag is not sym.QTAG_NONE: - # FIXME - raise NotImplementedError("quadrature") - + def opposite_face_connection(self): from meshmode.discretization.connection import \ make_opposite_face_connection - return make_opposite_face_connection( - self._interior_faces_connection(quadrature_tag)) + return make_opposite_face_connection(self._interior_faces_connection()) # }}} # {{{ all-faces @memoize_method - def _all_faces_volume_connection(self, quadrature_tag=None): + def _all_faces_volume_connection(self): from meshmode.discretization.connection import ( make_face_restriction, FACE_RESTR_ALL) return make_face_restriction( self._volume_discr, - self.group_factory_for_quadrature_tag(quadrature_tag), + self.group_factory_for_quadrature_tag(sym.QTAG_NONE), FACE_RESTR_ALL, # FIXME: This will need to change as soon as we support @@ -280,30 +272,6 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): # types. per_face_groups=False) - def _all_faces_discr(self, quadrature_tag=None): - return self._all_faces_volume_connection(quadrature_tag).to_discr - - @memoize_method - def _all_faces_connection(self, boundary_tag, to_qtag=None): - """Return a - :class:`meshmode.discretization.connection.DiscretizationConnection` - that goes from either - :meth:`_interior_faces_discr` (if *boundary_tag* is None) - or - :meth:`_boundary_discr` (if *boundary_tag* is not None) - to a discretization containing all the faces of the volume - discretization. - """ - from meshmode.discretization.connection import \ - make_face_to_all_faces_embedding - - if boundary_tag is None: - faces_conn = self._interior_faces_connection() - else: - faces_conn = self._boundary_connection(boundary_tag) - - return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr(to_qtag)) - # }}} @property -- GitLab From 874cc7ce59c923c7ba2a5743ae8b60f66da2ad0a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 4 Dec 2017 22:04:13 -0600 Subject: [PATCH 05/30] Remove qtag also from invocation of opposite_face_connection --- grudge/execution.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index 40dcb19a..4a51477c 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -252,14 +252,10 @@ class ExecutionMapper(mappers.Evaluator, return conn(self.queue, self.rec(field_expr)).with_queue(self.queue) def map_opposite_interior_face_swap(self, op, field_expr): - dd = op.dd_in + if op.dd_in.quadrature_tag is not sym.QTAG_NONE: + raise ValueError("face swap unavailable on quadrature meshes") - qtag = dd.quadrature_tag - if qtag is None: - # FIXME: Remove once proper quadrature support arrives - qtag = sym.QTAG_NONE - - return self.discrwb.opposite_face_connection(qtag)( + return self.discrwb.opposite_face_connection()( self.queue, self.rec(field_expr)).with_queue(self.queue) # {{{ face mass operator -- GitLab From e11a4c0373bebe0d4cfb389cf967999da5d9dc87 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 5 Dec 2017 00:00:08 -0600 Subject: [PATCH 06/30] Fixed mass matrix crash --- grudge/symbolic/operators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 91153604..3426b8e3 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -389,7 +389,7 @@ class RefMassOperator(RefMassOperatorBase): def matrix(out_element_group, in_element_group): #FIXME: Is this normal? Will bugs surface if I remove this? if out_element_group == in_element_group: - return element_group.mass_matrix() + return in_element_group.mass_matrix() n_quad_nodes = in_element_group.nunit_nodes n_reg_nodes = out_element_group.nunit_nodes -- GitLab From 1d446c2dbdfa759686d591eb69e9526a7ce53f7e Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 6 Dec 2017 12:22:16 -0600 Subject: [PATCH 07/30] Got quadrature working --- examples/advection/var-velocity.py | 17 ++++++++++------- grudge/execution.py | 7 +++---- grudge/models/advection.py | 27 +++++++++++++++++++++------ grudge/symbolic/operators.py | 15 +++++++++------ grudge/symbolic/primitives.py | 13 ++++++++++--- 5 files changed, 53 insertions(+), 26 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 0a69967e..97e793fd 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -30,6 +30,8 @@ logger = logging.getLogger(__name__) from grudge import sym, bind, DGDiscretizationWithBoundaries from pytools.obj_array import join_fields +import meshmode.discretization.poly_element as poly_element + def main(write_output=True, order=4): cl_ctx = cl.create_some_context() @@ -37,15 +39,15 @@ def main(write_output=True, order=4): dim = 2 - n_divs = 15 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=(n_divs, n_divs), order=order) + n=(10, 10), order=order) dt_factor = 5 - h = 1/n_divs + h = 1/10 + final_time = 5 + - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) sym_x = sym.nodes(2) @@ -69,10 +71,11 @@ def main(write_output=True, order=4): from grudge.models.advection import VariableCoefficientAdvectionOperator - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2*order}) op = VariableCoefficientAdvectionOperator(2, advec_v, - inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)), + u_analytic(sym.nodes(dim, sym.BTAG_ALL)),quad_tag="product", flux_type=flux_type) + #op = sym.interp("vol", sym.DOFDesc("vol", "product"))(sym.var("u")) bound_op = bind(discr, op.sym_operator()) @@ -81,7 +84,6 @@ def main(write_output=True, order=4): def rhs(t, u): return bound_op(queue, t=t, u=u) - final_time = 2 dt = dt_factor * h/order**2 nsteps = (final_time // dt) + 1 @@ -99,6 +101,7 @@ def main(write_output=True, order=4): if isinstance(event, dt_stepper.StateComputed): step += 1 + print(step) #print(step, event.t, norm(queue, u=event.state_component[0])) vis.write_vtk_file("fld-%04d.vtu" % step, diff --git a/grudge/execution.py b/grudge/execution.py index 4a51477c..b950a961 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -252,8 +252,8 @@ class ExecutionMapper(mappers.Evaluator, return conn(self.queue, self.rec(field_expr)).with_queue(self.queue) def map_opposite_interior_face_swap(self, op, field_expr): - if op.dd_in.quadrature_tag is not sym.QTAG_NONE: - raise ValueError("face swap unavailable on quadrature meshes") + #if op.dd_in.quadrature_tag is not sym.QTAG_NONE: + #raise ValueError("face swap unavailable on quadrature meshes") return self.discrwb.opposite_face_connection()( self.queue, self.rec(field_expr)).with_queue(self.queue) @@ -291,8 +291,7 @@ class ExecutionMapper(mappers.Evaluator, assert len(all_faces_discr.groups) == len(vol_discr.groups) - for cgrp, afgrp, volgrp in zip(all_faces_conn.groups, - all_faces_discr.groups, vol_discr.groups): + for afgrp, volgrp in zip(all_faces_discr.groups, vol_discr.groups): cache_key = "face_mass", afgrp, op, field.dtype nfaces = volgrp.mesh_el_group.nfaces diff --git a/grudge/models/advection.py b/grudge/models/advection.py index e0203378..4d678832 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -132,13 +132,18 @@ class WeakAdvectionOperator(AdvectionOperatorBase): # {{{ variable-coefficient advection class VariableCoefficientAdvectionOperator(HyperbolicOperator): - def __init__(self, dim, v, inflow_u, flux_type="central"): + def __init__(self, dim, v, inflow_u, flux_type="central", quad_tag="product"): self.ambient_dim = dim self.v = v self.inflow_u = inflow_u self.flux_type = flux_type + if quad_tag is None: + self.quad_tag = sym.QTAG_NONE + else: + self.quad_tag = quad_tag def flux(self, u): + normal = sym.normal(u. dd, self.ambient_dim) surf_v = sym.interp("vol", u.dd)(self.v) @@ -167,15 +172,25 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): # boundary conditions ------------------------------------------------- bc_in = self.inflow_u + all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) def flux(pair): - return sym.interp(pair.dd, "all_faces")( + return sym.interp(pair.dd, all_faces_dd)( self.flux(pair)) + quad_dd = sym.DOFDesc("vol", self.quad_tag) + + if self.quad_tag == sym.QTAG_NONE: + to_quad = lambda x : x + else: + to_quad = sym.interp("vol", quad_dd) + + from_quad_stiffness_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol") + return sym.InverseMassOperator()( - np.dot(sym.stiffness_t(self.ambient_dim), sym.cse(self.v*u)) - - sym.FaceMassOperator()( - flux(sym.int_tpair(u)) - #+ flux(sym.bv_tpair(sym.BTAG_ALL, u, bc_in)) + np.dot(from_quad_stiffness_t, to_quad(self.v)*to_quad(u)) + - sym.FaceMassOperator(all_faces_dd, "vol")( + flux(sym.int_tpair(u, self.quad_tag)) + + flux(sym.bv_tpair(sym.DOFDesc(sym.BTAG_ALL, self.quad_tag), u, bc_in)) # FIXME: Add back support for inflow/outflow tags #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in)) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 3426b8e3..ff909f7f 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -97,12 +97,15 @@ class ElementwiseLinearOperator(Operator): class InterpolationOperator(Operator): - #def __init__(self, dd_in, dd_out): - #print("ayy") - #print(dd_in) - #print(dd_out) - - #super().__init__(dd_in, dd_out) + def __init__(self, dd_in, dd_out): + official_dd_in = _sym().as_dofdesc(dd_in) + official_dd_out = _sym().as_dofdesc(dd_out) + if official_dd_in == official_dd_out: + from warnings import warn + raise ValueError("Interpolating from {} to {}" + " does not do anything.".format(official_dd_in, official_dd_out)) + + super().__init__(dd_in, dd_out) mapper_method = intern("map_interpolation") diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 91122a89..9234c225 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -654,10 +654,17 @@ class TracePair: def int_tpair(expression, qtag=None): - dd = _sym().DOFDesc("int_faces", qtag) - i = cse(_sym().interp("vol", dd)(expression)) + i = _sym().interp("vol", "int_faces")(expression) e = cse(_sym().OppositeInteriorFaceSwap()(i)) - return TracePair(dd, i, e) + + if qtag is not None and qtag != _sym().QTAG_NONE: + q_dd = _sym().DOFDesc("int_faces", qtag) + i = cse(_sym().interp("int_faces", q_dd)(i)) + e = cse(_sym().interp("int_faces", q_dd)(e)) + else: + q_dd = "int_faces" + + return TracePair(q_dd, i, e) #i = cse(_sym().interp("vol", "int_faces")(expression)) #e = cse(_sym().OppositeInteriorFaceSwap()(i)) -- GitLab From 7c2c6d6413bccc2a9a411b7e9b242971435f9185 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 6 Dec 2017 19:42:30 -0600 Subject: [PATCH 08/30] Cleaned things, placated flake --- examples/advection/var-velocity.py | 20 +++--- examples/quad/flux_quad.py | 91 --------------------------- examples/quad/mass_quad.py | 90 --------------------------- examples/quad/quad.py | 80 ------------------------ examples/quad/stiffness_quad.py | 99 ------------------------------ grudge/discretization.py | 2 +- grudge/execution.py | 20 +++--- grudge/models/advection.py | 15 +++-- grudge/symbolic/operators.py | 46 ++++---------- grudge/symbolic/primitives.py | 6 +- 10 files changed, 40 insertions(+), 429 deletions(-) delete mode 100644 examples/quad/flux_quad.py delete mode 100644 examples/quad/mass_quad.py delete mode 100644 examples/quad/quad.py delete mode 100644 examples/quad/stiffness_quad.py diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 97e793fd..80d656c4 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -30,8 +30,6 @@ logger = logging.getLogger(__name__) from grudge import sym, bind, DGDiscretizationWithBoundaries from pytools.obj_array import join_fields -import meshmode.discretization.poly_element as poly_element - def main(write_output=True, order=4): cl_ctx = cl.create_some_context() @@ -45,15 +43,12 @@ def main(write_output=True, order=4): dt_factor = 5 h = 1/10 - final_time = 5 - - sym_x = sym.nodes(2) advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2 - flux_type = "central" + flux_type = "upwind" source_center = np.array([0.1, 0.1]) source_width = 0.05 @@ -71,11 +66,12 @@ def main(write_output=True, order=4): from grudge.models.advection import VariableCoefficientAdvectionOperator - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2*order}) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, + order=order, quad_min_degrees={"product": 2*order}) + op = VariableCoefficientAdvectionOperator(2, advec_v, - u_analytic(sym.nodes(dim, sym.BTAG_ALL)),quad_tag="product", + u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product", flux_type=flux_type) - #op = sym.interp("vol", sym.DOFDesc("vol", "product"))(sym.var("u")) bound_op = bind(discr, op.sym_operator()) @@ -84,7 +80,7 @@ def main(write_output=True, order=4): def rhs(t, u): return bound_op(queue, t=t, u=u) - + final_time = 5 dt = dt_factor * h/order**2 nsteps = (final_time // dt) + 1 dt = final_time/nsteps + 1e-15 @@ -101,9 +97,9 @@ def main(write_output=True, order=4): if isinstance(event, dt_stepper.StateComputed): step += 1 - print(step) + if step % 10 == 0: + print(step) - #print(step, event.t, norm(queue, u=event.state_component[0])) vis.write_vtk_file("fld-%04d.vtu" % step, [("u", event.state_component)]) diff --git a/examples/quad/flux_quad.py b/examples/quad/flux_quad.py deleted file mode 100644 index 3d7f9a32..00000000 --- a/examples/quad/flux_quad.py +++ /dev/null @@ -1,91 +0,0 @@ -# Copyright (C) 2007 Andreas Kloeckner -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - - -import numpy as np -import pyopencl as cl # noqa -import pyopencl.array # noqa -import pyopencl.clmath # noqa - -import pytest # noqa - -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) - -import logging -logger = logging.getLogger(__name__) - -from grudge import sym, bind, DGDiscretizationWithBoundaries - -import numpy.linalg as la - - - - -def main(order=1): - cl_ctx = cl.create_some_context() - queue = cl.CommandQueue(cl_ctx) - - dim = 2 - - - from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh - #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) - mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), - n=(3, 3), order=7) - - - - - pquad_dd = sym.DOFDesc("vol", "product") - - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) - - - magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL, "product") - #magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL) - def fluxy_stuff(u): - return sym.interp(u.dd, magic_thing)(u.int) - - - ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0) - nds = sym.nodes(dim) - op = sym.FaceMassOperator(magic_thing, "vol")(fluxy_stuff(sym.int_tpair(nds, "product"))) - #op = (sym.int_tpair(nds, "product").int) - #op = sym.FaceMassOperator()(fluxy_stuff(sym.int_tpair(nds))) - print(op) - - print(sym.pretty(op)) - u = bind(discr, op)(queue, t=0) - - - print(u) - - - - - #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ]) - #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ]) - - - - - - - -if __name__ == "__main__": - main(4) - - diff --git a/examples/quad/mass_quad.py b/examples/quad/mass_quad.py deleted file mode 100644 index 05c306bd..00000000 --- a/examples/quad/mass_quad.py +++ /dev/null @@ -1,90 +0,0 @@ -# Copyright (C) 2007 Andreas Kloeckner -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - - -import numpy as np -import pyopencl as cl # noqa -import pyopencl.array # noqa -import pyopencl.clmath # noqa - -import pytest # noqa - -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) - -import logging -logger = logging.getLogger(__name__) - -from grudge import sym, bind, DGDiscretizationWithBoundaries - -import numpy.linalg as la - - - - -def main(order=1): - cl_ctx = cl.create_some_context() - queue = cl.CommandQueue(cl_ctx) - - dim = 2 - - - from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh - #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) - mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), - n=(3, 3), order=7) - - - - - pquad_dd = sym.DOFDesc("vol", "product") - to_pquad = sym.interp("vol", pquad_dd) - - def f(x): - x = x[0] / x[0] - #return sym.MassOperator()(x) - return sym.MassOperator(pquad_dd, "vol")(to_pquad(x)) - #return np.dot(from_pquad_stiffness_t, to_pquad(x)) - #return sym.interp("vol", pquad_dd)(sym.sin(3*x[0])) - - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) - - - - - - - - ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0) - u = bind(discr, f(sym.nodes(dim)))(queue, t=0) - - print(np.dot(ones.T, u)) - - - - - #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ]) - #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ]) - - - - - - - -if __name__ == "__main__": - main(4) - - diff --git a/examples/quad/quad.py b/examples/quad/quad.py deleted file mode 100644 index 67c41fa2..00000000 --- a/examples/quad/quad.py +++ /dev/null @@ -1,80 +0,0 @@ -# Copyright (C) 2007 Andreas Kloeckner -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - - -import numpy as np -import pyopencl as cl # noqa -import pyopencl.array # noqa -import pyopencl.clmath # noqa - -import pytest # noqa - -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) - -import logging -logger = logging.getLogger(__name__) - -from grudge import sym, bind, DGDiscretizationWithBoundaries - -import numpy.linalg as la - - - - -def main(order=1): - cl_ctx = cl.create_some_context() - queue = cl.CommandQueue(cl_ctx) - - dim = 2 - - - from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh - #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) - mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), - n=(3, 3), order=7) - - - - - pquad_dd = sym.DOFDesc("vol", "product") - - def f(x): - #return sym.sin(3*x[0]) - return sym.interp("vol", pquad_dd)(sym.sin(3*x[0])) - - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) - - u = bind(discr, f(sym.nodes(dim)))(queue, t=0) - - - - from meshmode.discretization.visualization import make_visualizer - #vis = make_visualizer(queue, discr.quad_volume_discr['product'], 2*order) - #vis = make_visualizer(queue, discr.volume_discr, order) - - #vis.write_vtk_file("fld-000.vtu", [ ("u", u) ]) - print(u) - - - - - - - -if __name__ == "__main__": - main(4) - - diff --git a/examples/quad/stiffness_quad.py b/examples/quad/stiffness_quad.py deleted file mode 100644 index ff334c4e..00000000 --- a/examples/quad/stiffness_quad.py +++ /dev/null @@ -1,99 +0,0 @@ -# Copyright (C) 2007 Andreas Kloeckner -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - - -import numpy as np -import pyopencl as cl # noqa -import pyopencl.array # noqa -import pyopencl.clmath # noqa - -import pytest # noqa - -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) - -import logging -logger = logging.getLogger(__name__) - -from grudge import sym, bind, DGDiscretizationWithBoundaries - -import numpy.linalg as la - - - - -def main(order=1): - cl_ctx = cl.create_some_context() - queue = cl.CommandQueue(cl_ctx) - - dim = 2 - - - from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh - #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9])) - mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), - n=(3, 3), order=7) - - - - - pquad_dd = sym.DOFDesc("vol", "product") - to_pquad = sym.interp("vol", pquad_dd) - from_pquad_stiffness_t = sym.stiffness_t(dim, pquad_dd, "vol") - from_pquad_mass = sym.MassOperator() - - def lin(x): - return x[0] - - def non_lin(x): - return x[1] * x[0] * x[0] + 3* x[0] * x[0] - 4 - - def f(x, fun): - #return sym.MassOperator()(x) - #return np.dot(sym.stiffness_t(2), to_pquad(fun(x))) - return np.dot(from_pquad_stiffness_t, fun(to_pquad(x))) - - def reg(x, fun): - #return sym.MassOperator()(x) - #return np.dot(sym.stiffness_t(2), to_pquad(fun(x))) - return np.dot(sym.stiffness_t(2), fun(x)) - - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} ) - - - #ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0) - nds = sym.nodes(dim) - #u = bind(discr, to_pquad(nds))(queue, t=0) - #u = bind(discr, reg(nds, lin))(queue, t=0) - u = bind(discr, f(nds, lin))(queue, t=0) - - print(u) - - - - - #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ]) - #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ]) - - - - - - - -if __name__ == "__main__": - main(4) - - diff --git a/grudge/discretization.py b/grudge/discretization.py index 2babc3a4..6f5cae13 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -1,6 +1,6 @@ from __future__ import division, absolute_import -__copyright__ = "Copyright (C) 2015 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy diff --git a/grudge/execution.py b/grudge/execution.py index b950a961..3fd47c3c 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -1,6 +1,6 @@ from __future__ import division, absolute_import -__copyright__ = "Copyright (C) 2015 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -206,7 +206,6 @@ class ExecutionMapper(mappers.Evaluator, knl = lp.split_iname(knl, "i", 16, inner_tag="l.0") return lp.tag_inames(knl, dict(k="g.0")) - in_discr = self.discrwb.discr_from_dd(op.dd_in) out_discr = self.discrwb.discr_from_dd(op.dd_out) @@ -217,24 +216,22 @@ class ExecutionMapper(mappers.Evaluator, for i in range(len(out_discr.groups)): in_grp = in_discr.groups[i] out_grp = out_discr.groups[i] - # FIXME: This cache thing + # FIXME: This cache thing #cache_key = "elwise_linear", grp, op, field.dtype #try: - #matrix = self.bound_op.operator_data_cache[cache_key] + #matrix = self.bound_op.operator_data_cache[cache_key] #except KeyError: matrix = ( - cl.array.to_device( - self.queue, - np.asarray(op.matrix(out_grp, in_grp), dtype=field.dtype)) - .with_queue(None)) + cl.array.to_device( + self.queue, + np.asarray(op.matrix(out_grp, in_grp), dtype=field.dtype)) + .with_queue(None)) - #self.bound_op.operator_data_cache[cache_key] = matrix + #self.bound_op.operator_data_cache[cache_key] = matrix knl()(self.queue, mat=matrix, result=out_grp.view(result), vec=in_grp.view(field)) - - return result def map_elementwise_max(self, op, field_expr): @@ -366,7 +363,6 @@ class ExecutionMapper(mappers.Evaluator, # This should be unified with map_elementwise_linear, which should # be extended to support batching. - # FIXME: Enable # assert repr_op.dd_in == repr_op.dd_out assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 4d678832..8d18b032 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -3,7 +3,7 @@ from __future__ import division, absolute_import -__copyright__ = "Copyright (C) 2009 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2009-2017 Andreas Kloeckner, Bogdan Enache" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -137,13 +137,14 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): self.v = v self.inflow_u = inflow_u self.flux_type = flux_type + if quad_tag is None: self.quad_tag = sym.QTAG_NONE else: self.quad_tag = quad_tag def flux(self, u): - + normal = sym.normal(u. dd, self.ambient_dim) surf_v = sym.interp("vol", u.dd)(self.v) @@ -173,6 +174,8 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): bc_in = self.inflow_u all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) + boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) + def flux(pair): return sym.interp(pair.dd, all_faces_dd)( self.flux(pair)) @@ -180,17 +183,19 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): quad_dd = sym.DOFDesc("vol", self.quad_tag) if self.quad_tag == sym.QTAG_NONE: - to_quad = lambda x : x + def to_quad(expr): + return expr else: to_quad = sym.interp("vol", quad_dd) from_quad_stiffness_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol") + quad_prod = (to_quad(self.v)*to_quad(u)) return sym.InverseMassOperator()( - np.dot(from_quad_stiffness_t, to_quad(self.v)*to_quad(u)) + np.dot(from_quad_stiffness_t, quad_prod) - sym.FaceMassOperator(all_faces_dd, "vol")( flux(sym.int_tpair(u, self.quad_tag)) - + flux(sym.bv_tpair(sym.DOFDesc(sym.BTAG_ALL, self.quad_tag), u, bc_in)) + + flux(sym.bv_tpair(boundary_dd, u, bc_in)) # FIXME: Add back support for inflow/outflow tags #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in)) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index ff909f7f..0ef8d722 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -2,7 +2,7 @@ from __future__ import division, absolute_import -__copyright__ = "Copyright (C) 2008 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -101,7 +101,6 @@ class InterpolationOperator(Operator): official_dd_in = _sym().as_dofdesc(dd_in) official_dd_out = _sym().as_dofdesc(dd_out) if official_dd_in == official_dd_out: - from warnings import warn raise ValueError("Interpolating from {} to {}" " does not do anything.".format(official_dd_in, official_dd_out)) @@ -152,7 +151,7 @@ class DiffOperatorBase(Operator): dd_in = _sym().DD_VOLUME if dd_out is None: dd_out = dd_in.with_qtag(_sym().QTAG_NONE) - if dd_out == "vol": #FIXME I'm pretty sure doing it this way is wrong + if dd_out == "vol": # FIXME I'm pretty sure doing it this way is wrong dd_out = _sym().DD_VOLUME if dd_out.uses_quadrature(): raise ValueError("differentiation outputs are not on " @@ -235,31 +234,21 @@ class RefStiffnessTOperator(RefDiffOperatorBase): mapper_method = intern("map_ref_stiffness_t") @staticmethod - def matrices(out_element_group, in_element_group): + def matrices(out_elem_grp, in_elem_grp): #FIXME: Is this right? Do we expose weird bugs if we delete this? - if in_element_group == out_element_group: - assert in_element_group.is_orthogonal_basis() - mmat = in_element_group.mass_matrix() - return [dmat.T.dot(mmat.T) for dmat in in_element_group.diff_matrices()] - - n_quad_nodes = in_element_group.nunit_nodes - n_reg_nodes = out_element_group.nunit_nodes + if in_elem_grp == out_elem_grp: + assert in_elem_grp.is_orthogonal_basis() + mmat = in_elem_grp.mass_matrix() + return [dmat.T.dot(mmat.T) for dmat in in_elem_grp.diff_matrices()] from modepy import vandermonde - vand = vandermonde(out_element_group.basis(), out_element_group.unit_nodes) - grad_vand = vandermonde(out_element_group.grad_basis(), in_element_group.unit_nodes) + vand = vandermonde(out_elem_grp.basis(), out_elem_grp.unit_nodes) + grad_vand = vandermonde(out_elem_grp.grad_basis(), in_elem_grp.unit_nodes) vand_inv_T = np.linalg.inv(vand).T - grad_basis = out_element_group.grad_basis() - weights = in_element_group.weights - nodes = in_element_group.unit_nodes + weights = in_elem_grp.weights - f_res = [np.zeros((n_reg_nodes, n_quad_nodes)) for d in range(out_element_group.dim)] - for d in range(out_element_group.dim): - for i in range(n_reg_nodes): - for j in range(n_quad_nodes): - f_res[d][i,j] = weights[j] * np.dot(vand_inv_T[i,:], grad_vand[d][j,:]) - return f_res + return np.einsum('c,bz,acz->abc', weights, vand_inv_T, grad_vand) # }}} @@ -394,25 +383,14 @@ class RefMassOperator(RefMassOperatorBase): if out_element_group == in_element_group: return in_element_group.mass_matrix() - n_quad_nodes = in_element_group.nunit_nodes - n_reg_nodes = out_element_group.nunit_nodes - from modepy import vandermonde vand = vandermonde(out_element_group.basis(), out_element_group.unit_nodes) o_vand = vandermonde(out_element_group.basis(), in_element_group.unit_nodes) vand_inv_T = np.linalg.inv(vand).T - grad_basis = out_element_group.grad_basis() weights = in_element_group.weights - nodes = in_element_group.unit_nodes - - f_res = np.zeros((n_reg_nodes, n_quad_nodes)) - for i in range(n_reg_nodes): - for j in range(n_quad_nodes): - f_res[i,j] = weights[j] * np.dot(vand_inv_T[i,:], o_vand[j,:]) - - return f_res + return np.einsum('c,bz,acz->abc', weights, vand_inv_T, o_vand) mapper_method = intern("map_ref_mass") diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 9234c225..00f7b8d3 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -2,7 +2,7 @@ from __future__ import division, absolute_import -__copyright__ = "Copyright (C) 2008 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -666,10 +666,6 @@ def int_tpair(expression, qtag=None): return TracePair(q_dd, i, e) - #i = cse(_sym().interp("vol", "int_faces")(expression)) - #e = cse(_sym().OppositeInteriorFaceSwap()(i)) - #return TracePair("int_faces", i, e) - def bdry_tpair(dd, interior, exterior): """ -- GitLab From a5292ed997d436958c8cc2764d375471734a453d Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 6 Dec 2017 19:49:53 -0600 Subject: [PATCH 09/30] Removed uppercase T in variable name --- grudge/symbolic/operators.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 0ef8d722..2f8642fa 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -244,11 +244,11 @@ class RefStiffnessTOperator(RefDiffOperatorBase): from modepy import vandermonde vand = vandermonde(out_elem_grp.basis(), out_elem_grp.unit_nodes) grad_vand = vandermonde(out_elem_grp.grad_basis(), in_elem_grp.unit_nodes) - vand_inv_T = np.linalg.inv(vand).T + vand_inv_t = np.linalg.inv(vand).T weights = in_elem_grp.weights - return np.einsum('c,bz,acz->abc', weights, vand_inv_T, grad_vand) + return np.einsum('c,bz,acz->abc', weights, vand_inv_t, grad_vand) # }}} @@ -386,11 +386,11 @@ class RefMassOperator(RefMassOperatorBase): from modepy import vandermonde vand = vandermonde(out_element_group.basis(), out_element_group.unit_nodes) o_vand = vandermonde(out_element_group.basis(), in_element_group.unit_nodes) - vand_inv_T = np.linalg.inv(vand).T + vand_inv_t = np.linalg.inv(vand).T weights = in_element_group.weights - return np.einsum('c,bz,acz->abc', weights, vand_inv_T, o_vand) + return np.einsum('c,bz,acz->abc', weights, vand_inv_t, o_vand) mapper_method = intern("map_ref_mass") -- GitLab From c9008e9df10e77e26a5dffd8f40ebd820dce07e4 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 6 Dec 2017 20:07:42 -0600 Subject: [PATCH 10/30] Fixed call to super() --- grudge/symbolic/operators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 2f8642fa..eece5cbf 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -104,7 +104,7 @@ class InterpolationOperator(Operator): raise ValueError("Interpolating from {} to {}" " does not do anything.".format(official_dd_in, official_dd_out)) - super().__init__(dd_in, dd_out) + super(InterpolationOperator, self).__init__(dd_in, dd_out) mapper_method = intern("map_interpolation") -- GitLab From 291c1a25d7708dbb9c029510c7a542ba4ec2aa6a Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Thu, 7 Dec 2017 18:33:12 -0600 Subject: [PATCH 11/30] Fixed warning on var-velocity, added a convergance test file, made things use 4*order --- examples/advection/test.py | 111 +++++++++++++++++++++++++++++ examples/advection/var-velocity.py | 2 +- grudge/models/advection.py | 7 +- 3 files changed, 116 insertions(+), 4 deletions(-) create mode 100644 examples/advection/test.py diff --git a/examples/advection/test.py b/examples/advection/test.py new file mode 100644 index 00000000..5553de06 --- /dev/null +++ b/examples/advection/test.py @@ -0,0 +1,111 @@ +import pyopencl as cl +import numpy as np + +from grudge import sym, bind, DGDiscretizationWithBoundaries +from pytools.obj_array import join_fields + + +def test_convergence_maxwell(order, visualize=False): + """Test whether 3D maxwells actually converges""" + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + dims = 2 + ns = [8, 10, 12] + for n in ns: + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*dims, + b=(0.5,)*dims, + n=(n,)*dims, + order=order) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, + order=order, quad_min_degrees={"product": 4*order}) + + sym_nds = sym.nodes(dims) + advec_v = join_fields(-1*sym_nds[1], sym_nds[0]) + + def angular_vel(radius): + return 1/radius + + def gaussian_mode(): + time = sym.var("t") + source_center = np.array([sym.cos(time), sym.sin(time)]) * 0.1 + source_width = 0.05 + + sym_x = sym.nodes(2) + sym_source_center_dist = sym_x - source_center + + return sym.exp( + -np.dot(sym_source_center_dist, sym_source_center_dist) + / source_width**2) + + def u_analytic(x): + return 0 + + from grudge.models.advection import VariableCoefficientAdvectionOperator + op = VariableCoefficientAdvectionOperator(2, advec_v, + u_analytic(sym.nodes(dims, sym.BTAG_ALL)), quad_tag="product", + flux_type="central") + + nsteps = 40 + dt = 0.015 + + bound_op = bind(discr, op.sym_operator()) + analytic_sol = bind(discr, gaussian_mode()) + fields = analytic_sol(queue, t=0) + + def rhs(t, u): + return bound_op(queue, t=t, u=u) + + from grudge.shortcuts import set_up_rk4 + dt_stepper = set_up_rk4("u", dt, fields, rhs) + + print("dt=%g nsteps=%d" % (dt, nsteps)) + + norm = bind(discr, sym.norm(2, sym.var("u"))) + + #from grudge.shortcuts import make_visualizer + #vis = make_visualizer(discr, vis_order=order) + #vis.write_vtk_file("fld-start.vtu", + #[("u", fields)]) + + step = 0 + final_t = nsteps * dt + for event in dt_stepper.run(t_end=final_t): + if isinstance(event, dt_stepper.StateComputed): + assert event.component_id == "u" + esc = event.state_component + + step += 1 + #vis.write_vtk_file("fld-%04d.vtu" % step, + #[("u", esc)]) + + if step % 10 == 0: + print(step) + + #vis.write_vtk_file("fld-end.vtu", + #[("u", analytic_sol(queue, t=step * dt))]) + + #vis.write_vtk_file("fld-actual-end.vtu", + #[("u", esc)]) + + sol = analytic_sol(queue, t=step * dt) + vals = [norm(queue, u=(esc - sol)) / norm(queue, u=sol)] # noqa E501 + total_error = sum(vals) + print(total_error) + eoc_rec.add_data_point(1.0/n, total_error) + + print(eoc_rec.pretty_print(abscissa_label="h", + error_label="L2 Error")) + + #assert eoc_rec.order_estimate() > order + + +if __name__ == "__main__": + test_convergence_maxwell(4) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 80d656c4..625b559b 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -67,7 +67,7 @@ def main(write_output=True, order=4): from grudge.models.advection import VariableCoefficientAdvectionOperator discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, - order=order, quad_min_degrees={"product": 2*order}) + order=order, quad_min_degrees={"product": 4*order}) op = VariableCoefficientAdvectionOperator(2, advec_v, u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product", diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 8d18b032..cf11bd15 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -188,11 +188,12 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): else: to_quad = sym.interp("vol", quad_dd) - from_quad_stiffness_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol") - quad_prod = (to_quad(self.v)*to_quad(u)) + stiff_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol") + quad_v = to_quad(self.v) + quad_u = to_quad(u) return sym.InverseMassOperator()( - np.dot(from_quad_stiffness_t, quad_prod) + (stiff_t[0](quad_u * quad_v[0]) + stiff_t[1](quad_u * quad_v[1])) - sym.FaceMassOperator(all_faces_dd, "vol")( flux(sym.int_tpair(u, self.quad_tag)) + flux(sym.bv_tpair(boundary_dd, u, bc_in)) -- GitLab From 91f258417b39f828a7bd0f5b3cf7a7ab52a286bb Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Fri, 8 Dec 2017 13:54:26 -0600 Subject: [PATCH 12/30] Removed some dead code, make things cleaner, re-added a operator cache --- grudge/execution.py | 27 +++++++++------------------ grudge/symbolic/operators.py | 6 ++++-- 2 files changed, 13 insertions(+), 20 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index 3fd47c3c..e128d741 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -213,21 +213,19 @@ class ExecutionMapper(mappers.Evaluator, queue=self.queue, dtype=field.dtype, allocator=self.bound_op.allocator) - for i in range(len(out_discr.groups)): - in_grp = in_discr.groups[i] - out_grp = out_discr.groups[i] - # FIXME: This cache thing - #cache_key = "elwise_linear", grp, op, field.dtype - #try: - #matrix = self.bound_op.operator_data_cache[cache_key] - #except KeyError: - matrix = ( + for in_grp, out_grp in zip(in_discr.groups, out_discr.groups): + + cache_key = "elwise_linear", in_grp, out_grp, op, field.dtype + try: + matrix = self.bound_op.operator_data_cache[cache_key] + except KeyError: + matrix = ( cl.array.to_device( self.queue, np.asarray(op.matrix(out_grp, in_grp), dtype=field.dtype)) .with_queue(None)) - #self.bound_op.operator_data_cache[cache_key] = matrix + self.bound_op.operator_data_cache[cache_key] = matrix knl()(self.queue, mat=matrix, result=out_grp.view(result), vec=in_grp.view(field)) @@ -249,9 +247,6 @@ class ExecutionMapper(mappers.Evaluator, return conn(self.queue, self.rec(field_expr)).with_queue(self.queue) def map_opposite_interior_face_swap(self, op, field_expr): - #if op.dd_in.quadrature_tag is not sym.QTAG_NONE: - #raise ValueError("face swap unavailable on quadrature meshes") - return self.discrwb.opposite_face_connection()( self.queue, self.rec(field_expr)).with_queue(self.queue) @@ -363,8 +358,6 @@ class ExecutionMapper(mappers.Evaluator, # This should be unified with map_elementwise_linear, which should # be extended to support batching. - # FIXME: Enable - # assert repr_op.dd_in == repr_op.dd_out assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag @memoize_in(self.discrwb, "reference_derivative_knl") @@ -394,9 +387,7 @@ class ExecutionMapper(mappers.Evaluator, dtype=field.dtype, extra_dims=(noperators,), allocator=self.bound_op.allocator) - for i in range(len(out_discr.groups)): - in_grp = in_discr.groups[i] - out_grp = out_discr.groups[i] + for in_grp, out_grp in zip(in_discr.groups, out_discr.groups): if in_grp.nelements == 0: continue diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index eece5cbf..0d579a32 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -149,10 +149,12 @@ class DiffOperatorBase(Operator): def __init__(self, xyz_axis, dd_in=None, dd_out=None): if dd_in is None: dd_in = _sym().DD_VOLUME + if dd_out is None: dd_out = dd_in.with_qtag(_sym().QTAG_NONE) - if dd_out == "vol": # FIXME I'm pretty sure doing it this way is wrong - dd_out = _sym().DD_VOLUME + else: + dd_out = _sym().as_dofdesc(dd_out) + if dd_out.uses_quadrature(): raise ValueError("differentiation outputs are not on " "quadrature grids") -- GitLab From d3b7322fb214b74d88042965de803bae27131782 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Fri, 8 Dec 2017 14:34:12 -0600 Subject: [PATCH 13/30] Removed FIXMEs that were fixed --- grudge/symbolic/operators.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 0d579a32..217fcf02 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -237,7 +237,6 @@ class RefStiffnessTOperator(RefDiffOperatorBase): @staticmethod def matrices(out_elem_grp, in_elem_grp): - #FIXME: Is this right? Do we expose weird bugs if we delete this? if in_elem_grp == out_elem_grp: assert in_elem_grp.is_orthogonal_basis() mmat = in_elem_grp.mass_matrix() @@ -381,7 +380,6 @@ class RefMassOperatorBase(ElementwiseLinearOperator): class RefMassOperator(RefMassOperatorBase): @staticmethod def matrix(out_element_group, in_element_group): - #FIXME: Is this normal? Will bugs surface if I remove this? if out_element_group == in_element_group: return in_element_group.mass_matrix() -- GitLab From d78fe7276c38c8389398d32be9945cdaa8fe4394 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sat, 13 Jan 2018 11:59:31 -0600 Subject: [PATCH 14/30] Moved test into pytest framework --- examples/advection/test.py | 111 ------------------------------------- test/test_grudge.py | 65 +++++++++++++++++++++- 2 files changed, 63 insertions(+), 113 deletions(-) delete mode 100644 examples/advection/test.py diff --git a/examples/advection/test.py b/examples/advection/test.py deleted file mode 100644 index 5553de06..00000000 --- a/examples/advection/test.py +++ /dev/null @@ -1,111 +0,0 @@ -import pyopencl as cl -import numpy as np - -from grudge import sym, bind, DGDiscretizationWithBoundaries -from pytools.obj_array import join_fields - - -def test_convergence_maxwell(order, visualize=False): - """Test whether 3D maxwells actually converges""" - - cl_ctx = cl.create_some_context() - queue = cl.CommandQueue(cl_ctx) - - from pytools.convergence import EOCRecorder - eoc_rec = EOCRecorder() - - dims = 2 - ns = [8, 10, 12] - for n in ns: - from meshmode.mesh.generation import generate_regular_rect_mesh - mesh = generate_regular_rect_mesh( - a=(-0.5,)*dims, - b=(0.5,)*dims, - n=(n,)*dims, - order=order) - - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, - order=order, quad_min_degrees={"product": 4*order}) - - sym_nds = sym.nodes(dims) - advec_v = join_fields(-1*sym_nds[1], sym_nds[0]) - - def angular_vel(radius): - return 1/radius - - def gaussian_mode(): - time = sym.var("t") - source_center = np.array([sym.cos(time), sym.sin(time)]) * 0.1 - source_width = 0.05 - - sym_x = sym.nodes(2) - sym_source_center_dist = sym_x - source_center - - return sym.exp( - -np.dot(sym_source_center_dist, sym_source_center_dist) - / source_width**2) - - def u_analytic(x): - return 0 - - from grudge.models.advection import VariableCoefficientAdvectionOperator - op = VariableCoefficientAdvectionOperator(2, advec_v, - u_analytic(sym.nodes(dims, sym.BTAG_ALL)), quad_tag="product", - flux_type="central") - - nsteps = 40 - dt = 0.015 - - bound_op = bind(discr, op.sym_operator()) - analytic_sol = bind(discr, gaussian_mode()) - fields = analytic_sol(queue, t=0) - - def rhs(t, u): - return bound_op(queue, t=t, u=u) - - from grudge.shortcuts import set_up_rk4 - dt_stepper = set_up_rk4("u", dt, fields, rhs) - - print("dt=%g nsteps=%d" % (dt, nsteps)) - - norm = bind(discr, sym.norm(2, sym.var("u"))) - - #from grudge.shortcuts import make_visualizer - #vis = make_visualizer(discr, vis_order=order) - #vis.write_vtk_file("fld-start.vtu", - #[("u", fields)]) - - step = 0 - final_t = nsteps * dt - for event in dt_stepper.run(t_end=final_t): - if isinstance(event, dt_stepper.StateComputed): - assert event.component_id == "u" - esc = event.state_component - - step += 1 - #vis.write_vtk_file("fld-%04d.vtu" % step, - #[("u", esc)]) - - if step % 10 == 0: - print(step) - - #vis.write_vtk_file("fld-end.vtu", - #[("u", analytic_sol(queue, t=step * dt))]) - - #vis.write_vtk_file("fld-actual-end.vtu", - #[("u", esc)]) - - sol = analytic_sol(queue, t=step * dt) - vals = [norm(queue, u=(esc - sol)) / norm(queue, u=sol)] # noqa E501 - total_error = sum(vals) - print(total_error) - eoc_rec.add_data_point(1.0/n, total_error) - - print(eoc_rec.pretty_print(abscissa_label="h", - error_label="L2 Error")) - - #assert eoc_rec.order_estimate() > order - - -if __name__ == "__main__": - test_convergence_maxwell(4) diff --git a/test/test_grudge.py b/test/test_grudge.py index 84c46168..66c786b3 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -334,8 +334,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type @pytest.mark.parametrize("order", [3, 4, 5]) -# test: 'test_convergence_advec(cl._csc, "disk", [0.1, 0.05], "strong", "upwind", 3)' -def test_convergence_maxwell(ctx_factory, order, visualize=False): +def test_convergence_maxwell(ctx_factory, order): """Test whether 3D maxwells actually converges""" cl_ctx = cl.create_some_context() @@ -403,6 +402,68 @@ def test_convergence_maxwell(ctx_factory, order, visualize=False): assert eoc_rec.order_estimate() > order +@pytest.mark.parametrize("order", [2, 3, 4]) +def test_improvement_quadrature(ctx_factory, order): + """Test whether quadrature improves things and converges""" + from meshmode.mesh.generation import generate_regular_rect_mesh + from grudge.models.advection import VariableCoefficientAdvectionOperator + from pytools.convergence import EOCRecorder + from pytools.obj_array import join_fields + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dims = 2 + sym_nds = sym.nodes(dims) + advec_v = join_fields(-1*sym_nds[1], sym_nds[0]) + + flux = "upwind" + op = VariableCoefficientAdvectionOperator(2, advec_v, + 0, quad_tag=None, + flux_type=flux) + q_op = VariableCoefficientAdvectionOperator(2, advec_v, + 0, quad_tag="product", + flux_type=flux) + + def gaussian_mode(): + source_width = 0.1 + sym_x = sym.nodes(2) + return sym.exp(-np.dot(sym_x, sym_x) / source_width**2) + + def conv_test(cl_ctx, queue, order, p_op): + + eoc_rec = EOCRecorder() + + ns = [20, 25] + for n in ns: + mesh = generate_regular_rect_mesh( + a=(-0.5,)*dims, + b=(0.5,)*dims, + n=(n,)*dims, + order=order) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, + order=order, quad_min_degrees={"product": 4*order}) + + bound_op = bind(discr, p_op.sym_operator()) + fields = bind(discr, gaussian_mode())(queue, t=0) + norm = bind(discr, sym.norm(2, sym.var("u"))) + + esc = bound_op(queue, u=fields) + total_error = norm(queue, u=esc) + eoc_rec.add_data_point(1.0/n, total_error) + + print(eoc_rec.pretty_print(abscissa_label="h", error_label="LInf Error")) + + return eoc_rec.order_estimate(), np.array([x[1] for x in eoc_rec.history]) + + eoc, errs = conv_test(cl_ctx, queue, order, op) + q_eoc, q_errs = conv_test(cl_ctx, queue, order, q_op) + assert q_eoc > eoc + assert (q_errs < errs).all() + assert q_eoc > order + + def test_foreign_points(ctx_factory): pytest.importorskip("sumpy") import sumpy.point_calculus as pc -- GitLab From 9ae6b17327d95feaba244c25458dcc4ea8a3ab29 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sun, 14 Jan 2018 15:02:07 -0600 Subject: [PATCH 15/30] changed quad min degrees to be more accurate --- examples/advection/var-velocity.py | 7 +++++-- grudge/discretization.py | 14 ++++++-------- test/test_grudge.py | 7 +++++-- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 625b559b..d01f4fd8 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -65,9 +65,12 @@ def main(write_output=True, order=4): return 0 from grudge.models.advection import VariableCoefficientAdvectionOperator + from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, - order=order, quad_min_degrees={"product": 4*order}) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, + quad_group_factory={ + "product": QuadratureSimplexGroupFactory(order=4*order) + }) op = VariableCoefficientAdvectionOperator(2, advec_v, u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product", diff --git a/grudge/discretization.py b/grudge/discretization.py index 6f5cae13..ccd9a5d1 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -47,17 +47,17 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): .. automethod :: zeros """ - def __init__(self, cl_ctx, mesh, order, quad_min_degrees=None): + def __init__(self, cl_ctx, mesh, order, quad_group_factory=None): """ :param quad_min_degrees: A mapping from quadrature tags to the degrees to which the desired quadrature is supposed to be exact. """ - if quad_min_degrees is None: - quad_min_degrees = {} + if quad_group_factory is None: + quad_group_factory = {} self.order = order - self.quad_min_degrees = quad_min_degrees + self.quad_group_factory = quad_group_factory from meshmode.discretization import Discretization @@ -203,12 +203,10 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): quadrature_tag = sym.QTAG_NONE from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory, \ - QuadratureSimplexGroupFactory + PolynomialWarpAndBlendGroupFactory if quadrature_tag is not sym.QTAG_NONE: - return QuadratureSimplexGroupFactory( - order=self.quad_min_degrees[quadrature_tag]) + return self.quad_group_factory[quadrature_tag] else: return PolynomialWarpAndBlendGroupFactory(order=self.order) diff --git a/test/test_grudge.py b/test/test_grudge.py index 66c786b3..79a88df3 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -409,6 +409,7 @@ def test_improvement_quadrature(ctx_factory, order): from grudge.models.advection import VariableCoefficientAdvectionOperator from pytools.convergence import EOCRecorder from pytools.obj_array import join_fields + from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) @@ -442,8 +443,10 @@ def test_improvement_quadrature(ctx_factory, order): n=(n,)*dims, order=order) - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, - order=order, quad_min_degrees={"product": 4*order}) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, + quad_group_factory={ + "product": QuadratureSimplexGroupFactory(order=4*order) + }) bound_op = bind(discr, p_op.sym_operator()) fields = bind(discr, gaussian_mode())(queue, t=0) -- GitLab From 3bcf19cf572922fda8409a874e1bc8e8b782d384 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 16 Jan 2018 00:28:21 -0600 Subject: [PATCH 16/30] Rename quad_group_factory to quad_tag_to_group_factory --- examples/advection/var-velocity.py | 2 +- grudge/discretization.py | 10 +++++----- test/test_grudge.py | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index d01f4fd8..e19d6aeb 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -68,7 +68,7 @@ def main(write_output=True, order=4): from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, - quad_group_factory={ + quad_tag_to_group_factory={ "product": QuadratureSimplexGroupFactory(order=4*order) }) diff --git a/grudge/discretization.py b/grudge/discretization.py index ccd9a5d1..8d80894f 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -47,17 +47,17 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): .. automethod :: zeros """ - def __init__(self, cl_ctx, mesh, order, quad_group_factory=None): + def __init__(self, cl_ctx, mesh, order, quad_tag_to_group_factory=None): """ :param quad_min_degrees: A mapping from quadrature tags to the degrees to which the desired quadrature is supposed to be exact. """ - if quad_group_factory is None: - quad_group_factory = {} + if quad_tag_to_group_factory is None: + quad_tag_to_group_factory = {} self.order = order - self.quad_group_factory = quad_group_factory + self.quad_tag_to_group_factory = quad_tag_to_group_factory from meshmode.discretization import Discretization @@ -206,7 +206,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): PolynomialWarpAndBlendGroupFactory if quadrature_tag is not sym.QTAG_NONE: - return self.quad_group_factory[quadrature_tag] + return self.quad_tag_to_group_factory[quadrature_tag] else: return PolynomialWarpAndBlendGroupFactory(order=self.order) diff --git a/test/test_grudge.py b/test/test_grudge.py index 79a88df3..64863ac6 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -444,7 +444,7 @@ def test_improvement_quadrature(ctx_factory, order): order=order) discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, - quad_group_factory={ + quad_tag_to_group_factory={ "product": QuadratureSimplexGroupFactory(order=4*order) }) -- GitLab From 3fe692428c44127a0a3fae3536eb2c966616735d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 24 Jan 2018 12:26:14 -0600 Subject: [PATCH 17/30] Remove an extra fold marker --- grudge/symbolic/operators.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 7009cdd8..194b0716 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -346,8 +346,6 @@ class VandermondeOperator(ElementwiseLinearOperator): # }}} -# }}} - # {{{ mass operators -- GitLab From 76772747fdadcc372d01b21cc16442c66da3ef2c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 24 Jan 2018 12:27:18 -0600 Subject: [PATCH 18/30] Change dim variable names in loopy matvecs --- grudge/execution.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index e128d741..2915ee85 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -198,8 +198,8 @@ class ExecutionMapper(mappers.Evaluator, knl = lp.make_kernel( """{[k,i,j]: 0<=k Date: Wed, 24 Jan 2018 12:27:47 -0600 Subject: [PATCH 19/30] Add a missing fold section name --- grudge/symbolic/compiler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index f6e23808..c7bd6c22 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -792,7 +792,7 @@ def aggregate_assignments(inf_mapper, instructions, result, # }}} -# {{{ +# {{{ to-loopy mapper def set_once(d, k, v): try: -- GitLab From 670eb218276298a1063f52ed3ed14f2585192c74 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 24 Jan 2018 12:28:12 -0600 Subject: [PATCH 20/30] Fix a broken debug import --- grudge/symbolic/compiler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index c7bd6c22..c555cea0 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -338,7 +338,7 @@ class Code(object): self.static_schedule_attempts = 5 def dump_dataflow_graph(self): - from grudge.tools import open_unique_debug_file + from pytools.debug import open_unique_debug_file open_unique_debug_file("dataflow", ".dot")\ .write(dot_dataflow_graph(self, max_node_label_length=None)) -- GitLab From 3fc59a848e1a1e879e07e6809d295750ac994c13 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 24 Jan 2018 12:30:38 -0600 Subject: [PATCH 21/30] Fix symbolic operator dumping for debugging --- examples/advection/var-velocity.py | 5 ++++- grudge/execution.py | 15 ++++++++------- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index e19d6aeb..366fd1e4 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -76,7 +76,10 @@ def main(write_output=True, order=4): u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product", flux_type=flux_type) - bound_op = bind(discr, op.sym_operator()) + bound_op = bind( + discr, op.sym_operator(), + #debug_flags=["dump_sym_operator_stages"] + ) u = bind(discr, f(sym.nodes(dim)))(queue, t=0) diff --git a/grudge/execution.py b/grudge/execution.py index 2915ee85..04482c86 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -559,19 +559,20 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x, stage = [0] - def dump_optemplate(name, sym_operator): - if "dump_optemplate_stages" in debug_flags: - from grudge.tools import open_unique_debug_file - from grudge.optemplate import pretty - open_unique_debug_file("%02d-%s" % (stage[0], name), ".txt").write( - pretty(sym_operator)) + def dump_sym_operator(name, sym_operator): + if "dump_sym_operator_stages" in debug_flags: + from pytools.debug import open_unique_debug_file + outf, name = open_unique_debug_file("%02d-%s" % (stage[0], name), ".txt") + with outf: + outf.write(sym.pretty(sym_operator)) + stage[0] += 1 sym_operator = process_sym_operator( discr, sym_operator, post_bind_mapper=post_bind_mapper, - dumper=dump_optemplate) + dumper=dump_sym_operator) from grudge.symbolic.compiler import OperatorCompiler discr_code, eval_code = OperatorCompiler(discr)(sym_operator) -- GitLab From 053247fb60925bab3943680cbda8abb825416349 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 24 Jan 2018 12:31:13 -0600 Subject: [PATCH 22/30] Fix license header in var-velocity example --- examples/advection/var-velocity.py | 37 +++++++++++++++++++----------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 366fd1e4..3477d2c7 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -1,17 +1,26 @@ -# 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, absolute_import + +__copyright__ = "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. +""" import numpy as np -- GitLab From e742a62fba3114777651f071a36b052696d98153 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 24 Jan 2018 12:31:38 -0600 Subject: [PATCH 23/30] Clean up stale comments in sym operator processing --- grudge/execution.py | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index 04482c86..cb765eb2 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -505,35 +505,9 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, complex_type=discrwb.complex_dtype.type, )(sym_operator) - # Ordering restriction: - # - # - Must run constant fold before first type inference pass, because zeros, - # while allowed, violate typing constraints (because they can't be assigned - # a unique type), and need to be killed before the type inferrer sees them. - - # FIXME: Reenable type inference - - # from grudge.symbolic.mappers.type_inference import TypeInferrer - # dumper("before-specializer", sym_operator) - # sym_operator = mappers.OperatorSpecializer( - # TypeInferrer()(sym_operator) - # )(sym_operator) - - # Ordering restriction: - # - # - Must run OperatorSpecializer before performing the GlobalToReferenceMapper, - # because otherwise it won't differentiate the type of grids (node or quadrature - # grids) that the operators will apply on. - dumper("before-global-to-reference", sym_operator) sym_operator = mappers.GlobalToReferenceMapper(discrwb.ambient_dim)(sym_operator) - # Ordering restriction: - # - # - Must specialize quadrature operators before performing inverse mass - # contraction, because there are no inverse-mass-contracted variants of the - # quadrature operators. - dumper("before-imass", sym_operator) sym_operator = mappers.InverseMassContractor()(sym_operator) -- GitLab From 0257e912e13b3c762fdc51c8d08c6db495af49e6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 24 Jan 2018 12:33:44 -0600 Subject: [PATCH 24/30] Implement quad removal mapper --- grudge/execution.py | 4 + grudge/symbolic/mappers/__init__.py | 117 +++++++++++----------------- 2 files changed, 48 insertions(+), 73 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index cb765eb2..2da2e508 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -490,6 +490,10 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, dumper("before-cfold", sym_operator) sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator) + dumper("before-qcheck", sym_operator) + sym_operator = mappers.QuadratureCheckerAndRemover( + discrwb.quad_tag_to_group_factory)(sym_operator) + # Work around https://github.com/numpy/numpy/issues/9438 # # The idea is that we need 1j as an expression to survive diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 405cb5f0..ca11cbf0 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -717,96 +717,67 @@ class NoCSEStringifyMapper(StringifyMapper): # {{{ quadrature support -class QuadratureUpsamplerRemover(CSECachingMapperMixin, IdentityMapper): - def __init__(self, quad_min_degrees, do_warn=True): +class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper): + """Checks whether all quadratu + """ + + def __init__(self, quad_tag_to_group_factory): IdentityMapper.__init__(self) CSECachingMapperMixin.__init__(self) - self.quad_min_degrees = quad_min_degrees - self.do_warn = do_warn + self.quad_tag_to_group_factory = quad_tag_to_group_factory map_common_subexpression_uncached = \ IdentityMapper.map_common_subexpression - def map_operator_binding(self, expr): - if isinstance(expr.op, (op.QuadratureGridUpsampler, - op.QuadratureInteriorFacesGridUpsampler, - op.QuadratureBoundaryGridUpsampler)): - try: - min_degree = self.quad_min_degrees[expr.op.quadrature_tag] - except KeyError: - if self.do_warn: - from warnings import warn - warn("No minimum degree for quadrature tag '%s' specified--" - "falling back to nodal evaluation" - % expr.op.quadrature_tag) - return self.rec(expr.field) - else: - if min_degree is None: - return self.rec(expr.field) - else: - return expr.op(self.rec(expr.field)) - else: - return IdentityMapper.map_operator_binding(self, expr) + def _process_dd(self, dd, location_descr): + from grudge.symbolic.primitives import DOFDesc, QTAG_NONE + if dd.quadrature_tag is not QTAG_NONE: + if dd.quadrature_tag not in self.quad_tag_to_group_factory: + raise ValueError("found unknown quadrature tag '%s' in '%s'" + % (dd.quadrature_tag, location_descr)) + grp_factory = self.quad_tag_to_group_factory[dd.quadrature_tag] + if grp_factory is None: + dd = DOFDesc(dd.domain_tag, QTAG_NONE) -class QuadratureDetector(CSECachingMapperMixin, CombineMapper): - """For a given expression, this mapper returns the upsampling - operator in effect at the root of the expression, or *None* - if there isn't one. - """ - class QuadStatusNotKnown: - pass - - map_common_subexpression_uncached = \ - CombineMapper.map_common_subexpression + return dd - def combine(self, values): - from pytools import single_valued - return single_valued([ - v for v in values if v is not self.QuadStatusNotKnown]) + def map_operator_binding(self, expr): + dd_in_orig = dd_in = expr.op.dd_in + dd_out_orig = dd_out = expr.op.dd_out - def map_variable(self, expr): - return None + dd_in = self._process_dd(dd_in, "dd_in of %s" % type(expr.op).__name__) + dd_out = self._process_dd(dd_out, "dd_out of %s" % type(expr.op).__name__) - def map_constant(self, expr): - return self.QuadStatusNotKnown + if dd_in_orig == dd_in and dd_out_orig == dd_out: + # unchanged + return IdentityMapper.map_operator_binding(self, expr) - def map_operator_binding(self, expr): - if isinstance(expr.op, ( - op.DiffOperatorBase, op.FluxOperatorBase, - op.MassOperatorBase)): - return None - elif isinstance(expr.op, (op.QuadratureGridUpsampler, - op.QuadratureInteriorFacesGridUpsampler)): - return expr.op - else: - return CombineMapper.map_operator_binding(self, expr) + import grudge.symbolic.operators as op + # changed + if dd_in == dd_out and isinstance(expr.op, op.InterpolationOperator): + # This used to be to-quad interpolation and has become a no-op. + # Remove it. + return self.rec(expr.field) -class QuadratureUpsamplerChanger(CSECachingMapperMixin, IdentityMapper): - """This mapper descends the expression tree, down to each - quadrature-consuming operator (diff, mass) along each branch. - It then change - """ - def __init__(self, desired_quad_op): - IdentityMapper.__init__(self) - CSECachingMapperMixin.__init__(self) + if isinstance(expr.op, op.StiffnessTOperator): + new_op = type(expr.op)(expr.op.xyz_axis, dd_in, dd_out) + elif isinstance(expr.op, (op.FaceMassOperator, op.InterpolationOperator)): + new_op = type(expr.op)(dd_in, dd_out) + else: + raise NotImplementedError("do not know how to modify dd_in and dd_out " + "in %s" % type(expr.op).__name__) - self.desired_quad_op = desired_quad_op + return new_op(self.rec(expr.field)) - map_common_subexpression_uncached = \ - IdentityMapper.map_common_subexpression + def map_ones(self, expr): + dd = self._process_dd(expr.dd, location_descr=type(expr).__name__) + return type(expr)(dd) - def map_operator_binding(self, expr): - if isinstance(expr.op, ( - op.DiffOperatorBase, op.FluxOperatorBase, - op.MassOperatorBase)): - return expr - elif isinstance(expr.op, (op.QuadratureGridUpsampler, - op.QuadratureInteriorFacesGridUpsampler)): - return self.desired_quad_op(expr.field) - else: - return IdentityMapper.map_operator_binding(self, expr) + def map_node_coordinate_component(self, expr): + dd = self._process_dd(expr.dd, location_descr=type(expr).__name__) + return type(expr)(expr.axis, dd) # }}} -- GitLab From 53e576e64b558261197061f7b9b2cea0ad5b349f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 24 Jan 2018 12:34:07 -0600 Subject: [PATCH 25/30] Remove no-quad accomodations from var-vel advection model --- grudge/models/advection.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index cf11bd15..091b9870 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -138,10 +138,7 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): self.inflow_u = inflow_u self.flux_type = flux_type - if quad_tag is None: - self.quad_tag = sym.QTAG_NONE - else: - self.quad_tag = quad_tag + self.quad_tag = quad_tag def flux(self, u): @@ -182,11 +179,7 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): quad_dd = sym.DOFDesc("vol", self.quad_tag) - if self.quad_tag == sym.QTAG_NONE: - def to_quad(expr): - return expr - else: - to_quad = sym.interp("vol", quad_dd) + to_quad = sym.interp("vol", quad_dd) stiff_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol") quad_v = to_quad(self.v) -- GitLab From 9e0bbdf718406e7cea93a488b8ae330575f073ee Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 25 Jan 2018 00:39:55 -0600 Subject: [PATCH 26/30] Fix quad_tag_to_group_factory docstring --- grudge/discretization.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/grudge/discretization.py b/grudge/discretization.py index 8d80894f..04cf4330 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -49,8 +49,12 @@ class DGDiscretizationWithBoundaries(DiscretizationBase): def __init__(self, cl_ctx, mesh, order, quad_tag_to_group_factory=None): """ - :param quad_min_degrees: A mapping from quadrature tags to the degrees - to which the desired quadrature is supposed to be exact. + :param quad_min_degrees: A mapping from quadrature tags (typically + strings--but may be any hashable/comparable object) to a + :class:`meshmode.discretization.ElementGroupFactory` indicating with + which quadrature discretization the operations are to be carried out, + or *None* to indicate that operations with this quadrature tag should + be carried out with the standard volume discretization. """ if quad_tag_to_group_factory is None: -- GitLab From 1ec29819862a33d59299f6754f74db9c78d0187c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 25 Jan 2018 00:40:28 -0600 Subject: [PATCH 27/30] Tweak var-velocity example --- examples/advection/var-velocity.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 3477d2c7..caa5ff35 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -46,12 +46,13 @@ def main(write_output=True, order=4): dim = 2 + resolution = 10 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=(10, 10), order=order) + n=(resolution, resolution), order=order) dt_factor = 5 - h = 1/10 + h = 1/resolution sym_x = sym.nodes(2) @@ -74,10 +75,11 @@ def main(write_output=True, order=4): return 0 from grudge.models.advection import VariableCoefficientAdvectionOperator - from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory + from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory # noqa discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_tag_to_group_factory={ + #"product": None, "product": QuadratureSimplexGroupFactory(order=4*order) }) @@ -95,7 +97,7 @@ def main(write_output=True, order=4): def rhs(t, u): return bound_op(queue, t=t, u=u) - final_time = 5 + final_time = 50 dt = dt_factor * h/order**2 nsteps = (final_time // dt) + 1 dt = final_time/nsteps + 1e-15 @@ -115,8 +117,8 @@ def main(write_output=True, order=4): if step % 10 == 0: print(step) - vis.write_vtk_file("fld-%04d.vtu" % step, - [("u", event.state_component)]) + vis.write_vtk_file("fld-%04d.vtu" % step, + [("u", event.state_component)]) if __name__ == "__main__": -- GitLab From 96f01bc65212906e91b4a95ad490751e04ad4ecd Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 25 Jan 2018 11:26:37 -0600 Subject: [PATCH 28/30] Fix test_improvement_quadrature --- test/test_grudge.py | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 64863ac6..8f72cdbf 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -419,20 +419,17 @@ def test_improvement_quadrature(ctx_factory, order): advec_v = join_fields(-1*sym_nds[1], sym_nds[0]) flux = "upwind" - op = VariableCoefficientAdvectionOperator(2, advec_v, - 0, quad_tag=None, - flux_type=flux) - q_op = VariableCoefficientAdvectionOperator(2, advec_v, - 0, quad_tag="product", - flux_type=flux) + op = VariableCoefficientAdvectionOperator(2, advec_v, 0, flux_type=flux) def gaussian_mode(): source_width = 0.1 sym_x = sym.nodes(2) return sym.exp(-np.dot(sym_x, sym_x) / source_width**2) - def conv_test(cl_ctx, queue, order, p_op): - + def conv_test(descr, use_quad): + print("-"*75) + print(descr) + print("-"*75) eoc_rec = EOCRecorder() ns = [20, 25] @@ -443,12 +440,17 @@ def test_improvement_quadrature(ctx_factory, order): n=(n,)*dims, order=order) + if use_quad: + quad_tag_to_group_factory = { + "product": QuadratureSimplexGroupFactory(order=4*order) + } + else: + quad_tag_to_group_factory = {"product": None} + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, - quad_tag_to_group_factory={ - "product": QuadratureSimplexGroupFactory(order=4*order) - }) + quad_tag_to_group_factory=quad_tag_to_group_factory) - bound_op = bind(discr, p_op.sym_operator()) + bound_op = bind(discr, op.sym_operator()) fields = bind(discr, gaussian_mode())(queue, t=0) norm = bind(discr, sym.norm(2, sym.var("u"))) @@ -456,12 +458,15 @@ def test_improvement_quadrature(ctx_factory, order): total_error = norm(queue, u=esc) eoc_rec.add_data_point(1.0/n, total_error) + print("-"*75) + print(descr) print(eoc_rec.pretty_print(abscissa_label="h", error_label="LInf Error")) + print("-"*75) return eoc_rec.order_estimate(), np.array([x[1] for x in eoc_rec.history]) - eoc, errs = conv_test(cl_ctx, queue, order, op) - q_eoc, q_errs = conv_test(cl_ctx, queue, order, q_op) + eoc, errs = conv_test("no quadrature", False) + q_eoc, q_errs = conv_test("with quadrature", True) assert q_eoc > eoc assert (q_errs < errs).all() assert q_eoc > order -- GitLab From dccba991e9992f9f7221e1465604d6d10b87918f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 25 Jan 2018 14:17:55 -0600 Subject: [PATCH 29/30] Tweak test_improvement_quadrature output --- test/test_grudge.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index 8f72cdbf..34532bfc 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -458,10 +458,7 @@ def test_improvement_quadrature(ctx_factory, order): total_error = norm(queue, u=esc) eoc_rec.add_data_point(1.0/n, total_error) - print("-"*75) - print(descr) print(eoc_rec.pretty_print(abscissa_label="h", error_label="LInf Error")) - print("-"*75) return eoc_rec.order_estimate(), np.array([x[1] for x in eoc_rec.history]) -- GitLab From fbf8c7e3428039705f5e7b9a040c26a19b10ae3e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 25 Jan 2018 15:51:55 -0600 Subject: [PATCH 30/30] Toy around with var-velocity example --- examples/advection/var-velocity.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index caa5ff35..00e4f0c8 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -46,7 +46,7 @@ def main(write_output=True, order=4): dim = 2 - resolution = 10 + resolution = 15 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=(resolution, resolution), order=order) @@ -66,11 +66,19 @@ def main(write_output=True, order=4): sym_x = sym.nodes(2) sym_source_center_dist = sym_x - source_center - def f(x): + def f_gaussian(x): return sym.exp( -np.dot(sym_source_center_dist, sym_source_center_dist) / source_width**2) + def f_step(x): + return sym.If( + sym.Comparison( + np.dot(sym_source_center_dist, sym_source_center_dist), + "<", + (4*source_width)**2), + 1, 0) + def u_analytic(x): return 0 @@ -92,7 +100,7 @@ def main(write_output=True, order=4): #debug_flags=["dump_sym_operator_stages"] ) - u = bind(discr, f(sym.nodes(dim)))(queue, t=0) + u = bind(discr, f_gaussian(sym.nodes(dim)))(queue, t=0) def rhs(t, u): return bound_op(queue, t=t, u=u) @@ -106,7 +114,7 @@ def main(write_output=True, order=4): dt_stepper = set_up_rk4("u", dt, u, rhs) from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr, vis_order=order) + vis = make_visualizer(discr, vis_order=2*order) step = 0 @@ -114,7 +122,7 @@ def main(write_output=True, order=4): if isinstance(event, dt_stepper.StateComputed): step += 1 - if step % 10 == 0: + if step % 30 == 0: print(step) vis.write_vtk_file("fld-%04d.vtu" % step, -- GitLab