From 3e06ba6b0bdc0a2f73f6256e4a0e8cc1a7129a1f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 20 Apr 2020 19:25:15 -0500 Subject: [PATCH 1/9] clean up advection example and doing it in 1d --- examples/advection/weak.py | 127 +++++++++++++++++++++--------------- grudge/models/advection.py | 4 +- grudge/symbolic/compiler.py | 5 ++ 3 files changed, 81 insertions(+), 55 deletions(-) diff --git a/examples/advection/weak.py b/examples/advection/weak.py index bd9d61ca..4d3736b1 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -13,106 +13,127 @@ # 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 numpy.linalg as la -import pytest # noqa +import pyopencl as cl +import pyopencl.array +import pyopencl.clmath -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from grudge import bind, sym import logging + logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO) -from grudge import sym, bind, DGDiscretizationWithBoundaries - -import numpy.linalg as la +def main(ctx_factory, dim=1, order=4, visualize=True): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + # {{{ parameters + # domain side [-d/2, d/2]^dim + d = 1.0 + # number of points in each dimension + npoints = 20 + # grid spacing + h = d / npoints + # cfl? + dt_factor = 4 + # final time + final_time = 1.0 -def main(write_output=True, order=4): - cl_ctx = cl.create_some_context() - queue = cl.CommandQueue(cl_ctx) + # velocity field + c = np.array([0.1] * dim) + norm_c = la.norm(c) + # flux + flux_type = "central" - dim = 2 + # compute number of steps + dt = dt_factor * h / order**2 + nsteps = int(final_time // dt) + 1 + dt = final_time/nsteps + 1.0e-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=(20, 20), order=order) + # {{{ discretization - dt_factor = 4 - h = 1/20 + from meshmode.mesh.generation import generate_box_mesh + mesh = generate_box_mesh( + [np.linspace(-d/2, d/2, npoints) for _ in range(dim)], + order=order) + from grudge import DGDiscretizationWithBoundaries discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) - c = np.array([0.1,0.1]) - norm_c = la.norm(c) - + # }}} - flux_type = "central" - + # {{{ solve advection def f(x): - return sym.sin(3*x) + return sym.sin(3 * x) def u_analytic(x): - return f(-np.dot(c, x)/norm_c+sym.var("t", sym.DD_SCALAR)*norm_c) + t = sym.var("t", sym.DD_SCALAR) + return f(-np.dot(c, x) / norm_c + t * norm_c) from grudge.models.advection import WeakAdvectionOperator - from meshmode.mesh import BTAG_ALL, BTAG_NONE - - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) op = WeakAdvectionOperator(c, inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)), flux_type=flux_type) bound_op = bind(discr, op.sym_operator()) - u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0) + return + def rhs(t, u): return bound_op(queue, t=t, u=u) - final_time = 0.3 - - dt = dt_factor * h/order**2 - nsteps = (final_time // dt) + 1 - dt = final_time/nsteps + 1e-15 - - from grudge.shortcuts import set_up_rk4 dt_stepper = set_up_rk4("u", dt, u, rhs) - last_u = None - - from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr, vis_order=order) + if dim == 1: + import matplotlib.pyplot as pt + pt.figure(figsize=(8, 8), dpi=300) + else: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) step = 0 - norm = bind(discr, sym.norm(2, sym.var("u"))) - for event in dt_stepper.run(t_end=final_time): - if isinstance(event, dt_stepper.StateComputed): - - step += 1 - - #print(step, event.t, norm(queue, u=event.state_component[0])) - vis.write_vtk_file("fld-weak-%04d.vtu" % step, - [ ("u", event.state_component) ]) - - + if not isinstance(event, dt_stepper.StateComputed): + continue + step += 1 + norm_u = norm(queue, u=event.state_component) + logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) + if dim == 1: + x = discr.discr_from_dd(sym.DTAG_VOLUME_ALL).nodes().get(queue) + u = event.state_component.get(queue) + pt.plot(x, u) + pt.xlabel("$x$") + pt.ylabel("$u$") + pt.title("t = {:.2f}".format(event.t)) + pt.savefig("fld-weak-%04d.png" % step) + pt.clf() + else: + vis.write_vtk_file("fld-weak-%04d.vtu" % step, [ + ("u", event.state_component) + ], overwrite=True) if __name__ == "__main__": - main() + import argparse + parser = argparse.ArgumentParser() + parser.add_argument("--dim", default=2, type=int) + args = parser.parse_args() + main(cl.create_some_context, + dim=args.dim) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 091b9870..c3b142e4 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -74,7 +74,7 @@ class AdvectionOperatorBase(HyperbolicOperator): class StrongAdvectionOperator(AdvectionOperatorBase): def flux(self, u): - normal = sym.normal(u. dd, self.ambient_dim) + normal = sym.normal(u.dd, self.ambient_dim) v_dot_normal = sym.cse(self.v.dot(normal), "v_dot_normal") return u.int * v_dot_normal - self.weak_flux(u) @@ -142,7 +142,7 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): def flux(self, u): - normal = sym.normal(u. dd, self.ambient_dim) + normal = sym.normal(u.dd, self.ambient_dim) surf_v = sym.interp("vol", u.dd)(self.v) diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index 51eb5202..4e13caaf 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -1047,6 +1047,11 @@ class ToLoopyInstructionMapper(object): self.insn_count += 1 + for expr in insn.exprs: + print(expr) + print(self.dd_inference_mapper(expr)) + print() + from pytools import single_valued governing_dd = single_valued( self.dd_inference_mapper(expr) -- GitLab From f2e31c38fb3a4fab9147e9e5c59b1d90158824a2 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 22 Apr 2020 14:25:03 -0500 Subject: [PATCH 2/9] add special cases for 1D/0D geometries --- examples/advection/weak.py | 48 +++++++++++++++++++++-------------- grudge/symbolic/compiler.py | 5 ---- grudge/symbolic/operators.py | 12 ++++++--- grudge/symbolic/primitives.py | 7 +++-- 4 files changed, 43 insertions(+), 29 deletions(-) diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 4d3736b1..68eeafe0 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -41,12 +41,12 @@ def main(ctx_factory, dim=1, order=4, visualize=True): # grid spacing h = d / npoints # cfl? - dt_factor = 4 + dt_factor = 2.0 # final time final_time = 1.0 # velocity field - c = np.array([0.1] * dim) + c = np.array([0.5] * dim) norm_c = la.norm(c) # flux flux_type = "central" @@ -68,6 +68,9 @@ def main(ctx_factory, dim=1, order=4, visualize=True): from grudge import DGDiscretizationWithBoundaries discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + faces_discr = discr.discr_from_dd(sym.FACE_RESTR_INTERIOR) + # }}} # {{{ solve advection @@ -75,8 +78,9 @@ def main(ctx_factory, dim=1, order=4, visualize=True): def f(x): return sym.sin(3 * x) - def u_analytic(x): - t = sym.var("t", sym.DD_SCALAR) + def u_analytic(x, t=None): + if t is None: + t = sym.var("t", sym.DD_SCALAR) return f(-np.dot(c, x) / norm_c + t * norm_c) from grudge.models.advection import WeakAdvectionOperator @@ -87,8 +91,6 @@ def main(ctx_factory, dim=1, order=4, visualize=True): bound_op = bind(discr, op.sym_operator()) u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0) - return - def rhs(t, u): return bound_op(queue, t=t, u=u) @@ -98,35 +100,43 @@ def main(ctx_factory, dim=1, order=4, visualize=True): if dim == 1: import matplotlib.pyplot as pt pt.figure(figsize=(8, 8), dpi=300) + + volume_x = volume_discr.nodes().get(queue) else: from grudge.shortcuts import make_visualizer vis = make_visualizer(discr, vis_order=order) - step = 0 - norm = bind(discr, sym.norm(2, sym.var("u"))) - for event in dt_stepper.run(t_end=final_time): - if not isinstance(event, dt_stepper.StateComputed): - continue - - step += 1 - norm_u = norm(queue, u=event.state_component) - logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) + def plot_solution(evt): + if not visualize: + return if dim == 1: - x = discr.discr_from_dd(sym.DTAG_VOLUME_ALL).nodes().get(queue) u = event.state_component.get(queue) + u_ = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=evt.t).get(queue) - pt.plot(x, u) + pt.plot(volume_x[0], u, 'o') + pt.plot(volume_x[0], u_, "k.") pt.xlabel("$x$") pt.ylabel("$u$") - pt.title("t = {:.2f}".format(event.t)) + pt.title("t = {:.2f}".format(evt.t)) pt.savefig("fld-weak-%04d.png" % step) pt.clf() else: vis.write_vtk_file("fld-weak-%04d.vtu" % step, [ - ("u", event.state_component) + ("u", evt.state_component) ], overwrite=True) + step = 0 + norm = bind(discr, sym.norm(2, sym.var("u"))) + for event in dt_stepper.run(t_end=final_time): + if not isinstance(event, dt_stepper.StateComputed): + continue + + step += 1 + norm_u = norm(queue, u=event.state_component) + logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) + plot_solution(event) + if __name__ == "__main__": import argparse diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index 4e13caaf..51eb5202 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -1047,11 +1047,6 @@ class ToLoopyInstructionMapper(object): self.insn_count += 1 - for expr in insn.exprs: - print(expr) - print(self.dd_inference_mapper(expr)) - print() - from pytools import single_valued governing_dd = single_valued( self.dd_inference_mapper(expr) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index a1899ed2..0cbb6eda 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -536,6 +536,13 @@ class RefFaceMassOperator(ElementwiseLinearOperator): volgrp.order, face_vertices) + if afgrp.dim == 0 and volgrp.dim == 1: + # NOTE: This complements a choice in `parametrization_derivative` + # where we don't really get a sign for the 0-dim case, so we add + # signs here in the hope that it'll only be used for fluxes where + # this makes sense + matrix[0, 0, 0] = -matrix[0, 0, 0] + # np.set_printoptions(linewidth=200, precision=3) # matrix[np.abs(matrix) < 1e-13] = 0 # print(matrix) @@ -602,8 +609,7 @@ def norm(p, arg, dd=None): if p == 2: norm_squared = sym.NodalSum(dd_in=dd)( - sym.FunctionSymbol("fabs")( - arg * sym.MassOperator()(arg))) + sym.fabs(arg * sym.MassOperator()(arg))) if isinstance(norm_squared, np.ndarray): norm_squared = norm_squared.sum() @@ -611,7 +617,7 @@ def norm(p, arg, dd=None): return sym.FunctionSymbol("sqrt")(norm_squared) elif p == np.Inf: - result = sym.NodalMax(dd_in=dd)(sym.FunctionSymbol("fabs")(arg)) + result = sym.NodalMax(dd_in=dd)(sym.fabs(arg)) from pymbolic.primitives import Max if isinstance(result, np.ndarray): diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 0b9fac20..80dd8413 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -364,6 +364,7 @@ class FunctionSymbol(ExpressionBase, pymbolic.primitives.Variable): mapper_method = "map_function_symbol" +fabs = FunctionSymbol("fabs") sqrt = FunctionSymbol("sqrt") exp = FunctionSymbol("exp") sin = FunctionSymbol("sin") @@ -506,7 +507,8 @@ def parametrization_derivative(ambient_dim, dim=None, dd=None): dim = ambient_dim if dim == 0: - return MultiVector(np.array([1.])) + inner_dd = (DD_VOLUME if dd is None else dd).with_qtag(QTAG_NONE) + return MultiVector(np.array([Ones(dd=inner_dd)])) from pytools import product return product( @@ -612,9 +614,10 @@ def mv_normal(dd, ambient_dim, dim=None): pseudoscalar(ambient_dim, dim, dd=dd) / # noqa: W504 area_element(ambient_dim, dim, dd=dd)) + return cse( # Dorst Section 3.7.2 - pder << pder.I.inv(), + pder if dim == 0 else pder << pder.I.inv(), "normal", cse_scope.DISCRETIZATION) -- GitLab From 2f9e7dbc178f1a8b56dfd1a0c3dc63f51403c0c8 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 22 Apr 2020 15:21:18 -0500 Subject: [PATCH 3/9] add 1d case to test_convergence_advec --- examples/advection/weak.py | 2 +- requirements.txt | 3 ++- test/test_grudge.py | 16 ++++++++++++---- 3 files changed, 15 insertions(+), 6 deletions(-) diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 68eeafe0..85545303 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -28,7 +28,7 @@ logger = logging.getLogger(__name__) logging.basicConfig(level=logging.INFO) -def main(ctx_factory, dim=1, order=4, visualize=True): +def main(ctx_factory, dim=2, order=4, visualize=True): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) diff --git a/requirements.txt b/requirements.txt index deb09394..4997171b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,4 +7,5 @@ git+https://gitlab.tiker.net/inducer/dagrt.git git+https://gitlab.tiker.net/inducer/leap.git git+https://github.com/inducer/meshpy.git git+https://github.com/inducer/modepy.git -git+https://github.com/inducer/meshmode.git +# git+https://github.com/inducer/meshmode.git +git+https://gitlab.tiker.net/fikl2/meshmode.git@face-connection-1d diff --git a/test/test_grudge.py b/test/test_grudge.py index e373e133..35f72a8d 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -204,18 +204,18 @@ def test_2d_gauss_theorem(ctx_factory): @pytest.mark.parametrize(("mesh_name", "mesh_pars"), [ + ("segment", [8, 16, 32]), ("disk", [0.1, 0.05]), ("rect2", [4, 8]), ("rect3", [4, 6]), ]) @pytest.mark.parametrize("op_type", ["strong", "weak"]) -@pytest.mark.parametrize("flux_type", ["upwind"]) +@pytest.mark.parametrize("flux_type", ["central"]) @pytest.mark.parametrize("order", [3, 4, 5]) # test: 'test_convergence_advec(cl._csc, "disk", [0.1, 0.05], "strong", "upwind", 3)' def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type, order, visualize=False): """Test whether 2D advection actually converges""" - cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) @@ -223,7 +223,16 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type eoc_rec = EOCRecorder() for mesh_par in mesh_pars: - if mesh_name == "disk": + if mesh_name == "segment": + from meshmode.mesh.generation import generate_box_mesh + mesh = generate_box_mesh( + [np.linspace(-1.0, 1.0, mesh_par)], + order=order) + + h = 2.0 / mesh_par + dim = 1 + dt_factor = 1.0 + elif mesh_name == "disk": pytest.importorskip("meshpy") from meshpy.geometry import make_circle, GeometryBuilder @@ -241,7 +250,6 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type h = np.sqrt(mesh_par) dim = 2 dt_factor = 4 - elif mesh_name.startswith("rect"): dim = int(mesh_name[4:]) from meshmode.mesh.generation import generate_regular_rect_mesh -- GitLab From 2284f61ae080cc0e09f176fcac8f6e721154e7fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Thu, 23 Apr 2020 00:02:23 +0200 Subject: [PATCH 4/9] Add CI runs based on Intel CL --- .gitlab-ci.yml | 52 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 57caec74..68a91853 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -32,6 +32,24 @@ Python 3 POCL: reports: junit: test/pytest.xml +Python 3 Intel: + script: + - export PY_EXE=python3 + - export EXTRA_INSTALL="pybind11 numpy mako" + - source /opt/enable-intel-cl.sh + - export PYOPENCL_TEST=intel(r):pu + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh + - ". ./build-and-test-py-project.sh" + tags: + - python3 + - pocl + - mpi + except: + - tags + artifacts: + reports: + junit: test/pytest.xml + Python 2.7 POCL MPI: script: - export PY_EXE=python2.7 @@ -68,6 +86,25 @@ Python 3 POCL MPI: reports: junit: test/pytest.xml +Python 3 Intel MPI: + script: + - export PY_EXE=python3 + - source /opt/enable-intel-cl.sh + - export PYOPENCL_TEST=intel(r):pu + - export EXTRA_INSTALL="pybind11 numpy mako mpi4py pymetis" + - export PYTEST_ADDOPTS="-k mpi" + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh + - ". ./build-and-test-py-project.sh" + tags: + - python3 + - pocl + - mpi + except: + - tags + artifacts: + reports: + junit: test/pytest.xml + Python 3 POCL Examples: script: - export PY_EXE=python3 @@ -82,6 +119,21 @@ Python 3 POCL Examples: except: - tags +Python 3 Intel Examples: + script: + - export PY_EXE=python3 + - source /opt/enable-intel-cl.sh + - export PYOPENCL_TEST=intel(r):pu + - export EXTRA_INSTALL="pybind11 numpy mako mpi4py pyvisfile pymetis" + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project-and-run-examples.sh + - ". ./build-py-project-and-run-examples.sh" + tags: + - python3 + - pocl + - large-node + except: + - tags + Documentation: script: - EXTRA_INSTALL="pybind11 numpy" -- GitLab From c6c954456f8703d03367c179ef83e2f8ca4993b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Thu, 23 Apr 2020 00:06:10 +0200 Subject: [PATCH 5/9] Add missing quoting for Intel CI jobs --- .gitlab-ci.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 68a91853..2de6d958 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -37,7 +37,7 @@ Python 3 Intel: - export PY_EXE=python3 - export EXTRA_INSTALL="pybind11 numpy mako" - source /opt/enable-intel-cl.sh - - export PYOPENCL_TEST=intel(r):pu + - export PYOPENCL_TEST="intel(r):pu" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" tags: @@ -90,7 +90,7 @@ Python 3 Intel MPI: script: - export PY_EXE=python3 - source /opt/enable-intel-cl.sh - - export PYOPENCL_TEST=intel(r):pu + - export PYOPENCL_TEST="intel(r):pu" - export EXTRA_INSTALL="pybind11 numpy mako mpi4py pymetis" - export PYTEST_ADDOPTS="-k mpi" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh @@ -123,7 +123,7 @@ Python 3 Intel Examples: script: - export PY_EXE=python3 - source /opt/enable-intel-cl.sh - - export PYOPENCL_TEST=intel(r):pu + - export PYOPENCL_TEST="intel(r):pu" - export EXTRA_INSTALL="pybind11 numpy mako mpi4py pyvisfile pymetis" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project-and-run-examples.sh - ". ./build-py-project-and-run-examples.sh" -- GitLab From 771c5e0c0c989d8c163f9b09edb04c08d4a3083c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 22 Apr 2020 17:14:04 -0500 Subject: [PATCH 6/9] update requirements.txt --- requirements.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index 4997171b..deb09394 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,5 +7,4 @@ git+https://gitlab.tiker.net/inducer/dagrt.git git+https://gitlab.tiker.net/inducer/leap.git git+https://github.com/inducer/meshpy.git git+https://github.com/inducer/modepy.git -# git+https://github.com/inducer/meshmode.git -git+https://gitlab.tiker.net/fikl2/meshmode.git@face-connection-1d +git+https://github.com/inducer/meshmode.git -- GitLab From 7ccdba2729038aeb59739dad1f99df9c56e5645a Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 23 Apr 2020 21:16:58 -0500 Subject: [PATCH 7/9] another attempt at making 1D work nicely The problem is that in 1D faces are 0D and so carry exactly no useful information about which side of the element they may be on. To hack around that, this commit introduces a hacky operator that hacks it. Less ambiguously, we just use the connection to figure out which face we're on on set the normal to `+1` or `-1` based on that --- examples/advection/weak.py | 2 -- examples/dagrt-fusion.py | 2 +- grudge/execution.py | 15 +++++++++ grudge/symbolic/dofdesc_inference.py | 1 + grudge/symbolic/mappers/__init__.py | 3 ++ grudge/symbolic/operators.py | 7 ----- grudge/symbolic/primitives.py | 46 +++++++++++++++++++++------- 7 files changed, 55 insertions(+), 21 deletions(-) diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 85545303..6e30094d 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -67,9 +67,7 @@ def main(ctx_factory, dim=2, order=4, visualize=True): from grudge import DGDiscretizationWithBoundaries discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) - volume_discr = discr.discr_from_dd(sym.DD_VOLUME) - faces_discr = discr.discr_from_dd(sym.FACE_RESTR_INTERIOR) # }}} diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py index 27e35abd..bfb6f761 100755 --- a/examples/dagrt-fusion.py +++ b/examples/dagrt-fusion.py @@ -159,7 +159,7 @@ def transcribe_phase(dag, field_var_name, field_components, phase_name, "
": sym.var("input_dt", sym.DD_SCALAR), f"{field_var_name}": sym.make_sym_array( f"input_{field_var_name}", field_components), - f"

residual": sym.make_sym_array( + "

residual": sym.make_sym_array( "input_residual", field_components), } diff --git a/grudge/execution.py b/grudge/execution.py index 948be131..fade7fca 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -283,6 +283,21 @@ class ExecutionMapper(mappers.Evaluator, # }}} + def map_signed_face_ones(self, expr): + face_discr = self.discrwb.discr_from_dd(expr.dd) + assert face_discr.dim == 0 + + all_faces_conn = self.discrwb.connection_from_dds("vol", expr.dd) + field = face_discr.empty(self.queue, allocator=self.bound_op.allocator) + field.fill(1) + + for grp in all_faces_conn.groups: + for batch in grp.batches: + i = batch.to_element_indices.with_queue(self.queue) + field[i] = (2.0 * (batch.to_element_face % 2) - 1.0) * field[i] + + return field + # }}} # {{{ instruction execution functions diff --git a/grudge/symbolic/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py index c44a7940..38c77512 100644 --- a/grudge/symbolic/dofdesc_inference.py +++ b/grudge/symbolic/dofdesc_inference.py @@ -183,6 +183,7 @@ class DOFDescInferenceMapper(RecursiveMapper, CSECachingMapperMixin): return expr.dd map_node_coordinate_component = map_ones + map_signed_face_ones = map_ones def map_call(self, expr): arg_dds = [ diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 0e9b8e02..b277c4c9 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -213,6 +213,7 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin): map_function_symbol = map_grudge_variable map_ones = map_grudge_variable + map_signed_face_ones = map_grudge_variable map_node_coordinate_component = map_grudge_variable # }}} @@ -266,6 +267,7 @@ class DependencyMapper( return set() map_ones = _map_leaf + map_signed_face_ones = _map_leaf map_node_coordinate_component = _map_leaf @@ -1235,6 +1237,7 @@ class CollectorMixin(OperatorReducerMixin, LocalOpReducerMixin, FluxOpReducerMix map_function_symbol = map_constant map_ones = map_grudge_variable + map_signed_face_ones = map_grudge_variable map_node_coordinate_component = map_grudge_variable map_operator = map_grudge_variable diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 0cbb6eda..f1997b1c 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -536,13 +536,6 @@ class RefFaceMassOperator(ElementwiseLinearOperator): volgrp.order, face_vertices) - if afgrp.dim == 0 and volgrp.dim == 1: - # NOTE: This complements a choice in `parametrization_derivative` - # where we don't really get a sign for the 0-dim case, so we add - # signs here in the hope that it'll only be used for fluxes where - # this makes sense - matrix[0, 0, 0] = -matrix[0, 0, 0] - # np.set_printoptions(linewidth=200, precision=3) # matrix[np.abs(matrix) < 1e-13] = 0 # print(matrix) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 80dd8413..875a7bff 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -426,6 +426,20 @@ class Ones(HasDOFDesc, ExpressionBase): mapper_method = intern("map_ones") +class _SignedFaceOnes(HasDOFDesc, ExpressionBase): + """Sets the face to :math:`-1` if its odd and :math:`+1` if its even. + This is based on the face order of + :meth:`meshmode.mesh.MeshElementGroup.face_vertex_indices`. + """ + + def __init__(self, dd): + dd = as_dofdesc(dd) + assert dd.is_trace() + super(_SignedFaceOnes, self).__init__(dd) + + mapper_method = intern("map_signed_face_ones") + + # {{{ geometry data class DiscretizationProperty(HasDOFDesc, ExpressionBase): @@ -507,8 +521,7 @@ def parametrization_derivative(ambient_dim, dim=None, dd=None): dim = ambient_dim if dim == 0: - inner_dd = (DD_VOLUME if dd is None else dd).with_qtag(QTAG_NONE) - return MultiVector(np.array([Ones(dd=inner_dd)])) + return MultiVector(np.array([_SignedFaceOnes(dd)])) from pytools import product return product( @@ -607,18 +620,29 @@ def mv_normal(dd, ambient_dim, dim=None): if dim is None: dim = ambient_dim - 1 - # Don't be tempted to add a sign here. As it is, it produces + # NOTE: Don't be tempted to add a sign here. As it is, it produces # exterior normals for positively oriented curves. - pder = ( - pseudoscalar(ambient_dim, dim, dd=dd) - / # noqa: W504 - area_element(ambient_dim, dim, dd=dd)) + pder = pseudoscalar(ambient_dim, dim, dd=dd) \ + / area_element(ambient_dim, dim, dd=dd) - return cse( - # Dorst Section 3.7.2 - pder if dim == 0 else pder << pder.I.inv(), - "normal", cse_scope.DISCRETIZATION) + # Dorst Section 3.7.2 + mv = pder << pder.I.inv() + + if dim == 0: + # NOTE: when the mesh is 0D, we do not have a clear notion of + # `exterior normal`, so we just take the tangent of the 1D element, + # interpolate it to the element faces and make it signed + tangent = parametrization_derivative( + ambient_dim, dim=1, dd=DD_VOLUME) + tangent = tangent / sqrt(tangent.norm_squared()) + + interp = _sym().interp(DD_VOLUME, dd) + mv = MultiVector(np.array([ + mv.as_scalar() * interp(t) for t in tangent.as_vector() + ])) + + return cse(mv, "normal", cse_scope.DISCRETIZATION) def normal(dd, ambient_dim, dim=None, quadrature_tag=None): -- GitLab From b471be590b4f31291704ed5a9bc0f063eb7ad4cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Sun, 26 Apr 2020 20:10:15 +0200 Subject: [PATCH 8/9] Apply suggestion to grudge/symbolic/primitives.py --- grudge/symbolic/primitives.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 875a7bff..c2397c47 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -427,9 +427,17 @@ class Ones(HasDOFDesc, ExpressionBase): class _SignedFaceOnes(HasDOFDesc, ExpressionBase): - """Sets the face to :math:`-1` if its odd and :math:`+1` if its even. + """Produces DoFs on a face that are :math:`-1` if their corresponding + face number is odd and :math:`+1` if it is even. + *dd* must refer to a 0D (point-shaped) trace domain. This is based on the face order of :meth:`meshmode.mesh.MeshElementGroup.face_vertex_indices`. + + .. note:: + + This is used as a hack to generate normals with the correct orientation + in 1D problems, and so far only intended for this particular use cases. + (If you can think of a better way, please speak up!) """ def __init__(self, dd): -- GitLab From 100210e67dd1428bb4393a8b992f57cf327b11ab Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 26 Apr 2020 14:51:40 -0500 Subject: [PATCH 9/9] fix flake8 --- grudge/symbolic/primitives.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 25350e44..f8bff1bf 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -641,7 +641,8 @@ def mv_normal(dd, ambient_dim, dim=None): ambient_dim, dim=1, dd=DD_VOLUME) tangent = tangent / sqrt(tangent.norm_squared()) - interp = _sym().interp(DD_VOLUME, dd) + from grudge.symbolic.operators import interp + interp = interp(DD_VOLUME, dd) mv = MultiVector(np.array([ mv.as_scalar() * interp(t) for t in tangent.as_vector() ])) -- GitLab