diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0f0f630b023b09f998fb9fe2902939cb9a32ba9f..25e262d074551a616c5d65ad5874871afb936137 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -52,6 +52,21 @@ jobs:
build_py_project_in_conda_env
run_examples
+ docs:
+ name: Documentation
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v2
+ - name: "Main Script"
+ run: |
+ sudo apt-get update
+ sudo apt-get install openmpi-bin libopenmpi-dev
+ CONDA_ENVIRONMENT=.test-conda-env-py3.yml
+ curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/ci-support.sh
+ . ci-support.sh
+ build_py_project_in_conda_env
+ build_docs
+
# vim: sw=4
diff --git a/.gitignore b/.gitignore
index 94648fab1309424c88e2836ab4034eab05049871..14be95d194bd112b5596328cf5bafcfc092e0cce 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,3 +33,5 @@ run-debug-*
.cache
.pytest_cache
+
+*.json
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 66fc4e491c73a020bba5baee77f85048805696be..ebe4262476c11ae643380f4111a45f0d083f304c 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -128,8 +128,6 @@ Documentation:
- ". ./build-docs.sh"
tags:
- python3
- only:
- - master
Flake8:
script:
diff --git a/doc/conf.py b/doc/conf.py
index d75d2aa2ac850052f2c0ad196affaea5b0ec9e98..dc8f0e244b6a4212b817e4751f54e348bbe9661b 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -170,7 +170,7 @@ html_sidebars = {
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
+#html_static_path = ['_static']
# Add any extra paths that contain custom files (such as robots.txt or
# .htaccess) here, relative to this directory. These files are copied
@@ -327,4 +327,4 @@ intersphinx_mapping = {
'https://documen.tician.de/meshmode/': None,
'http://documen.tician.de/loopy/': None,
}
-autoclass_content = "both"
+autoclass_content = "class"
diff --git a/doc/discretization.rst b/doc/discretization.rst
new file mode 100644
index 0000000000000000000000000000000000000000..414111348b5de7a8261b088bdf5b5a2c28680a24
--- /dev/null
+++ b/doc/discretization.rst
@@ -0,0 +1,12 @@
+Discretization
+==============
+
+.. module:: grudge
+
+.. automodule:: grudge.discretization
+
+Discretization with Eager Evaluation
+====================================
+
+.. automodule:: grudge.eager
+
diff --git a/doc/eager.rst b/doc/eager.rst
deleted file mode 100644
index 5d2bdd5b071a15623134a1f549589ce6d0178f13..0000000000000000000000000000000000000000
--- a/doc/eager.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-Eagerly-evaluated Operation
-===========================
-
-.. automodule:: grudge.eager
-
diff --git a/doc/index.rst b/doc/index.rst
index 6d7782bc950b12431961f3bd0b8ee30deba4d050..7b6ce79d0a108cdcc8df3d05c00af34ee7213190 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -6,9 +6,9 @@ Contents:
.. toctree::
:maxdepth: 2
- misc
+ discretization
symbolic
- eager
+ misc
Indices and Tables
diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py
index 637da92e270d16ecb34881bab56521f054c1d430..bab82d545df9a7ba9bf4f6b41937dc10f45191df 100755
--- a/examples/dagrt-fusion.py
+++ b/examples/dagrt-fusion.py
@@ -29,6 +29,25 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+# FIXME:
+# Results before https://github.com/inducer/grudge/pull/15 were better:
+#
+# Operator | \parbox{1in}{\centering \% Memory Ops. Due to Scalar Assignments}
+# -------------+-------------------------------------------------------------------
+# 2D: Baseline | 51.1
+# 2D: Inlined | 48.9
+# 3D: Baseline | 50.1
+# 3D: Inlined | 48.6
+# INFO:__main__:Wrote '<stdout>'
+# ==== Scalar Assigment Inlining Impact ====
+# Operator | Bytes Read | Bytes Written | Total | \% of Baseline
+# -------------+------------+---------------+------------+----------------
+# 2D: Baseline | 9489600 | 3348000 | 12837600 | 100
+# 2D: Inlined | 8949600 | 2808000 | 11757600 | 91.6
+# 3D: Baseline | 1745280000 | 505440000 | 2250720000 | 100
+# 3D: Inlined | 1680480000 | 440640000 | 2121120000 | 94.2
+# INFO:__main__:Wrote '<stdout>'
+
import contextlib
import logging
@@ -413,66 +432,23 @@ class FusedRK4TimeStepper(RK4TimeStepperBase):
# {{{ problem setup code
-def get_strong_wave_op_with_discr(actx, dims=2, order=4):
- from meshmode.mesh.generation import generate_regular_rect_mesh
- mesh = generate_regular_rect_mesh(
- a=(-0.5,)*dims,
- b=(0.5,)*dims,
- n=(16,)*dims)
-
- logger.debug("%d elements", mesh.nelements)
-
- discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
-
+def _get_source_term(dims):
source_center = np.array([0.1, 0.22, 0.33])[:dims]
source_width = 0.05
source_omega = 3
- sym_x = sym.nodes(mesh.dim)
+ sym_x = sym.nodes(dims)
sym_source_center_dist = sym_x - source_center
sym_t = sym.ScalarVariable("t")
- from grudge.models.wave import StrongWaveOperator
- from meshmode.mesh import BTAG_ALL, BTAG_NONE
- op = StrongWaveOperator(-0.1, dims,
- source_f=(
+ return (
sym.sin(source_omega*sym_t)
* sym.exp(
-np.dot(sym_source_center_dist, sym_source_center_dist)
- / source_width**2)),
- dirichlet_tag=BTAG_NONE,
- neumann_tag=BTAG_NONE,
- radiation_tag=BTAG_ALL,
- flux_type="upwind")
-
- op.check_bc_coverage(mesh)
-
- return (op.sym_operator(), discr)
-
-
-def dg_flux(c, tpair):
- u = tpair[0]
- v = tpair[1:]
-
- dims = len(v)
-
- normal = sym.normal(tpair.dd, dims)
- flux_weak = flat_obj_array(
- np.dot(v.avg, normal),
- u.avg * normal)
-
- flux_weak -= (1 if c > 0 else -1)*flat_obj_array(
- 0.5*(u.int-u.ext),
- 0.5*(normal * np.dot(normal, v.int-v.ext)))
-
- flux_strong = flat_obj_array(
- np.dot(v.int, normal),
- u.int * normal) - flux_weak
-
- return sym.project(tpair.dd, "all_faces")(c*flux_strong)
+ / source_width**2))
-def get_strong_wave_op_with_discr_direct(actx, dims=2, order=4):
+def get_wave_op_with_discr(actx, dims=2, order=4):
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(
a=(-0.5,)*dims,
@@ -483,54 +459,21 @@ def get_strong_wave_op_with_discr_direct(actx, dims=2, order=4):
discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
- source_center = np.array([0.1, 0.22, 0.33])[:dims]
- source_width = 0.05
- source_omega = 3
-
- sym_x = sym.nodes(mesh.dim)
- sym_source_center_dist = sym_x - source_center
- sym_t = sym.ScalarVariable("t")
-
- from meshmode.mesh import BTAG_ALL
-
- c = -0.1
- sign = -1
-
- w = sym.make_sym_array("w", dims+1)
- u = w[0]
- v = w[1:]
-
- source_f = (
- sym.sin(source_omega*sym_t)
- * sym.exp(
- -np.dot(sym_source_center_dist, sym_source_center_dist)
- / source_width**2))
-
- rad_normal = sym.normal(BTAG_ALL, dims)
-
- rad_u = sym.cse(sym.project("vol", BTAG_ALL)(u))
- rad_v = sym.cse(sym.project("vol", BTAG_ALL)(v))
-
- rad_bc = sym.cse(flat_obj_array(
- 0.5*(rad_u - sign*np.dot(rad_normal, rad_v)),
- 0.5*rad_normal*(np.dot(rad_normal, rad_v) - sign*rad_u)
- ), "rad_bc")
+ from grudge.models.wave import WeakWaveOperator
+ from meshmode.mesh import BTAG_ALL, BTAG_NONE
+ op = WeakWaveOperator(0.1, dims,
+ source_f=_get_source_term(dims),
+ dirichlet_tag=BTAG_NONE,
+ neumann_tag=BTAG_NONE,
+ radiation_tag=BTAG_ALL,
+ flux_type="upwind")
- sym_operator = (
- - flat_obj_array(
- -c*np.dot(sym.nabla(dims), v) - source_f,
- -c*(sym.nabla(dims)*u)
- )
- + sym.InverseMassOperator()(
- sym.FaceMassOperator()(
- dg_flux(c, sym.int_tpair(w))
- + dg_flux(c, sym.bv_tpair(BTAG_ALL, w, rad_bc))
- )))
+ op.check_bc_coverage(mesh)
- return (sym_operator, discr)
+ return (op.sym_operator(), discr)
-def get_strong_wave_component(state_component):
+def get_wave_component(state_component):
return (state_component[0], state_component[1:])
# }}}
@@ -545,10 +488,10 @@ def test_stepper_equivalence(ctx_factory, order=4):
dims = 2
- sym_operator, _ = get_strong_wave_op_with_discr(
- actx, dims=dims, order=order)
- sym_operator_direct, discr = get_strong_wave_op_with_discr_direct(
+ sym_operator, discr = get_wave_op_with_discr(
actx, dims=dims, order=order)
+ #sym_operator_direct, discr = get_wave_op_with_discr_direct(
+ # actx, dims=dims, order=order)
if dims == 2:
dt = 0.04
@@ -561,11 +504,11 @@ def test_stepper_equivalence(ctx_factory, order=4):
bound_op = bind(discr, sym_operator)
stepper = RK4TimeStepper(
- discr, "w", bound_op, 1 + discr.dim, get_strong_wave_component)
+ discr, "w", bound_op, 1 + discr.dim, get_wave_component)
fused_stepper = FusedRK4TimeStepper(
- discr, "w", sym_operator_direct, 1 + discr.dim,
- get_strong_wave_component)
+ discr, "w", sym_operator, 1 + discr.dim,
+ get_wave_component)
t_start = 0
t_end = 0.5
@@ -808,7 +751,7 @@ def test_assignment_memory_model(ctx_factory):
queue = cl.CommandQueue(cl_ctx)
actx = PyOpenCLArrayContext(queue)
- _, discr = get_strong_wave_op_with_discr_direct(actx, dims=2, order=3)
+ _, discr = get_wave_op_with_discr(actx, dims=2, order=3)
# Assignment instruction
bound_op = bind(
@@ -838,7 +781,7 @@ def test_stepper_mem_ops(ctx_factory, use_fusion):
dims = 2
- sym_operator, discr = get_strong_wave_op_with_discr_direct(
+ sym_operator, discr = get_wave_op_with_discr(
actx, dims=dims, order=3)
t_start = 0
@@ -855,13 +798,13 @@ def test_stepper_mem_ops(ctx_factory, use_fusion):
stepper = RK4TimeStepper(
discr, "w", bound_op, 1 + discr.dim,
- get_strong_wave_component,
+ get_wave_component,
exec_mapper_factory=ExecutionMapperWithMemOpCounting)
else:
stepper = FusedRK4TimeStepper(
discr, "w", sym_operator, 1 + discr.dim,
- get_strong_wave_component,
+ get_wave_component,
exec_mapper_factory=ExecutionMapperWithMemOpCounting)
step = 0
@@ -1009,7 +952,7 @@ def test_stepper_timing(ctx_factory, use_fusion):
dims = 3
- sym_operator, discr = get_strong_wave_op_with_discr_direct(
+ sym_operator, discr = get_wave_op_with_discr(
actx, dims=dims, order=3)
t_start = 0
@@ -1026,13 +969,13 @@ def test_stepper_timing(ctx_factory, use_fusion):
stepper = RK4TimeStepper(
discr, "w", bound_op, 1 + discr.dim,
- get_strong_wave_component,
+ get_wave_component,
exec_mapper_factory=ExecutionMapperWithTiming)
else:
stepper = FusedRK4TimeStepper(
discr, "w", sym_operator, 1 + discr.dim,
- get_strong_wave_component,
+ get_wave_component,
exec_mapper_factory=ExecutionMapperWithTiming)
step = 0
@@ -1060,7 +1003,7 @@ def test_stepper_timing(ctx_factory, use_fusion):
def get_example_stepper(actx, dims=2, order=3, use_fusion=True,
exec_mapper_factory=ExecutionMapper,
return_ic=False):
- sym_operator, discr = get_strong_wave_op_with_discr_direct(
+ sym_operator, discr = get_wave_op_with_discr(
actx, dims=dims, order=3)
if not use_fusion:
@@ -1070,13 +1013,13 @@ def get_example_stepper(actx, dims=2, order=3, use_fusion=True,
stepper = RK4TimeStepper(
discr, "w", bound_op, 1 + discr.dim,
- get_strong_wave_component,
+ get_wave_component,
exec_mapper_factory=exec_mapper_factory)
else:
stepper = FusedRK4TimeStepper(
discr, "w", sym_operator, 1 + discr.dim,
- get_strong_wave_component,
+ get_wave_component,
exec_mapper_factory=exec_mapper_factory)
if return_ic:
@@ -1131,7 +1074,7 @@ def problem_stats(order=3):
actx = PyOpenCLArrayContext(queue)
with open_output_file("grudge-problem-stats.txt") as outf:
- _, dg_discr_2d = get_strong_wave_op_with_discr_direct(
+ _, dg_discr_2d = get_wave_op_with_discr(
actx, dims=2, order=order)
print("Number of 2D elements:", dg_discr_2d.mesh.nelements, file=outf)
vol_discr_2d = dg_discr_2d.discr_from_dd("vol")
@@ -1139,7 +1082,7 @@ def problem_stats(order=3):
from pytools import one
print("Number of DOFs per 2D element:", one(dofs_2d), file=outf)
- _, dg_discr_3d = get_strong_wave_op_with_discr_direct(
+ _, dg_discr_3d = get_wave_op_with_discr(
actx, dims=3, order=order)
print("Number of 3D elements:", dg_discr_3d.mesh.nelements, file=outf)
vol_discr_3d = dg_discr_3d.discr_from_dd("vol")
diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py
index f392b420d4e06c6fb547c6461c0290726388cf29..858e72549413b8120962a8052ebd7174622fff18 100644
--- a/examples/wave/var-propagation-speed.py
+++ b/examples/wave/var-propagation-speed.py
@@ -56,7 +56,7 @@ def main(write_output=True, order=4):
sym_t = sym.ScalarVariable("t")
c = sym.If(sym.Comparison(
np.dot(sym_x, sym_x), "<", 0.15),
- np.float32(-0.1), np.float32(-0.2))
+ np.float32(0.1), np.float32(0.2))
from grudge.models.wave import VariableCoefficientWeakWaveOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
diff --git a/examples/wave/wave-eager-mpi.py b/examples/wave/wave-eager-mpi.py
index 40ddbb7b1fb2049f54e173e0c5696a8a5a2d19e3..a6ca4381c9c84345d5d7eec64e470c43e05a20e1 100644
--- a/examples/wave/wave-eager-mpi.py
+++ b/examples/wave/wave-eager-mpi.py
@@ -59,9 +59,9 @@ def wave_flux(discr, c, w_tpair):
)
# upwind
- v_jump = np.dot(normal, v.int-v.ext)
+ v_jump = np.dot(normal, v.ext-v.int)
flux_weak += flat_obj_array(
- 0.5*(u.int-u.ext),
+ 0.5*(u.ext-u.int),
0.5*normal*scalar(v_jump),
)
@@ -87,7 +87,7 @@ def wave_operator(discr, c, w):
discr.face_mass(
wave_flux(discr, c=c, w_tpair=interior_trace_pair(discr, w))
+ wave_flux(discr, c=c, w_tpair=TracePair(
- BTAG_ALL, dir_bval, dir_bc))
+ BTAG_ALL, interior=dir_bval, exterior=dir_bc))
+ sum(
wave_flux(discr, c=c, w_tpair=tpair)
for tpair in cross_rank_trace_pairs(discr, w))
@@ -180,8 +180,6 @@ def main():
def rhs(t, w):
return wave_operator(discr, c=1, w=w)
- rank = comm.Get_rank()
-
t = 0
t_final = 3
istep = 0
@@ -190,7 +188,9 @@ def main():
if istep % 10 == 0:
print(istep, t, discr.norm(fields[0]))
- vis.write_vtk_file("fld-wave-eager-mpi-%03d-%04d.vtu" % (rank, istep),
+ vis.write_parallel_vtk_file(
+ comm,
+ f"fld-wave-eager-mpi-{{rank:03d}}-{istep:04d}.vtu",
[
("u", fields[0]),
("v", fields[1:]),
diff --git a/examples/wave/wave-eager-var-velocity.py b/examples/wave/wave-eager-var-velocity.py
index 9ac8a4c4f2a4a47f46d0a62f7c9271a7dda4f140..5164562cd03453f4a544fd2feb724bc44f8ce1d6 100644
--- a/examples/wave/wave-eager-var-velocity.py
+++ b/examples/wave/wave-eager-var-velocity.py
@@ -60,10 +60,9 @@ def wave_flux(discr, c, w_tpair):
)
# upwind
- v_jump = np.dot(normal, v.int-v.ext)
flux_weak += flat_obj_array(
- 0.5*(u.int-u.ext),
- 0.5*normal*scalar(v_jump),
+ 0.5*(u.ext-u.int),
+ 0.5*normal*scalar(np.dot(normal, v.ext-v.int)),
)
# FIXME this flux is only correct for continuous c
@@ -102,7 +101,7 @@ def wave_operator(discr, c, w):
dd_allfaces_quad,
wave_flux(discr, c=c, w_tpair=interior_trace_pair(discr, w))
+ wave_flux(discr, c=c, w_tpair=TracePair(
- BTAG_ALL, dir_bval, dir_bc))
+ BTAG_ALL, interior=dir_bval, exterior=dir_bc))
))
)
diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py
index 52ab08244d480ccfe6d58a4476e440fe4ad0e100..21c6d73a79d32adaae116e2183016ad88e0a82b5 100644
--- a/examples/wave/wave-eager.py
+++ b/examples/wave/wave-eager.py
@@ -57,10 +57,9 @@ def wave_flux(discr, c, w_tpair):
)
# upwind
- v_jump = np.dot(normal, v.int-v.ext)
flux_weak += flat_obj_array(
- 0.5*(u.int-u.ext),
- 0.5*normal*scalar(v_jump),
+ 0.5*(u.ext-u.int),
+ 0.5*normal*scalar(np.dot(normal, v.ext-v.int)),
)
return discr.project(w_tpair.dd, "all_faces", c*flux_weak)
@@ -85,7 +84,7 @@ def wave_operator(discr, c, w):
discr.face_mass(
wave_flux(discr, c=c, w_tpair=interior_trace_pair(discr, w))
+ wave_flux(discr, c=c, w_tpair=TracePair(
- BTAG_ALL, dir_bval, dir_bc))
+ BTAG_ALL, interior=dir_bval, exterior=dir_bc))
))
)
diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py
index bca2610199ed1cf50ae0c178d28caf1579730e3c..3e0063d866942bcc921f583c9464ee6d78f33bd7 100644
--- a/examples/wave/wave-min-mpi.py
+++ b/examples/wave/wave-min-mpi.py
@@ -79,9 +79,9 @@ def main(write_output=True, order=4):
sym_source_center_dist = sym_x - source_center
sym_t = sym.ScalarVariable("t")
- from grudge.models.wave import StrongWaveOperator
+ from grudge.models.wave import WeakWaveOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
- op = StrongWaveOperator(-0.1, discr.dim,
+ op = WeakWaveOperator(0.1, discr.dim,
source_f=(
sym.sin(source_omega*sym_t)
* sym.exp(
@@ -124,8 +124,6 @@ def main(write_output=True, order=4):
from time import time
t_last_step = time()
- rank = comm.Get_rank()
-
for event in dt_stepper.run(t_end=final_t):
if isinstance(event, dt_stepper.StateComputed):
assert event.component_id == "w"
@@ -135,11 +133,9 @@ def main(write_output=True, order=4):
print(step, event.t, norm(u=event.state_component[0]),
time()-t_last_step)
if step % 10 == 0:
- vis.write_vtk_file(
- "fld-wave-min-mpi-%03d-%04d.vtu" % (
- rank,
- step,
- ),
+ vis.write_parallel_vtk_file(
+ comm,
+ f"fld-wave-min-mpi-{{rank:03d}}-{step:04d}.vtu",
[
("u", event.state_component[0]),
("v", event.state_component[1:]),
diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py
index de6c9f92e75d449e6a42de80f7b6fa7270df8a15..e567251aa9631b1fe1da7e648208d91802c42ad7 100644
--- a/examples/wave/wave-min.py
+++ b/examples/wave/wave-min.py
@@ -61,9 +61,9 @@ def main(write_output=True, order=4):
sym_source_center_dist = sym_x - source_center
sym_t = sym.ScalarVariable("t")
- from grudge.models.wave import StrongWaveOperator
+ from grudge.models.wave import WeakWaveOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
- op = StrongWaveOperator(-0.1, discr.dim,
+ op = WeakWaveOperator(0.1, discr.dim,
source_f=(
sym.sin(source_omega*sym_t)
* sym.exp(
diff --git a/grudge/discretization.py b/grudge/discretization.py
index 7c0573c2e5ec3515843b6f8c377ad10301256aaa..a5452a24b9e690ec14cd0b5bf455f711f3b7f2e7 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -30,12 +30,12 @@ import numpy as np # noqa: F401
from meshmode.array_context import ArrayContext
-# FIXME Naming not ideal
-class DiscretizationBase(object):
- pass
+__doc__ = """
+.. autoclass:: DGDiscretizationWithBoundaries
+"""
-class DGDiscretizationWithBoundaries(DiscretizationBase):
+class DGDiscretizationWithBoundaries(object):
"""
.. automethod :: discr_from_dd
.. automethod :: connection_from_dds
@@ -53,10 +53,11 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
"""
:param quad_tag_to_group_factory: 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.
+ :class:`~meshmode.discretization.poly_element.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.
"""
self._setup_actx = array_context
@@ -390,7 +391,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
def order(self):
from warnings import warn
warn("DGDiscretizationWithBoundaries.order is deprecated, "
- "consider the orders of element groups instead. "
+ "consider using the orders of element groups instead. "
"'order' will go away in 2021.",
DeprecationWarning, stacklevel=2)
diff --git a/grudge/eager.py b/grudge/eager.py
index 87368161e5045b2a7c1b508b89cf60ebb336053a..38efbc02739560344b1c8084761dce5cd9bbc6f1 100644
--- a/grudge/eager.py
+++ b/grudge/eager.py
@@ -45,15 +45,28 @@ __doc__ = """
class EagerDGDiscretization(DGDiscretizationWithBoundaries):
"""
+ Inherits from :class:`~grudge.discretization.DGDiscretizationWithBoundaries`.
+
+ .. automethod:: __init__
.. automethod:: project
.. automethod:: nodes
+
.. automethod:: grad
+ .. automethod:: d_dx
.. automethod:: div
+
.. automethod:: weak_grad
+ .. automethod:: weak_d_dx
.. automethod:: weak_div
+
.. automethod:: normal
.. automethod:: inverse_mass
.. automethod:: face_mass
+
+ .. automethod:: norm
+ .. automethod:: nodal_sum
+ .. automethod:: nodal_min
+ .. automethod:: nodal_max
"""
def interp(self, src, tgt, vec):
@@ -64,6 +77,14 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
return self.project(src, tgt, vec)
def project(self, src, tgt, vec):
+ """Project from one discretization to another, e.g. from the
+ volume to the boundary, or from the base to the an overintegrated
+ quadrature discretization.
+
+ :arg src: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one
+ :arg tgt: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+ """
if (isinstance(vec, np.ndarray)
and vec.dtype.char == "O"
and not isinstance(vec, DOFArray)):
@@ -72,19 +93,80 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
return self.connection_from_dds(src, tgt)(vec)
- def nodes(self):
- return self._volume_discr.nodes()
+ def nodes(self, dd=None):
+ r"""Return the nodes of a discretization.
+
+ :arg dd: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one.
+ Defaults to the base volume discretization.
+ :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s
+ """
+ if dd is None:
+ return self._volume_discr.nodes()
+ else:
+ return self.discr_from_dd(dd).nodes()
+
+ # {{{ derivatives
@memoize_method
def _bound_grad(self):
return bind(self, sym.nabla(self.dim) * sym.Variable("u"), local_only=True)
def grad(self, vec):
+ r"""Return the gradient of the volume function represented by *vec*.
+
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+ :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s
+ """
return self._bound_grad()(u=vec)
+ @memoize_method
+ def _bound_d_dx(self, xyz_axis):
+ return bind(self, sym.nabla(self.dim)[xyz_axis] * sym.Variable("u"),
+ local_only=True)
+
+ def d_dx(self, xyz_axis, vec):
+ r"""Return the derivative along axis *xyz_axis* of the volume function
+ represented by *vec*.
+
+ :arg xyz_axis: an integer indicating the axis along which the derivative
+ is taken
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+ :returns: a :class:`~meshmode.dof_array.DOFArray`\ s
+ """
+ return self._bound_d_dx(xyz_axis)(u=vec)
+
+ def _div_helper(self, diff_func, vecs):
+ if not isinstance(vecs, np.ndarray):
+ raise TypeError("argument must be an object array")
+ assert vecs.dtype == np.object
+
+ if vecs.shape[-1] != self.ambient_dim:
+ raise ValueError("last dimension of *vecs* argument must match "
+ "ambient dimension")
+
+ if len(vecs.shape) == 1:
+ return sum(diff_func(i, vec_i) for i, vec_i in enumerate(vecs))
+ else:
+ result = np.zeros(vecs.shape[:-1], dtype=object)
+ for idx in np.ndindex(vecs.shape[:-1]):
+ result[idx] = sum(
+ diff_func(i, vec_i) for i, vec_i in enumerate(vecs[idx]))
+ return result
+
def div(self, vecs):
- return sum(
- self.grad(vec_i)[i] for i, vec_i in enumerate(vecs))
+ r"""Return the divergence of the vector volume function
+ represented by *vecs*.
+
+ :arg vec: an object array of
+ a :class:`~meshmode.dof_array.DOFArray`\ s,
+ where the last axis of the array must have length
+ matching the volume dimension.
+ :returns: a :class:`~meshmode.dof_array.DOFArray`
+ """
+
+ return self._div_helper(
+ lambda i, subvec: self.d_dx(i, subvec),
+ vecs)
@memoize_method
def _bound_weak_grad(self, dd):
@@ -93,6 +175,16 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
local_only=True)
def weak_grad(self, *args):
+ r"""Return the "weak gradient" of the volume function represented by
+ *vec*.
+
+ May be called with ``(vecs)`` or ``(dd, vecs)``.
+
+ :arg dd: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one.
+ Defaults to the base volume discretization if not provided.
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+ :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s
+ """
if len(args) == 1:
vec, = args
dd = sym.DOFDesc("vol", sym.QTAG_NONE)
@@ -103,7 +195,48 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
return self._bound_weak_grad(dd)(u=vec)
+ @memoize_method
+ def _bound_weak_d_dx(self, dd, xyz_axis):
+ return bind(self,
+ sym.stiffness_t(self.dim, dd_in=dd)[xyz_axis]
+ * sym.Variable("u", dd),
+ local_only=True)
+
+ def weak_d_dx(self, *args):
+ r"""Return the derivative along axis *xyz_axis* of the volume function
+ represented by *vec*.
+
+ May be called with ``(xyz_axis, vecs)`` or ``(dd, xyz_axis, vecs)``.
+
+ :arg xyz_axis: an integer indicating the axis along which the derivative
+ is taken
+ :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+ :returns: a :class:`~meshmode.dof_array.DOFArray`\ s
+ """
+ if len(args) == 2:
+ xyz_axis, vec = args
+ dd = sym.DOFDesc("vol", sym.QTAG_NONE)
+ elif len(args) == 3:
+ dd, xyz_axis, vec = args
+ else:
+ raise TypeError("invalid number of arguments")
+
+ return self._bound_weak_d_dx(dd, xyz_axis)(u=vec)
+
def weak_div(self, *args):
+ r"""Return the "weak divergence" of the vector volume function
+ represented by *vecs*.
+
+ May be called with ``(vecs)`` or ``(dd, vecs)``.
+
+ :arg dd: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one.
+ Defaults to the base volume discretization if not provided.
+ :arg vec: a object array of
+ a :class:`~meshmode.dof_array.DOFArray`\ s,
+ where the last axis of the array must have length
+ matching the volume dimension.
+ :returns: a :class:`~meshmode.dof_array.DOFArray`
+ """
if len(args) == 1:
vecs, = args
dd = sym.DOFDesc("vol", sym.QTAG_NONE)
@@ -112,8 +245,11 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
else:
raise TypeError("invalid number of arguments")
- return sum(
- self.weak_grad(dd, vec_i)[i] for i, vec_i in enumerate(vecs))
+ return self._div_helper(
+ lambda i, subvec: self.weak_d_dx(dd, i, subvec),
+ vecs)
+
+ # }}}
@memoize_method
def normal(self, dd):
@@ -218,9 +354,11 @@ def interior_trace_pair(discrwb, vec):
lambda el: discrwb.opposite_face_connection()(el),
i)
- return TracePair("int_faces", i, e)
+ return TracePair("int_faces", interior=i, exterior=e)
+# {{{ distributed-memory functionality
+
class _RankBoundaryCommunication:
base_tag = 1273
@@ -259,8 +397,9 @@ class _RankBoundaryCommunication:
self.send_req.Wait()
- return TracePair(self.remote_btag, self.local_dof_array,
- swapped_remote_dof_array)
+ return TracePair(self.remote_btag,
+ interior=self.local_dof_array,
+ exterior=swapped_remote_dof_array)
def _cross_rank_trace_pairs_scalar_field(discrwb, vec, tag=None):
@@ -295,5 +434,7 @@ def cross_rank_trace_pairs(discrwb, vec, tag=None):
else:
return _cross_rank_trace_pairs_scalar_field(discrwb, vec, tag=tag)
+# }}}
+
# vim: foldmethod=marker
diff --git a/grudge/models/em.py b/grudge/models/em.py
index 4810b4d0ebd8f023adc2eebf01dd62bc67957d6d..e21ad4994544f93a6cd135988f963ddec17c1d47 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -36,11 +36,9 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE
from grudge import sym
from pytools.obj_array import flat_obj_array, make_obj_array
-# TODO: Check PML
-
class MaxwellOperator(HyperbolicOperator):
- """A 3D Maxwell operator which supports fixed or variable
+ """A strong-form 3D Maxwell operator which supports fixed or variable
isotropic, non-dispersive, positive epsilon and mu.
Field order is [Ex Ey Ez Hx Hy Hz].
@@ -109,7 +107,7 @@ class MaxwellOperator(HyperbolicOperator):
self.incident_bc_data = incident_bc
def flux(self, w):
- """The template for the numerical flux for variable coefficients.
+ """The numerical flux for variable coefficients.
:param flux_type: can be in [0,1] for anything between central and upwind,
or "lf" for Lax-Friedrichs.
@@ -136,13 +134,13 @@ class MaxwellOperator(HyperbolicOperator):
return flat_obj_array(
# flux e,
1/2*(
- -self.space_cross_h(normal, h.int-h.ext)
+ -self.space_cross_h(normal, h.ext-h.int)
# multiplication by epsilon undoes material divisor below
#-max_c*(epsilon*e.int - epsilon*e.ext)
),
# flux h
1/2*(
- self.space_cross_e(normal, e.int-e.ext)
+ self.space_cross_e(normal, e.ext-e.int)
# multiplication by mu undoes material divisor below
#-max_c*(mu*h.int - mu*h.ext)
))
@@ -152,14 +150,14 @@ class MaxwellOperator(HyperbolicOperator):
# flux e,
(
-1/(Z_int+Z_ext)*self.space_cross_h(normal,
- Z_ext*(h.int-h.ext)
- - self.flux_type*self.space_cross_e(normal, e.int-e.ext))
+ Z_ext*(h.ext-h.int)
+ - self.flux_type*self.space_cross_e(normal, e.ext-e.int))
),
# flux h
(
1/(Y_int + Y_ext)*self.space_cross_e(normal,
- Y_ext*(e.int-e.ext)
- + self.flux_type*self.space_cross_h(normal, h.int-h.ext))
+ Y_ext*(e.ext-e.int)
+ + self.flux_type*self.space_cross_h(normal, h.ext-h.int))
),
)
else:
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 09a954e873573df374402c28762668573ae71882..863b0112052f71853fd746a80e29385b5ebf187f 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -35,153 +35,6 @@ from pytools.obj_array import flat_obj_array
# {{{ constant-velocity
-class StrongWaveOperator(HyperbolicOperator):
- r"""This operator discretizes the wave equation
- :math:`\partial_t^2 u = c^2 \Delta u`.
-
- To be precise, we discretize the hyperbolic system
-
- .. math::
-
- \partial_t u - c \\nabla \\cdot v = 0
-
- \partial_t v - c \\nabla u = 0
-
- The sign of :math:`v` determines whether we discretize the forward or the
- backward wave equation.
-
- :math:`c` is assumed to be constant across all space.
- """
-
- def __init__(self, c, ambient_dim, source_f=0,
- flux_type="upwind",
- dirichlet_tag=BTAG_ALL,
- dirichlet_bc_f=0,
- neumann_tag=BTAG_NONE,
- radiation_tag=BTAG_NONE):
- assert isinstance(ambient_dim, int)
-
- self.c = c
- self.ambient_dim = ambient_dim
- self.source_f = source_f
-
- if self.c > 0:
- self.sign = 1
- else:
- self.sign = -1
-
- self.dirichlet_tag = dirichlet_tag
- self.neumann_tag = neumann_tag
- self.radiation_tag = radiation_tag
-
- self.dirichlet_bc_f = dirichlet_bc_f
-
- self.flux_type = flux_type
-
- def flux(self, w):
- u = w[0]
- v = w[1:]
- normal = sym.normal(w.dd, self.ambient_dim)
-
- flux_weak = flat_obj_array(
- np.dot(v.avg, normal),
- u.avg * normal)
-
- if self.flux_type == "central":
- pass
- elif self.flux_type == "upwind":
- # see doc/notes/grudge-notes.tm
- flux_weak -= self.sign*flat_obj_array(
- 0.5*(u.int-u.ext),
- 0.5*(normal * np.dot(normal, v.int-v.ext)))
- else:
- raise ValueError("invalid flux type '%s'" % self.flux_type)
-
- flux_strong = flat_obj_array(
- np.dot(v.int, normal),
- u.int * normal) - flux_weak
-
- return self.c*flux_strong
-
- def sym_operator(self):
- d = self.ambient_dim
-
- w = sym.make_sym_array("w", d+1)
- u = w[0]
- v = w[1:]
-
- # boundary conditions -------------------------------------------------
-
- # dirichlet BCs -------------------------------------------------------
- dir_u = sym.cse(sym.project("vol", self.dirichlet_tag)(u))
- dir_v = sym.cse(sym.project("vol", self.dirichlet_tag)(v))
- if self.dirichlet_bc_f:
- # FIXME
- from warnings import warn
- warn("Inhomogeneous Dirichlet conditions on the wave equation "
- "are still having issues.")
-
- dir_g = sym.Field("dir_bc_u")
- dir_bc = flat_obj_array(2*dir_g - dir_u, dir_v)
- else:
- dir_bc = flat_obj_array(-dir_u, dir_v)
-
- dir_bc = sym.cse(dir_bc, "dir_bc")
-
- # neumann BCs ---------------------------------------------------------
- neu_u = sym.cse(sym.project("vol", self.neumann_tag)(u))
- neu_v = sym.cse(sym.project("vol", self.neumann_tag)(v))
- neu_bc = sym.cse(flat_obj_array(neu_u, -neu_v), "neu_bc")
-
- # radiation BCs -------------------------------------------------------
- rad_normal = sym.normal(self.radiation_tag, d)
-
- rad_u = sym.cse(sym.project("vol", self.radiation_tag)(u))
- rad_v = sym.cse(sym.project("vol", self.radiation_tag)(v))
-
- rad_bc = sym.cse(flat_obj_array(
- 0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
- 0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.sign*rad_u)
- ), "rad_bc")
-
- # entire operator -----------------------------------------------------
- nabla = sym.nabla(d)
-
- def flux(pair):
- return sym.project(pair.dd, "all_faces")(
- self.flux(pair))
-
- result = (
- - flat_obj_array(
- -self.c*np.dot(nabla, v),
- -self.c*(nabla*u)
- )
- + # noqa: W504
- sym.InverseMassOperator()(
- sym.FaceMassOperator()(
- flux(sym.int_tpair(w))
- + flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
- + flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
- + flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
- )
- )
- )
-
- result[0] += self.source_f
-
- return result
-
- def check_bc_coverage(self, mesh):
- from meshmode.mesh import check_bc_coverage
- check_bc_coverage(mesh, [
- self.dirichlet_tag,
- self.neumann_tag,
- self.radiation_tag])
-
- def max_eigenvalue(self, t, fields=None, discr=None):
- return abs(self.c)
-
-
class WeakWaveOperator(HyperbolicOperator):
r"""This operator discretizes the wave equation
:math:`\partial_t^2 u = c^2 \Delta u`.
@@ -230,16 +83,16 @@ class WeakWaveOperator(HyperbolicOperator):
v = w[1:]
normal = sym.normal(w.dd, self.ambient_dim)
- flux_weak = flat_obj_array(
+ central_flux_weak = -self.c*flat_obj_array(
np.dot(v.avg, normal),
u.avg * normal)
if self.flux_type == "central":
- return -self.c*flux_weak
+ return central_flux_weak
elif self.flux_type == "upwind":
- return -self.c*(flux_weak + self.sign*flat_obj_array(
- 0.5*(u.int-u.ext),
- 0.5*(normal * np.dot(normal, v.int-v.ext))))
+ return central_flux_weak - self.c*self.sign*flat_obj_array(
+ 0.5*(u.ext-u.int),
+ 0.5*(normal * np.dot(normal, v.ext-v.int)))
else:
raise ValueError("invalid flux type '%s'" % self.flux_type)
@@ -294,7 +147,6 @@ class WeakWaveOperator(HyperbolicOperator):
-self.c*(sym.stiffness_t(self.ambient_dim)*u)
)
-
- sym.FaceMassOperator()(flux(sym.int_tpair(w))
+ flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
+ flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
@@ -370,18 +222,18 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
v = w[2:]
normal = sym.normal(w.dd, self.ambient_dim)
+ flux_central_weak = -0.5 * flat_obj_array(
+ np.dot(v.int*c.int + v.ext*c.ext, normal),
+ (u.int * c.int + u.ext*c.ext) * normal)
+
if self.flux_type == "central":
- return -0.5 * flat_obj_array(
- np.dot(v.int*c.int + v.ext*c.ext, normal),
- (u.int * c.int + u.ext*c.ext) * normal)
+ return flux_central_weak
elif self.flux_type == "upwind":
- return -0.5 * flat_obj_array(
- np.dot(normal, c.ext * v.ext + c.int * v.int)
- + c.ext*u.ext - c.int * u.int,
+ return flux_central_weak - 0.5 * flat_obj_array(
+ c.ext*u.ext - c.int * u.int,
- normal * (np.dot(normal, c.ext * v.ext - c.int * v.int)
- + c.ext*u.ext + c.int * u.int))
+ normal * (np.dot(normal, c.ext * v.ext - c.int * v.int)))
else:
raise ValueError("invalid flux type '%s'" % self.flux_type)
@@ -443,7 +295,6 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
-self.c*(sym.stiffness_t(self.ambient_dim)*u)
)
-
- sym.FaceMassOperator()(flux(sym.int_tpair(flux_w))
+ flux(sym.bv_tpair(self.dirichlet_tag, flux_w, dir_bc))
+ flux(sym.bv_tpair(self.neumann_tag, flux_w, neu_bc))
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index a8ab366b52cca01bb6bedd060fe05adf7c1295df..6a5a982ce29cf3a1168b1c3851f1c9a7195f59ec 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -94,12 +94,12 @@ class Operator(pymbolic.primitives.Expression):
"""
.. attribute:: dd_in
- an instance of :class:`~grudge.symbolic.primitives.DOFDesc` describing the
+ an instance of :class:`~grudge.sym.DOFDesc` describing the
input discretization.
.. attribute:: dd_out
- an instance of :class:`~grudge.symbolic.primitives.DOFDesc` describing the
+ an instance of :class:`~grudge.sym.DOFDesc` describing the
output discretization.
"""
@@ -477,7 +477,10 @@ class VandermondeOperator(ElementwiseLinearOperator):
class MassOperatorBase(Operator):
"""
- :attr:`dd_in` and :attr:`dd_out` may be surface or volume discretizations.
+ Inherits from :class:`Operator`.
+
+ :attr:`~Operator.dd_in` and :attr:`~Operator.dd_out` may be surface or volume
+ discretizations.
"""
def __init__(self, dd_in=None, dd_out=None):
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index d706314cdd03b23fa46ff43a7df16234dce6cd70..8d7d81328db86b4f107e3efa1b915150811796c0 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -55,6 +55,8 @@ class ExpressionBase(prim.Expression):
__doc__ = """
+.. currentmodule:: grudge.sym
+
DOF description
^^^^^^^^^^^^^^^
@@ -863,7 +865,7 @@ class TracePair:
an expression representing the exterior value to
be used for the flux.
"""
- def __init__(self, dd, interior, exterior):
+ def __init__(self, dd, *, interior, exterior):
"""
"""
@@ -874,8 +876,8 @@ class TracePair:
def __getitem__(self, index):
return TracePair(
self.dd,
- self.exterior[index],
- self.interior[index])
+ interior=self.interior[index],
+ exterior=self.exterior[index])
def __len__(self):
assert len(self.exterior) == len(self.interior)
@@ -912,7 +914,7 @@ def int_tpair(expression, qtag=None, from_dd=None):
i = cse(project(trace_dd.with_qtag(None), trace_dd)(i))
e = cse(project(trace_dd.with_qtag(None), trace_dd)(e))
- return TracePair(trace_dd, i, e)
+ return TracePair(trace_dd, interior=i, exterior=e)
def bdry_tpair(dd, interior, exterior):
@@ -924,7 +926,7 @@ def bdry_tpair(dd, interior, exterior):
representing the exterior value to be used
for the flux.
"""
- return TracePair(dd, interior, exterior)
+ return TracePair(dd, interior=interior, exterior=exterior)
def bv_tpair(dd, interior, exterior):
@@ -938,7 +940,7 @@ def bv_tpair(dd, interior, exterior):
"""
from grudge.symbolic.operators import project
interior = cse(project("vol", dd)(interior))
- return TracePair(dd, interior, exterior)
+ return TracePair(dd, interior=interior, exterior=exterior)
# }}}
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 25707b217fed39acb7956e8a9ab6aefaccf9d358..faf4bc05c4ffb28a4cd98c0a518fca9be58c5c0b 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -34,8 +34,8 @@ from pytools.obj_array import flat_obj_array, make_obj_array
from grudge import sym, bind, DGDiscretizationWithBoundaries
import pytest
-from pyopencl.tools import ( # noqa
- pytest_generate_tests_for_pyopencl
+from meshmode.array_context import ( # noqa
+ pytest_generate_tests_for_pyopencl_array_context
as pytest_generate_tests)
import logging
@@ -46,10 +46,8 @@ logger = logging.getLogger(__name__)
# {{{ inverse metric
@pytest.mark.parametrize("dim", [2, 3])
-def test_inverse_metric(ctx_factory, dim):
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+def test_inverse_metric(actx_factory, dim):
+ actx = actx_factory()
from meshmode.mesh.generation import generate_regular_rect_mesh
mesh = generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim,
@@ -97,13 +95,11 @@ def test_inverse_metric(ctx_factory, dim):
@pytest.mark.parametrize("ambient_dim", [1, 2, 3])
@pytest.mark.parametrize("quad_tag", [sym.QTAG_NONE, "OVSMP"])
-def test_mass_mat_trig(ctx_factory, ambient_dim, quad_tag):
+def test_mass_mat_trig(actx_factory, ambient_dim, quad_tag):
"""Check the integral of some trig functions on an interval using the mass
matrix.
"""
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+ actx = actx_factory()
nelements = 17
order = 4
@@ -196,10 +192,8 @@ def _spheroid_surface_area(radius, aspect_ratio):
@pytest.mark.parametrize("name", [
"2-1-ellipse", "spheroid", "box2d", "box3d"
])
-def test_mass_surface_area(ctx_factory, name):
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+def test_mass_surface_area(actx_factory, name):
+ actx = actx_factory()
# {{{ cases
@@ -266,7 +260,7 @@ def test_mass_surface_area(ctx_factory, name):
# {{{ surface mass inverse
@pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"])
-def test_surface_mass_operator_inverse(ctx_factory, name):
+def test_surface_mass_operator_inverse(actx_factory, name):
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
actx = PyOpenCLArrayContext(queue)
@@ -339,12 +333,9 @@ def _avg_face_normal(x, side=-1):
@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"])
-def test_face_normal_surface(ctx_factory, mesh_name):
+def test_face_normal_surface(actx_factory, mesh_name):
"""Check that face normals are orthogonal to the surface normal"""
-
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+ actx = actx_factory()
# {{{ geometry
@@ -453,15 +444,13 @@ def test_face_normal_surface(ctx_factory, mesh_name):
# {{{ diff operator
@pytest.mark.parametrize("dim", [1, 2, 3])
-def test_tri_diff_mat(ctx_factory, dim, order=4):
+def test_tri_diff_mat(actx_factory, dim, order=4):
"""Check differentiation matrix along the coordinate axes on a disk
Uses sines as the function to differentiate.
"""
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+ actx = actx_factory()
from meshmode.mesh.generation import generate_regular_rect_mesh
@@ -497,7 +486,7 @@ def test_tri_diff_mat(ctx_factory, dim, order=4):
# {{{ divergence theorem
-def test_2d_gauss_theorem(ctx_factory):
+def test_2d_gauss_theorem(actx_factory):
"""Verify Gauss's theorem explicitly on a mesh"""
pytest.importorskip("meshpy")
@@ -515,9 +504,7 @@ def test_2d_gauss_theorem(ctx_factory):
from meshmode.mesh.io import from_meshpy
mesh = from_meshpy(mesh_info, order=1)
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+ actx = actx_factory()
discr = DGDiscretizationWithBoundaries(actx, mesh, order=2)
@@ -541,7 +528,7 @@ def test_2d_gauss_theorem(ctx_factory):
@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"])
-def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False):
+def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
r"""Check the surface divergence theorem.
.. math::
@@ -555,10 +542,7 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False):
face normal (which should be orthogonal to both the surface normal
and the face tangent).
"""
-
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+ actx = actx_factory()
# {{{ cases
@@ -709,13 +693,11 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False):
@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,
+def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_type,
order, visualize=False):
"""Test whether 2D advection actually converges"""
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+ actx = actx_factory()
from pytools.convergence import EOCRecorder
eoc_rec = EOCRecorder()
@@ -855,12 +837,10 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
# {{{ models: maxwell
@pytest.mark.parametrize("order", [3, 4, 5])
-def test_convergence_maxwell(ctx_factory, order):
+def test_convergence_maxwell(actx_factory, order):
"""Test whether 3D Maxwell's actually converges"""
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+ actx = actx_factory()
from pytools.convergence import EOCRecorder
eoc_rec = EOCRecorder()
@@ -930,16 +910,14 @@ def test_convergence_maxwell(ctx_factory, order):
# {{{ models: variable coefficient advection oversampling
@pytest.mark.parametrize("order", [2, 3, 4])
-def test_improvement_quadrature(ctx_factory, order):
+def test_improvement_quadrature(actx_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 meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+ actx = actx_factory()
dims = 2
sym_nds = sym.nodes(dims)
@@ -1032,10 +1010,8 @@ def test_op_collector_order_determinism():
# {{{ bessel
-def test_bessel(ctx_factory):
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+def test_bessel(actx_factory):
+ actx = actx_factory()
dims = 2
@@ -1066,10 +1042,8 @@ def test_bessel(ctx_factory):
# {{{ function symbol
-def test_external_call(ctx_factory):
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(queue)
+def test_external_call(actx_factory):
+ actx = actx_factory()
def double(queue, x):
return 2 * x
diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py
index a5109737106279f0a76d147a9d77a4e35d2d6d84..91656bda4d1e123cc6a69eafe0c2f1a9b4eb1fa4 100644
--- a/test/test_mpi_communication.py
+++ b/test/test_mpi_communication.py
@@ -140,9 +140,9 @@ def mpi_communication_entrypoint():
sym_source_center_dist = sym_x - source_center
sym_t = sym.ScalarVariable("t")
- from grudge.models.wave import StrongWaveOperator
+ from grudge.models.wave import WeakWaveOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
- op = StrongWaveOperator(-0.1, vol_discr.dim,
+ op = WeakWaveOperator(0.1, vol_discr.dim,
source_f=(
sym.sin(source_omega*sym_t)
* sym.exp(