diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py
index a493059e642e24867c2b6658722cad837189d3ae..fd37f27f867bd9db0bd195371cd768a80a243f25 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
@@ -807,7 +750,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(
@@ -837,7 +780,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
@@ -854,13 +797,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
@@ -1008,7 +951,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
@@ -1025,13 +968,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
@@ -1059,7 +1002,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:
@@ -1069,13 +1012,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:
@@ -1130,7 +1073,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")
@@ -1138,7 +1081,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/eager.py b/grudge/eager.py
index 4cff62c5efe9dda454cff20b06789f0db0412f44..38efbc02739560344b1c8084761dce5cd9bbc6f1 100644
--- a/grudge/eager.py
+++ b/grudge/eager.py
@@ -354,7 +354,7 @@ 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
@@ -397,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):
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/primitives.py b/grudge/symbolic/primitives.py
index 0a49034568f6ec45c25b45018e9451fb866fe361..5376865883c3e397a152f432037d9ac963adbdff 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -690,7 +690,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):
"""
"""
@@ -701,8 +701,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)
@@ -734,7 +734,7 @@ def int_tpair(expression, qtag=None):
else:
q_dd = "int_faces"
- return TracePair(q_dd, i, e)
+ return TracePair(q_dd, interior=i, exterior=e)
def bdry_tpair(dd, interior, exterior):
@@ -746,7 +746,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):
@@ -760,7 +760,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 03e09fb658108d1125f2bd5ce4297930f87f3e37..25fdbc409227ba01135cb97efb05f4c80945ee42 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -34,8 +34,10 @@ import pyopencl as cl
from meshmode.array_context import PyOpenCLArrayContext
from meshmode.dof_array import unflatten, flatten, flat_norm
-from pyopencl.tools import ( # noqa
- pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+from meshmode.array_context import ( # noqa
+ pytest_generate_tests_for_pyopencl_array_context
+ as pytest_generate_tests)
+
import logging
@@ -47,10 +49,8 @@ from grudge import sym, bind, DGDiscretizationWithBoundaries
@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,
@@ -94,13 +94,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
@@ -160,15 +158,13 @@ def test_mass_mat_trig(ctx_factory, ambient_dim, quad_tag):
@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
@@ -200,7 +196,7 @@ def test_tri_diff_mat(ctx_factory, dim, order=4):
assert eoc_rec.order_estimate() >= order
-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")
@@ -218,9 +214,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)
@@ -253,13 +247,11 @@ def test_2d_gauss_theorem(ctx_factory):
@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()
@@ -381,12 +373,10 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
@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()
@@ -452,16 +442,14 @@ def test_convergence_maxwell(ctx_factory, order):
@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)
@@ -546,10 +534,8 @@ def test_op_collector_order_determinism():
assert list(TestBoundOperatorCollector(TestOperator)(ob0 + ob1)) == [ob0, ob1]
-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
@@ -576,10 +562,8 @@ def test_bessel(ctx_factory):
assert z < 1e-15
-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(