From 99326552c0d1b516c073969a13c86304a6e12cc3 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sat, 1 Apr 2017 20:37:58 -0500 Subject: [PATCH] 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