diff --git a/grudge/execution.py b/grudge/execution.py index 66d2cd53d9cd3bddb37156d2613be55669f5f181..5199f8fbbf67550e0b54b0a70ad0a9fc6a9453d3 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -250,7 +250,7 @@ class ExecutionMapper(mappers.Evaluator, if dd_in.is_volume(): if dd_out.domain_tag is sym.FACE_RESTR_ALL: - conn = self.discr.all_faces_connection(qtag) + conn = self.discr.all_faces_volume_connection(qtag) elif dd_out.domain_tag is sym.FACE_RESTR_INTERIOR: conn = self.discr.interior_faces_connection(qtag) elif dd_out.is_boundary(): diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 14049676e008b42256b21724c12e768236a7ff65..d2ef5c66b60d587078857b74fcd01b81c9e54358 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -338,19 +338,6 @@ class OperatorBinder(CSECachingMapperMixin, IdentityMapper): class DistributedMapper(CSECachingMapperMixin, IdentityMapper): map_common_subexpression_uncached = IdentityMapper.map_common_subexpression - def __init__(self, connected_parts): - self.connected_parts = connected_parts - - def map_operator_binding(self, expr): - if isinstance(expr.op, op.RefFaceMassOperator): - return expr.op(RankCommunicationMapper(self.connected_parts)(expr.field)) - else: - return IdentityMapper.map_operator_binding(self, expr) - - -class RankCommunicationMapper(CSECachingMapperMixin, IdentityMapper): - map_common_subexpression_uncached = IdentityMapper.map_common_subexpression - def __init__(self, connected_parts): self.connected_parts = connected_parts @@ -368,11 +355,37 @@ class RankCommunicationMapper(CSECachingMapperMixin, IdentityMapper): distributed_work += op.InterpolationOperator(dd_in=btag_part, dd_out=expr.op.dd_out)(mapped_field) return expr + distributed_work - + # if isinstance(expr.op, op.RefFaceMassOperator): + # return expr.op(RankCommunicationMapper(self.connected_parts)(expr.field)) else: return IdentityMapper.map_operator_binding(self, expr) +# class RankCommunicationMapper(CSECachingMapperMixin, IdentityMapper): +# map_common_subexpression_uncached = IdentityMapper.map_common_subexpression +# +# def __init__(self, connected_parts): +# self.connected_parts = connected_parts +# +# def map_operator_binding(self, expr): +# from meshmode.mesh import BTAG_PARTITION +# from meshmode.discretization.connection import (FACE_RESTR_ALL, +# FACE_RESTR_INTERIOR) +# if (isinstance(expr.op, op.InterpolationOperator) +# and expr.op.dd_in.domain_tag is FACE_RESTR_INTERIOR +# and expr.op.dd_out.domain_tag is FACE_RESTR_ALL): +# distributed_work = 0 +# for i_remote_part in self.connected_parts: +# mapped_field = RankGeometryChanger(i_remote_part)(expr.field) +# btag_part = BTAG_PARTITION(i_remote_part) +# distributed_work += op.InterpolationOperator(dd_in=btag_part, +# dd_out=expr.op.dd_out)(mapped_field) +# return expr + distributed_work +# +# else: +# return IdentityMapper.map_operator_binding(self, expr) + + class RankGeometryChanger(CSECachingMapperMixin, IdentityMapper): map_common_subexpression_uncached = IdentityMapper.map_common_subexpression diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 3d9ddf16f0af207b2c991d2cfd4a97689ad4fe1e..7dc28669a405f21c0ed680ad2d44eafc0f481ed2 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -383,19 +383,19 @@ class OppositePartitionFaceSwap(Operator): def __init__(self, dd_in=None, dd_out=None): sym = _sym() - if dd_in is None and dd_in is None: + if dd_in is None and dd_out is None: raise ValueError("dd_in or dd_out must be specified") elif dd_in is None: dd_in = dd_out elif dd_out is None: dd_out = dd_in - if not isinstance(dd_in.domain_tag, sym.BTAG_PARTITION): + super(OppositePartitionFaceSwap, self).__init__(dd_in, dd_out) + if not isinstance(self.dd_in.domain_tag, sym.BTAG_PARTITION): raise ValueError("dd_in must be a partition boundary faces domain") - if dd_out != dd_in: + if self.dd_out != self.dd_in: raise ValueError("dd_out and dd_in must be identical") - super(OppositePartitionFaceSwap, self).__init__(dd_in, dd_out) self.i_remote_part = dd_in.domain_tag.part_nr mapper_method = intern("map_opposite_partition_face_swap") @@ -410,12 +410,12 @@ class OppositeInteriorFaceSwap(Operator): if dd_out is None: dd_out = dd_in - if dd_in.domain_tag is not sym.FACE_RESTR_INTERIOR: + super(OppositeInteriorFaceSwap, self).__init__(dd_in, dd_out) + if self.dd_in.domain_tag is not sym.FACE_RESTR_INTERIOR: raise ValueError("dd_in must be an interior faces domain") - if dd_out != dd_in: + if self.dd_out != self.dd_in: raise ValueError("dd_out and dd_in must be identical") - super(OppositeInteriorFaceSwap, self).__init__(dd_in, dd_out) mapper_method = intern("map_opposite_interior_face_swap") diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py index 29aab0d954a55c9cad89ea1009024453715dff3a..30710d677921a94c307c6f3549d63fcf713fd37b 100644 --- a/test/test_mpi_communication.py +++ b/test/test_mpi_communication.py @@ -36,6 +36,154 @@ from grudge import sym, bind, Discretization from grudge.shortcuts import set_up_rk4 +# TODO: Make new test +# Create a partitioned mesh and apply sin(2x + 3y) to its field +# If everything is working, the boundaries of the partitions should be continuous +# Look at int_tpair +# Interpolate volume to boundary, ask for the opposite partition at the boundary +# then compare +# def mpi_communication_entrypoint(): +# cl_ctx = cl.create_some_context() +# queue = cl.CommandQueue(cl_ctx) +# from meshmode.distributed import MPIMeshDistributor +# +# from mpi4py import MPI +# comm = MPI.COMM_WORLD +# rank = comm.Get_rank() +# num_parts = comm.Get_size() +# +# mesh_dist = MPIMeshDistributor(comm) +# +# dims = 2 +# dt = 0.04 +# order = 6 +# +# if mesh_dist.is_mananger_rank(): +# 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) +# +# from pymetis import part_graph +# _, p = part_graph(num_parts, +# xadj=mesh.nodal_adjacency.neighbors_starts.tolist(), +# adjncy=mesh.nodal_adjacency.neighbors.tolist()) +# part_per_element = np.array(p) +# +# local_mesh = mesh_dist.send_mesh_parts(mesh, part_per_element, num_parts) +# else: +# local_mesh = mesh_dist.receive_mesh_part() +# +# vol_discr = Discretization(cl_ctx, local_mesh, order=order) +# +# if 0: +# sym_x = sym.nodes(local_mesh.dim) +# myfunc_symb = sym.sin(np.dot(sym_x, [2, 3])) +# myfunc = bind(vol_discr, myfunc_symb)(queue) +# +# sym_all_faces_func = sym.cse( +# sym.interp("vol", "all_faces")(sym.var("myfunc"))) +# sym_int_faces_func = sym.cse( +# sym.interp("vol", "int_faces")(sym.var("myfunc"))) +# sym_bdry_faces_func = sym.cse( +# sym.interp(sym.BTAG_ALL, "all_faces")( +# sym.interp("vol", sym.BTAG_ALL)(sym.var("myfunc")))) +# +# bound_face_swap = bind(vol_discr, +# sym.interp("int_faces", "all_faces")( +# sym.OppositeInteriorFaceSwap("int_faces")( +# sym_int_faces_func) +# ) - (sym_all_faces_func - sym_bdry_faces_func) +# ) +# +# hopefully_zero = bound_face_swap(queue, myfunc=myfunc) +# np.set_printoptions(threshold=100000000, suppress=True) +# print(hopefully_zero) +# +# import numpy.linalg as la +# print(la.norm(hopefully_zero.get())) +# else: +# sym_x = sym.nodes(local_mesh.dim) +# myfunc_symb = sym.sin(np.dot(sym_x, [2, 3])) +# myfunc = bind(vol_discr, myfunc_symb)(queue) +# +# sym_all_faces_func = sym.cse( +# sym.interp("vol", "all_faces")(sym.var("myfunc")) +# - sym.interp(sym.BTAG_ALL, "all_faces")( +# sym.interp("vol", sym.BTAG_ALL)(sym.var("myfunc"))) +# ) +# sym_int_faces_func = sym.cse( +# sym.interp("vol", "int_faces")(sym.var("myfunc"))) +# +# swapped = bind(vol_discr, +# sym.interp("int_faces", "all_faces")( +# sym.OppositeInteriorFaceSwap("int_faces")( +# sym_int_faces_func) +# ))(queue, myfunc=myfunc) +# unswapped = bind(vol_discr, sym_all_faces_func)(queue, myfunc=myfunc) +# +# together = np.zeros((3,)+swapped.shape) +# print(together.shape) +# together[0] = swapped.get() +# together[1] = unswapped.get() +# together[2] = together[1]-together[0] +# +# np.set_printoptions(threshold=100000000, suppress=True, linewidth=150) +# print(together.T) +# +# import numpy.linalg as la +# print(la.norm(hopefully_zero.get())) +# 1/0 +# +# w = sym.make_sym_array("w", vol_discr.dim+1) +# operator = sym.InverseMassOperator()( +# sym.FaceMassOperator()(sym.int_tpair(w))) +# +# # print(sym.pretty(operator) +# bound_op = bind(vol_discr, operator) +# # print(bound_op) +# # 1/0 +# +# def rhs(t, w): +# return bound_op(queue, t=t, w=w) +# +# from pytools.obj_array import join_fields +# fields = join_fields(vol_discr.zeros(queue), +# [vol_discr.zeros(queue) for i in range(vol_discr.dim)]) +# +# dt_stepper = set_up_rk4("w", dt, fields, rhs) +# +# final_t = 10 +# nsteps = int(final_t/dt) +# print("rank=%d dt=%g nsteps=%d" % (rank, dt, nsteps)) +# +# from grudge.shortcuts import make_visualizer +# vis = make_visualizer(vol_discr, vis_order=order) +# +# step = 0 +# +# norm = bind(vol_discr, sym.norm(2, sym.var("u"))) +# +# from time import time +# t_last_step = time() +# +# for event in dt_stepper.run(t_end=final_t): +# if isinstance(event, dt_stepper.StateComputed): +# assert event.component_id == "w" +# +# step += 1 +# +# print(step, event.t, norm(queue, u=event.state_component[0]), +# time()-t_last_step) +# if step % 10 == 0: +# vis.write_vtk_file("rank%d-fld-%04d.vtu" % (rank, step), +# [ +# ("u", event.state_component[0]), +# ("v", event.state_component[1:]), +# ]) +# t_last_step = time() +# logger.debug("Rank %d exiting", rank) + def mpi_communication_entrypoint(): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) @@ -112,9 +260,9 @@ def mpi_communication_entrypoint(): dt_stepper = set_up_rk4("w", dt, fields, rhs) - final_t = 1 + final_t = 10 nsteps = int(final_t/dt) - print("dt=%g nsteps=%d" % (dt, nsteps)) + print("rank=%d dt=%g nsteps=%d" % (rank, dt, nsteps)) from grudge.shortcuts import make_visualizer vis = make_visualizer(vol_discr, vis_order=order)