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 <http://www.gnu.org/licenses/>.
-
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,
"<dt>": sym.var("input_dt", sym.DD_SCALAR),
f"<state>{field_var_name}": sym.make_sym_array(
f"input_{field_var_name}", field_components),
- f"<p>residual": sym.make_sym_array(
+ "<p>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