diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 2ab36d516416796c6760023000638ead0004fbda..3c9304ae68ca7fa2ff47aa79f0c8c8c364c87f94 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -148,8 +148,8 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): quad_tag_to_group_factory[product_tag] = \ QuadratureSimplexGroupFactory(order=4*order) - from grudge import DGDiscretizationWithBoundaries - discr = DGDiscretizationWithBoundaries(actx, mesh, + from grudge import DiscretizationCollection + discr = DiscretizationCollection(actx, mesh, quad_tag_to_group_factory=quad_tag_to_group_factory) volume_discr = discr.discr_from_dd(sym.DD_VOLUME) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 9aa4d5b5faf5535fd0b0b1dbd4c958a690fbfec0..35905d7a14611eca799ae918ef036fef92d29ea6 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -144,8 +144,8 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=False): else: quad_tag_to_group_factory = {} - from grudge import DGDiscretizationWithBoundaries - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order, + from grudge import DiscretizationCollection + discr = DiscretizationCollection(actx, mesh, order=order, quad_tag_to_group_factory=quad_tag_to_group_factory) # }}} diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 6e4edd0170175f4fb705d235317c6d80f997c752..0fbf6eda064e334c4b79b79991b3f10bcb5d6777 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -127,8 +127,8 @@ def main(ctx_factory, dim=2, order=4, visualize=False): [np.linspace(-d/2, d/2, npoints) for _ in range(dim)], order=order) - from grudge import DGDiscretizationWithBoundaries - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order) + from grudge import DiscretizationCollection + discr = DiscretizationCollection(actx, mesh, order=order) # }}} diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py index 51033dba8fdaf00aff280c6c16c732ba35fc50b1..779a783cc991f5a5be1979142c1d45a2bd27789b 100755 --- a/examples/dagrt-fusion.py +++ b/examples/dagrt-fusion.py @@ -72,7 +72,7 @@ from pymbolic.mapper.evaluator import EvaluationMapper \ from pytools import memoize from pytools.obj_array import flat_obj_array -from grudge import sym, bind, DGDiscretizationWithBoundaries +from grudge import sym, bind, DiscretizationCollection from leap.rk import LSRK4MethodBuilder from pyopencl.tools import ( # noqa @@ -457,7 +457,7 @@ def get_wave_op_with_discr(actx, dims=2, order=4): logger.debug("%d elements", mesh.nelements) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order) + discr = DiscretizationCollection(actx, mesh, order=order) from grudge.models.wave import WeakWaveOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE @@ -634,7 +634,7 @@ class ExecutionMapperWithMemOpCounting(ExecutionMapperWrapper): def map_insn_loopy_kernel(self, insn, profile_data): kdescr = insn.kernel_descriptor - discr = self.inner_mapper.discrwb.discr_from_dd(kdescr.governing_dd) + discr = self.inner_mapper.dcoll.discr_from_dd(kdescr.governing_dd) dof_array_kwargs = {} other_kwargs = {} @@ -705,7 +705,7 @@ class ExecutionMapperWithMemOpCounting(ExecutionMapperWrapper): logger.debug("assignment not profiled: %s <- %s", name, expr) inner_mapper = self.inner_mapper value = inner_mapper.rec(expr) - inner_mapper.discrwb._discr_scoped_subexpr_name_to_value[name] = value + inner_mapper.dcoll._discr_scoped_subexpr_name_to_value[name] = value assignments.append((name, value)) return assignments, [] @@ -714,7 +714,7 @@ class ExecutionMapperWithMemOpCounting(ExecutionMapperWrapper): return [( insn.name, self.inner_mapper. - discrwb._discr_scoped_subexpr_name_to_value[insn.name])], [] + dcoll._discr_scoped_subexpr_name_to_value[insn.name])], [] def map_insn_rank_data_swap(self, insn, profile_data): raise NotImplementedError("no profiling for instruction: %s" % insn) diff --git a/examples/geometry.py b/examples/geometry.py index 0affdc8ecf5ba887a409758152c58420609f74db..60c7c828bc6053bbc7d9fae868ab5139695c0f26 100644 --- a/examples/geometry.py +++ b/examples/geometry.py @@ -25,7 +25,7 @@ THE SOFTWARE. import numpy as np # noqa import pyopencl as cl -from grudge import sym, bind, DGDiscretizationWithBoundaries, shortcuts +from grudge import sym, bind, DiscretizationCollection, shortcuts from meshmode.array_context import PyOpenCLArrayContext @@ -38,7 +38,7 @@ def main(write_output=True): from meshmode.mesh.generation import generate_warped_rect_mesh mesh = generate_warped_rect_mesh(dim=2, order=4, n=6) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=4) + discr = DiscretizationCollection(actx, mesh, order=4) sym_op = sym.normal(sym.BTAG_ALL, mesh.dim) #sym_op = sym.nodes(mesh.dim, where=sym.BTAG_ALL) diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index d31e08219f2fb944e980d93d5f63b9aee754d6d2..fd6dcd0845353fea9610e1f42b49dacc1d320ee6 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -29,7 +29,7 @@ import pyopencl as cl from meshmode.array_context import PyOpenCLArrayContext from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, DGDiscretizationWithBoundaries +from grudge import sym, bind, DiscretizationCollection from grudge.models.em import get_rectangular_cavity_mode @@ -48,7 +48,7 @@ def main(dims, write_output=True, order=4): b=(1.0,)*dims, n=(5,)*dims) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order) + discr = DiscretizationCollection(actx, mesh, order=order) if 0: epsilon0 = 8.8541878176e-12 # C**2 / (N m**2) diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index 6a075650184f94e7ee820806c7706fb337e7278c..1389af3eef335ff8b6220f3a4fb36a8d81f696fb 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -26,7 +26,7 @@ THE SOFTWARE. import numpy as np import pyopencl as cl from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, DGDiscretizationWithBoundaries +from grudge import sym, bind, DiscretizationCollection from meshmode.array_context import PyOpenCLArrayContext @@ -43,7 +43,7 @@ def main(write_output=True, order=4): b=(0.5,)*dims, n=(20,)*dims) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order) + discr = DiscretizationCollection(actx, mesh, order=order) source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim] source_width = 0.05 diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py index a65582b7dd1227868a3ae81f67ca08b1c37b4237..354adc3a376f76cb34cae8f6d415889096e1f1ce 100644 --- a/examples/wave/wave-min-mpi.py +++ b/examples/wave/wave-min-mpi.py @@ -27,7 +27,7 @@ import numpy as np import pyopencl as cl from meshmode.array_context import PyOpenCLArrayContext from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, DGDiscretizationWithBoundaries +from grudge import sym, bind, DiscretizationCollection from mpi4py import MPI @@ -61,7 +61,7 @@ def main(write_output=True, order=4): else: local_mesh = mesh_dist.receive_mesh_part() - discr = DGDiscretizationWithBoundaries(actx, local_mesh, order=order, + discr = DiscretizationCollection(actx, local_mesh, order=order, mpi_communicator=comm) if local_mesh.dim == 2: diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py index c5d0e6017274f60996e9db838d0b4707e00ca4e2..05dbd758044dc0e269dca2c9def33a1f6c5aac28 100644 --- a/examples/wave/wave-min.py +++ b/examples/wave/wave-min.py @@ -28,7 +28,7 @@ import numpy as np import pyopencl as cl from meshmode.array_context import PyOpenCLArrayContext from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, DGDiscretizationWithBoundaries +from grudge import sym, bind, DiscretizationCollection def main(write_output=True, order=4): @@ -50,7 +50,7 @@ def main(write_output=True, order=4): print("%d elements" % mesh.nelements) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order) + discr = DiscretizationCollection(actx, mesh, order=order) source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim] source_width = 0.05 diff --git a/examples/wave/wave-eager-mpi.py b/examples/wave/wave-op-mpi.py similarity index 77% rename from examples/wave/wave-eager-mpi.py rename to examples/wave/wave-op-mpi.py index b6ca9ae386874607fdaf31533952a013dbbe8c1e..2157927a15847503f054bdced305c8ffcacd6743 100644 --- a/examples/wave/wave-eager-mpi.py +++ b/examples/wave/wave-op-mpi.py @@ -32,8 +32,8 @@ from meshmode.dof_array import thaw from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import ( - EagerDGDiscretization, interior_trace_pair, cross_rank_trace_pairs) +from grudge.discretization import DiscretizationCollection +import grudge.op as op from grudge.shortcuts import make_visualizer from grudge.symbolic.primitives import TracePair from mpi4py import MPI @@ -41,11 +41,11 @@ from mpi4py import MPI # {{{ wave equation bits -def wave_flux(discr, c, w_tpair): +def wave_flux(dcoll, c, w_tpair): u = w_tpair[0] v = w_tpair[1:] - normal = thaw(u.int.array_context, discr.normal(w_tpair.dd)) + normal = thaw(u.int.array_context, op.normal(dcoll, w_tpair.dd)) flux_weak = flat_obj_array( np.dot(v.avg, normal), @@ -59,32 +59,32 @@ def wave_flux(discr, c, w_tpair): 0.5*normal*v_jump, ) - return discr.project(w_tpair.dd, "all_faces", c*flux_weak) + return op.project(dcoll, w_tpair.dd, "all_faces", c*flux_weak) -def wave_operator(discr, c, w): +def wave_operator(dcoll, c, w): u = w[0] v = w[1:] - dir_u = discr.project("vol", BTAG_ALL, u) - dir_v = discr.project("vol", BTAG_ALL, v) + dir_u = op.project(dcoll, "vol", BTAG_ALL, u) + dir_v = op.project(dcoll, "vol", BTAG_ALL, v) dir_bval = flat_obj_array(dir_u, dir_v) dir_bc = flat_obj_array(-dir_u, dir_v) return ( - discr.inverse_mass( + op.inverse_mass(dcoll, flat_obj_array( - -c*discr.weak_div(v), - -c*discr.weak_grad(u) + -c*op.weak_local_div(dcoll, v), + -c*op.weak_local_grad(dcoll, u) ) + # noqa: W504 - discr.face_mass( - wave_flux(discr, c=c, w_tpair=interior_trace_pair(discr, w)) - + wave_flux(discr, c=c, w_tpair=TracePair( + op.face_mass(dcoll, + wave_flux(dcoll, c=c, w_tpair=op.interior_trace_pair(dcoll, w)) + + wave_flux(dcoll, c=c, w_tpair=TracePair( 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)) + wave_flux(dcoll, c=c, w_tpair=tpair) + for tpair in op.cross_rank_trace_pairs(dcoll, w)) ) ) ) @@ -100,15 +100,15 @@ def rk4_step(y, t, h, f): return y + h/6*(k1 + 2*k2 + 2*k3 + k4) -def bump(actx, discr, t=0): - source_center = np.array([0.2, 0.35, 0.1])[:discr.dim] +def bump(actx, dcoll, t=0): + source_center = np.array([0.2, 0.35, 0.1])[:dcoll.dim] source_width = 0.05 source_omega = 3 - nodes = thaw(actx, discr.nodes()) + nodes = thaw(actx, op.nodes(dcoll)) center_dist = flat_obj_array([ nodes[i] - source_center[i] - for i in range(discr.dim) + for i in range(dcoll.dim) ]) return ( @@ -152,7 +152,7 @@ def main(): order = 3 - discr = EagerDGDiscretization(actx, local_mesh, order=order, + dcoll = DiscretizationCollection(actx, local_mesh, order=order, mpi_communicator=comm) if dim == 2: @@ -165,14 +165,14 @@ def main(): raise ValueError("don't have a stable time step guesstimate") fields = flat_obj_array( - bump(actx, discr), - [discr.zeros(actx) for i in range(discr.dim)] + bump(actx, dcoll), + [dcoll.zeros(actx) for i in range(dcoll.dim)] ) - vis = make_visualizer(discr, order+3 if dim == 2 else order) + vis = make_visualizer(dcoll, order+3 if dim == 2 else order) def rhs(t, w): - return wave_operator(discr, c=1, w=w) + return wave_operator(dcoll, c=1, w=w) t = 0 t_final = 3 @@ -181,7 +181,7 @@ def main(): fields = rk4_step(fields, t, dt, rhs) if istep % 10 == 0: - print(istep, t, discr.norm(fields[0])) + print(istep, t, op.norm(dcoll, fields[0], p=2)) vis.write_parallel_vtk_file( comm, f"fld-wave-eager-mpi-{{rank:03d}}-{istep:04d}.vtu", diff --git a/examples/wave/wave-eager-var-velocity.py b/examples/wave/wave-op-var-velocity.py similarity index 75% rename from examples/wave/wave-eager-var-velocity.py rename to examples/wave/wave-op-var-velocity.py index 194a1d654e602d5add77ea6d5722849001ca87b0..b531c939e5b71da5e76c39c92579f2abb9759a2c 100644 --- a/examples/wave/wave-eager-var-velocity.py +++ b/examples/wave/wave-op-var-velocity.py @@ -32,21 +32,22 @@ from meshmode.dof_array import thaw from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization, interior_trace_pair +from grudge.discretization import DiscretizationCollection +import grudge.op as op from grudge.shortcuts import make_visualizer from grudge.symbolic.primitives import TracePair, QTAG_NONE, DOFDesc # {{{ wave equation bits -def wave_flux(discr, c, w_tpair): +def wave_flux(dcoll, c, w_tpair): dd = w_tpair.dd dd_quad = dd.with_qtag("vel_prod") u = w_tpair[0] v = w_tpair[1:] - normal = thaw(u.int.array_context, discr.normal(dd)) + normal = thaw(u.int.array_context, op.normal(dcoll, dd)) flux_weak = flat_obj_array( np.dot(v.avg, normal), @@ -61,40 +62,40 @@ def wave_flux(discr, c, w_tpair): # FIXME this flux is only correct for continuous c dd_allfaces_quad = dd_quad.with_dtag("all_faces") - c_quad = discr.project("vol", dd_quad, c) - flux_quad = discr.project(dd, dd_quad, flux_weak) + c_quad = op.project(dcoll, "vol", dd_quad, c) + flux_quad = op.project(dcoll, dd, dd_quad, flux_weak) - return discr.project(dd_quad, dd_allfaces_quad, c_quad*flux_quad) + return op.project(dcoll, dd_quad, dd_allfaces_quad, c_quad*flux_quad) -def wave_operator(discr, c, w): +def wave_operator(dcoll, c, w): u = w[0] v = w[1:] - dir_u = discr.project("vol", BTAG_ALL, u) - dir_v = discr.project("vol", BTAG_ALL, v) + dir_u = op.project(dcoll, "vol", BTAG_ALL, u) + dir_v = op.project(dcoll, "vol", BTAG_ALL, v) dir_bval = flat_obj_array(dir_u, dir_v) dir_bc = flat_obj_array(-dir_u, dir_v) dd_quad = DOFDesc("vol", "vel_prod") - c_quad = discr.project("vol", dd_quad, c) - w_quad = discr.project("vol", dd_quad, w) + c_quad = op.project(dcoll, "vol", dd_quad, c) + w_quad = op.project(dcoll, "vol", dd_quad, w) u_quad = w_quad[0] v_quad = w_quad[1:] dd_allfaces_quad = DOFDesc("all_faces", "vel_prod") return ( - discr.inverse_mass( + op.inverse_mass(dcoll, flat_obj_array( - -discr.weak_div(dd_quad, c_quad*v_quad), - -discr.weak_grad(dd_quad, c_quad*u_quad) + -op.weak_local_div(dcoll, dd_quad, c_quad*v_quad), + -op.weak_local_grad(dcoll, dd_quad, c_quad*u_quad) ) + # noqa: W504 - discr.face_mass( + op.face_mass(dcoll, dd_allfaces_quad, - wave_flux(discr, c=c, w_tpair=interior_trace_pair(discr, w)) - + wave_flux(discr, c=c, w_tpair=TracePair( + wave_flux(dcoll, c=c, w_tpair=op.interior_trace_pair(dcoll, w)) + + wave_flux(dcoll, c=c, w_tpair=TracePair( BTAG_ALL, interior=dir_bval, exterior=dir_bc)) )) ) @@ -110,17 +111,17 @@ def rk4_step(y, t, h, f): return y + h/6*(k1 + 2*k2 + 2*k3 + k4) -def bump(actx, discr, t=0, width=0.05, center=None): +def bump(actx, dcoll, t=0, width=0.05, center=None): if center is None: center = np.array([0.2, 0.35, 0.1]) - center = center[:discr.dim] + center = center[:dcoll.dim] source_omega = 3 - nodes = thaw(actx, discr.nodes()) + nodes = thaw(actx, op.nodes(dcoll)) center_dist = flat_obj_array([ nodes[i] - center[i] - for i in range(discr.dim) + for i in range(dcoll.dim) ]) return ( @@ -159,24 +160,24 @@ def main(): from meshmode.discretization.poly_element import \ QuadratureSimplexGroupFactory, \ PolynomialWarpAndBlendGroupFactory - discr = EagerDGDiscretization(actx, mesh, + dcoll = DiscretizationCollection(actx, mesh, quad_tag_to_group_factory={ QTAG_NONE: PolynomialWarpAndBlendGroupFactory(order), "vel_prod": QuadratureSimplexGroupFactory(3*order), }) # bounded above by 1 - c = 0.2 + 0.8*bump(actx, discr, center=np.zeros(3), width=0.5) + c = 0.2 + 0.8*bump(actx, dcoll, center=np.zeros(3), width=0.5) fields = flat_obj_array( - bump(actx, discr, ), - [discr.zeros(actx) for i in range(discr.dim)] + bump(actx, dcoll, ), + [dcoll.zeros(actx) for i in range(dcoll.dim)] ) - vis = make_visualizer(discr, order+3 if dim == 2 else order) + vis = make_visualizer(dcoll, order+3 if dim == 2 else order) def rhs(t, w): - return wave_operator(discr, c=c, w=w) + return wave_operator(dcoll, c=c, w=w) t = 0 t_final = 3 @@ -185,7 +186,7 @@ def main(): fields = rk4_step(fields, t, dt, rhs) if istep % 10 == 0: - print(istep, t, discr.norm(fields[0])) + print(istep, t, op.norm(dcoll, fields[0], p=2)) vis.write_vtk_file("fld-wave-eager-var-velocity-%04d.vtu" % istep, [ ("c", c), diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-op.py similarity index 74% rename from examples/wave/wave-eager.py rename to examples/wave/wave-op.py index a8c4438217f5a09c03ddd9c00f3e1209b2edd49a..50a1f40102cad3eed3cca2ecaffe9c9d212f5e20 100644 --- a/examples/wave/wave-eager.py +++ b/examples/wave/wave-op.py @@ -32,18 +32,19 @@ from meshmode.dof_array import thaw from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization, interior_trace_pair +from grudge.discretization import DiscretizationCollection +import grudge.op as op from grudge.shortcuts import make_visualizer from grudge.symbolic.primitives import TracePair # {{{ wave equation bits -def wave_flux(discr, c, w_tpair): +def wave_flux(dcoll, c, w_tpair): u = w_tpair[0] v = w_tpair[1:] - normal = thaw(u.int.array_context, discr.normal(w_tpair.dd)) + normal = thaw(u.int.array_context, op.normal(dcoll, w_tpair.dd)) flux_weak = flat_obj_array( np.dot(v.avg, normal), @@ -56,28 +57,28 @@ def wave_flux(discr, c, w_tpair): 0.5*normal*np.dot(normal, v.ext-v.int), ) - return discr.project(w_tpair.dd, "all_faces", c*flux_weak) + return op.project(dcoll, w_tpair.dd, "all_faces", c*flux_weak) -def wave_operator(discr, c, w): +def wave_operator(dcoll, c, w): u = w[0] v = w[1:] - dir_u = discr.project("vol", BTAG_ALL, u) - dir_v = discr.project("vol", BTAG_ALL, v) + dir_u = op.project(dcoll, "vol", BTAG_ALL, u) + dir_v = op.project(dcoll, "vol", BTAG_ALL, v) dir_bval = flat_obj_array(dir_u, dir_v) dir_bc = flat_obj_array(-dir_u, dir_v) return ( - discr.inverse_mass( + op.inverse_mass(dcoll, flat_obj_array( - -c*discr.weak_div(v), - -c*discr.weak_grad(u) + -c*op.weak_local_div(dcoll, v), + -c*op.weak_local_grad(dcoll, u) ) + # noqa: W504 - discr.face_mass( - wave_flux(discr, c=c, w_tpair=interior_trace_pair(discr, w)) - + wave_flux(discr, c=c, w_tpair=TracePair( + op.face_mass(dcoll, + wave_flux(dcoll, c=c, w_tpair=op.interior_trace_pair(dcoll, w)) + + wave_flux(dcoll, c=c, w_tpair=TracePair( BTAG_ALL, interior=dir_bval, exterior=dir_bc)) )) ) @@ -93,15 +94,15 @@ def rk4_step(y, t, h, f): return y + h/6*(k1 + 2*k2 + 2*k3 + k4) -def bump(actx, discr, t=0): - source_center = np.array([0.2, 0.35, 0.1])[:discr.dim] +def bump(actx, dcoll, t=0): + source_center = np.array([0.2, 0.35, 0.1])[:dcoll.dim] source_width = 0.05 source_omega = 3 - nodes = thaw(actx, discr.nodes()) + nodes = thaw(actx, op.nodes(dcoll)) center_dist = flat_obj_array([ nodes[i] - source_center[i] - for i in range(discr.dim) + for i in range(dcoll.dim) ]) return ( @@ -137,17 +138,17 @@ def main(): print("%d elements" % mesh.nelements) - discr = EagerDGDiscretization(actx, mesh, order=order) + dcoll = DiscretizationCollection(actx, mesh, order=order) fields = flat_obj_array( - bump(actx, discr), - [discr.zeros(actx) for i in range(discr.dim)] + bump(actx, dcoll), + [dcoll.zeros(actx) for i in range(dcoll.dim)] ) - vis = make_visualizer(discr, order+3 if dim == 2 else order) + vis = make_visualizer(dcoll, order+3 if dim == 2 else order) def rhs(t, w): - return wave_operator(discr, c=1, w=w) + return wave_operator(dcoll, c=1, w=w) t = 0 t_final = 3 @@ -156,8 +157,8 @@ def main(): fields = rk4_step(fields, t, dt, rhs) if istep % 10 == 0: - print(f"step: {istep} t: {t} L2: {discr.norm(fields[0])} " - f"sol max: {discr.nodal_max('vol', fields[0])}") + print(f"step: {istep} t: {t} L2: {op.norm(dcoll, fields[0], 2)} " + f"sol max: {op.nodal_max(dcoll, 'vol', fields[0])}") vis.write_vtk_file("fld-wave-eager-%04d.vtu" % istep, [ ("u", fields[0]), diff --git a/grudge/__init__.py b/grudge/__init__.py index cbc7de30a8847e908ff013038d69d3f022e5061d..04cccfb1392d573828770bc6a87beb7749c1c791 100644 --- a/grudge/__init__.py +++ b/grudge/__init__.py @@ -23,8 +23,8 @@ THE SOFTWARE. import grudge.symbolic as sym from grudge.execution import bind -from grudge.discretization import DGDiscretizationWithBoundaries +from grudge.discretization import DiscretizationCollection __all__ = [ - "sym", "bind", "DGDiscretizationWithBoundaries" + "sym", "bind", "DiscretizationCollection" ] diff --git a/grudge/discretization.py b/grudge/discretization.py index 54bad00c7159f1d6cd1b5f330abd0b0dfd413c88..01951a6a6d4a723659b594f1b587c015b7317015 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -27,11 +27,11 @@ from meshmode.array_context import ArrayContext __doc__ = """ -.. autoclass:: DGDiscretizationWithBoundaries +.. autoclass:: DiscretizationCollection """ -class DGDiscretizationWithBoundaries: +class DiscretizationCollection: """ .. automethod :: discr_from_dd .. automethod :: connection_from_dds @@ -390,7 +390,7 @@ class DGDiscretizationWithBoundaries: @property def order(self): from warnings import warn - warn("DGDiscretizationWithBoundaries.order is deprecated, " + warn("DiscretizationCollection.order is deprecated, " "consider using the orders of element groups instead. " "'order' will go away in 2021.", DeprecationWarning, stacklevel=2) @@ -399,4 +399,11 @@ class DGDiscretizationWithBoundaries: return single_valued(egrp.order for egrp in self._volume_discr.groups) +class DGDiscretizationWithBoundaries(DiscretizationCollection): + def __init__(self, *args, **kwargs): + from warnings import warn + warn("DGDiscretizationWithBoundaries is deprecated and will go away " + "in 2022. Use DiscretizationCollection instead.", + DeprecationWarning, stacklevel=2) + # vim: foldmethod=marker diff --git a/grudge/eager.py b/grudge/eager.py index 3314c942ba7fa04e106da3426ddf00f4a52050bd..cb84df6ee0db9a5e80500c7432cd9e11f535d01c 100644 --- a/grudge/eager.py +++ b/grudge/eager.py @@ -20,452 +20,77 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +import grudge.op as op +from grudge.discretization import DiscretizationCollection -import numpy as np # noqa -from pytools import memoize_method -from pytools.obj_array import obj_array_vectorize, make_obj_array -import pyopencl.array as cla # noqa -from grudge import sym, bind -from meshmode.mesh import BTAG_ALL, BTAG_NONE, BTAG_PARTITION # noqa -from meshmode.dof_array import freeze, flatten, unflatten - -from grudge.discretization import DGDiscretizationWithBoundaries -from grudge.symbolic.primitives import TracePair - -from numbers import Number - -__doc__ = """ -.. autoclass:: EagerDGDiscretization -.. autofunction:: interior_trace_pair -.. autofunction:: cross_rank_trace_pairs -""" - - -class EagerDGDiscretization(DGDiscretizationWithBoundaries): +class EagerDGDiscretization(DiscretizationCollection): """ - Inherits from :class:`~grudge.discretization.DGDiscretizationWithBoundaries`. + Inherits from :class:`~grudge.discretization.DiscretizationCollection`. .. 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:: mass - .. automethod:: inverse_mass - .. automethod:: face_mass - - .. automethod:: norm - .. automethod:: nodal_sum - .. automethod:: nodal_min - .. automethod:: nodal_max """ - def interp(self, src, tgt, vec): + def __init__(self, *args, **kwargs): from warnings import warn - warn("using 'interp' is deprecated, use 'project' instead.", + warn("EagerDGDiscretization is deprecated and will go away in 2022. " + "Use the base DiscretizationCollection with grudge.op " + "instead.", DeprecationWarning, stacklevel=2) - return self.project(src, tgt, vec) + super().__init__(*args, **kwargs) 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` - """ - src = sym.as_dofdesc(src) - tgt = sym.as_dofdesc(tgt) - if src == tgt: - return vec - - if isinstance(vec, np.ndarray): - return obj_array_vectorize( - lambda el: self.project(src, tgt, el), vec) - - if isinstance(vec, Number): - return vec - - return self.connection_from_dds(src, tgt)(vec) + return op.project(self, src, tgt, vec) 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 + return op.nodes(self, dd) - @memoize_method - def _bound_grad(self): - return bind(self, sym.nabla(self.dim) * sym.Variable("u"), local_only=True) + def normal(self, dd): + return op.normal(self, dd) 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) + return op.local_grad(self, vec) 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 == 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 + return op.local_d_dx(self, xyz_axis, vec) def div(self, 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): - return bind(self, - sym.stiffness_t(self.dim, dd_in=dd) * sym.Variable("u", dd), - local_only=True) + return op.local_d_dx(self, vecs) 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) - elif len(args) == 2: - dd, vec = args - else: - raise TypeError("invalid number of arguments") - - 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) + return op.weak_local_grad(self, *args) 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) + return op.weak_local_d_dx(self, *args) 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) - elif len(args) == 2: - dd, vecs = args - else: - raise TypeError("invalid number of arguments") - - return self._div_helper( - lambda i, subvec: self.weak_d_dx(dd, i, subvec), - vecs) - - # }}} - - @memoize_method - def normal(self, dd): - """Get unit normal to specified surface discretization, *dd*. - - :arg dd: a :class:`~grudge.sym.DOFDesc` as the surface discretization. - :returns: an object array of :class:`~meshmode.dof_array.DOFArray`. - """ - surface_discr = self.discr_from_dd(dd) - actx = surface_discr._setup_actx - return freeze( - bind(self, - sym.normal(dd, surface_discr.ambient_dim, surface_discr.dim), - local_only=True) - (array_context=actx)) - - @memoize_method - def _bound_mass(self, dd): - return bind(self, sym.MassOperator(dd_in=dd)(sym.Variable("u", dd)), - local_only=True) + return op.weak_local_div(self, *args) def mass(self, *args): - if len(args) == 1: - vec, = args - dd = sym.DOFDesc("vol", sym.QTAG_NONE) - elif len(args) == 2: - dd, vec = args - else: - raise TypeError("invalid number of arguments") - - if isinstance(vec, np.ndarray): - return obj_array_vectorize( - lambda el: self.mass(dd, el), vec) - - return self._bound_mass(dd)(u=vec) - - @memoize_method - def _bound_inverse_mass(self): - return bind(self, sym.InverseMassOperator()(sym.Variable("u")), - local_only=True) + return op.mass(self, *args) def inverse_mass(self, vec): - if isinstance(vec, np.ndarray): - return obj_array_vectorize( - lambda el: self.inverse_mass(el), vec) - - return self._bound_inverse_mass()(u=vec) - - @memoize_method - def _bound_face_mass(self, dd): - u = sym.Variable("u", dd=dd) - return bind(self, sym.FaceMassOperator(dd_in=dd)(u), local_only=True) + return op.inverse_mass(self, vec) def face_mass(self, *args): - if len(args) == 1: - vec, = args - dd = sym.DOFDesc("all_faces", sym.QTAG_NONE) - elif len(args) == 2: - dd, vec = args - else: - raise TypeError("invalid number of arguments") - - if isinstance(vec, np.ndarray): - return obj_array_vectorize( - lambda el: self.face_mass(dd, el), vec) - - return self._bound_face_mass(dd)(u=vec) - - @memoize_method - def _norm(self, p, dd): - return bind(self, - sym.norm(p, sym.var("arg", dd=dd), dd=dd), - local_only=True) + return op.face_mass(self, *args) def norm(self, vec, p=2, dd=None): - if dd is None: - dd = "vol" - - dd = sym.as_dofdesc(dd) - - if isinstance(vec, np.ndarray): - if p == 2: - return sum( - self.norm(vec[idx], dd=dd)**2 - for idx in np.ndindex(vec.shape))**0.5 - elif p == np.inf: - return max( - self.norm(vec[idx], np.inf, dd=dd) - for idx in np.ndindex(vec.shape)) - else: - raise ValueError("unsupported norm order") - - return self._norm(p, dd)(arg=vec) - - @memoize_method - def _nodal_reduction(self, operator, dd): - return bind(self, operator(dd)(sym.var("arg")), local_only=True) + return op.norm(self, vec, p, dd) def nodal_sum(self, dd, vec): - return self._nodal_reduction(sym.NodalSum, dd)(arg=vec) + return op.nodal_sum(self, dd, vec) def nodal_min(self, dd, vec): - return self._nodal_reduction(sym.NodalMin, dd)(arg=vec) + return op.nodal_min(self, dd, vec) def nodal_max(self, dd, vec): - return self._nodal_reduction(sym.NodalMax, dd)(arg=vec) - - @memoize_method - def connected_ranks(self): - from meshmode.distributed import get_connected_partitions - return get_connected_partitions(self._volume_discr.mesh) - - -def interior_trace_pair(discrwb, vec): - """Return a :class:`grudge.sym.TracePair` for the interior faces of - *discrwb*. - """ - i = discrwb.project("vol", "int_faces", vec) - - def get_opposite_face(el): - if isinstance(el, Number): - return el - else: - return discrwb.opposite_face_connection()(el) - - e = obj_array_vectorize(get_opposite_face, i) - - return TracePair("int_faces", interior=i, exterior=e) - - -# {{{ distributed-memory functionality - -class _RankBoundaryCommunication: - base_tag = 1273 - - def __init__(self, discrwb, remote_rank, vol_field, tag=None): - self.tag = self.base_tag - if tag is not None: - self.tag += tag - - self.discrwb = discrwb - self.array_context = vol_field.array_context - self.remote_btag = BTAG_PARTITION(remote_rank) - - self.bdry_discr = discrwb.discr_from_dd(self.remote_btag) - self.local_dof_array = discrwb.project("vol", self.remote_btag, vol_field) - - local_data = self.array_context.to_numpy(flatten(self.local_dof_array)) - - comm = self.discrwb.mpi_communicator - - self.send_req = comm.Isend( - local_data, remote_rank, tag=self.tag) - - self.remote_data_host = np.empty_like(local_data) - self.recv_req = comm.Irecv(self.remote_data_host, remote_rank, self.tag) - - def finish(self): - self.recv_req.Wait() - - actx = self.array_context - remote_dof_array = unflatten(self.array_context, self.bdry_discr, - actx.from_numpy(self.remote_data_host)) - - bdry_conn = self.discrwb.get_distributed_boundary_swap_connection( - sym.as_dofdesc(sym.DTAG_BOUNDARY(self.remote_btag))) - swapped_remote_dof_array = bdry_conn(remote_dof_array) - - self.send_req.Wait() - - 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): - if isinstance(vec, Number): - return [TracePair(BTAG_PARTITION(remote_rank), interior=vec, exterior=vec) - for remote_rank in discrwb.connected_ranks()] - else: - rbcomms = [_RankBoundaryCommunication(discrwb, remote_rank, vec, tag=tag) - for remote_rank in discrwb.connected_ranks()] - return [rbcomm.finish() for rbcomm in rbcomms] - - -def cross_rank_trace_pairs(discrwb, vec, tag=None): - if isinstance(vec, np.ndarray): - - n, = vec.shape - result = {} - for ivec in range(n): - for rank_tpair in _cross_rank_trace_pairs_scalar_field( - discrwb, vec[ivec]): - assert isinstance(rank_tpair.dd.domain_tag, sym.DTAG_BOUNDARY) - assert isinstance(rank_tpair.dd.domain_tag.tag, BTAG_PARTITION) - result[rank_tpair.dd.domain_tag.tag.part_nr, ivec] = rank_tpair - - return [ - TracePair( - dd=sym.as_dofdesc(sym.DTAG_BOUNDARY(BTAG_PARTITION(remote_rank))), - interior=make_obj_array([ - result[remote_rank, i].int for i in range(n)]), - exterior=make_obj_array([ - result[remote_rank, i].ext for i in range(n)]) - ) - for remote_rank in discrwb.connected_ranks()] - else: - return _cross_rank_trace_pairs_scalar_field(discrwb, vec, tag=tag) + return op.nodal_max(self, dd, vec) -# }}} +connected_ranks = op.connected_ranks +interior_trace_pair = op.interior_trace_pair +cross_rank_trace_pairs = op.cross_rank_trace_pairs # vim: foldmethod=marker diff --git a/grudge/execution.py b/grudge/execution.py index 3e99556bdaa98e84312539db2dbedae3226f5999..ee33a989e267ffee625bd21475345801e5fb7ff7 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -53,7 +53,7 @@ class ExecutionMapper(mappers.Evaluator, mappers.LocalOpReducerMixin): def __init__(self, array_context, context, bound_op): super().__init__(context) - self.discrwb = bound_op.discrwb + self.dcoll = bound_op.dcoll self.bound_op = bound_op self.function_registry = bound_op.function_registry self.array_context = array_context @@ -64,7 +64,7 @@ class ExecutionMapper(mappers.Evaluator, if expr.dd.is_scalar(): return 1 - discr = self.discrwb.discr_from_dd(expr.dd) + discr = self.dcoll.discr_from_dd(expr.dd) result = discr.empty(self.array_context) for grp_ary in result: @@ -72,7 +72,7 @@ class ExecutionMapper(mappers.Evaluator, return result def map_node_coordinate_component(self, expr): - discr = self.discrwb.discr_from_dd(expr.dd) + discr = self.dcoll.discr_from_dd(expr.dd) return thaw(self.array_context, discr.nodes()[expr.axis]) def map_grudge_variable(self, expr): @@ -80,7 +80,7 @@ class ExecutionMapper(mappers.Evaluator, value = self.context[expr.name] if not expr.dd.is_scalar() and isinstance(value, Number): - discr = self.discrwb.discr_from_dd(expr.dd) + discr = self.dcoll.discr_from_dd(expr.dd) ary = discr.empty(self.array_context) for grp_ary in ary: grp_ary.fill(value) @@ -96,7 +96,7 @@ class ExecutionMapper(mappers.Evaluator, from numbers import Number if not dd.is_scalar() and isinstance(value, Number): - discr = self.discrwb.discr_from_dd(dd) + discr = self.dcoll.discr_from_dd(dd) ary = discr.empty(self.array_context) for grp_ary in ary: grp_ary.fill(value) @@ -123,7 +123,7 @@ class ExecutionMapper(mappers.Evaluator, name="grudge_elementwise_%s" % op_name) field = self.rec(field_expr) - discr = self.discrwb.discr_from_dd(dd) + discr = self.dcoll.discr_from_dd(dd) assert field.shape == (len(discr.groups),) result = discr.empty(self.array_context, dtype=field.entry_dtype) @@ -276,8 +276,8 @@ class ExecutionMapper(mappers.Evaluator, if is_zero(field): return 0 - in_discr = self.discrwb.discr_from_dd(op.dd_in) - out_discr = self.discrwb.discr_from_dd(op.dd_out) + in_discr = self.dcoll.discr_from_dd(op.dd_in) + out_discr = self.dcoll.discr_from_dd(op.dd_out) result = out_discr.empty(self.array_context, dtype=field.entry_dtype) @@ -304,17 +304,17 @@ class ExecutionMapper(mappers.Evaluator, return result def map_projection(self, op, field_expr): - conn = self.discrwb.connection_from_dds(op.dd_in, op.dd_out) + conn = self.dcoll.connection_from_dds(op.dd_in, op.dd_out) return conn(self.rec(field_expr)) def map_opposite_partition_face_swap(self, op, field_expr): assert op.dd_in == op.dd_out - bdry_conn = self.discrwb.get_distributed_boundary_swap_connection(op.dd_in) + bdry_conn = self.dcoll.get_distributed_boundary_swap_connection(op.dd_in) remote_bdry_vec = self.rec(field_expr) # swapped by RankDataSwapAssign return bdry_conn(remote_bdry_vec) def map_opposite_interior_face_swap(self, op, field_expr): - return self.discrwb.opposite_face_connection()(self.rec(field_expr)) + return self.dcoll.opposite_face_connection()(self.rec(field_expr)) # }}} @@ -340,7 +340,7 @@ class ExecutionMapper(mappers.Evaluator, """, name="face_mass") - all_faces_conn = self.discrwb.connection_from_dds("vol", op.dd_in) + all_faces_conn = self.dcoll.connection_from_dds("vol", op.dd_in) all_faces_discr = all_faces_conn.to_discr vol_discr = all_faces_conn.from_discr @@ -374,16 +374,16 @@ class ExecutionMapper(mappers.Evaluator, def map_signed_face_ones(self, expr): assert expr.dd.is_trace() - face_discr = self.discrwb.discr_from_dd(expr.dd) + face_discr = self.dcoll.discr_from_dd(expr.dd) assert face_discr.dim == 0 # NOTE: ignore quadrature_tags on expr.dd, since we only care about # the face_id here - all_faces_conn = self.discrwb.connection_from_dds( + all_faces_conn = self.dcoll.connection_from_dds( sym.DD_VOLUME, sym.DOFDesc(expr.dd.domain_tag)) - field = face_discr.empty(self.array_context, dtype=self.discrwb.real_dtype) + field = face_discr.empty(self.array_context, dtype=self.dcoll.real_dtype) for grp_ary in field: grp_ary.fill(1) @@ -402,7 +402,7 @@ class ExecutionMapper(mappers.Evaluator, def map_insn_rank_data_swap(self, insn, profile_data=None): local_data = self.array_context.to_numpy(flatten(self.rec(insn.field))) - comm = self.discrwb.mpi_communicator + comm = self.dcoll.mpi_communicator # print("Sending data to rank %d with tag %d" # % (insn.i_remote_rank, insn.send_tag)) @@ -414,7 +414,7 @@ class ExecutionMapper(mappers.Evaluator, return [], [ MPIRecvFuture( array_context=self.array_context, - bdry_discr=self.discrwb.discr_from_dd(insn.dd_out), + bdry_discr=self.dcoll.discr_from_dd(insn.dd_out), recv_req=recv_req, insn_name=insn.name, remote_data_host=remote_data_host), @@ -422,7 +422,7 @@ class ExecutionMapper(mappers.Evaluator, def map_insn_loopy_kernel(self, insn, profile_data=None): kdescr = insn.kernel_descriptor - discr = self.discrwb.discr_from_dd(kdescr.governing_dd) + discr = self.dcoll.discr_from_dd(kdescr.governing_dd) dof_array_kwargs = {} other_kwargs = {} @@ -475,14 +475,14 @@ class ExecutionMapper(mappers.Evaluator, assignments = [] for name, expr in zip(insn.names, insn.exprs): value = self.rec(expr) - self.discrwb._discr_scoped_subexpr_name_to_value[name] = value + self.dcoll._discr_scoped_subexpr_name_to_value[name] = value assignments.append((name, value)) return assignments, [] def map_insn_assign_from_discr_scoped(self, insn, profile_data=None): return [(insn.name, - self.discrwb._discr_scoped_subexpr_name_to_value[insn.name])], [] + self.dcoll._discr_scoped_subexpr_name_to_value[insn.name])], [] def map_insn_diff_batch_assign(self, insn, profile_data=None): field = self.rec(insn.field) @@ -490,8 +490,8 @@ class ExecutionMapper(mappers.Evaluator, assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag - in_discr = self.discrwb.discr_from_dd(repr_op.dd_in) - out_discr = self.discrwb.discr_from_dd(repr_op.dd_out) + in_discr = self.dcoll.discr_from_dd(repr_op.dd_in) + out_discr = self.dcoll.discr_from_dd(repr_op.dd_out) prg = self._elwise_linear_loopy_prg() @@ -565,9 +565,9 @@ class MPISendFuture: # {{{ bound operator class BoundOperator: - def __init__(self, discrwb, discr_code, eval_code, debug_flags, + def __init__(self, dcoll, discr_code, eval_code, debug_flags, function_registry, exec_mapper_factory): - self.discrwb = discrwb + self.dcoll = dcoll self.discr_code = discr_code self.eval_code = eval_code self.operator_data_cache = {} @@ -632,17 +632,17 @@ class BoundOperator: # }}} - # {{{ discrwb-scope evaluation + # {{{ dcoll-scope evaluation if any( (result_var.name not in - self.discrwb._discr_scoped_subexpr_name_to_value) + self.dcoll._discr_scoped_subexpr_name_to_value) for result_var in self.discr_code.result): - # need to do discrwb-scope evaluation - discrwb_eval_context: Dict[str, ResultType] = {} + # need to do dcoll-scope evaluation + dcoll_eval_context: Dict[str, ResultType] = {} self.discr_code.execute( self.exec_mapper_factory( - array_context, discrwb_eval_context, self)) + array_context, dcoll_eval_context, self)) # }}} @@ -656,7 +656,7 @@ class BoundOperator: # {{{ process_sym_operator function -def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, dumper=None, +def process_sym_operator(dcoll, sym_operator, post_bind_mapper=None, dumper=None, local_only=None): if local_only is None: local_only = False @@ -671,7 +671,7 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, dumper=No dumper("before-bind", sym_operator) sym_operator = mappers.OperatorBinder()(sym_operator) - mappers.ErrorChecker(discrwb.mesh)(sym_operator) + mappers.ErrorChecker(dcoll.mesh)(sym_operator) sym_operator = \ mappers.OppositeInteriorFaceSwapUniqueIDAssigner()(sym_operator) @@ -681,17 +681,17 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, dumper=No # also make sure all ranks had same orig_sym_operator - if discrwb.mpi_communicator is not None: + if dcoll.mpi_communicator is not None: (mgmt_rank_orig_sym_operator, mgmt_rank_sym_operator) = \ - discrwb.mpi_communicator.bcast( + dcoll.mpi_communicator.bcast( (orig_sym_operator, sym_operator), - discrwb.get_management_rank_index()) + dcoll.get_management_rank_index()) if not np.array_equal(mgmt_rank_orig_sym_operator, orig_sym_operator): raise ValueError("rank %d received a different symbolic " "operator to bind from rank %d" - % (discrwb.mpi_communicator.Get_rank(), - discrwb.get_management_rank_index())) + % (dcoll.mpi_communicator.Get_rank(), + dcoll.get_management_rank_index())) sym_operator = mgmt_rank_sym_operator @@ -702,14 +702,14 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, dumper=No sym_operator = post_bind_mapper(sym_operator) dumper("before-empty-flux-killer", sym_operator) - sym_operator = mappers.EmptyFluxKiller(discrwb.mesh)(sym_operator) + sym_operator = mappers.EmptyFluxKiller(dcoll.mesh)(sym_operator) dumper("before-cfold", sym_operator) sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator) dumper("before-qcheck", sym_operator) sym_operator = mappers.QuadratureCheckerAndRemover( - discrwb.quad_tag_to_group_factory)(sym_operator) + dcoll.quad_tag_to_group_factory)(sym_operator) # Work around https://github.com/numpy/numpy/issues/9438 # @@ -722,17 +722,17 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, dumper=No dumper("before-csize", sym_operator) sym_operator = mappers.ConstantToNumpyConversionMapper( - real_type=discrwb.real_dtype.type, - complex_type=discrwb.complex_dtype.type, + real_type=dcoll.real_dtype.type, + complex_type=dcoll.complex_dtype.type, )(sym_operator) dumper("before-global-to-reference", sym_operator) - sym_operator = mappers.GlobalToReferenceMapper(discrwb)(sym_operator) + sym_operator = mappers.GlobalToReferenceMapper(dcoll)(sym_operator) dumper("before-distributed", sym_operator) if not local_only: - volume_mesh = discrwb.discr_from_dd("vol").mesh + volume_mesh = dcoll.discr_from_dd("vol").mesh from meshmode.distributed import get_connected_partitions connected_parts = get_connected_partitions(volume_mesh) diff --git a/grudge/op.py b/grudge/op.py new file mode 100644 index 0000000000000000000000000000000000000000..cf372f02d3c55cb6d44f492962b1f4e53ea9c5c1 --- /dev/null +++ b/grudge/op.py @@ -0,0 +1,507 @@ +""" +.. autofunction:: project +.. autofunction:: nodes + +.. autofunction:: local_grad +.. autofunction:: local_d_dx +.. autofunction:: local_div + +.. autofunction:: weak_local_grad +.. autofunction:: weak_local_d_dx +.. autofunction:: weak_local_div + +.. autofunction:: normal +.. autofunction:: mass +.. autofunction:: inverse_mass +.. autofunction:: face_mass + +.. autofunction:: norm +.. autofunction:: nodal_sum +.. autofunction:: nodal_min +.. autofunction:: nodal_max + +.. autofunction:: interior_trace_pair +.. autofunction:: cross_rank_trace_pairs +""" + +__copyright__ = "Copyright (C) 2021 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + +from numbers import Number +from pytools import memoize_on_first_arg + +import numpy as np # noqa +from pytools.obj_array import obj_array_vectorize, make_obj_array +import pyopencl.array as cla # noqa +from grudge import sym, bind + +from meshmode.mesh import BTAG_ALL, BTAG_NONE, BTAG_PARTITION # noqa +from meshmode.dof_array import freeze, flatten, unflatten + +from grudge.symbolic.primitives import TracePair + + +# def interp(dcoll, src, tgt, vec): +# from warnings import warn +# warn("using 'interp' is deprecated, use 'project' instead.", +# DeprecationWarning, stacklevel=2) +# +# return dcoll.project(src, tgt, vec) + + +def project(dcoll, 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` + """ + src = sym.as_dofdesc(src) + tgt = sym.as_dofdesc(tgt) + if src == tgt: + return vec + + if isinstance(vec, np.ndarray): + return obj_array_vectorize( + lambda el: project(dcoll, src, tgt, el), vec) + + if isinstance(vec, Number): + return vec + + return dcoll.connection_from_dds(src, tgt)(vec) + + +# {{{ geometric properties + +def nodes(dcoll, 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 dcoll._volume_discr.nodes() + else: + return dcoll.discr_from_dd(dd).nodes() + + +@memoize_on_first_arg +def normal(dcoll, dd): + """Get unit normal to specified surface discretization, *dd*. + + :arg dd: a :class:`~grudge.sym.DOFDesc` as the surface discretization. + :returns: an object array of :class:`~meshmode.dof_array.DOFArray`. + """ + surface_discr = dcoll.discr_from_dd(dd) + actx = surface_discr._setup_actx + return freeze( + bind(dcoll, + sym.normal(dd, surface_discr.ambient_dim, surface_discr.dim), + local_only=True) + (array_context=actx)) + +# }}} + + +# {{{ derivatives + +@memoize_on_first_arg +def _bound_grad(dcoll): + return bind(dcoll, sym.nabla(dcoll.dim) * sym.Variable("u"), local_only=True) + + +def local_grad(dcoll, vec): + r"""Return the element-local 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 _bound_grad(dcoll)(u=vec) + + +@memoize_on_first_arg +def _bound_d_dx(dcoll, xyz_axis): + return bind(dcoll, sym.nabla(dcoll.dim)[xyz_axis] * sym.Variable("u"), + local_only=True) + + +def local_d_dx(dcoll, xyz_axis, vec): + r"""Return the element-local 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 _bound_d_dx(dcoll, xyz_axis)(u=vec) + + +def _div_helper(dcoll, diff_func, vecs): + if not isinstance(vecs, np.ndarray): + raise TypeError("argument must be an object array") + assert vecs.dtype == object + + if vecs.shape[-1] != dcoll.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 local_div(dcoll, vecs): + r"""Return the element-local 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 _div_helper(dcoll, + lambda i, subvec: local_d_dx(dcoll, i, subvec), + vecs) + + +@memoize_on_first_arg +def _bound_weak_grad(dcoll, dd): + return bind(dcoll, + sym.stiffness_t(dcoll.dim, dd_in=dd) * sym.Variable("u", dd), + local_only=True) + + +def weak_local_grad(dcoll, *args): + r"""Return the element-local 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) + elif len(args) == 2: + dd, vec = args + else: + raise TypeError("invalid number of arguments") + + return _bound_weak_grad(dcoll, dd)(u=vec) + + +@memoize_on_first_arg +def _bound_weak_d_dx(dcoll, dd, xyz_axis): + return bind(dcoll, + sym.stiffness_t(dcoll.dim, dd_in=dd)[xyz_axis] + * sym.Variable("u", dd), + local_only=True) + + +def weak_local_d_dx(dcoll, *args): + r"""Return the element-local weak 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 _bound_weak_d_dx(dcoll, dd, xyz_axis)(u=vec) + + +def weak_local_div(dcoll, *args): + r"""Return the element-local 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) + elif len(args) == 2: + dd, vecs = args + else: + raise TypeError("invalid number of arguments") + + return _div_helper(dcoll, + lambda i, subvec: weak_local_d_dx(dcoll, dd, i, subvec), + vecs) + +# }}} + + +# {{{ mass-like + +@memoize_on_first_arg +def _bound_mass(dcoll, dd): + return bind(dcoll, sym.MassOperator(dd_in=dd)(sym.Variable("u", dd)), + local_only=True) + + +def mass(dcoll, *args): + if len(args) == 1: + vec, = args + dd = sym.DOFDesc("vol", sym.QTAG_NONE) + elif len(args) == 2: + dd, vec = args + else: + raise TypeError("invalid number of arguments") + + if isinstance(vec, np.ndarray): + return obj_array_vectorize( + lambda el: mass(dcoll, dd, el), vec) + + return _bound_mass(dcoll, dd)(u=vec) + + +@memoize_on_first_arg +def _bound_inverse_mass(dcoll): + return bind(dcoll, sym.InverseMassOperator()(sym.Variable("u")), + local_only=True) + + +def inverse_mass(dcoll, vec): + if isinstance(vec, np.ndarray): + return obj_array_vectorize( + lambda el: inverse_mass(dcoll, el), vec) + + return _bound_inverse_mass(dcoll)(u=vec) + + +@memoize_on_first_arg +def _bound_face_mass(dcoll, dd): + u = sym.Variable("u", dd=dd) + return bind(dcoll, sym.FaceMassOperator(dd_in=dd)(u), local_only=True) + + +def face_mass(dcoll, *args): + if len(args) == 1: + vec, = args + dd = sym.DOFDesc("all_faces", sym.QTAG_NONE) + elif len(args) == 2: + dd, vec = args + else: + raise TypeError("invalid number of arguments") + + if isinstance(vec, np.ndarray): + return obj_array_vectorize( + lambda el: face_mass(dcoll, dd, el), vec) + + return _bound_face_mass(dcoll, dd)(u=vec) + +# }}} + + +# {{{ reductions + +@memoize_on_first_arg +def _norm(dcoll, p, dd): + return bind(dcoll, + sym.norm(p, sym.var("arg", dd=dd), dd=dd), + local_only=True) + + +def norm(dcoll, vec, p, dd=None): + if dd is None: + dd = "vol" + + dd = sym.as_dofdesc(dd) + + if isinstance(vec, np.ndarray): + if p == 2: + return sum( + norm(dcoll, vec[idx], dd=dd)**2 + for idx in np.ndindex(vec.shape))**0.5 + elif p == np.inf: + return max( + norm(dcoll, vec[idx], np.inf, dd=dd) + for idx in np.ndindex(vec.shape)) + else: + raise ValueError("unsupported norm order") + + return _norm(dcoll, p, dd)(arg=vec) + + +@memoize_on_first_arg +def _nodal_reduction(dcoll, operator, dd): + return bind(dcoll, operator(dd)(sym.var("arg")), local_only=True) + + +def nodal_sum(dcoll, dd, vec): + return _nodal_reduction(dcoll, sym.NodalSum, dd)(arg=vec) + + +def nodal_min(dcoll, dd, vec): + return _nodal_reduction(dcoll, sym.NodalMin, dd)(arg=vec) + + +def nodal_max(dcoll, dd, vec): + return _nodal_reduction(dcoll, sym.NodalMax, dd)(arg=vec) + +# }}} + + +@memoize_on_first_arg +def connected_ranks(dcoll): + from meshmode.distributed import get_connected_partitions + return get_connected_partitions(dcoll._volume_discr.mesh) + + +# {{{ interior_trace_pair + +def interior_trace_pair(dcoll, vec): + """Return a :class:`grudge.sym.TracePair` for the interior faces of + *dcoll*. + """ + i = project(dcoll, "vol", "int_faces", vec) + + def get_opposite_face(el): + if isinstance(el, Number): + return el + else: + return dcoll.opposite_face_connection()(el) + + e = obj_array_vectorize(get_opposite_face, i) + + return TracePair("int_faces", interior=i, exterior=e) + +# }}} + + +# {{{ distributed-memory functionality + +class _RankBoundaryCommunication: + base_tag = 1273 + + def __init__(self, dcoll, remote_rank, vol_field, tag=None): + self.tag = self.base_tag + if tag is not None: + self.tag += tag + + self.dcoll = dcoll + self.array_context = vol_field.array_context + self.remote_btag = BTAG_PARTITION(remote_rank) + + self.bdry_discr = dcoll.discr_from_dd(self.remote_btag) + self.local_dof_array = project(dcoll, "vol", self.remote_btag, vol_field) + + local_data = self.array_context.to_numpy(flatten(self.local_dof_array)) + + comm = self.dcoll.mpi_communicator + + self.send_req = comm.Isend( + local_data, remote_rank, tag=self.tag) + + self.remote_data_host = np.empty_like(local_data) + self.recv_req = comm.Irecv(self.remote_data_host, remote_rank, self.tag) + + def finish(self): + self.recv_req.Wait() + + actx = self.array_context + remote_dof_array = unflatten(self.array_context, self.bdry_discr, + actx.from_numpy(self.remote_data_host)) + + bdry_conn = self.dcoll.get_distributed_boundary_swap_connection( + sym.as_dofdesc(sym.DTAG_BOUNDARY(self.remote_btag))) + swapped_remote_dof_array = bdry_conn(remote_dof_array) + + self.send_req.Wait() + + return TracePair(self.remote_btag, + interior=self.local_dof_array, + exterior=swapped_remote_dof_array) + + +def _cross_rank_trace_pairs_scalar_field(dcoll, vec, tag=None): + if isinstance(vec, Number): + return [TracePair(BTAG_PARTITION(remote_rank), interior=vec, exterior=vec) + for remote_rank in connected_ranks(dcoll)] + else: + rbcomms = [_RankBoundaryCommunication(dcoll, remote_rank, vec, tag=tag) + for remote_rank in connected_ranks(dcoll)] + return [rbcomm.finish() for rbcomm in rbcomms] + + +def cross_rank_trace_pairs(dcoll, vec, tag=None): + if isinstance(vec, np.ndarray): + + n, = vec.shape + result = {} + for ivec in range(n): + for rank_tpair in _cross_rank_trace_pairs_scalar_field( + dcoll, vec[ivec]): + assert isinstance(rank_tpair.dd.domain_tag, sym.DTAG_BOUNDARY) + assert isinstance(rank_tpair.dd.domain_tag.tag, BTAG_PARTITION) + result[rank_tpair.dd.domain_tag.tag.part_nr, ivec] = rank_tpair + + return [ + TracePair( + dd=sym.as_dofdesc(sym.DTAG_BOUNDARY(BTAG_PARTITION(remote_rank))), + interior=make_obj_array([ + result[remote_rank, i].int for i in range(n)]), + exterior=make_obj_array([ + result[remote_rank, i].ext for i in range(n)]) + ) + for remote_rank in connected_ranks(dcoll)] + else: + return _cross_rank_trace_pairs_scalar_field(dcoll, vec, tag=tag) + +# }}} + + +# vim: foldmethod=marker diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py index a4de4bb7a13d6dc7462f17698a23f858f7af6a2b..d8b1c53d05712eaf95dc764d832c2a8014818b2e 100644 --- a/grudge/shortcuts.py +++ b/grudge/shortcuts.py @@ -39,16 +39,16 @@ def set_up_rk4(field_var_name, dt, fields, rhs, t_start=0): return dt_stepper -def make_visualizer(discrwb, vis_order): +def make_visualizer(dcoll, vis_order): from meshmode.discretization.visualization import make_visualizer return make_visualizer( - discrwb._setup_actx, - discrwb.discr_from_dd("vol"), vis_order) + dcoll._setup_actx, + dcoll.discr_from_dd("vol"), vis_order) -def make_boundary_visualizer(discrwb, vis_order): +def make_boundary_visualizer(dcoll, vis_order): from meshmode.discretization.visualization import make_visualizer from grudge import sym return make_visualizer( - discrwb._setup_actx, discrwb.discr_from_dd(sym.BTAG_ALL), + dcoll._setup_actx, dcoll.discr_from_dd(sym.BTAG_ALL), vis_order) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 555b5ffe8358233bc1d10331a0d2f40f66c55c92..fcbb89ea3e0645b3e6b30fd6a5b83089862fa89b 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -582,14 +582,14 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): reference elements, together with explicit multiplication by geometric factors. """ - def __init__(self, discrwb): + def __init__(self, dcoll): CSECachingMapperMixin.__init__(self) IdentityMapper.__init__(self) - self.ambient_dim = discrwb.ambient_dim - self.dim = discrwb.dim + self.ambient_dim = dcoll.ambient_dim + self.dim = dcoll.dim - volume_discr = discrwb.discr_from_dd(sym.DD_VOLUME) + volume_discr = dcoll.discr_from_dd(sym.DD_VOLUME) self.use_wadg = not all(grp.is_affine for grp in volume_discr.groups) map_common_subexpression_uncached = \ diff --git a/test/test_grudge.py b/test/test_grudge.py index d8c3046eab4413036b8f9e1c62a49a300699a362..edbdc10362b300a4bfc14a3d894a02a276d3f996 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -29,7 +29,7 @@ import meshmode.mesh.generation as mgen from pytools.obj_array import flat_obj_array, make_obj_array -from grudge import sym, bind, DGDiscretizationWithBoundaries +from grudge import sym, bind, DiscretizationCollection import pytest from meshmode.array_context import ( # noqa @@ -65,7 +65,7 @@ def test_inverse_metric(actx_factory, dim): from meshmode.mesh.processing import map_mesh mesh = map_mesh(mesh, m) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=4) + discr = DiscretizationCollection(actx, mesh, order=4) sym_op = ( sym.forward_metric_derivative_mat(mesh.dim) @@ -117,7 +117,7 @@ def test_mass_mat_trig(actx_factory, ambient_dim, quad_tag): mesh = mgen.generate_regular_rect_mesh( a=(a,)*ambient_dim, b=(b,)*ambient_dim, n=(nelements,)*ambient_dim, order=1) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order, + discr = DiscretizationCollection(actx, mesh, order=order, quad_tag_to_group_factory=quad_tag_to_group_factory) def _get_variables_on(dd): @@ -221,7 +221,7 @@ def test_mass_surface_area(actx_factory, name): for resolution in builder.resolutions: mesh = builder.get_mesh(resolution, builder.mesh_order) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=builder.order) + discr = DiscretizationCollection(actx, mesh, order=builder.order) volume_discr = discr.discr_from_dd(sym.DD_VOLUME) logger.info("ndofs: %d", volume_discr.ndofs) @@ -279,7 +279,7 @@ def test_surface_mass_operator_inverse(actx_factory, name): for resolution in builder.resolutions: mesh = builder.get_mesh(resolution, builder.mesh_order) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=builder.order) + discr = DiscretizationCollection(actx, mesh, order=builder.order) volume_discr = discr.discr_from_dd(sym.DD_VOLUME) logger.info("ndofs: %d", volume_discr.ndofs) @@ -335,7 +335,7 @@ def test_face_normal_surface(actx_factory, mesh_name): raise ValueError("unknown mesh name: %s" % mesh_name) mesh = builder.get_mesh(builder.resolutions[0], builder.mesh_order) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=builder.order) + discr = DiscretizationCollection(actx, mesh, order=builder.order) volume_discr = discr.discr_from_dd(sym.DD_VOLUME) logger.info("ndofs: %d", volume_discr.ndofs) @@ -420,7 +420,7 @@ def test_tri_diff_mat(actx_factory, dim, order=4): mesh = mgen.generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim, n=(n,)*dim, order=4) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=4) + discr = DiscretizationCollection(actx, mesh, order=4) nabla = sym.nabla(dim) for axis in range(dim): @@ -465,7 +465,7 @@ def test_2d_gauss_theorem(actx_factory): actx = actx_factory() - discr = DGDiscretizationWithBoundaries(actx, mesh, order=2) + discr = DiscretizationCollection(actx, mesh, order=2) def f(x): return flat_obj_array( @@ -561,7 +561,7 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False): from meshmode.discretization.poly_element import \ QuadratureSimplexGroupFactory - discr = DGDiscretizationWithBoundaries(actx, mesh, order=builder.order, + discr = DiscretizationCollection(actx, mesh, order=builder.order, quad_tag_to_group_factory={ "product": QuadratureSimplexGroupFactory(2 * builder.order) }) @@ -723,7 +723,7 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ from grudge.models.advection import ( StrongAdvectionOperator, WeakAdvectionOperator) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order) + discr = DiscretizationCollection(actx, mesh, order=order) op_class = { "strong": StrongAdvectionOperator, "weak": WeakAdvectionOperator, @@ -809,7 +809,7 @@ def test_convergence_maxwell(actx_factory, order): b=(1.0,)*dims, n=(n,)*dims) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order) + discr = DiscretizationCollection(actx, mesh, order=order) epsilon = 1 mu = 1 @@ -906,7 +906,7 @@ def test_improvement_quadrature(actx_factory, order): else: quad_tag_to_group_factory = {"product": None} - discr = DGDiscretizationWithBoundaries(actx, mesh, order=order, + discr = DiscretizationCollection(actx, mesh, order=order, quad_tag_to_group_factory=quad_tag_to_group_factory) bound_op = bind(discr, op.sym_operator()) @@ -974,7 +974,7 @@ def test_bessel(actx_factory): b=(1.0,)*dims, n=(8,)*dims) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=3) + discr = DiscretizationCollection(actx, mesh, order=3) nodes = sym.nodes(dims) r = sym.cse(sym.sqrt(nodes[0]**2 + nodes[1]**2)) @@ -1005,7 +1005,7 @@ def test_external_call(actx_factory): mesh = mgen.generate_regular_rect_mesh( a=(0,) * dims, b=(1,) * dims, n=(4,) * dims) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=1) + discr = DiscretizationCollection(actx, mesh, order=1) ones = sym.Ones(sym.DD_VOLUME) op = ( @@ -1037,7 +1037,7 @@ def test_function_symbol_array(actx_factory, array_type): mesh = mgen.generate_regular_rect_mesh( a=(-0.5,)*dim, b=(0.5,)*dim, n=(8,)*dim, order=4) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=4) + discr = DiscretizationCollection(actx, mesh, order=4) volume_discr = discr.discr_from_dd(sym.DD_VOLUME) if array_type == "scalar": @@ -1065,7 +1065,7 @@ def test_norm_obj_array(actx_factory, p): mesh = mgen.generate_regular_rect_mesh( a=(-0.5,)*dim, b=(0.5,)*dim, n=(8,)*dim, order=1) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=4) + discr = DiscretizationCollection(actx, mesh, order=4) w = make_obj_array([1.0, 2.0, 3.0])[:dim] @@ -1103,7 +1103,7 @@ def test_map_if(actx_factory): mesh = mgen.generate_regular_rect_mesh( a=(-0.5,)*dim, b=(0.5,)*dim, n=(8,)*dim, order=4) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=4) + discr = DiscretizationCollection(actx, mesh, order=4) sym_if = sym.If(sym.Comparison(2.0, "<", 1.0e-14), 1.0, 2.0) bind(discr, sym_if)(actx) @@ -1118,7 +1118,7 @@ def test_empty_boundary(actx_factory): mesh = mgen.generate_regular_rect_mesh( a=(-0.5,)*dim, b=(0.5,)*dim, n=(8,)*dim, order=4) - discr = DGDiscretizationWithBoundaries(actx, mesh, order=4) + discr = DiscretizationCollection(actx, mesh, order=4) normal = bind(discr, sym.normal(sym.BTAG_NONE, dim, dim=dim - 1))(actx) from meshmode.dof_array import DOFArray diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py index dedc35d16116aa31973c84f245f44dbf479c458c..39d42632d85f4186c1a4b5523daa0262f498c444 100644 --- a/test/test_mpi_communication.py +++ b/test/test_mpi_communication.py @@ -35,7 +35,7 @@ logger = logging.getLogger(__name__) logging.basicConfig() logger.setLevel(logging.INFO) -from grudge import sym, bind, DGDiscretizationWithBoundaries +from grudge import sym, bind, DiscretizationCollection from grudge.shortcuts import set_up_rk4 @@ -64,7 +64,7 @@ def simple_mpi_communication_entrypoint(): else: local_mesh = mesh_dist.receive_mesh_part() - vol_discr = DGDiscretizationWithBoundaries(actx, local_mesh, order=5, + vol_discr = DiscretizationCollection(actx, local_mesh, order=5, mpi_communicator=comm) sym_x = sym.nodes(local_mesh.dim) @@ -126,7 +126,7 @@ def mpi_communication_entrypoint(): else: local_mesh = mesh_dist.receive_mesh_part() - vol_discr = DGDiscretizationWithBoundaries(actx, local_mesh, order=order, + vol_discr = DiscretizationCollection(actx, local_mesh, order=order, mpi_communicator=comm) source_center = np.array([0.1, 0.22, 0.33])[:local_mesh.dim] diff --git a/unported-examples/benchmark_grudge/benchmark_mpi.py b/unported-examples/benchmark_grudge/benchmark_mpi.py index 2d47b6a32231f591e92ae6319607e9efb7b00997..8c60ea3cf867b00432e3f992c7815572ce45f33d 100644 --- a/unported-examples/benchmark_grudge/benchmark_mpi.py +++ b/unported-examples/benchmark_grudge/benchmark_mpi.py @@ -2,7 +2,7 @@ import os import numpy as np import pyopencl as cl -from grudge import sym, bind, DGDiscretizationWithBoundaries +from grudge import sym, bind, DiscretizationCollection from grudge.shortcuts import set_up_rk4 @@ -35,7 +35,7 @@ def simple_wave_entrypoint(dim=2, num_elems=256, order=4, num_steps=30, else: local_mesh = mesh_dist.receive_mesh_part() - vol_discr = DGDiscretizationWithBoundaries(cl_ctx, local_mesh, order=order, + vol_discr = DiscretizationCollection(cl_ctx, local_mesh, order=order, mpi_communicator=comm) source_center = np.array([0.1, 0.22, 0.33])[:local_mesh.dim]