From 625db070a9bbb957bbc4e69965c82ff714a6110c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 9 May 2020 20:39:49 -0500 Subject: [PATCH 1/3] examples: clean up variable velocity example --- examples/advection/var-velocity.py | 165 +++++++++++++++++------------ examples/advection/weak.py | 137 ++++++++++++++---------- grudge/models/advection.py | 113 +++++++++----------- 3 files changed, 224 insertions(+), 191 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index bc720731..c8791efd 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -22,114 +22,141 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ - import numpy as np -import pyopencl as cl # noqa -import pyopencl.array # noqa -import pyopencl.clmath # noqa -import pytest # noqa +import pyopencl as cl +import pyopencl.array +import pyopencl.clmath -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from grudge import bind, sym +from pytools.obj_array import join_fields import logging logger = logging.getLogger(__name__) -from grudge import sym, bind, DGDiscretizationWithBoundaries -from pytools.obj_array import join_fields +def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) -FINAL_TIME = 5 + # {{{ parameters + # domain [0, d]^dim + d = 1.0 + # number of points in each dimension + npoints = 25 + # grid spacing + h = d / npoints -def main(write_output=True, order=4): - cl_ctx = cl.create_some_context() - queue = cl.CommandQueue(cl_ctx) + # cfl + dt_factor = 1.0 + # finale time + final_time = 0.5 + # time steps + dt = dt_factor * h/order**2 + nsteps = int(final_time // dt) + 1 + dt = final_time/nsteps + 1.0e-15 - dim = 2 + # flux + flux_type = "upwind" + # velocity field + sym_x = sym.nodes(dim) + if dim == 1: + c = sym_x + else: + # solid body rotation + c = join_fields( + np.pi * (d/2 - sym_x[1]), + np.pi * (sym_x[0] - d/2), + 0)[:dim] - resolution = 15 - from meshmode.mesh.generation import generate_regular_rect_mesh - mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), - n=(resolution, resolution), order=order) + # }}} - dt_factor = 5 - h = 1/resolution + # {{{ discretization - sym_x = sym.nodes(2) + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(0,)*dim, b=(d,)*dim, + n=(npoints,)*dim, + order=order) + + from meshmode.discretization.poly_element import \ + QuadratureSimplexGroupFactory + + if product_tag: + quad_tag_to_group_factory = { + product_tag: QuadratureSimplexGroupFactory(order=4*order) + } + else: + quad_tag_to_group_factory = {} + + from grudge import DGDiscretizationWithBoundaries + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, + quad_tag_to_group_factory=quad_tag_to_group_factory) - advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2 + # }}} - flux_type = "upwind" + # {{{ symbolic operators - source_center = np.array([0.1, 0.1]) + # gaussian parameters + source_center = np.array([0.5, 0.75, 0.0])[:dim] source_width = 0.05 - - sym_x = sym.nodes(2) - sym_source_center_dist = sym_x - source_center + dist_squared = np.dot(sym_x - source_center, sym_x - source_center) def f_gaussian(x): - return sym.exp( - -np.dot(sym_source_center_dist, sym_source_center_dist) - / source_width**2) - - def f_step(x): - return sym.If( - sym.Comparison( - np.dot(sym_source_center_dist, sym_source_center_dist), - "<", - (4*source_width)**2), - 1, 0) + return sym.exp(-dist_squared / source_width**2) - def u_analytic(x): - return 0 + def u_bc(x): + return 0.0 from grudge.models.advection import VariableCoefficientAdvectionOperator - from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory # noqa - - discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, - quad_tag_to_group_factory={ - #"product": None, - "product": QuadratureSimplexGroupFactory(order=4*order) - }) - - op = VariableCoefficientAdvectionOperator(2, advec_v, - u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product", - flux_type=flux_type) - - bound_op = bind( - discr, op.sym_operator(), - #debug_flags=["dump_sym_operator_stages"] - ) + op = VariableCoefficientAdvectionOperator( + c, + u_bc(sym.nodes(dim, sym.BTAG_ALL)), + quad_tag=product_tag, + flux_type=flux_type) + bound_op = bind(discr, op.sym_operator()) u = bind(discr, f_gaussian(sym.nodes(dim)))(queue, t=0) def rhs(t, u): return bound_op(queue, t=t, u=u) - dt = dt_factor * h/order**2 - nsteps = (FINAL_TIME // dt) + 1 - dt = FINAL_TIME/nsteps + 1e-15 + # }}} + + # {{{ time stepping from grudge.shortcuts import set_up_rk4 dt_stepper = set_up_rk4("u", dt, u, rhs) - from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr, vis_order=2*order) + from weak import Plotter + plot = Plotter(queue, discr, order, visualize=visualize) 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 - for event in dt_stepper.run(t_end=FINAL_TIME): - if isinstance(event, dt_stepper.StateComputed): + if step % 5 == 0: + norm_u = norm(queue, u=event.state_component) - step += 1 - if step % 30 == 0: - print(step) + step += 1 + logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) + plot(event, "fld-var-velocity-%04d" % step) - vis.write_vtk_file("fld-var-velocity-%04d.vtu" % step, - [("u", event.state_component)]) + # }}} if __name__ == "__main__": - main() + import argparse + + parser = argparse.ArgumentParser() + parser.add_argument("--dim", default=2, type=int) + parser.add_argument("--qtag", default="product") + args = parser.parse_args() + + logging.basicConfig(level=logging.INFO) + main(cl.create_some_context, + dim=args.dim, + product_tag=args.qtag) diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 6e30094d..770e9717 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -1,17 +1,21 @@ -# Copyright (C) 2007 Andreas Kloeckner -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2007 Andreas Kloeckner" + +__license__ = """ +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . +""" import numpy as np import numpy.linalg as la @@ -23,9 +27,49 @@ import pyopencl.clmath from grudge import bind, sym import logging - logger = logging.getLogger(__name__) -logging.basicConfig(level=logging.INFO) + + +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, 'o') + 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, visualize=True): @@ -34,16 +78,21 @@ def main(ctx_factory, dim=2, order=4, visualize=True): # {{{ parameters - # domain side [-d/2, d/2]^dim + # domain [-d/2, d/2]^dim d = 1.0 # number of points in each dimension npoints = 20 # grid spacing h = d / npoints - # cfl? + + # cfl dt_factor = 2.0 # final time final_time = 1.0 + # compute number of steps + dt = dt_factor * h/order**2 + nsteps = int(final_time // dt) + 1 + dt = final_time/nsteps + 1.0e-15 # velocity field c = np.array([0.5] * dim) @@ -51,11 +100,6 @@ def main(ctx_factory, dim=2, order=4, visualize=True): # flux flux_type = "central" - # compute number of steps - dt = dt_factor * h / order**2 - nsteps = int(final_time // dt) + 1 - dt = final_time/nsteps + 1.0e-15 - # }}} # {{{ discretization @@ -67,18 +111,16 @@ def main(ctx_factory, dim=2, order=4, visualize=True): from grudge import DGDiscretizationWithBoundaries discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) - volume_discr = discr.discr_from_dd(sym.DD_VOLUME) # }}} - # {{{ solve advection + # {{{ symbolic operators def f(x): return sym.sin(3 * x) - def u_analytic(x, t=None): - if t is None: - t = sym.var("t", sym.DD_SCALAR) + def u_analytic(x): + t = sym.var("t", sym.DD_SCALAR) return f(-np.dot(c, x) / norm_c + t * norm_c) from grudge.models.advection import WeakAdvectionOperator @@ -92,37 +134,13 @@ def main(ctx_factory, dim=2, order=4, visualize=True): def rhs(t, u): return bound_op(queue, t=t, u=u) - from grudge.shortcuts import set_up_rk4 - dt_stepper = set_up_rk4("u", dt, u, rhs) - - 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) + # }}} - def plot_solution(evt): - if not visualize: - return + # {{{ time stepping - if dim == 1: - u = event.state_component.get(queue) - u_ = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=evt.t).get(queue) - - 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(evt.t)) - pt.savefig("fld-weak-%04d.png" % step) - pt.clf() - else: - vis.write_vtk_file("fld-weak-%04d.vtu" % step, [ - ("u", evt.state_component) - ], overwrite=True) + from grudge.shortcuts import set_up_rk4 + dt_stepper = set_up_rk4("u", dt, u, rhs) + plot = Plotter(queue, discr, order, visualize=visualize) step = 0 norm = bind(discr, sym.norm(2, sym.var("u"))) @@ -133,7 +151,9 @@ def main(ctx_factory, dim=2, order=4, visualize=True): 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) + plot(event, "fld-weak-%04d" % step) + + # }}} if __name__ == "__main__": @@ -143,5 +163,6 @@ if __name__ == "__main__": parser.add_argument("--dim", default=2, type=int) args = parser.parse_args() + logging.basicConfig(level=logging.INFO) main(cl.create_some_context, dim=args.dim) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index c3b142e4..36ce5c15 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -33,6 +33,31 @@ from grudge.models import HyperbolicOperator from grudge import sym +# {{{ fluxes + +def advection_weak_flux(flux_type, u, velocity): + normal = sym.normal(u.dd, len(velocity)) + v_dot_n = sym.cse(velocity.dot(normal), "v_dot_normal") + + flux_type = flux_type.lower() + if flux_type == "central": + return u.avg * v_dot_n + elif flux_type == "lf": + norm_v = sym.sqrt((velocity**2).sum()) + return u.avg * v_dot_n + 0.5 * norm_v * (u.int - u.ext) + elif flux_type == "upwind": + u_upwind = sym.If( + sym.Comparison(v_dot_n, ">", 0), + u.int, # outflow + u.ext # inflow + ) + return u_upwind * v_dot_n + else: + raise ValueError("flux `{}` is not implemented".format(flux_type)) + +# }}} + + # {{{ constant-coefficient advection class AdvectionOperatorBase(HyperbolicOperator): @@ -48,25 +73,11 @@ class AdvectionOperatorBase(HyperbolicOperator): self.inflow_u = inflow_u self.flux_type = flux_type - def weak_flux(self, u): - normal = sym.normal(u. dd, self.ambient_dim) + if flux_type not in self.flux_types: + raise ValueError("unknown flux type: {}".format(flux_type)) - v_dot_normal = sym.cse(self.v.dot(normal), "v_dot_normal") - norm_v = sym.sqrt((self.v**2).sum()) - - if self.flux_type == "central": - return u.avg*v_dot_normal - elif self.flux_type == "lf": - return u.avg*v_dot_normal + 0.5*norm_v*(u.int - u.ext) - elif self.flux_type == "upwind": - return ( - v_dot_normal * sym.If( - sym.Comparison(v_dot_normal, ">", 0), - u.int, # outflow - u.ext, # inflow - )) - else: - raise ValueError("invalid flux type") + def weak_flux(self, u): + return advection_weak_flux(self.flux_type, u, self.v) def max_eigenvalue(self, t=None, fields=None, discr=None): return la.norm(self.v) @@ -106,14 +117,13 @@ class WeakAdvectionOperator(AdvectionOperatorBase): def sym_operator(self): u = sym.var("u") - # boundary conditions ------------------------------------------------- - bc_in = self.inflow_u - # bc_out = sym.interp("vol", self.outflow_tag)(u) - def flux(pair): return sym.interp(pair.dd, "all_faces")( self.flux(pair)) + bc_in = self.inflow_u + # bc_out = sym.interp("vol", self.outflow_tag)(u) + return sym.InverseMassOperator()( np.dot( self.v, sym.stiffness_t(self.ambient_dim)*u) @@ -131,65 +141,40 @@ class WeakAdvectionOperator(AdvectionOperatorBase): # {{{ variable-coefficient advection -class VariableCoefficientAdvectionOperator(HyperbolicOperator): - def __init__(self, dim, v, inflow_u, flux_type="central", quad_tag="product"): - self.ambient_dim = dim - self.v = v - self.inflow_u = inflow_u - self.flux_type = flux_type +class VariableCoefficientAdvectionOperator(AdvectionOperatorBase): + def __init__(self, v, inflow_u, flux_type="central", quad_tag="product"): + super(VariableCoefficientAdvectionOperator, self).__init__( + v, inflow_u, flux_type=flux_type) self.quad_tag = quad_tag def flux(self, u): - - normal = sym.normal(u.dd, self.ambient_dim) - - surf_v = sym.interp("vol", u.dd)(self.v) - - v_dot_normal = sym.cse(np.dot(surf_v, normal), "v_dot_normal") - norm_v = sym.sqrt(np.sum(self.v**2)) - - if self.flux_type == "central": - return u.avg*v_dot_normal - - elif self.flux_type == "lf": - return u.avg*v_dot_normal + 0.5*norm_v*(u.int - u.ext) - elif self.flux_type == "upwind": - return ( - v_dot_normal * sym.If( - sym.Comparison(v_dot_normal, ">", 0), - u.int, # outflow - u.ext, # inflow - )) - else: - raise ValueError("invalid flux type") + surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v) + return advection_weak_flux(self.flux_type, u, surf_v) def sym_operator(self): u = sym.var("u") - # boundary conditions ------------------------------------------------- - bc_in = self.inflow_u - - all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) - boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) - def flux(pair): - return sym.interp(pair.dd, all_faces_dd)( - self.flux(pair)) + return sym.interp(pair.dd, face_dd)(self.flux(pair)) - quad_dd = sym.DOFDesc("vol", self.quad_tag) + face_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) + boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) + quad_dd = sym.DOFDesc(sym.DTAG_VOLUME_ALL, self.quad_tag) - to_quad = sym.interp("vol", quad_dd) + to_quad = sym.interp(sym.DD_VOLUME, quad_dd) + stiff_t_op = sym.stiffness_t(self.ambient_dim, + dd_in=quad_dd, dd_out=sym.DD_VOLUME) - stiff_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol") quad_v = to_quad(self.v) quad_u = to_quad(u) return sym.InverseMassOperator()( - (stiff_t[0](quad_u * quad_v[0]) + stiff_t[1](quad_u * quad_v[1])) - - sym.FaceMassOperator(all_faces_dd, "vol")( + sum(stiff_t_op[n](quad_u * quad_v[n]) + for n in range(self.ambient_dim)) + - sym.FaceMassOperator(face_dd, sym.DD_VOLUME)( flux(sym.int_tpair(u, self.quad_tag)) - + flux(sym.bv_tpair(boundary_dd, u, bc_in)) + + flux(sym.bv_tpair(boundary_dd, u, self.inflow_u)) # FIXME: Add back support for inflow/outflow tags #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in)) -- GitLab From b7fbbdfa91148aedd56f3ad9f322ceab44c90ae4 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 9 May 2020 21:00:56 -0500 Subject: [PATCH 2/3] fix RefStiffnessTOp for the 1D oversampled case --- examples/advection/var-velocity.py | 58 +++++++++++++++++++++++++++--- examples/advection/weak.py | 7 +++- grudge/symbolic/operators.py | 5 ++- test/test_grudge.py | 2 +- 4 files changed, 65 insertions(+), 7 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index c8791efd..7d4c91f5 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 770e9717..e780f7f8 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 8c2f7fd0..1950d3f7 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 12b7b116..f068d34d 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 -- GitLab From 59038c352efe870f5ebe274c0e4c4814f39cb10d Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 14 May 2020 11:20:21 -0500 Subject: [PATCH 3/3] examples: add ylim to 1d plots --- examples/advection/var-velocity.py | 32 +++++++++++++--------- examples/advection/weak.py | 43 +++++++++++++++++++----------- 2 files changed, 48 insertions(+), 27 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 7d4c91f5..b6beb7a2 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -22,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +import os import numpy as np import pyopencl as cl @@ -38,11 +39,9 @@ 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) - + def __init__(self, queue, discr, order, visualize=True, ylim=None): self.queue = queue - self.x = volume_discr.nodes().get(self.queue) + self.dim = discr.ambient_dim self.visualize = visualize if not self.visualize: @@ -51,33 +50,41 @@ class Plotter: if self.dim == 1: import matplotlib.pyplot as pt self.fig = pt.figure(figsize=(8, 8), dpi=300) + self.ylim = ylim + + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + self.x = volume_discr.nodes().get(self.queue) 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): + def __call__(self, evt, basename, overwrite=True): if not self.visualize: return if self.dim == 1: u = evt.state_component.get(self.queue) + filename = "%s.png" % basename + if not overwrite and os.path.exists(filename): + from meshmode import FileExistsError + raise FileExistsError("output file '%s' already exists" % filename) + ax = self.fig.gca() ax.plot(self.x[0], u, '-') ax.plot(self.x[0], u, 'k.') + if self.ylim is not None: + ax.set_ylim(self.ylim) + ax.set_xlabel("$x$") ax.set_ylabel("$u$") ax.set_title("t = {:.2f}".format(evt.t)) - self.fig.savefig("%s.png" % basename) + self.fig.savefig(filename) self.fig.clf() else: self.vis.write_vtk_file("%s.vtu" % basename, [ ("u", evt.state_component) - ], overwrite=True) + ], overwrite=overwrite) # }}} @@ -180,7 +187,8 @@ 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) - plot = Plotter(queue, discr, order, visualize=visualize) + plot = Plotter(queue, discr, order, visualize=visualize, + ylim=[-0.1, 1.1]) step = 0 norm = bind(discr, sym.norm(2, sym.var("u"))) diff --git a/examples/advection/weak.py b/examples/advection/weak.py index e780f7f8..14ed08e1 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -17,6 +17,7 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . """ +import os import numpy as np import numpy.linalg as la @@ -33,11 +34,9 @@ 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) - + def __init__(self, queue, discr, order, visualize=True, ylim=None): self.queue = queue - self.x = volume_discr.nodes().get(self.queue) + self.dim = discr.ambient_dim self.visualize = visualize if not self.visualize: @@ -46,33 +45,42 @@ class Plotter: if self.dim == 1: import matplotlib.pyplot as pt self.fig = pt.figure(figsize=(8, 8), dpi=300) + self.ylim = ylim + + volume_discr = discr.discr_from_dd(sym.DD_VOLUME) + self.x = volume_discr.nodes().get(self.queue) 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): + def __call__(self, evt, basename, overwrite=True): if not self.visualize: return if self.dim == 1: u = evt.state_component.get(self.queue) + filename = "%s.png" % basename + if not overwrite and os.path.exists(filename): + from meshmode import FileExistsError + raise FileExistsError("output file '%s' already exists" % filename) + ax = self.fig.gca() ax.plot(self.x[0], u, '-') ax.plot(self.x[0], u, 'k.') + if self.ylim is not None: + ax.set_ylim(self.ylim) + ax.set_xlabel("$x$") ax.set_ylabel("$u$") ax.set_title("t = {:.2f}".format(evt.t)) - self.fig.savefig("%s.png" % basename) + + self.fig.savefig(filename) self.fig.clf() else: self.vis.write_vtk_file("%s.vtu" % basename, [ ("u", evt.state_component) - ], overwrite=True) + ], overwrite=overwrite) # }}} @@ -145,18 +153,23 @@ def main(ctx_factory, dim=2, order=4, visualize=True): from grudge.shortcuts import set_up_rk4 dt_stepper = set_up_rk4("u", dt, u, rhs) - plot = Plotter(queue, discr, order, visualize=visualize) + plot = Plotter(queue, discr, order, visualize=visualize, + ylim=[-1.1, 1.1]) - step = 0 norm = bind(discr, sym.norm(2, sym.var("u"))) + + step = 0 + norm_u = 0.0 for event in dt_stepper.run(t_end=final_time): if not isinstance(event, dt_stepper.StateComputed): continue + if step % 10 == 0: + norm_u = norm(queue, u=event.state_component) + plot(event, "fld-weak-%04d" % step) + step += 1 - norm_u = norm(queue, u=event.state_component) logger.info("[%04d] t = %.5f |u| = %.5e", step, event.t, norm_u) - plot(event, "fld-weak-%04d" % step) # }}} -- GitLab