diff --git a/examples/wave/simple-var-prop.py b/examples/wave/simple-var-prop.py new file mode 100644 index 0000000000000000000000000000000000000000..e2fd97958379b118ab22efc12c66754ee5cdc7af --- /dev/null +++ b/examples/wave/simple-var-prop.py @@ -0,0 +1,129 @@ +"""Minimal example of a grudge driver.""" + +from __future__ import division, print_function + +__copyright__ = "Copyright (C) 2015 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +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 +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 + + 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 VariableCoefficientWeakWaveOperator + from meshmode.mesh import BTAG_ALL, BTAG_NONE + op = VariableCoefficientWeakWaveOperator(c, + discr.dim, + dirichlet_tag=BTAG_NONE, + neumann_tag=BTAG_NONE, + radiation_tag=BTAG_ALL, + flux_type="upwind") + + queue = cl.CommandQueue(discr.cl_context) + from pytools.obj_array import join_fields + def f(x): + return sym.exp( + -np.dot(sym_source_center_dist, sym_source_center_dist) + / source_width**2) + + fields = join_fields(bind(discr, f(sym.nodes(dims)))(queue, t=0), + [discr.zeros(queue) for i in range(dims)]) + + + + + op.check_bc_coverage(mesh) + + c_eval = bind(discr, c)(queue) + + bound_op = bind(discr, op.sym_operator()) + + 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:]), + ("c", c_eval), + ]) + t_last_step = time() + + +if __name__ == "__main__": + main() diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index 6e1a7ea7368198871e2852e1c7ff43c51ddb75c0..8f7961383210910ae42285a8539e0e4bd58696bc 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -35,12 +35,12 @@ def main(write_output=True, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - dims = 3 + 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=(20,)*dims) + n=(8,)*dims) discr = Discretization(cl_ctx, mesh, order=order) @@ -84,9 +84,9 @@ def main(write_output=True, order=4): return bound_op(queue, t=t, w=w) if mesh.dim == 2: - dt = 0.04 * 0.3 + dt = 0.04 elif mesh.dim == 3: - dt = 0.02 * 0.1 + dt = 0.02 dt_stepper = set_up_rk4("w", dt, fields, rhs) diff --git a/grudge/models/wave.py b/grudge/models/wave.py index 19229abe81cb3ed567cb4346174f3f2eff8e577c..28523ee8000642257a12f4b704173ddb895582d3 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -350,9 +350,9 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): self.ambient_dim = ambient_dim self.source_f = source_f - self.sign = sym.If(sym.Comparison( + self.sign = sym.cse(sym.If(sym.Comparison( self.c, ">", 0), - np.int32(1), np.int32(-1)) + np.int32(1), np.int32(-1))) self.dirichlet_tag = dirichlet_tag self.neumann_tag = neumann_tag @@ -367,6 +367,11 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): u = w[1] v = w[2:] normal = sym.normal(w.dd, self.ambient_dim) + sign = sym.cse(sym.If(sym.Comparison( + c.int, ">", 0), + np.int32(1), np.int32(-1))) + cabsi = sym.CFunction("fabs")(c.int) + cabse = sym.CFunction("fabs")(c.ext) if self.flux_type == "central": return -0.5 * join_fields( @@ -375,8 +380,8 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): 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) + + cabse*u.ext - cabsi * u.int, + normal * (np.dot(normal, cabse * v.ext - cabsi * v.int) + c.ext*u.ext + c.int * u.int)) else: