diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 4d3736b1012c00ce4f21cd650bd477e075fbb491..68eeafe0765ce499ec7d131de249bc9715a24ee5 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -41,12 +41,12 @@ def main(ctx_factory, dim=1, order=4, visualize=True): # grid spacing h = d / npoints # cfl? - dt_factor = 4 + dt_factor = 2.0 # final time final_time = 1.0 # velocity field - c = np.array([0.1] * dim) + c = np.array([0.5] * dim) norm_c = la.norm(c) # flux flux_type = "central" @@ -68,6 +68,9 @@ def main(ctx_factory, dim=1, order=4, visualize=True): from grudge import DGDiscretizationWithBoundaries discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + faces_discr = discr.discr_from_dd(sym.FACE_RESTR_INTERIOR) + # }}} # {{{ solve advection @@ -75,8 +78,9 @@ def main(ctx_factory, dim=1, order=4, visualize=True): def f(x): return sym.sin(3 * x) - def u_analytic(x): - t = sym.var("t", sym.DD_SCALAR) + def u_analytic(x, t=None): + if t is None: + t = sym.var("t", sym.DD_SCALAR) return f(-np.dot(c, x) / norm_c + t * norm_c) from grudge.models.advection import WeakAdvectionOperator @@ -87,8 +91,6 @@ def main(ctx_factory, dim=1, order=4, visualize=True): bound_op = bind(discr, op.sym_operator()) u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0) - return - def rhs(t, u): return bound_op(queue, t=t, u=u) @@ -98,35 +100,43 @@ def main(ctx_factory, dim=1, order=4, visualize=True): if dim == 1: import matplotlib.pyplot as pt pt.figure(figsize=(8, 8), dpi=300) + + volume_x = volume_discr.nodes().get(queue) else: from grudge.shortcuts import make_visualizer vis = make_visualizer(discr, vis_order=order) - step = 0 - norm = bind(discr, sym.norm(2, sym.var("u"))) - for event in dt_stepper.run(t_end=final_time): - if not isinstance(event, dt_stepper.StateComputed): - continue - - step += 1 - norm_u = norm(queue, u=event.state_component) - logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) + def plot_solution(evt): + if not visualize: + return if dim == 1: - x = discr.discr_from_dd(sym.DTAG_VOLUME_ALL).nodes().get(queue) u = event.state_component.get(queue) + u_ = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=evt.t).get(queue) - pt.plot(x, u) + pt.plot(volume_x[0], u, 'o') + pt.plot(volume_x[0], u_, "k.") pt.xlabel("$x$") pt.ylabel("$u$") - pt.title("t = {:.2f}".format(event.t)) + pt.title("t = {:.2f}".format(evt.t)) pt.savefig("fld-weak-%04d.png" % step) pt.clf() else: vis.write_vtk_file("fld-weak-%04d.vtu" % step, [ - ("u", event.state_component) + ("u", evt.state_component) ], overwrite=True) + step = 0 + norm = bind(discr, sym.norm(2, sym.var("u"))) + for event in dt_stepper.run(t_end=final_time): + if not isinstance(event, dt_stepper.StateComputed): + continue + + step += 1 + norm_u = norm(queue, u=event.state_component) + logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) + plot_solution(event) + if __name__ == "__main__": import argparse diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index 4e13caaf00354df11276d26648078857c8773b6a..51eb5202acaa1d1ffa6c4cbc8b293eef4a25e54e 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -1047,11 +1047,6 @@ class ToLoopyInstructionMapper(object): self.insn_count += 1 - for expr in insn.exprs: - print(expr) - print(self.dd_inference_mapper(expr)) - print() - from pytools import single_valued governing_dd = single_valued( self.dd_inference_mapper(expr) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index a1899ed21ef6431f7b57a81fcf143680433b7b53..0cbb6eda7c675d027677ba0f16fe5d11ecfe4a7c 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -536,6 +536,13 @@ class RefFaceMassOperator(ElementwiseLinearOperator): volgrp.order, face_vertices) + if afgrp.dim == 0 and volgrp.dim == 1: + # NOTE: This complements a choice in `parametrization_derivative` + # where we don't really get a sign for the 0-dim case, so we add + # signs here in the hope that it'll only be used for fluxes where + # this makes sense + matrix[0, 0, 0] = -matrix[0, 0, 0] + # np.set_printoptions(linewidth=200, precision=3) # matrix[np.abs(matrix) < 1e-13] = 0 # print(matrix) @@ -602,8 +609,7 @@ def norm(p, arg, dd=None): if p == 2: norm_squared = sym.NodalSum(dd_in=dd)( - sym.FunctionSymbol("fabs")( - arg * sym.MassOperator()(arg))) + sym.fabs(arg * sym.MassOperator()(arg))) if isinstance(norm_squared, np.ndarray): norm_squared = norm_squared.sum() @@ -611,7 +617,7 @@ def norm(p, arg, dd=None): return sym.FunctionSymbol("sqrt")(norm_squared) elif p == np.Inf: - result = sym.NodalMax(dd_in=dd)(sym.FunctionSymbol("fabs")(arg)) + result = sym.NodalMax(dd_in=dd)(sym.fabs(arg)) from pymbolic.primitives import Max if isinstance(result, np.ndarray): diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 0b9fac20332e505b038ee2ed9bda94485048d8d2..80dd8413f01dbb3a6961a4d740e13c14d31de11d 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -364,6 +364,7 @@ class FunctionSymbol(ExpressionBase, pymbolic.primitives.Variable): mapper_method = "map_function_symbol" +fabs = FunctionSymbol("fabs") sqrt = FunctionSymbol("sqrt") exp = FunctionSymbol("exp") sin = FunctionSymbol("sin") @@ -506,7 +507,8 @@ def parametrization_derivative(ambient_dim, dim=None, dd=None): dim = ambient_dim if dim == 0: - return MultiVector(np.array([1.])) + inner_dd = (DD_VOLUME if dd is None else dd).with_qtag(QTAG_NONE) + return MultiVector(np.array([Ones(dd=inner_dd)])) from pytools import product return product( @@ -612,9 +614,10 @@ def mv_normal(dd, ambient_dim, dim=None): pseudoscalar(ambient_dim, dim, dd=dd) / # noqa: W504 area_element(ambient_dim, dim, dd=dd)) + return cse( # Dorst Section 3.7.2 - pder << pder.I.inv(), + pder if dim == 0 else pder << pder.I.inv(), "normal", cse_scope.DISCRETIZATION)