diff --git a/contrib/maxima/wave.mac b/contrib/maxima/wave.mac index cc812f5bc8725e21c24dbebe77e470e146881fc8..9ced321a7e9e95531c39faa11da862e911a4aa38 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/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index 3f4b9927c0436896d94e95e78f2d1325ec4a2045..03dcdc8e931062a5faf793720d95b5ea9c3457b0 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,98 @@ 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 +import pyopencl as cl +from grudge.shortcuts import set_up_rk4 +from grudge import sym, bind, Discretization + + +def main(write_output=True, order=4): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + 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) + + 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 - import grudge.symbolic as sym - sym_x = sym.nodes(2) + 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.15), + np.float32(-0.1), np.float32(-0.2)) - 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() + 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="upwind") - fields = stepper(fields, t, dt, rhs) + 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)]) - assert discr.norm(fields) < 1 - finally: - if write_output: - vis.close() + op.check_bc_coverage(mesh) - logmgr.close() - discr.close() + bound_op = bind(discr, op.sym_operator()) - # }}} + def rhs(t, w): + return bound_op(queue, t=t, w=w) -if __name__ == "__main__": - main(flux_type_arg="upwind") + 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)) -# entry points for py.test ---------------------------------------------------- -def test_var_velocity_wave(): - from pytools.test import mark_test - mark_long = mark_test.long + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) - 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) + step = 0 -# vim: foldmethod=marker + 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() + + +if __name__ == "__main__": + main() diff --git a/grudge/execution.py b/grudge/execution.py index 0241c82b532c9167c60362917a08c73123c385df..e3646b7361d8b1606a1a234fbe6089a51127f8b7 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -119,12 +119,39 @@ 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 + var = p.Variable + + i = var("i") + if isinstance(then, pyopencl.array.Array): + sym_then = var("a")[i] + elif isinstance(then, np.number): + sym_then = var("a") + else: + raise TypeError( + "Expected parameter to be of type np.number or pyopencl.array.Array") - return result + if isinstance(else_, pyopencl.array.Array): + sym_else = var("b")[i] + elif isinstance(then, np.number): + sym_else = var("b") + 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: + 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) - - c = sym.RestrictToBoundary(w.tag)(c_vol) + normal = sym.normal(w.dd, self.ambient_dim) - 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)) + flux_weak = join_fields( + np.dot(v.avg, normal), + u.avg * normal) 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)))) else: raise ValueError("invalid flux type '%s'" % self.flux_type) - return flux + def sym_operator(self): + d = self.ambient_dim - # }}} + w = sym.make_sym_array("w", d+1) + u = w[0] + v = w[1:] - def bind_characteristic_velocity(self, discr): - from grudge.symbolic.operators import ElementwiseMaxOperator + # boundary conditions ------------------------------------------------- - compiled = discr.compile(ElementwiseMaxOperator()(self.c)) + # 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.") - def do(t, w, **extra_context): - return compiled(t=t, w=w, **extra_context) + 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) - return do + dir_bc = sym.cse(dir_bc, "dir_bc") - def sym_operator(self, with_sensor=False): - d = self.ambient_dim + # 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") - w = sym.make_sym_array("w", d+1) - u = w[0] - v = w[1:] + # radiation BCs ------------------------------------------------------- + rad_normal = sym.normal(self.radiation_tag, d) - flux_w = join_fields(self.c, w) + 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 - 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 ----------------------------------------------------- + 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[0] += self.source_f + + 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]) - # {{{ boundary conditions + def max_eigenvalue(self, t, fields=None, discr=None): + return abs(self.c) - # 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) +# }}} - # 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 - neu_bc = join_fields(neu_c, neu_u, -neu_v) +# {{{ variable-velocity - # Radiation - rad_normal = sym.make_normal(self.radiation_tag, d) +class VariableCoefficientWeakWaveOperator(HyperbolicOperator): + """This operator discretizes the wave equation + :math:`\\partial_t^2 u = c^2 \\Delta u`. - 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 + To be precise, we discretize the hyperbolic system - 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) - ) + .. math:: - # }}} + \partial_t u - c \\nabla \\cdot v = 0 - # {{{ diffusion ------------------------------------------------------- - from pytools.obj_array import with_object_array_or_scalar + \partial_t v - c \\nabla u = 0 - 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 + The sign of :math:`v` determines whether we discretize the forward or the + backward wave equation. - if with_sensor: - diffusion_coeff += sym.Field("sensor") + :math:`c` is assumed to be constant across all space. + """ - from grudge.second_order import SecondDerivativeTarget + 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) - # strong_form here allows the reuse the value of grad u. - grad_tgt = SecondDerivativeTarget( - self.ambient_dim, strong_form=True, - operand=arg) + self.c = c + self.ambient_dim = ambient_dim + self.source_f = source_f - self.diffusion_scheme.grad(grad_tgt, bc_getter=None, - dirichlet_tags=[], neumann_tags=[]) + self.sign = sym.If(sym.Comparison( + self.c, ">", 0), + np.int32(1), np.int32(-1)) - div_tgt = SecondDerivativeTarget( - self.ambient_dim, strong_form=False, - operand=diffusion_coeff*grad_tgt.minv_all) + self.dirichlet_tag = dirichlet_tag + self.neumann_tag = neumann_tag + self.radiation_tag = radiation_tag - self.diffusion_scheme.div(div_tgt, - bc_getter=None, - dirichlet_tags=[], neumann_tags=[]) + self.dirichlet_bc_f = dirichlet_bc_f - return div_tgt.minv_all - else: - return 0 + self.flux_type = flux_type - # }}} + def flux(self, w): + c = w[0] + u = w[1] + v = w[2:] + normal = sym.normal(w.dd, self.ambient_dim) - # entire operator ----------------------------------------------------- - nabla = sym.nabla(d) - flux_op = sym.get_flux_operator(self.flux()) + if self.flux_type == "central": + 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) - return ( - - join_fields( - - self.time_sign*self.c*np.dot(nabla, v) - make_diffusion(u) - + self.source, + elif self.flux_type == "upwind": + 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)) - -self.time_sign*self.c*(nabla*u) - with_object_array_or_scalar( - make_diffusion, v) + else: + raise ValueError("invalid flux type '%s'" % self.flux_type) + + def sym_operator(self): + d = self.ambient_dim + + 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: + # 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(dir_c, 2*dir_g - dir_u, dir_v) + else: + 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_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_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) + - sym.interp("vol", self.radiation_tag)(self.sign)*rad_u) + ), "rad_bc") + + # entire operator ----------------------------------------------------- + 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.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(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)) + )) + result[0] += self.source_f + + return result + def check_bc_coverage(self, mesh): from meshmode.mesh import check_bc_coverage check_bc_coverage(mesh, [ @@ -369,8 +458,7 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator): self.neumann_tag, self.radiation_tag]) - def max_eigenvalue_expr(self): - import grudge.symbolic as sym + def max_eigenvalue(self, t, fields=None, discr=None): return sym.NodalMax()(sym.CFunction("fabs")(self.c)) # }}}