From c487a10e5e16f94c6ae3be86b435c68b81d2b249 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 28 Mar 2017 05:58:36 -0500 Subject: [PATCH 1/5] Implemented variable coefficient waves, and used loopy to do symbolic if statements --- examples/wave/var-propagation-speed.py | 263 +++++++---------- grudge/execution.py | 35 ++- grudge/models/wave.py | 378 +++++++++++++++++-------- 3 files changed, 385 insertions(+), 291 deletions(-) diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index 3f4b9927..24719059 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -1,11 +1,8 @@ -"""Variable-coefficient wave propagation.""" +"""Minimal example of a grudge driver.""" -from __future__ import division -from __future__ import absolute_import -from __future__ import print_function -from six.moves import range +from __future__ import division, print_function -__copyright__ = "Copyright (C) 2009 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2015 Andreas Kloeckner" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -29,171 +26,105 @@ THE SOFTWARE. import numpy as np -from grudge.mesh import BTAG_ALL, BTAG_NONE - - -def main(write_output=True, - dir_tag=BTAG_NONE, - neu_tag=BTAG_NONE, - rad_tag=BTAG_ALL, - flux_type_arg="upwind"): - from math import sin, cos, pi, exp, sqrt # noqa - - from grudge.backends import guess_run_context - rcon = guess_run_context() - - dim = 2 - - if dim == 1: - if rcon.is_head_rank: - from grudge.mesh.generator import make_uniform_1d_mesh - mesh = make_uniform_1d_mesh(-10, 10, 500) - elif dim == 2: - from grudge.mesh.generator import make_rect_mesh - if rcon.is_head_rank: - mesh = make_rect_mesh(a=(-1, -1), b=(1, 1), max_area=0.003) - elif dim == 3: - if rcon.is_head_rank: - from grudge.mesh.generator import make_ball_mesh - mesh = make_ball_mesh(max_volume=0.0005) - else: - raise RuntimeError("bad number of dimensions") - - if rcon.is_head_rank: - print("%d elements" % len(mesh.elements)) - mesh_data = rcon.distribute_mesh(mesh) - else: - mesh_data = rcon.receive_mesh() - - discr = rcon.make_discretization(mesh_data, order=4) - - from grudge.timestep.runge_kutta import LSRK4TimeStepper - stepper = LSRK4TimeStepper() - - from grudge.visualization import VtkVisualizer - if write_output: - vis = VtkVisualizer(discr, rcon, "fld") - - source_center = np.array([0.7, 0.4]) - source_width = 1/16 - source_omega = 3 - - import grudge.symbolic as sym - sym_x = sym.nodes(2) - sym_source_center_dist = sym_x - source_center +import pyopencl as cl +from grudge.shortcuts import set_up_rk4 +from grudge import sym, bind, Discretization - from grudge.models.wave import VariableVelocityStrongWaveOperator - op = VariableVelocityStrongWaveOperator( - c=sym.If(sym.Comparison( - np.dot(sym_x, sym_x), "<", 0.4**2), - 1, 0.5), - dimensions=discr.dimensions, - source= - sym.CFunction("sin")(source_omega*sym.ScalarParameter("t")) - * sym.CFunction("exp")( - -np.dot(sym_source_center_dist, sym_source_center_dist) - / source_width**2), - dirichlet_tag=dir_tag, - neumann_tag=neu_tag, - radiation_tag=rad_tag, - flux_type=flux_type_arg - ) - - from grudge.tools import join_fields - fields = join_fields(discr.volume_zeros(), - [discr.volume_zeros() for i in range(discr.dimensions)]) - - # {{{ diagnostics setup - - from pytools.log import LogManager, \ - add_general_quantities, \ - add_simulation_quantities, \ - add_run_info - - if write_output: - log_file_name = "wave.dat" - else: - log_file_name = None - - logmgr = LogManager(log_file_name, "w", rcon.communicator) - add_run_info(logmgr) - add_general_quantities(logmgr) - add_simulation_quantities(logmgr) - discr.add_instrumentation(logmgr) - - from pytools.log import IntervalTimer - vis_timer = IntervalTimer("t_vis", "Time spent visualizing") - logmgr.add_quantity(vis_timer) - stepper.add_instrumentation(logmgr) - - from grudge.log import LpNorm - u_getter = lambda: fields[0] - logmgr.add_quantity(LpNorm(u_getter, discr, 1, name="l1_u")) - logmgr.add_quantity(LpNorm(u_getter, discr, name="l2_u")) - - logmgr.add_watches(["step.max", "t_sim.max", "l2_u", "t_step.max"]) - - # }}} - - # {{{ timestep loop - - rhs = op.bind(discr) - try: - from grudge.timestep.stability import \ - approximate_rk4_relative_imag_stability_region - max_dt = ( - 1/discr.compile(op.max_eigenvalue_expr())() - * discr.dt_non_geometric_factor() - * discr.dt_geometric_factor() - * approximate_rk4_relative_imag_stability_region(stepper)) - if flux_type_arg == "central": - max_dt *= 0.25 - - from grudge.timestep import times_and_steps - step_it = times_and_steps(final_time=3, logmgr=logmgr, - max_dt_getter=lambda t: max_dt) - - for step, t, dt in step_it: - if step % 10 == 0 and write_output: - visf = vis.make_file("fld-%04d" % step) - - vis.add_data(visf, - [ - ("u", fields[0]), - ("v", fields[1:]), - ], - time=t, - step=step) - visf.close() - - fields = stepper(fields, t, dt, rhs) - assert discr.norm(fields) < 1 - finally: - if write_output: - vis.close() +def main(write_output=True, order=4): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) - logmgr.close() - discr.close() + dims = 2 + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*dims, + b=(0.5,)*dims, + n=(8,)*dims) - # }}} - -if __name__ == "__main__": - main(flux_type_arg="upwind") + discr = Discretization(cl_ctx, mesh, order=order) + source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim] + source_width = 0.05 + source_omega = 3 -# entry points for py.test ---------------------------------------------------- -def test_var_velocity_wave(): - from pytools.test import mark_test - mark_long = mark_test.long + 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.2), + -0.1, -0.2) + #c = np.dot(sym_x,sym_x) + + 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="central") + + queue = cl.CommandQueue(discr.cl_context) + from pytools.obj_array import join_fields + fields = join_fields(discr.zeros(queue), + [discr.zeros(queue) for i in range(discr.dim)]) + + # FIXME + #dt = op.estimate_rk4_timestep(discr, fields=fields) + + op.check_bc_coverage(mesh) + + # print(sym.pretty(op.sym_operator())) + bound_op = bind(discr, op.sym_operator()) + print(bound_op) + # 1/0 + + def rhs(t, w): + return bound_op(queue, t=t, w=w) + + if mesh.dim == 2: + dt = 0.04 + elif mesh.dim == 3: + dt = 0.02 + + dt_stepper = set_up_rk4("w", dt, fields, rhs) + + final_t = 10 + nsteps = int(final_t/dt) + print("dt=%g nsteps=%d" % (dt, nsteps)) + + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + + step = 0 + + norm = bind(discr, sym.norm(2, sym.var("u"))) + + from time import time + t_last_step = time() + + for event in dt_stepper.run(t_end=final_t): + if isinstance(event, dt_stepper.StateComputed): + assert event.component_id == "w" + + step += 1 + + print(step, event.t, norm(queue, u=event.state_component[0]), + time()-t_last_step) + if step % 10 == 0: + vis.write_vtk_file("fld-%04d.vtu" % step, + [ + ("u", event.state_component[0]), + ("v", event.state_component[1:]), + ]) + t_last_step = time() - for flux_type in ["upwind", "central"]: - yield ("dirichlet var-v wave equation with %s flux" % flux_type, - mark_long(main), - False, BTAG_ALL, BTAG_NONE, TAG_NONE, flux_type) - yield ("neumann var-v wave equation", mark_long(main), - False, BTAG_NONE, BTAG_ALL, TAG_NONE) - yield ("radiation-bc var-v wave equation", mark_long(main), - False, BTAG_NONE, TAG_NONE, BTAG_ALL) -# vim: foldmethod=marker +if __name__ == "__main__": + main() diff --git a/grudge/execution.py b/grudge/execution.py index 0241c82b..f78c345f 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -119,12 +119,37 @@ class ExecutionMapper(mappers.Evaluator, then = self.rec(expr.then) else_ = self.rec(expr.else_) - result = cl.array.empty_like(then, queue=self.queue, - allocator=self.bound_op.allocator) - cl.array.if_positive(bool_crit, then, else_, out=result, - queue=self.queue) + import pymbolic.primitives as p - return result + var = p.Variable + + i = var("i") #sym.Variable("i") + if isinstance(then, pyopencl.array.Array): + sym_then = var("a")[i] #sym.Variable("a")[i] + else: + sym_then = var("a") # sym.Variable("a") + then = np.float64(then) + + if isinstance(else_, pyopencl.array.Array): + sym_else = var("b")[i] # sym.Variable("b")[i] + else: + sym_else = var("b") # sym.Variable("b") + else_ = np.float64(else_) + + @memoize_in(self.bound_op, "map_if_knl") + def knl(): + knl = lp.make_kernel( + "{[i]: 0<=i 0: + self.sign = 1 + else: + self.sign = -1 self.dirichlet_tag = dirichlet_tag self.neumann_tag = neumann_tag self.radiation_tag = radiation_tag - self.flux_type = flux_type + self.dirichlet_bc_f = dirichlet_bc_f - self.diffusion_coeff = diffusion_coeff - self.diffusion_scheme = diffusion_scheme + self.flux_type = flux_type - # {{{ flux ---------------------------------------------------------------- - def flux(self, w, c_vol): + def flux(self, w): u = w[0] v = w[1:] - normal = sym.normal(w.tag, self.ambient_dim) + normal = sym.normal(w.dd, self.ambient_dim) - c = sym.RestrictToBoundary(w.tag)(c_vol) + flux_weak = join_fields( + np.dot(v.avg, normal), + u.avg * normal) - flux = self.time_sign*1/2*join_fields( - c.ext * np.dot(v.ext, normal) - - c.int * np.dot(v.int, normal), - normal*(c.ext*u.ext - c.int*u.int)) if self.flux_type == "central": - pass + return -self.c*flux_weak elif self.flux_type == "upwind": - flux += join_fields( - c.ext*u.ext - c.int*u.int, - c.ext*normal*np.dot(normal, v.ext) - - c.int*normal*np.dot(normal, v.int) - ) + return -self.c*(flux_weak + self.sign*join_fields( + 0.5*(u.int-u.ext), + 0.5*(normal * np.dot(normal, v.int-v.ext)))) + # see doc/notes/grudge-notes.tm + # THis is terrible + #flux_weak -= self.sign*join_fields( + #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( + #np.dot(v.int, normal), + #u.int * normal) - flux_weak + + #return self.c*flux_strong + + def sym_operator(self): + d = self.ambient_dim + + w = sym.make_sym_array("w", d+1) + u = w[0] + v = w[1:] + + # boundary conditions ------------------------------------------------- + + # dirichlet BCs ------------------------------------------------------- + dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u)) + dir_v = sym.cse(sym.interp("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.Field("dir_bc_u") + dir_bc = join_fields(2*dir_g - dir_u, dir_v) else: - raise ValueError("invalid flux type '%s'" % self.flux_type) + dir_bc = join_fields(-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") + + # radiation BCs ------------------------------------------------------- + rad_normal = sym.normal(self.radiation_tag, d) + + rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u)) + rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v)) - return flux + rad_bc = sym.cse(join_fields( + 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") - # }}} + # entire operator ----------------------------------------------------- + nabla = sym.nabla(d) - def bind_characteristic_velocity(self, discr): - from grudge.symbolic.operators import ElementwiseMaxOperator + def flux(pair): + return sym.interp(pair.dd, "all_faces")(self.flux(pair)) + + + #result = ( + #- join_fields( + #-self.c*np.dot(nabla, v), + #-self.c*(nabla*u) + #) + #+ + #sym.InverseMassOperator()( + #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)) + #))) + + result = sym.InverseMassOperator()( + join_fields( + -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v), + -self.c*(sym.stiffness_t(self.ambient_dim)*u) + ) - compiled = discr.compile(ElementwiseMaxOperator()(self.c)) - def do(t, w, **extra_context): - return compiled(t=t, w=w, **extra_context) + - 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)) - return do + ) ) - def sym_operator(self, with_sensor=False): - d = self.ambient_dim + result[0] += self.source_f - w = sym.make_sym_array("w", d+1) - u = w[0] - v = w[1:] + return result + + def check_bc_coverage(self, mesh): + from meshmode.mesh import check_bc_coverage + check_bc_coverage(mesh, [ + self.dirichlet_tag, + self.neumann_tag, + self.radiation_tag]) + + def max_eigenvalue(self, t, fields=None, discr=None): + return abs(self.c) - flux_w = join_fields(self.c, w) - # {{{ boundary conditions +# }}} - # Dirichlet - dir_c = sym.RestrictToBoundary(self.dirichlet_tag) * self.c - dir_u = sym.RestrictToBoundary(self.dirichlet_tag) * u - dir_v = sym.RestrictToBoundary(self.dirichlet_tag) * v - dir_bc = join_fields(dir_c, -dir_u, dir_v) +# {{{ variable-velocity - # Neumann - neu_c = sym.RestrictToBoundary(self.neumann_tag) * self.c - neu_u = sym.RestrictToBoundary(self.neumann_tag) * u - neu_v = sym.RestrictToBoundary(self.neumann_tag) * v +class VariableCoefficientWeakWaveOperator(HyperbolicOperator): + """This operator discretizes the wave equation + :math:`\\partial_t^2 u = c^2 \\Delta u`. - neu_bc = join_fields(neu_c, neu_u, -neu_v) + To be precise, we discretize the hyperbolic system - # Radiation - rad_normal = sym.make_normal(self.radiation_tag, d) + .. math:: - rad_c = sym.RestrictToBoundary(self.radiation_tag) * self.c - rad_u = sym.RestrictToBoundary(self.radiation_tag) * u - rad_v = sym.RestrictToBoundary(self.radiation_tag) * v + \partial_t u - c \\nabla \\cdot v = 0 - rad_bc = join_fields( - rad_c, - 0.5*(rad_u - self.time_sign*np.dot(rad_normal, rad_v)), - 0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.time_sign*rad_u) - ) + \partial_t v - c \\nabla u = 0 - # }}} + The sign of :math:`v` determines whether we discretize the forward or the + backward wave equation. - # {{{ diffusion ------------------------------------------------------- - from pytools.obj_array import with_object_array_or_scalar + :math:`c` is assumed to be constant across all space. + """ - def make_diffusion(arg): - if with_sensor or ( - self.diffusion_coeff is not None and self.diffusion_coeff != 0): - if self.diffusion_coeff is None: - diffusion_coeff = 0 - else: - diffusion_coeff = self.diffusion_coeff + def __init__(self, c, ambient_dim, source_f=0, + 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 with_sensor: - diffusion_coeff += sym.Field("sensor") + self.c = c + self.ambient_dim = ambient_dim + self.source_f = source_f - from grudge.second_order import SecondDerivativeTarget + self.sign = sym.If(sym.Comparison( + self.c, ">", 0), + 1, -1) - # strong_form here allows the reuse the value of grad u. - grad_tgt = SecondDerivativeTarget( - self.ambient_dim, strong_form=True, - operand=arg) + self.dirichlet_tag = dirichlet_tag + self.neumann_tag = neumann_tag + self.radiation_tag = radiation_tag - self.diffusion_scheme.grad(grad_tgt, bc_getter=None, - dirichlet_tags=[], neumann_tags=[]) + self.dirichlet_bc_f = dirichlet_bc_f + + self.flux_type = flux_type + + def flux(self, w): + u = w[0] + v = w[1:] + normal = sym.normal(w.dd, self.ambient_dim) + + surf_c = sym.interp("vol", u.dd)(self.c) + + flux_weak = join_fields( + np.dot(v.avg, normal), + u.avg * normal) + return -surf_c*flux_weak + + + if self.flux_type == "central": + return -surf_c*flux_weak + #elif self.flux_type == "upwind": + # see doc/notes/grudge-notes.tm + #flux_weak -= self.sign*join_fields( + #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( + #np.dot(v.int, normal), + #u.int * normal) - flux_weak - div_tgt = SecondDerivativeTarget( - self.ambient_dim, strong_form=False, - operand=diffusion_coeff*grad_tgt.minv_all) + #return self.c*flux_strong - self.diffusion_scheme.div(div_tgt, - bc_getter=None, - dirichlet_tags=[], neumann_tags=[]) + def sym_operator(self): + d = self.ambient_dim - return div_tgt.minv_all - else: - return 0 + w = sym.make_sym_array("w", d+1) + u = w[0] + v = w[1:] - # }}} + # boundary conditions ------------------------------------------------- + + # dirichlet BCs ------------------------------------------------------- + dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u)) + dir_v = sym.cse(sym.interp("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.Field("dir_bc_u") + dir_bc = join_fields(2*dir_g - dir_u, dir_v) + else: + dir_bc = join_fields(-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") + + # radiation BCs ------------------------------------------------------- + rad_normal = sym.normal(self.radiation_tag, d) + + 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( + 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) -sym.interp("vol",self.radiation_tag)(self.sign)*rad_u) + ), "rad_bc") # entire operator ----------------------------------------------------- nabla = sym.nabla(d) - flux_op = sym.get_flux_operator(self.flux()) - return ( - - join_fields( - - self.time_sign*self.c*np.dot(nabla, v) - make_diffusion(u) - + self.source, + def flux(pair): + return sym.interp(pair.dd, "all_faces")(self.flux(pair)) + + + #result = sym.InverseMassOperator()( + #join_fields( + #-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)) - -self.time_sign*self.c*(nabla*u) - with_object_array_or_scalar( - make_diffusion, v) + #)# ) + + result = sym.InverseMassOperator()( + join_fields( + -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v), + -self.c*(sym.stiffness_t(self.ambient_dim)*u) ) - + - sym.InverseMassOperator() * ( - flux_op(flux_w) - + flux_op(sym.BoundaryPair(flux_w, dir_bc, self.dirichlet_tag)) - + flux_op(sym.BoundaryPair(flux_w, neu_bc, self.neumann_tag)) - + flux_op(sym.BoundaryPair(flux_w, rad_bc, self.radiation_tag)) - )) + + + - 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)) + + ) ) + + result[0] += self.source_f + + return result def check_bc_coverage(self, mesh): from meshmode.mesh import check_bc_coverage @@ -369,9 +507,9 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator): self.neumann_tag, self.radiation_tag]) - def max_eigenvalue_expr(self): - import grudge.symbolic as sym - return sym.NodalMax()(sym.CFunction("fabs")(self.c)) + def max_eigenvalue(self, t, fields=None, discr=None): + return abs(self.c) + # }}} -- GitLab From 99326552c0d1b516c073969a13c86304a6e12cc3 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sat, 1 Apr 2017 20:37:58 -0500 Subject: [PATCH 2/5] Cleaning up code to satisfy flake --- examples/wave/var-propagation-speed.py | 13 +--- grudge/execution.py | 26 +++---- grudge/models/wave.py | 93 ++++++-------------------- 3 files changed, 37 insertions(+), 95 deletions(-) diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index 24719059..cb5d3e9e 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -51,10 +51,9 @@ def main(write_output=True, order=4): 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.2), - -0.1, -0.2) - #c = np.dot(sym_x,sym_x) + c = sym.If(sym.Comparison( + np.dot(sym_x, sym_x), "<", 0.15), + np.float32(-0.1), np.float32(-0.2)) from grudge.models.wave import VariableCoefficientWeakWaveOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE @@ -75,15 +74,9 @@ def main(write_output=True, order=4): fields = join_fields(discr.zeros(queue), [discr.zeros(queue) for i in range(discr.dim)]) - # FIXME - #dt = op.estimate_rk4_timestep(discr, fields=fields) - op.check_bc_coverage(mesh) - # print(sym.pretty(op.sym_operator())) bound_op = bind(discr, op.sym_operator()) - print(bound_op) - # 1/0 def rhs(t, w): return bound_op(queue, t=t, w=w) diff --git a/grudge/execution.py b/grudge/execution.py index f78c345f..e3646b73 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -120,30 +120,32 @@ class ExecutionMapper(mappers.Evaluator, else_ = self.rec(expr.else_) import pymbolic.primitives as p - var = p.Variable - i = var("i") #sym.Variable("i") + i = var("i") if isinstance(then, pyopencl.array.Array): - sym_then = var("a")[i] #sym.Variable("a")[i] + sym_then = var("a")[i] + elif isinstance(then, np.number): + sym_then = var("a") else: - sym_then = var("a") # sym.Variable("a") - then = np.float64(then) + raise TypeError( + "Expected parameter to be of type np.number or pyopencl.array.Array") if isinstance(else_, pyopencl.array.Array): - sym_else = var("b")[i] # sym.Variable("b")[i] + sym_else = var("b")[i] + elif isinstance(then, np.number): + sym_else = var("b") else: - sym_else = var("b") # sym.Variable("b") - else_ = np.float64(else_) + raise TypeError( + "Expected parameter to be of type np.number or pyopencl.array.Array") @memoize_in(self.bound_op, "map_if_knl") def knl(): knl = lp.make_kernel( - "{[i]: 0<=i", 0), - 1, -1) + np.int32(1), np.int32(-1)) self.dirichlet_tag = dirichlet_tag self.neumann_tag = neumann_tag @@ -401,24 +372,15 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): flux_weak = join_fields( np.dot(v.avg, normal), u.avg * normal) - return -surf_c*flux_weak - if self.flux_type == "central": return -surf_c*flux_weak - #elif self.flux_type == "upwind": - # see doc/notes/grudge-notes.tm - #flux_weak -= self.sign*join_fields( - #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( - #np.dot(v.int, normal), - #u.int * normal) - flux_weak - - #return self.c*flux_strong + elif self.flux_type == "upwind": + return -self.c*(flux_weak + self.sign*join_fields( + 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) def sym_operator(self): d = self.ambient_dim @@ -457,44 +419,29 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v)) rad_bc = sym.cse(join_fields( - 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) -sym.interp("vol",self.radiation_tag)(self.sign)*rad_u) + 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) + - sym.interp("vol", self.radiation_tag)(self.sign)*rad_u) ), "rad_bc") # entire operator ----------------------------------------------------- - nabla = sym.nabla(d) - def flux(pair): return sym.interp(pair.dd, "all_faces")(self.flux(pair)) - - #result = sym.InverseMassOperator()( - #join_fields( - #-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)) - - #)# ) - - result = sym.InverseMassOperator()( + result = sym.InverseMassOperator()( join_fields( -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)) + - 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)) - ) ) + )) result[0] += self.source_f -- GitLab From a6da4eab18b8069c9c44b3fa2cfb94ef3dc4b23e Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 10 Apr 2017 17:34:24 -0500 Subject: [PATCH 3/5] There was an attempt at the right way of doing upwind fluxes --- examples/wave/var-propagation-speed.py | 2 +- grudge/models/wave.py | 55 +++++++++++++++++--------- 2 files changed, 38 insertions(+), 19 deletions(-) diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index cb5d3e9e..03dcdc8e 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -67,7 +67,7 @@ def main(write_output=True, order=4): dirichlet_tag=BTAG_NONE, neumann_tag=BTAG_NONE, radiation_tag=BTAG_ALL, - flux_type="central") + flux_type="upwind") queue = cl.CommandQueue(discr.cl_context) from pytools.obj_array import join_fields diff --git a/grudge/models/wave.py b/grudge/models/wave.py index 01969125..02aea3a5 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -363,22 +363,37 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): self.flux_type = flux_type def flux(self, w): - u = w[0] - v = w[1:] + c = w[0] + u = w[1] + v = w[2:] normal = sym.normal(w.dd, self.ambient_dim) - surf_c = sym.interp("vol", u.dd)(self.c) - flux_weak = join_fields( - np.dot(v.avg, normal), - u.avg * normal) + np.dot(v.int*c.int + v.ext*c.ext, normal), + (u.int * c.int + u.ext*c.ext) * normal) * -0.5 + + flux = sym.interp("vol",u.dd)(self.sign) * join_fields( + c.ext * np.dot(v.ext, normal) - c.int*np.dot(v.int,normal), + normal * (c.ext*u.ext - c.int* u.int)) * -0.5 + + flux_strong = join_fields( + np.dot(v.int, normal), + u.int * normal) * -c.int if self.flux_type == "central": - return -surf_c*flux_weak + return flux_strong - flux elif self.flux_type == "upwind": - return -self.c*(flux_weak + self.sign*join_fields( - 0.5*(u.int-u.ext), - 0.5*(normal * np.dot(normal, v.int-v.ext)))) + flux += join_fields( + c.ext*u.ext - c.int*u.int, + c.ext*normal*np.dot(normal, v.ext) + - c.int*normal*np.dot(normal, v.int) + ) + return flux_strong - flux + return (flux_weak + join_fields( + c.ext*u.ext - c.int*u.int, + c.ext*normal*np.dot(normal, v.ext) + - c.int*normal*np.dot(normal, v.int) + )) else: raise ValueError("invalid flux type '%s'" % self.flux_type) @@ -388,10 +403,12 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): w = sym.make_sym_array("w", d+1) u = w[0] v = w[1:] + flux_w = join_fields(self.c,w) # boundary conditions ------------------------------------------------- # dirichlet BCs ------------------------------------------------------- + dir_c = sym.cse(sym.interp("vol", self.dirichlet_tag)(self.c)) dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u)) dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v)) if self.dirichlet_bc_f: @@ -401,24 +418,26 @@ class VariableCoefficientWeakWaveOperator(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 = join_fields(dir_c,2*dir_g - dir_u, dir_v) else: - dir_bc = join_fields(-dir_u, dir_v) + dir_bc = join_fields(dir_c,-dir_u, dir_v) dir_bc = sym.cse(dir_bc, "dir_bc") # neumann BCs --------------------------------------------------------- + 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_u, -neu_v), "neu_bc") + neu_bc = sym.cse(join_fields(neu_c,neu_u, -neu_v), "neu_bc") # radiation BCs ------------------------------------------------------- rad_normal = sym.normal(self.radiation_tag, d) + rad_c = sym.cse(sym.interp("vol", self.radiation_tag)(self.c)) 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(join_fields(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,10 +455,10 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): ) - - 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)) + - 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)) )) -- GitLab From 255f2659d4982c908c8887ceceaff4dc644e84d5 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 17 Apr 2017 20:28:02 -0500 Subject: [PATCH 4/5] Did variable coefficient upwind flux the right way, maybe? --- contrib/maxima/wave.mac | 4 ++-- grudge/models/wave.py | 30 ++++++------------------------ 2 files changed, 8 insertions(+), 26 deletions(-) diff --git a/contrib/maxima/wave.mac b/contrib/maxima/wave.mac index cc812f5b..9ced321a 100644 --- a/contrib/maxima/wave.mac +++ b/contrib/maxima/wave.mac @@ -79,11 +79,11 @@ print("Wave equation flux in terms of characteristic variables:"); print(wave_sflux); -print("Wave equation flux in terms of physical variables:"); -print(wave_wflux); print("Weak flux divided by (-c), as implemented in StrongWaveOperator:"); print(hypsimp(ev(wave_wflux, cp=c, cm=c)/(-c))); */ +print("Wave equation weak flux in terms of physical variables:"); +print(wave_wflux); print("Strong-form wave equation flux in terms of physical variables:"); /*print(wave_strongwflux);*/ diff --git a/grudge/models/wave.py b/grudge/models/wave.py index 02aea3a5..6db93661 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -368,32 +368,15 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): v = w[2:] normal = sym.normal(w.dd, self.ambient_dim) - flux_weak = join_fields( + if self.flux_type == "central": + return join_fields( np.dot(v.int*c.int + v.ext*c.ext, normal), (u.int * c.int + u.ext*c.ext) * normal) * -0.5 - flux = sym.interp("vol",u.dd)(self.sign) * join_fields( - c.ext * np.dot(v.ext, normal) - c.int*np.dot(v.int,normal), - normal * (c.ext*u.ext - c.int* u.int)) * -0.5 - - flux_strong = join_fields( - np.dot(v.int, normal), - u.int * normal) * -c.int - - if self.flux_type == "central": - return flux_strong - flux elif self.flux_type == "upwind": - flux += join_fields( - c.ext*u.ext - c.int*u.int, - c.ext*normal*np.dot(normal, v.ext) - - c.int*normal*np.dot(normal, v.int) - ) - return flux_strong - flux - return (flux_weak + join_fields( - c.ext*u.ext - c.int*u.int, - c.ext*normal*np.dot(normal, v.ext) - - c.int*normal*np.dot(normal, v.int) - )) + return join_fields(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) + c.ext*u.ext + c.int * u.int ) ) * -0.5 + else: raise ValueError("invalid flux type '%s'" % self.flux_type) @@ -474,8 +457,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): self.radiation_tag]) def max_eigenvalue(self, t, fields=None, discr=None): - return abs(self.c) - + return sym.NodalMax()(sym.CFunction("fabs")(self.c)) # }}} -- GitLab From 2a553e5af02222c8dd35fd836eec446212a5915e Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 17 Apr 2017 20:34:28 -0500 Subject: [PATCH 5/5] Forgot to lint --- grudge/models/wave.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/grudge/models/wave.py b/grudge/models/wave.py index 6db93661..19229abe 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -369,13 +369,15 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): normal = sym.normal(w.dd, self.ambient_dim) if self.flux_type == "central": - return join_fields( + return -0.5 * join_fields( np.dot(v.int*c.int + v.ext*c.ext, normal), - (u.int * c.int + u.ext*c.ext) * normal) * -0.5 + (u.int * c.int + u.ext*c.ext) * normal) elif self.flux_type == "upwind": - return join_fields(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) + c.ext*u.ext + c.int * u.int ) ) * -0.5 + return -0.5 * join_fields(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) + + c.ext*u.ext + c.int * u.int)) else: raise ValueError("invalid flux type '%s'" % self.flux_type) @@ -386,7 +388,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 = join_fields(self.c, w) # boundary conditions ------------------------------------------------- @@ -401,9 +403,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 = join_fields(dir_c, 2*dir_g - dir_u, dir_v) else: - dir_bc = join_fields(dir_c,-dir_u, dir_v) + dir_bc = join_fields(dir_c, -dir_u, dir_v) dir_bc = sym.cse(dir_bc, "dir_bc") @@ -411,7 +413,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(join_fields(neu_c, neu_u, -neu_v), "neu_bc") # radiation BCs ------------------------------------------------------- rad_normal = sym.normal(self.radiation_tag, d) -- GitLab