diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 57caec74f47a8cc99b8c94ba1c6dfd460097715f..2de6d958a74cc452ee8e130ec7fe534d99e5f8a2 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" diff --git a/examples/advection/weak.py b/examples/advection/weak.py index bd9d61ca3e2bbfd416a271ffbbabca22f366e3fd..6e30094d080033f225a9dbc45227e0d0c7665ad0 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -13,106 +13,135 @@ # 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__) - -from grudge import sym, bind, DGDiscretizationWithBoundaries - -import numpy.linalg as la - +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO) -def main(write_output=True, order=4): - cl_ctx = cl.create_some_context() +def main(ctx_factory, dim=2, order=4, visualize=True): + cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) - dim = 2 + # {{{ 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 = 2.0 + # final time + final_time = 1.0 + + # velocity field + c = np.array([0.5] * dim) + norm_c = la.norm(c) + # flux + flux_type = "central" + # 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) + # }}} - dt_factor = 4 - h = 1/20 + # {{{ discretization - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + 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) - c = np.array([0.1,0.1]) - norm_c = la.norm(c) + from grudge import DGDiscretizationWithBoundaries + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + # }}} - 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) + 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 - 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) 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) + + volume_x = volume_discr.nodes().get(queue) + else: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + + def plot_solution(evt): + if not visualize: + return + + if dim == 1: + u = event.state_component.get(queue) + u_ = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=evt.t).get(queue) + + 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(evt.t)) + pt.savefig("fld-weak-%04d.png" % step) + pt.clf() + else: + vis.write_vtk_file("fld-weak-%04d.vtu" % step, [ + ("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 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) + plot_solution(event) 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/examples/dagrt-fusion.py b/examples/dagrt-fusion.py index b9ef892e675ec70b9df309541013ebc89dd36878..d4b977b129de7ef7019d4b237023c73095d91406 100755 --- a/examples/dagrt-fusion.py +++ b/examples/dagrt-fusion.py @@ -160,7 +160,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 948be131205034a2d0ddb597a11d2c8ad8f42327..fade7fca37f92d728c2573c82f7af9a8c8cbcbdf 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/models/advection.py b/grudge/models/advection.py index 091b9870c4b82a3fde2fc21b84bf92e5ec62ab17..c3b142e411cd5fd7b34b17bce69e529ff565d316 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/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py index c44a7940b9ec142171a0b0a67885c041ec9d43cb..38c775127488b79624264439778c334fac29fcd3 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 0e9b8e02d67b6f2b86d217d8c230263281d49a19..b277c4c99e97280e4f82be63cac8bbcbb00bee8e 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/primitives.py b/grudge/symbolic/primitives.py index ea9e89f6bb6a607ec96944b45c3f7f8c5f8a3193..f8bff1bf2af1403a652f37a32995fd3a70b88f6a 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -420,6 +420,28 @@ class Ones(HasDOFDesc, ExpressionBase): mapper_method = intern("map_ones") +class _SignedFaceOnes(HasDOFDesc, ExpressionBase): + """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): + 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): @@ -503,7 +525,7 @@ def parametrization_derivative(ambient_dim, dim=None, dd=None): dim = ambient_dim if dim == 0: - return MultiVector(np.array([1.])) + return MultiVector(np.array([_SignedFaceOnes(dd)])) from pytools import product return product( @@ -602,17 +624,30 @@ 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)) - return cse( - # Dorst Section 3.7.2 - pder << pder.I.inv(), - "normal", cse_scope.DISCRETIZATION) + pder = pseudoscalar(ambient_dim, dim, dd=dd) \ + / area_element(ambient_dim, dim, dd=dd) + + # 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()) + + 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() + ])) + + return cse(mv, "normal", cse_scope.DISCRETIZATION) def normal(dd, ambient_dim, dim=None, quadrature_tag=None): diff --git a/test/test_grudge.py b/test/test_grudge.py index 1d2732a9acfb9b51c0aadf4534ddb85d49d46995..65cd9d8975d49ebdd46c9029fd803aafb8089766 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -207,18 +207,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) @@ -226,7 +226,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 @@ -244,7 +253,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