diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index 5fb31617a2787de8f8334f709cf49f7a75f7d709..ffdb3b84bf03075e3520b7917af4446375886e5f 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -25,13 +25,19 @@ THE SOFTWARE. import numpy as np import pyopencl as cl -from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, DiscretizationCollection from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw + +from grudge.shortcuts import set_up_rk4 +from grudge import DiscretizationCollection + +from pytools.obj_array import flat_obj_array + +import grudge.op as op -def main(write_output=True, order=4): +def main(write_output=False, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) @@ -43,45 +49,52 @@ def main(write_output=True, order=4): b=(0.5,)*dims, nelements_per_axis=(20,)*dims) - discr = DiscretizationCollection(actx, mesh, order=order) - - source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim] - 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") - c = sym.If(sym.Comparison( - np.dot(sym_x, sym_x), "<", 0.15), - np.float32(0.1), np.float32(0.2)) + dcoll = DiscretizationCollection(actx, mesh, order=order) + + def source_f(actx, dcoll, t=0): + source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim] + source_width = 0.05 + source_omega = 3 + nodes = thaw(actx, op.nodes(dcoll)) + source_center_dist = flat_obj_array( + [nodes[i] - source_center[i] for i in range(dcoll.dim)] + ) + return ( + np.sin(source_omega*t) + * actx.np.exp( + -np.dot(source_center_dist, source_center_dist) + / source_width**2 + ) + ) + + x = thaw(actx, op.nodes(dcoll)) + ones = dcoll.zeros(actx) + 1 + c = actx.np.where(np.dot(x, x) < 0.15, + np.float32(0.1)*ones, + np.float32(0.2)*ones) from grudge.models.wave import VariableCoefficientWeakWaveOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE - op = VariableCoefficientWeakWaveOperator(c, - discr.dim, - source_f=( - 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") - from pytools.obj_array import flat_obj_array - fields = flat_obj_array(discr.zeros(actx), - [discr.zeros(actx) for i in range(discr.dim)]) + wave_op = VariableCoefficientWeakWaveOperator( + dcoll, + c, + source_f=source_f, + dirichlet_tag=BTAG_NONE, + neumann_tag=BTAG_NONE, + radiation_tag=BTAG_ALL, + flux_type="upwind" + ) - op.check_bc_coverage(mesh) + fields = flat_obj_array( + dcoll.zeros(actx), + [dcoll.zeros(actx) for i in range(dcoll.dim)] + ) - c_eval = bind(discr, c)(actx) - - bound_op = bind(discr, op.sym_operator()) + wave_op.check_bc_coverage(mesh) def rhs(t, w): - return bound_op(t=t, w=w) + return wave_op.operator(t, w) if mesh.dim == 2: dt = 0.04 * 0.3 @@ -95,15 +108,28 @@ def main(write_output=True, order=4): print("dt=%g nsteps=%d" % (dt, nsteps)) from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr) + vis = make_visualizer(dcoll) step = 0 - norm = bind(discr, sym.norm(2, sym.var("u"))) + def norm(u): + return op.norm(dcoll, u, 2) from time import time t_last_step = time() + if write_output: + u = fields[0] + v = fields[1:] + vis.write_vtk_file( + "fld-var-propogation-speed-%04d.vtu" % step, + [ + ("u", u), + ("v", v), + ("c", c), + ] + ) + for event in dt_stepper.run(t_end=final_t): if isinstance(event, dt_stepper.StateComputed): assert event.component_id == "w" @@ -113,12 +139,15 @@ def main(write_output=True, order=4): if step % 10 == 0: print(f"step: {step} t: {time()-t_last_step} " f"L2: {norm(u=event.state_component[0])}") - vis.write_vtk_file("fld-var-propogation-speed-%04d.vtu" % step, + if write_output: + vis.write_vtk_file( + "fld-var-propogation-speed-%04d.vtu" % step, [ ("u", event.state_component[0]), ("v", event.state_component[1:]), - ("c", c_eval), - ]) + ("c", c), + ] + ) t_last_step = time() assert norm(u=event.state_component[0]) < 1 diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py index 29f59dfff8342227638c6bdc503f7ad79d3251ed..b375f514b620e2865d6410049a2a1c9380040697 100644 --- a/examples/wave/wave-min-mpi.py +++ b/examples/wave/wave-min-mpi.py @@ -25,13 +25,21 @@ THE SOFTWARE. import numpy as np import pyopencl as cl + from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw + from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, DiscretizationCollection +from grudge import DiscretizationCollection + from mpi4py import MPI +from pytools.obj_array import flat_obj_array -def main(write_output=True, order=4): +import grudge.op as op + + +def main(write_output=False, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) @@ -61,7 +69,7 @@ def main(write_output=True, order=4): else: local_mesh = mesh_dist.receive_mesh_part() - discr = DiscretizationCollection(actx, local_mesh, order=order, + dcoll = DiscretizationCollection(actx, local_mesh, order=order, mpi_communicator=comm) if local_mesh.dim == 2: @@ -69,59 +77,79 @@ def main(write_output=True, order=4): elif local_mesh.dim == 3: dt = 0.02 - source_center = np.array([0.1, 0.22, 0.33])[:local_mesh.dim] - source_width = 0.05 - source_omega = 3 - - sym_x = sym.nodes(local_mesh.dim) - sym_source_center_dist = sym_x - source_center - sym_t = sym.ScalarVariable("t") + def source_f(actx, dcoll, t=0): + source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim] + source_width = 0.05 + source_omega = 3 + nodes = thaw(actx, op.nodes(dcoll)) + source_center_dist = flat_obj_array( + [nodes[i] - source_center[i] for i in range(dcoll.dim)] + ) + return ( + np.sin(source_omega*t) + * actx.np.exp( + -np.dot(source_center_dist, source_center_dist) + / source_width**2 + ) + ) from grudge.models.wave import WeakWaveOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE - op = WeakWaveOperator(0.1, discr.dim, - source_f=( - 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") - - from pytools.obj_array import flat_obj_array + + wave_op = WeakWaveOperator( + dcoll, + 0.1, + source_f=source_f, + dirichlet_tag=BTAG_NONE, + neumann_tag=BTAG_NONE, + radiation_tag=BTAG_ALL, + flux_type="upwind" + ) + fields = flat_obj_array( - discr.zeros(actx), - [discr.zeros(actx) for i in range(discr.dim)]) + dcoll.zeros(actx), + [dcoll.zeros(actx) for i in range(dcoll.dim)] + ) # FIXME - #dt = op.estimate_rk4_timestep(discr, fields=fields) - - op.check_bc_coverage(local_mesh) + # dt = wave_op.estimate_rk4_timestep(dcoll, fields=fields) - # print(sym.pretty(op.sym_operator())) - bound_op = bind(discr, op.sym_operator()) + wave_op.check_bc_coverage(local_mesh) def rhs(t, w): - return bound_op(t=t, w=w) + return wave_op.operator(t, w) dt_stepper = set_up_rk4("w", dt, fields, rhs) final_t = 10 nsteps = int(final_t/dt) - print("dt=%g nsteps=%d" % (dt, nsteps)) + + if comm.rank == 0: + print("dt=%g nsteps=%d" % (dt, nsteps)) from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr) + vis = make_visualizer(dcoll) step = 0 - norm = bind(discr, sym.norm(2, sym.var("u"))) + def norm(u): + return op.norm(dcoll, u, 2) from time import time t_last_step = time() + if write_output: + u = fields[0] + v = fields[1:] + vis.write_parallel_vtk_file( + comm, + f"fld-wave-min-mpi-{{rank:03d}}-{step:04d}.vtu", + [ + ("u", u), + ("v", v), + ] + ) + for event in dt_stepper.run(t_end=final_t): if isinstance(event, dt_stepper.StateComputed): assert event.component_id == "w" @@ -129,16 +157,17 @@ def main(write_output=True, order=4): step += 1 if step % 10 == 0: - if comm.rank == 0: - print(f"step: {step} t: {time()-t_last_step} " - f"L2: {norm(u=event.state_component[0])}") - vis.write_parallel_vtk_file( + print(f"step: {step} t: {time()-t_last_step} " + f"L2: {norm(u=event.state_component[0])}") + if write_output: + 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:]), - ]) + ] + ) t_last_step = time() assert norm(u=event.state_component[0]) < 1 diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py index bd306b0fbed822cf428504912b8e024965888150..e995f5ff48cba88838fba6570ac85a9b40cbb081 100644 --- a/examples/wave/wave-min.py +++ b/examples/wave/wave-min.py @@ -26,12 +26,19 @@ THE SOFTWARE. import numpy as np import pyopencl as cl + from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw + from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, DiscretizationCollection +from grudge import DiscretizationCollection + +from pytools.obj_array import flat_obj_array +import grudge.op as op -def main(write_output=True, order=4): + +def main(write_output=False, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) @@ -50,43 +57,49 @@ def main(write_output=True, order=4): print("%d elements" % mesh.nelements) - discr = DiscretizationCollection(actx, mesh, order=order) - - source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim] - 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") + dcoll = DiscretizationCollection(actx, mesh, order=order) + + def source_f(actx, dcoll, t=0): + source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim] + source_width = 0.05 + source_omega = 3 + nodes = thaw(actx, op.nodes(dcoll)) + source_center_dist = flat_obj_array( + [nodes[i] - source_center[i] for i in range(dcoll.dim)] + ) + return ( + np.sin(source_omega*t) + * actx.np.exp( + -np.dot(source_center_dist, source_center_dist) + / source_width**2 + ) + ) from grudge.models.wave import WeakWaveOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE - op = WeakWaveOperator(0.1, discr.dim, - source_f=( - 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") - - from pytools.obj_array import flat_obj_array - fields = flat_obj_array(discr.zeros(actx), - [discr.zeros(actx) for i in range(discr.dim)]) - # FIXME - #dt = op.estimate_rk4_timestep(discr, fields=fields) + wave_op = WeakWaveOperator( + dcoll, + 0.1, + source_f=source_f, + dirichlet_tag=BTAG_NONE, + neumann_tag=BTAG_NONE, + radiation_tag=BTAG_ALL, + flux_type="upwind" + ) + + fields = flat_obj_array( + dcoll.zeros(actx), + [dcoll.zeros(actx) for i in range(dcoll.dim)] + ) - op.check_bc_coverage(mesh) + # FIXME + # dt = wave_op.estimate_rk4_timestep(dcoll, fields=fields) - # print(sym.pretty(op.sym_operator())) - bound_op = bind(discr, op.sym_operator()) + wave_op.check_bc_coverage(mesh) def rhs(t, w): - return bound_op(t=t, w=w) + return wave_op.operator(t, w) dt_stepper = set_up_rk4("w", dt, fields, rhs) @@ -95,15 +108,27 @@ def main(write_output=True, order=4): print("dt=%g nsteps=%d" % (dt, nsteps)) from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr) + vis = make_visualizer(dcoll) step = 0 - norm = bind(discr, sym.norm(2, sym.var("u"))) + def norm(u): + return op.norm(dcoll, u, 2) from time import time t_last_step = time() + if write_output: + u = fields[0] + v = fields[1:] + vis.write_vtk_file( + "fld-wave-min-%04d.vtu" % step, + [ + ("u", u), + ("v", v), + ] + ) + for event in dt_stepper.run(t_end=final_t): if isinstance(event, dt_stepper.StateComputed): assert event.component_id == "w" @@ -113,11 +138,14 @@ def main(write_output=True, order=4): if step % 10 == 0: print(f"step: {step} t: {time()-t_last_step} " f"L2: {norm(u=event.state_component[0])}") - vis.write_vtk_file("fld-wave-min-%04d.vtu" % step, + if write_output: + vis.write_vtk_file( + "fld-wave-min-%04d.vtu" % step, [ ("u", event.state_component[0]), ("v", event.state_component[1:]), - ]) + ] + ) t_last_step = time() assert norm(u=event.state_component[0]) < 1 diff --git a/grudge/models/wave.py b/grudge/models/wave.py index 5fc65a70e2235ea03d33ef95ad7e91e302f313c4..f470a959224ee04787fd736f830c0cb1cf17dd7f 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -23,11 +23,17 @@ THE SOFTWARE. """ import numpy as np + from grudge.models import HyperbolicOperator +from grudge.symbolic.primitives import TracePair + +from meshmode.dof_array import thaw from meshmode.mesh import BTAG_ALL, BTAG_NONE -from grudge import sym + from pytools.obj_array import flat_obj_array +import grudge.op as op + # {{{ constant-velocity @@ -49,16 +55,18 @@ class WeakWaveOperator(HyperbolicOperator): :math:`c` is assumed to be constant across all space. """ - def __init__(self, c, ambient_dim, source_f=0, + def __init__(self, dcoll, c, source_f=None, flux_type="upwind", dirichlet_tag=BTAG_ALL, dirichlet_bc_f=0, neumann_tag=BTAG_NONE, radiation_tag=BTAG_NONE): - assert isinstance(ambient_dim, int) + if source_f is None: + source_f = lambda actx, dcoll, t: dcoll.zeros(actx) # noqa: E731 + + self.dcoll = dcoll self.c = c - self.ambient_dim = ambient_dim self.source_f = source_f if self.c > 0: @@ -74,10 +82,11 @@ class WeakWaveOperator(HyperbolicOperator): self.flux_type = flux_type - def flux(self, w): - u = w[0] - v = w[1:] - normal = sym.normal(w.dd, self.ambient_dim) + def flux(self, wtpair): + u = wtpair[0] + v = wtpair[1:] + actx = u.int.array_context + normal = thaw(actx, op.normal(self.dcoll, wtpair.dd)) central_flux_weak = -self.c*flat_obj_array( np.dot(v.avg, normal), @@ -92,65 +101,74 @@ class WeakWaveOperator(HyperbolicOperator): else: raise ValueError("invalid flux type '%s'" % self.flux_type) - def sym_operator(self): - d = self.ambient_dim - - w = sym.make_sym_array("w", d+1) + def operator(self, t, w): + dcoll = self.dcoll u = w[0] v = w[1:] + actx = u.array_context # 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)) + dir_u = op.project(dcoll, "vol", self.dirichlet_tag, u) + dir_v = op.project(dcoll, "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.var("dir_bc_u") + dir_g = self.dirichlet_bc_f 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") + neu_u = op.project(dcoll, "vol", self.neumann_tag, u) + neu_v = op.project(dcoll, "vol", self.neumann_tag, v) + neu_bc = flat_obj_array(neu_u, -neu_v) # radiation BCs ------------------------------------------------------- - rad_normal = sym.normal(self.radiation_tag, d) + rad_normal = thaw(actx, op.normal(dcoll, dd=self.radiation_tag)) - rad_u = sym.cse(sym.project("vol", self.radiation_tag)(u)) - rad_v = sym.cse(sym.project("vol", self.radiation_tag)(v)) + rad_u = op.project(dcoll, "vol", self.radiation_tag, u) + rad_v = op.project(dcoll, "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") + rad_bc = 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) + ) # entire operator ----------------------------------------------------- - def flux(pair): - return sym.project(pair.dd, "all_faces")(self.flux(pair)) - - result = sym.InverseMassOperator()( - flat_obj_array( - -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v), - -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)) - + flux(sym.bv_tpair(self.radiation_tag, w, rad_bc)) + def flux(tpair): + return op.project(dcoll, tpair.dd, "all_faces", self.flux(tpair)) - )) + def bv_tpair(tag, w, bc): + return TracePair(tag, interior=op.project(dcoll, "vol", tag, w), + exterior=bc) - result[0] += self.source_f + result = ( + op.inverse_mass( + dcoll, + flat_obj_array( + -self.c*op.weak_local_div(dcoll, v), + -self.c*op.weak_local_grad(dcoll, u) + ) + - op.face_mass( + dcoll, + flux(op.interior_trace_pair(dcoll, w)) + # communication of interface fluxes between + # parallel boundaries + + sum(flux(tpair) + for tpair in op.cross_rank_trace_pairs(dcoll, w)) + + flux(bv_tpair(self.dirichlet_tag, w, dir_bc)) + + flux(bv_tpair(self.neumann_tag, w, neu_bc)) + + flux(bv_tpair(self.radiation_tag, w, rad_bc)) + ) + ) + ) + + result[0] += self.source_f(actx, dcoll, t) return result @@ -188,21 +206,23 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): :math:`c` is assumed to be constant across all space. """ - def __init__(self, c, ambient_dim, source_f=0, + def __init__(self, dcoll, c, source_f=None, flux_type="upwind", dirichlet_tag=BTAG_ALL, dirichlet_bc_f=0, neumann_tag=BTAG_NONE, radiation_tag=BTAG_NONE): - assert isinstance(ambient_dim, int) + if source_f is None: + source_f = lambda actx, dcoll, t: dcoll.zeros(actx) # noqa: E731 + + self.dcoll = dcoll self.c = c - self.ambient_dim = ambient_dim self.source_f = source_f - self.sign = sym.If(sym.Comparison( - self.c, ">", 0), - np.int32(1), np.int32(-1)) + actx = dcoll._setup_actx + ones = dcoll.zeros(actx) + 1 + self.sign = actx.np.where(self.c > 0, ones, -ones) self.dirichlet_tag = dirichlet_tag self.neumann_tag = neumann_tag @@ -212,11 +232,12 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): self.flux_type = flux_type - def flux(self, w): - c = w[0] - u = w[1] - v = w[2:] - normal = sym.normal(w.dd, self.ambient_dim) + def flux(self, wtpair): + c = wtpair[0] + u = wtpair[1] + v = wtpair[2:] + actx = u.int.array_context + normal = thaw(actx, op.normal(self.dcoll, wtpair.dd)) flux_central_weak = -0.5 * flat_obj_array( np.dot(v.int*c.int + v.ext*c.ext, normal), @@ -234,71 +255,80 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): else: raise ValueError("invalid flux type '%s'" % self.flux_type) - def sym_operator(self): - d = self.ambient_dim - - w = sym.make_sym_array("w", d+1) + def operator(self, t, w): + dcoll = self.dcoll u = w[0] v = w[1:] flux_w = flat_obj_array(self.c, w) + actx = u.array_context # boundary conditions ------------------------------------------------- # dirichlet BCs ------------------------------------------------------- - dir_c = sym.cse(sym.project("vol", self.dirichlet_tag)(self.c)) - dir_u = sym.cse(sym.project("vol", self.dirichlet_tag)(u)) - dir_v = sym.cse(sym.project("vol", self.dirichlet_tag)(v)) + dir_c = op.project(dcoll, "vol", self.dirichlet_tag, self.c) + dir_u = op.project(dcoll, "vol", self.dirichlet_tag, u) + dir_v = op.project(dcoll, "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.var("dir_bc_u") + dir_g = self.dirichlet_bc_f dir_bc = flat_obj_array(dir_c, 2*dir_g - dir_u, dir_v) else: dir_bc = flat_obj_array(dir_c, -dir_u, dir_v) - dir_bc = sym.cse(dir_bc, "dir_bc") - # neumann BCs --------------------------------------------------------- - neu_c = sym.cse(sym.project("vol", self.neumann_tag)(self.c)) - 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_c, neu_u, -neu_v), "neu_bc") + neu_c = op.project(dcoll, "vol", self.neumann_tag, self.c) + neu_u = op.project(dcoll, "vol", self.neumann_tag, u) + neu_v = op.project(dcoll, "vol", self.neumann_tag, v) + neu_bc = flat_obj_array(neu_c, neu_u, -neu_v) # radiation BCs ------------------------------------------------------- - rad_normal = sym.normal(self.radiation_tag, d) + rad_normal = thaw(actx, op.normal(dcoll, dd=self.radiation_tag)) - rad_c = sym.cse(sym.project("vol", self.radiation_tag)(self.c)) - rad_u = sym.cse(sym.project("vol", self.radiation_tag)(u)) - rad_v = sym.cse(sym.project("vol", self.radiation_tag)(v)) + rad_c = op.project(dcoll, "vol", self.radiation_tag, self.c) + rad_u = op.project(dcoll, "vol", self.radiation_tag, u) + rad_v = op.project(dcoll, "vol", self.radiation_tag, v) + rad_sign = op.project(dcoll, "vol", self.radiation_tag, self.sign) - rad_bc = sym.cse(flat_obj_array(rad_c, - 0.5*(rad_u - sym.project("vol", self.radiation_tag)(self.sign) - * np.dot(rad_normal, rad_v)), - 0.5*rad_normal*(np.dot(rad_normal, rad_v) - - sym.project("vol", self.radiation_tag)(self.sign)*rad_u) - ), "rad_bc") + rad_bc = flat_obj_array( + rad_c, + 0.5*(rad_u - rad_sign * np.dot(rad_normal, rad_v)), + 0.5*rad_normal*(np.dot(rad_normal, rad_v) - rad_sign*rad_u) + ) # entire operator ----------------------------------------------------- - def flux(pair): - return sym.project(pair.dd, "all_faces")(self.flux(pair)) - - result = sym.InverseMassOperator()( - flat_obj_array( - -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v), - -self.c*(sym.stiffness_t(self.ambient_dim)*u) - ) + def flux(tpair): + return op.project(dcoll, tpair.dd, "all_faces", self.flux(tpair)) - - 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)) - + flux(sym.bv_tpair(self.radiation_tag, flux_w, rad_bc)) + def bv_tpair(tag, w, bc): + return TracePair(tag, interior=op.project(dcoll, "vol", tag, w), + exterior=bc) - )) - - result[0] += self.source_f + result = ( + op.inverse_mass( + dcoll, + flat_obj_array( + -self.c*op.weak_local_div(dcoll, v), + -self.c*op.weak_local_grad(dcoll, u) + ) + - op.face_mass( + dcoll, + flux(op.interior_trace_pair(dcoll, flux_w)) + # communication of interface fluxes between + # parallel boundaries + + sum(flux(tpair) + for tpair in op.cross_rank_trace_pairs(dcoll, w)) + + flux(bv_tpair(self.dirichlet_tag, flux_w, dir_bc)) + + flux(bv_tpair(self.neumann_tag, flux_w, neu_bc)) + + flux(bv_tpair(self.radiation_tag, flux_w, rad_bc)) + ) + ) + ) + + result[0] += self.source_f(actx, dcoll, t) return result @@ -310,7 +340,8 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): self.radiation_tag]) def max_eigenvalue(self, t, fields=None, discr=None): - return sym.NodalMax("vol")(sym.fabs(self.c)) + actx = self.dcoll._setup_actx + return op.nodal_maximum(actx.np.fabs(self.c)) # }}}