From 3d322b523546d7e9b59611d13fd67cc91c0b5971 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 16 May 2017 09:47:59 -0500 Subject: [PATCH 1/2] Negative c values attempt --- examples/wave/simple-var-prop.py | 129 +++++++++++++++++++++++++ examples/wave/var-propagation-speed.py | 12 +-- grudge/models/wave.py | 17 ++-- 3 files changed, 145 insertions(+), 13 deletions(-) create mode 100644 examples/wave/simple-var-prop.py diff --git a/examples/wave/simple-var-prop.py b/examples/wave/simple-var-prop.py new file mode 100644 index 00000000..e2fd9795 --- /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 6e1a7ea7..3d3be21d 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) @@ -53,7 +53,7 @@ def main(write_output=True, order=4): 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)) + np.float32(0.1), np.float32(0.2)) from grudge.models.wave import VariableCoefficientWeakWaveOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE @@ -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="upwind") + flux_type="central") queue = cl.CommandQueue(discr.cl_context) from pytools.obj_array import join_fields @@ -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 19229abe..bf07fbd3 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,16 +367,19 @@ 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))) if self.flux_type == "central": - return -0.5 * join_fields( + return 0.5 * sign * join_fields( np.dot(v.int*c.int + v.ext*c.ext, normal), (u.int * c.int + u.ext*c.ext) * normal) elif self.flux_type == "upwind": - return -0.5 * join_fields(np.dot(normal, c.ext * v.ext + c.int * v.int) + return 0.5 * join_fields(sign * 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) + sign * normal * (np.dot(normal, c.ext * v.ext - c.int * v.int) + c.ext*u.ext + c.int * u.int)) else: @@ -435,8 +438,8 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): 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) + self.sign * self.c*np.dot(sym.stiffness_t(self.ambient_dim), v), + self.sign * self.c*(sym.stiffness_t(self.ambient_dim)*u) ) -- GitLab From f798d2aad5b63588af1576234f7541aa5dc001d5 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sat, 27 May 2017 15:17:00 -0500 Subject: [PATCH 2/2] An attempt was made at making positive coefficient waves work --- examples/wave/var-propagation-speed.py | 4 ++-- grudge/models/wave.py | 14 ++++++++------ 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index 3d3be21d..8f796138 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -53,7 +53,7 @@ def main(write_output=True, order=4): 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)) + np.float32(-0.1), np.float32(-0.2)) from grudge.models.wave import VariableCoefficientWeakWaveOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE @@ -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 bf07fbd3..28523ee8 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -370,16 +370,18 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): 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 * sign * 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) elif self.flux_type == "upwind": - return 0.5 * join_fields(sign * np.dot(normal, c.ext * v.ext + c.int * v.int) - + c.ext*u.ext - c.int * u.int, - sign * normal * (np.dot(normal, c.ext * v.ext - c.int * v.int) + return -0.5 * join_fields(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: @@ -438,8 +440,8 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): result = sym.InverseMassOperator()( join_fields( - self.sign * self.c*np.dot(sym.stiffness_t(self.ambient_dim), v), - self.sign * self.c*(sym.stiffness_t(self.ambient_dim)*u) + -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v), + -self.c*(sym.stiffness_t(self.ambient_dim)*u) ) -- GitLab