diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py index aa9e87ffa0a7885adf9cc8d67d82d1be15a44e93..5737dcd57ed41482ecaf96749e7923ae235d0b62 100644 --- a/examples/wave/wave-eager.py +++ b/examples/wave/wave-eager.py @@ -29,13 +29,13 @@ import pyopencl as cl import pyopencl.array as cla # noqa import pyopencl.clmath as clmath from pytools.obj_array import ( - join_fields, make_obj_array, with_object_array_or_scalar) from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.eager import EagerDGDiscretization, with_queue from grudge.symbolic.primitives import TracePair from grudge.shortcuts import make_visualizer +from pytools.obj_array import flat_obj_array, make_obj_array def interior_trace_pair(discr, vec): i = discr.interp("vol", "int_faces", vec) @@ -57,14 +57,14 @@ def wave_flux(discr, c, w_tpair): # workaround for object array behavior return make_obj_array([ni*scalar for ni in normal]) - flux_weak = join_fields( + flux_weak = flat_obj_array( np.dot(v.avg, normal), normal_times(u.avg), ) # upwind v_jump = np.dot(normal, v.int-v.ext) - flux_weak -= join_fields( + flux_weak -= flat_obj_array( 0.5*(u.int-u.ext), 0.5*normal_times(v_jump), ) @@ -78,12 +78,12 @@ def wave_operator(discr, c, w): dir_u = discr.interp("vol", BTAG_ALL, u) dir_v = discr.interp("vol", BTAG_ALL, v) - dir_bval = join_fields(dir_u, dir_v) - dir_bc = join_fields(-dir_u, dir_v) + dir_bval = flat_obj_array(dir_u, dir_v) + dir_bc = flat_obj_array(-dir_u, dir_v) return ( discr.inverse_mass( - join_fields( + flat_obj_array( c*discr.weak_div(v), c*discr.weak_grad(u) ) @@ -112,7 +112,7 @@ def bump(discr, queue, t=0): source_omega = 3 nodes = discr.nodes().with_queue(queue) - center_dist = join_fields([ + center_dist = flat_obj_array([ nodes[i] - source_center[i] for i in range(discr.dim) ]) @@ -151,7 +151,7 @@ def main(): discr = EagerDGDiscretization(cl_ctx, mesh, order=order) - fields = join_fields( + fields = flat_obj_array( bump(discr, queue), [discr.zeros(queue) for i in range(discr.dim)] ) diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py index 08705655cf874bf7701bb1c3f998c0d6ea7a6e06..9a8e70bf8b2f2ee76bca2c586089f7de30958f5a 100644 --- a/examples/wave/wave-min-mpi.py +++ b/examples/wave/wave-min-mpi.py @@ -91,8 +91,8 @@ def main(write_output=True, order=4): flux_type="upwind") queue = cl.CommandQueue(discr.cl_context) - from pytools.obj_array import join_fields - fields = join_fields(discr.zeros(queue), + from pytools.obj_array import flat_obj_array + fields = flat_obj_array(discr.zeros(queue), [discr.zeros(queue) for i in range(discr.dim)]) # FIXME diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py index d793e3a09716b883d91c90f1fddfb2f96bacf1a2..665ce71f254a61f105b6220d0eff07d91a6b2429 100644 --- a/examples/wave/wave-min.py +++ b/examples/wave/wave-min.py @@ -73,8 +73,8 @@ def main(write_output=True, order=4): flux_type="upwind") queue = cl.CommandQueue(discr.cl_context) - from pytools.obj_array import join_fields - fields = join_fields(discr.zeros(queue), + from pytools.obj_array import flat_obj_array + fields = flat_obj_array(discr.zeros(queue), [discr.zeros(queue) for i in range(discr.dim)]) # FIXME diff --git a/grudge/models/em.py b/grudge/models/em.py index 225b3b4a411beae96359df0548ff5366f00ea675..c3a651c4572ec7caf64ca8068d644826ebdc6cce 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -34,7 +34,7 @@ from pytools import memoize_method from grudge.models import HyperbolicOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE from grudge import sym -from pytools.obj_array import join_fields, make_obj_array +from pytools.obj_array import flat_obj_array, make_obj_array # TODO: Check PML @@ -133,7 +133,7 @@ class MaxwellOperator(HyperbolicOperator): # if self.fixed_material: # max_c = (self.epsilon*self.mu)**(-0.5) - return join_fields( + return flat_obj_array( # flux e, 1/2*( -self.space_cross_h(normal, h.int-h.ext) @@ -148,7 +148,7 @@ class MaxwellOperator(HyperbolicOperator): )) elif isinstance(self.flux_type, (int, float)): # see doc/maxima/maxwell.mac - return join_fields( + return flat_obj_array( # flux e, ( -1/(Z_int+Z_ext)*self.space_cross_h(normal, @@ -182,7 +182,7 @@ class MaxwellOperator(HyperbolicOperator): return self.space_cross_h(nabla, field) # in conservation form: u_t + A u_x = 0 - return join_fields( + return flat_obj_array( (self.current - h_curl(h)), e_curl(e) ) @@ -194,7 +194,7 @@ class MaxwellOperator(HyperbolicOperator): pec_e = sym.cse(sym.interp("vol", self.pec_tag)(e)) pec_h = sym.cse(sym.interp("vol", self.pec_tag)(h)) - return join_fields(-pec_e, pec_h) + return flat_obj_array(-pec_e, pec_h) def pmc_bc(self, w): "Construct part of the flux operator template for PMC boundary conditions" @@ -203,7 +203,7 @@ class MaxwellOperator(HyperbolicOperator): pmc_e = sym.cse(sym.interp("vol", self.pmc_tag)(e)) pmc_h = sym.cse(sym.interp("vol", self.pmc_tag)(h)) - return join_fields(pmc_e, -pmc_h) + return flat_obj_array(pmc_e, -pmc_h) def absorbing_bc(self, w): """Construct part of the flux operator template for 1st order @@ -224,7 +224,7 @@ class MaxwellOperator(HyperbolicOperator): absorb_e = sym.cse(sym.interp("vol", self.absorb_tag)(e)) absorb_h = sym.cse(sym.interp("vol", self.absorb_tag)(h)) - bc = join_fields( + bc = flat_obj_array( absorb_e + 1/2*(self.space_cross_h(absorb_normal, self.space_cross_e( absorb_normal, absorb_e)) - absorb_Z*self.space_cross_h(absorb_normal, absorb_h)), @@ -437,7 +437,7 @@ def get_rectangular_cavity_mode(E_0, mode_indices): # noqa: N803 if dims == 2: tfac = sym.ScalarVariable("t") * omega - result = join_fields( + result = flat_obj_array( 0, 0, sym.sin(kx * x) * sym.sin(ky * y) * sym.cos(tfac), # ez @@ -449,7 +449,7 @@ def get_rectangular_cavity_mode(E_0, mode_indices): # noqa: N803 tdep = sym.exp(-1j * omega * sym.ScalarVariable("t")) gamma_squared = ky**2 + kx**2 - result = join_fields( + result = flat_obj_array( -kx * kz * E_0*cx*sy*sz*tdep / gamma_squared, # ex -ky * kz * E_0*sx*cy*sz*tdep / gamma_squared, # ey E_0 * sx*sy*cz*tdep, # ez diff --git a/grudge/models/wave.py b/grudge/models/wave.py index 6defafc878f856fe83449b12dc055cfdf0b1f432..90f248d6556f6615733a361158ba573f4a6d9575 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -30,7 +30,7 @@ import numpy as np from grudge.models import HyperbolicOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE from grudge import sym -from pytools.obj_array import join_fields +from pytools.obj_array import flat_obj_array # {{{ constant-velocity @@ -83,7 +83,7 @@ class StrongWaveOperator(HyperbolicOperator): v = w[1:] normal = sym.normal(w.dd, self.ambient_dim) - flux_weak = join_fields( + flux_weak = flat_obj_array( np.dot(v.avg, normal), u.avg * normal) @@ -91,13 +91,13 @@ class StrongWaveOperator(HyperbolicOperator): pass elif self.flux_type == "upwind": # see doc/notes/grudge-notes.tm - flux_weak -= self.sign*join_fields( + flux_weak -= self.sign*flat_obj_array( 0.5*(u.int-u.ext), 0.5*(normal * np.dot(normal, v.int-v.ext))) else: raise ValueError("invalid flux type '%s'" % self.flux_type) - flux_strong = join_fields( + flux_strong = flat_obj_array( np.dot(v.int, normal), u.int * normal) - flux_weak @@ -122,16 +122,16 @@ class StrongWaveOperator(HyperbolicOperator): "are still having issues.") dir_g = sym.Field("dir_bc_u") - dir_bc = join_fields(2*dir_g - dir_u, dir_v) + dir_bc = flat_obj_array(2*dir_g - dir_u, dir_v) else: - dir_bc = join_fields(-dir_u, dir_v) + dir_bc = flat_obj_array(-dir_u, dir_v) dir_bc = sym.cse(dir_bc, "dir_bc") # neumann BCs --------------------------------------------------------- neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u)) neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v)) - neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc") + neu_bc = sym.cse(flat_obj_array(neu_u, -neu_v), "neu_bc") # radiation BCs ------------------------------------------------------- rad_normal = sym.normal(self.radiation_tag, d) @@ -139,7 +139,7 @@ class StrongWaveOperator(HyperbolicOperator): rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u)) rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v)) - rad_bc = sym.cse(join_fields( + 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") @@ -152,7 +152,7 @@ class StrongWaveOperator(HyperbolicOperator): self.flux(pair)) result = ( - - join_fields( + - flat_obj_array( -self.c*np.dot(nabla, v), -self.c*(nabla*u) ) @@ -230,14 +230,14 @@ class WeakWaveOperator(HyperbolicOperator): v = w[1:] normal = sym.normal(w.dd, self.ambient_dim) - flux_weak = join_fields( + flux_weak = flat_obj_array( np.dot(v.avg, normal), u.avg * normal) if self.flux_type == "central": return -self.c*flux_weak elif self.flux_type == "upwind": - return -self.c*(flux_weak + self.sign*join_fields( + return -self.c*(flux_weak + self.sign*flat_obj_array( 0.5*(u.int-u.ext), 0.5*(normal * np.dot(normal, v.int-v.ext)))) else: @@ -262,16 +262,16 @@ class WeakWaveOperator(HyperbolicOperator): "are still having issues.") dir_g = sym.Field("dir_bc_u") - dir_bc = join_fields(2*dir_g - dir_u, dir_v) + dir_bc = flat_obj_array(2*dir_g - dir_u, dir_v) else: - dir_bc = join_fields(-dir_u, dir_v) + dir_bc = flat_obj_array(-dir_u, dir_v) dir_bc = sym.cse(dir_bc, "dir_bc") # neumann BCs --------------------------------------------------------- neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u)) neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v)) - neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc") + neu_bc = sym.cse(flat_obj_array(neu_u, -neu_v), "neu_bc") # radiation BCs ------------------------------------------------------- rad_normal = sym.normal(self.radiation_tag, d) @@ -279,7 +279,7 @@ class WeakWaveOperator(HyperbolicOperator): rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u)) rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v)) - rad_bc = sym.cse(join_fields( + 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") @@ -289,7 +289,7 @@ class WeakWaveOperator(HyperbolicOperator): return sym.interp(pair.dd, "all_faces")(self.flux(pair)) result = sym.InverseMassOperator()( - join_fields( + flat_obj_array( -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v), -self.c*(sym.stiffness_t(self.ambient_dim)*u) ) @@ -371,14 +371,16 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): normal = sym.normal(w.dd, self.ambient_dim) if self.flux_type == "central": - return -0.5 * join_fields( + return -0.5 * flat_obj_array( np.dot(v.int*c.int + v.ext*c.ext, normal), (u.int * c.int + u.ext*c.ext) * normal) elif self.flux_type == "upwind": - return -0.5 * join_fields(np.dot(normal, c.ext * v.ext + c.int * v.int) + return -0.5 * flat_obj_array( + np.dot(normal, c.ext * v.ext + c.int * v.int) + c.ext*u.ext - c.int * u.int, - normal * (np.dot(normal, c.ext * v.ext - c.int * v.int) + + normal * (np.dot(normal, c.ext * v.ext - c.int * v.int) + c.ext*u.ext + c.int * u.int)) else: @@ -390,7 +392,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): w = sym.make_sym_array("w", d+1) u = w[0] v = w[1:] - flux_w = join_fields(self.c, w) + flux_w = flat_obj_array(self.c, w) # boundary conditions ------------------------------------------------- @@ -405,9 +407,9 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): "are still having issues.") dir_g = sym.Field("dir_bc_u") - dir_bc = join_fields(dir_c, 2*dir_g - dir_u, dir_v) + dir_bc = flat_obj_array(dir_c, 2*dir_g - dir_u, dir_v) else: - dir_bc = join_fields(dir_c, -dir_u, dir_v) + dir_bc = flat_obj_array(dir_c, -dir_u, dir_v) dir_bc = sym.cse(dir_bc, "dir_bc") @@ -415,7 +417,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): neu_c = sym.cse(sym.interp("vol", self.neumann_tag)(self.c)) neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u)) neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v)) - neu_bc = sym.cse(join_fields(neu_c, neu_u, -neu_v), "neu_bc") + neu_bc = sym.cse(flat_obj_array(neu_c, neu_u, -neu_v), "neu_bc") # radiation BCs ------------------------------------------------------- rad_normal = sym.normal(self.radiation_tag, d) @@ -424,7 +426,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u)) rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v)) - rad_bc = sym.cse(join_fields(rad_c, + rad_bc = sym.cse(flat_obj_array(rad_c, 0.5*(rad_u - sym.interp("vol", self.radiation_tag)(self.sign) * np.dot(rad_normal, rad_v)), 0.5*rad_normal*(np.dot(rad_normal, rad_v) @@ -436,7 +438,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): return sym.interp(pair.dd, "all_faces")(self.flux(pair)) result = sym.InverseMassOperator()( - join_fields( + flat_obj_array( -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v), -self.c*(sym.stiffness_t(self.ambient_dim)*u) ) diff --git a/grudge/tools.py b/grudge/tools.py index 8c5f57b6f0704bc3a5d872cc77ddac5fe366fcf8..67107e48048b4e78727a97e5d828f9a53cffdcad 100644 --- a/grudge/tools.py +++ b/grudge/tools.py @@ -92,11 +92,11 @@ class SubsettableCrossProduct: used in place of the product *sign*xj*yk*. Defaults to just this product if not given. """ - from pytools.obj_array import join_fields + from pytools.obj_array import flat_obj_array if three_mult is None: - return join_fields(*[f(x, y) for f in self.functions]) + return flat_obj_array(*[f(x, y) for f in self.functions]) else: - return join_fields( + return flat_obj_array( *[sum(three_mult(lc, x[j], y[k]) for lc, j, k in lcjk) for lcjk in self.component_lcjk]) diff --git a/test/test_grudge.py b/test/test_grudge.py index fc20811f9e126dea877fa61d4111b3b336db606f..26f274b764e727c2eae5f952c5f9b3a7c2727a46 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -29,8 +29,8 @@ import numpy.linalg as la import pyopencl as cl import pyopencl.array import pyopencl.clmath +from pytools.obj_array import flat_obj_array, make_obj_array -from pytools.obj_array import join_fields, make_obj_array import pytest # noqa @@ -221,7 +221,7 @@ def test_2d_gauss_theorem(ctx_factory): discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=2) def f(x): - return join_fields( + return flat_obj_array( sym.sin(3*x[0])+sym.cos(3*x[1]), sym.sin(2*x[0])+sym.cos(x[1])) @@ -458,7 +458,7 @@ def test_improvement_quadrature(ctx_factory, order): dims = 2 sym_nds = sym.nodes(dims) - advec_v = join_fields(-1*sym_nds[1], sym_nds[0]) + advec_v = flat_obj_array(-1*sym_nds[1], sym_nds[0]) flux = "upwind" op = VariableCoefficientAdvectionOperator(advec_v, 0, flux_type=flux) diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py index f9b12d3b4ec8922204cb7201c715adbb51f037b4..1047a71e0583e3de6efae7921fea86ad52440826 100644 --- a/test/test_mpi_communication.py +++ b/test/test_mpi_communication.py @@ -148,8 +148,8 @@ def mpi_communication_entrypoint(): radiation_tag=BTAG_ALL, flux_type="upwind") - from pytools.obj_array import join_fields - fields = join_fields(vol_discr.zeros(queue), + from pytools.obj_array import flat_obj_array + fields = flat_obj_array(vol_discr.zeros(queue), [vol_discr.zeros(queue) for i in range(vol_discr.dim)]) # FIXME