diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index c8791efd250fea80ecb8dbd94a1ba9b85f77590d..7d4c91f53c0aae94bd72929d9c9950e90acad5a1 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -35,6 +35,53 @@ import logging logger = logging.getLogger(__name__) +# {{{ plotting (keep in sync with `weak.py`) + +class Plotter: + def __init__(self, queue, discr, order, visualize=True): + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + + self.queue = queue + self.x = volume_discr.nodes().get(self.queue) + + self.visualize = visualize + if not self.visualize: + return + + if self.dim == 1: + import matplotlib.pyplot as pt + self.fig = pt.figure(figsize=(8, 8), dpi=300) + else: + from grudge.shortcuts import make_visualizer + self.vis = make_visualizer(discr, vis_order=order) + + @property + def dim(self): + return len(self.x) + + def __call__(self, evt, basename): + if not self.visualize: + return + + if self.dim == 1: + u = evt.state_component.get(self.queue) + + ax = self.fig.gca() + ax.plot(self.x[0], u, '-') + ax.plot(self.x[0], u, 'k.') + ax.set_xlabel("$x$") + ax.set_ylabel("$u$") + ax.set_title("t = {:.2f}".format(evt.t)) + self.fig.savefig("%s.png" % basename) + self.fig.clf() + else: + self.vis.write_vtk_file("%s.vtu" % basename, [ + ("u", evt.state_component) + ], overwrite=True) + +# }}} + + def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) @@ -106,6 +153,11 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): def f_gaussian(x): return sym.exp(-dist_squared / source_width**2) + def f_step(x): + return sym.If(sym.Comparison( + dist_squared, "<", (4*source_width) ** 2), + 1, 0) + def u_bc(x): return 0.0 @@ -128,8 +180,6 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): from grudge.shortcuts import set_up_rk4 dt_stepper = set_up_rk4("u", dt, u, rhs) - - from weak import Plotter plot = Plotter(queue, discr, order, visualize=visualize) step = 0 @@ -138,12 +188,12 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): if not isinstance(event, dt_stepper.StateComputed): continue - if step % 5 == 0: + if step % 10 == 0: norm_u = norm(queue, u=event.state_component) + plot(event, "fld-var-velocity-%04d" % step) step += 1 logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) - plot(event, "fld-var-velocity-%04d" % step) # }}} diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 770e9717fa55e53ad9ba75f58f5d6444b4effa74..e780f7f8b8491cd772aecac39a61072a8e9fa86e 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -30,6 +30,8 @@ import logging logger = logging.getLogger(__name__) +# {{{ plotting (keep in sync with `var-velocity.py`) + class Plotter: def __init__(self, queue, discr, order, visualize=True): volume_discr = discr.discr_from_dd(sym.DD_VOLUME) @@ -60,7 +62,8 @@ class Plotter: u = evt.state_component.get(self.queue) ax = self.fig.gca() - ax.plot(self.x[0], u, 'o') + ax.plot(self.x[0], u, '-') + ax.plot(self.x[0], u, 'k.') ax.set_xlabel("$x$") ax.set_ylabel("$u$") ax.set_title("t = {:.2f}".format(evt.t)) @@ -71,6 +74,8 @@ class Plotter: ("u", evt.state_component) ], overwrite=True) +# }}} + def main(ctx_factory, dim=2, order=4, visualize=True): cl_ctx = ctx_factory() diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 8c2f7fd0fecad5ef4c4dcae4e2d708d317edc3a7..1950d3f7248705af923d1e9148294d7e402a3028 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -307,8 +307,11 @@ class RefStiffnessTOperator(RefDiffOperatorBase): grad_vand = vandermonde(out_elem_grp.grad_basis(), in_elem_grp.unit_nodes) vand_inv_t = np.linalg.inv(vand).T - weights = in_elem_grp.weights + if not isinstance(grad_vand, tuple): + # NOTE: special case for 1d + grad_vand = (grad_vand,) + weights = in_elem_grp.weights return np.einsum('c,bz,acz->abc', weights, vand_inv_t, grad_vand) diff --git a/test/test_grudge.py b/test/test_grudge.py index 12b7b116e010006c957d34a210d3455e223c25c4..f068d34d8ea8b3de15053d20187ff13c0ed66223 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -460,7 +460,7 @@ def test_improvement_quadrature(ctx_factory, order): advec_v = join_fields(-1*sym_nds[1], sym_nds[0]) flux = "upwind" - op = VariableCoefficientAdvectionOperator(2, advec_v, 0, flux_type=flux) + op = VariableCoefficientAdvectionOperator(advec_v, 0, flux_type=flux) def gaussian_mode(): source_width = 0.1