hedge.pde

1 # -*- coding: utf8 -*- 2 """Canned operators for several PDEs, such as Maxwell's, heat, Poisson, etc. (deprecated module)""" 3 4 from __future__ import division 5 6 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" 7 8 __license__ = """ 9 This program is free software: you can redistribute it and/or modify 10 it under the terms of the GNU General Public License as published by 11 the Free Software Foundation, either version 3 of the License, or 12 (at your option) any later version. 13 14 This program is distributed in the hope that it will be useful, 15 but WITHOUT ANY WARRANTY; without even the implied warranty of 16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17 GNU General Public License for more details. 18 19 You should have received a copy of the GNU General Public License 20 along with this program. If not, see U{http://www.gnu.org/licenses/}. 21 """ 22 23 24 25 26 # THIS FILE HAS GROWN TOO LARGE AND IS BEING DEPRECATED. DON'T ADD STUFF HERE. 27 # ADD IT TO A NEW OR EXISTING FILE UNDER hedge.models INSTEAD. 28 29 30 31 32 33 import numpy 34 import numpy.linalg as la 35 import hedge.tools 36 import hedge.mesh 37 import hedge.data 38 from pytools import memoize_method 39 from hedge.models import Operator, TimeDependentOperator 40 41 42 43

44 -class StrongWaveOperator:

45 """This operator discretizes the Wave equation S{part}tt u = c^2 S{Delta} u. 46 47 To be precise, we discretize the hyperbolic system 48 49 - S{part}t u - c div v = 0 50 - S{part}t v - c grad u = 0 51 52 The sign of M{v} determines whether we discretize the forward or the 53 backward wave equation. 54 55 c is assumed to be constant across all space. 56 """ 57

58 - def __init__(self, c, dimensions, source_f=None, 59 flux_type="upwind", 60 dirichlet_tag=hedge.mesh.TAG_ALL, 61 neumann_tag=hedge.mesh.TAG_NONE, 62 radiation_tag=hedge.mesh.TAG_NONE):

63 assert isinstance(dimensions, int) 64 65 self.c = c 66 self.dimensions = dimensions 67 self.source_f = source_f 68 69 if self.c > 0: 70 self.sign = 1 71 else: 72 self.sign = -1 73 74 self.dirichlet_tag = dirichlet_tag 75 self.neumann_tag = neumann_tag 76 self.radiation_tag = radiation_tag 77 78 self.flux_type = flux_type

79

80 - def flux(self):

81 from hedge.flux import FluxVectorPlaceholder, make_normal 82 83 dim = self.dimensions 84 w = FluxVectorPlaceholder(1+dim) 85 u = w[0] 86 v = w[1:] 87 normal = make_normal(dim) 88 89 from hedge.tools import join_fields 90 flux_weak = join_fields( 91 numpy.dot(v.avg, normal), 92 u.avg * normal) 93 94 if self.flux_type == "central": 95 pass 96 elif self.flux_type == "upwind": 97 # see doc/notes/hedge-notes.tm 98 flux_weak -= self.sign*join_fields( 99 0.5*(u.int-u.ext), 100 0.5*(normal * numpy.dot(normal, v.int-v.ext))) 101 else: 102 raise ValueError, "invalid flux type '%s'" % self.flux_type 103 104 flux_strong = join_fields( 105 numpy.dot(v.int, normal), 106 u.int * normal) - flux_weak 107 108 return -self.c*flux_strong

109

110 - def op_template(self):

111 from hedge.optemplate import \ 112 make_vector_field, \ 113 pair_with_boundary, \ 114 get_flux_operator, \ 115 make_nabla, \ 116 InverseMassOperator, \ 117 BoundarizeOperator 118 119 d = self.dimensions 120 121 w = make_vector_field("w", d+1) 122 u = w[0] 123 v = w[1:] 124 125 # boundary conditions ------------------------------------------------- 126 from hedge.tools import join_fields 127 128 129 # dirichlet BC's ------------------------------------------------------ 130 dir_u = BoundarizeOperator(self.dirichlet_tag) * u 131 dir_v = BoundarizeOperator(self.dirichlet_tag) * v 132 dir_bc = join_fields(-dir_u, dir_v) 133 134 # neumann BC's -------------------------------------------------------- 135 neu_u = BoundarizeOperator(self.neumann_tag) * u 136 neu_v = BoundarizeOperator(self.neumann_tag) * v 137 neu_bc = join_fields(neu_u, -neu_v) 138 139 # radiation BC's ------------------------------------------------------ 140 from hedge.optemplate import make_normal 141 rad_normal = make_normal(self.radiation_tag, d) 142 143 rad_u = BoundarizeOperator(self.radiation_tag) * u 144 rad_v = BoundarizeOperator(self.radiation_tag) * v 145 146 rad_bc = join_fields( 147 0.5*(rad_u - self.sign*numpy.dot(rad_normal, rad_v)), 148 0.5*rad_normal*(numpy.dot(rad_normal, rad_v) - self.sign*rad_u) 149 ) 150 151 # entire operator ----------------------------------------------------- 152 nabla = make_nabla(d) 153 flux_op = get_flux_operator(self.flux()) 154 155 from hedge.tools import join_fields 156 return ( 157 - join_fields( 158 -self.c*numpy.dot(nabla, v), 159 -self.c*(nabla*u) 160 ) 161 + 162 InverseMassOperator() * ( 163 flux_op*w 164 + flux_op * pair_with_boundary(w, dir_bc, self.dirichlet_tag) 165 + flux_op * pair_with_boundary(w, neu_bc, self.neumann_tag) 166 + flux_op * pair_with_boundary(w, rad_bc, self.radiation_tag) 167 ))

168 169

170 - def bind(self, discr):

171 from hedge.mesh import check_bc_coverage 172 check_bc_coverage(discr.mesh, [ 173 self.dirichlet_tag, 174 self.neumann_tag, 175 self.radiation_tag]) 176 177 compiled_op_template = discr.compile(self.op_template()) 178 179 def rhs(t, w): 180 rhs = compiled_op_template(w=w) 181 182 if self.source_f is not None: 183 rhs[0] += self.source_f(t) 184 185 return rhs

186 187 return rhs

188

189 - def max_eigenvalue(self):

190 return abs(self.c)

191 192 193 194

195 -class VariableVelocityStrongWaveOperator:

196 """This operator discretizes the Wave equation S{part}tt u = c^2 S{Delta} u. 197 198 To be precise, we discretize the hyperbolic system 199 200 - S{part}t u - c div v = 0 201 - S{part}t v - c grad u = 0 202 """ 203

204 - def __init__(self, c, dimensions, source=None, 205 flux_type="upwind", 206 dirichlet_tag=hedge.mesh.TAG_ALL, 207 neumann_tag=hedge.mesh.TAG_NONE, 208 radiation_tag=hedge.mesh.TAG_NONE, 209 time_sign=1):

210 """`c` is assumed to be positive and conforms to the 211 `hedge.data.ITimeDependentGivenFunction` interface. 212 213 `source` also conforms to the 214 `hedge.data.ITimeDependentGivenFunction` interface. 215 """ 216 assert isinstance(dimensions, int) 217 218 self.c = c 219 self.time_sign = time_sign 220 self.dimensions = dimensions 221 self.source = source 222 223 self.dirichlet_tag = dirichlet_tag 224 self.neumann_tag = neumann_tag 225 self.radiation_tag = radiation_tag 226 227 self.flux_type = flux_type

228

229 - def flux(self):

230 from hedge.flux import FluxVectorPlaceholder, make_normal 231 232 dim = self.dimensions 233 w = FluxVectorPlaceholder(2+dim) 234 c = w[0] 235 u = w[1] 236 v = w[2:] 237 normal = make_normal(dim) 238 239 from hedge.tools import join_fields 240 flux = self.time_sign*1/2*join_fields( 241 c.ext * numpy.dot(v.ext, normal) 242 - c.int * numpy.dot(v.int, normal), 243 normal*(c.ext*u.ext - c.int*u.int)) 244 245 if self.flux_type == "central": 246 pass 247 elif self.flux_type == "upwind": 248 flux += join_fields( 249 c.ext*u.ext - c.int*u.int, 250 c.ext*normal*numpy.dot(normal, v.ext) 251 - c.int*normal*numpy.dot(normal, v.int) 252 ) 253 else: 254 raise ValueError, "invalid flux type '%s'" % self.flux_type 255 256 return flux

257

258 - def op_template(self):

259 from hedge.optemplate import \ 260 Field, \ 261 make_vector_field, \ 262 pair_with_boundary, \ 263 get_flux_operator, \ 264 make_nabla, \ 265 InverseMassOperator, \ 266 BoundarizeOperator 267 268 d = self.dimensions 269 270 w = make_vector_field("w", d+1) 271 u = w[0] 272 v = w[1:] 273 274 from hedge.tools import join_fields 275 c = Field("c") 276 flux_w = join_fields(c, w) 277 278 # boundary conditions ------------------------------------------------- 279 from hedge.flux import make_normal 280 normal = make_normal(d) 281 282 from hedge.tools import join_fields 283 # dirichlet BC's ------------------------------------------------------ 284 dir_c = BoundarizeOperator(self.dirichlet_tag) * c 285 dir_u = BoundarizeOperator(self.dirichlet_tag) * u 286 dir_v = BoundarizeOperator(self.dirichlet_tag) * v 287 288 dir_bc = join_fields(dir_c, -dir_u, dir_v) 289 290 # neumann BC's -------------------------------------------------------- 291 neu_c = BoundarizeOperator(self.neumann_tag) * c 292 neu_u = BoundarizeOperator(self.neumann_tag) * u 293 neu_v = BoundarizeOperator(self.neumann_tag) * v 294 295 neu_bc = join_fields(neu_c, neu_u, -neu_v) 296 297 # radiation BC's ------------------------------------------------------ 298 from hedge.optemplate import make_normal 299 rad_normal = make_normal(self.radiation_tag, d) 300 301 rad_c = BoundarizeOperator(self.radiation_tag) * c 302 rad_u = BoundarizeOperator(self.radiation_tag) * u 303 rad_v = BoundarizeOperator(self.radiation_tag) * v 304 305 rad_bc = join_fields( 306 rad_c, 307 0.5*(rad_u - self.time_sign*numpy.dot(rad_normal, rad_v)), 308 0.5*rad_normal*(numpy.dot(rad_normal, rad_v) - self.time_sign*rad_u) 309 ) 310 311 # entire operator ----------------------------------------------------- 312 nabla = make_nabla(d) 313 flux_op = get_flux_operator(self.flux()) 314 315 return ( 316 - join_fields( 317 -numpy.dot(nabla, self.time_sign*c*v), 318 -(nabla*(self.time_sign*c*u)) 319 ) 320 + 321 InverseMassOperator() * ( 322 flux_op*flux_w 323 + flux_op * pair_with_boundary(flux_w, dir_bc, self.dirichlet_tag) 324 + flux_op * pair_with_boundary(flux_w, neu_bc, self.neumann_tag) 325 + flux_op * pair_with_boundary(flux_w, rad_bc, self.radiation_tag) 326 ))

327 328

329 - def bind(self, discr):

330 from hedge.mesh import check_bc_coverage 331 check_bc_coverage(discr.mesh, [ 332 self.dirichlet_tag, 333 self.neumann_tag, 334 self.radiation_tag]) 335 336 compiled_op_template = discr.compile(self.op_template()) 337 338 def rhs(t, w): 339 rhs = compiled_op_template(w=w, 340 c=self.c.volume_interpolant(t, discr)) 341 342 if self.source is not None: 343 rhs[0] += self.source.volume_interpolant(t, discr) 344 345 return rhs

346 347 return rhs

348 349 #def max_eigenvalue(self): 350 #return abs(self.c) 351 352 353 354

355 -class EulerOperator(TimeDependentOperator):

356 """An nD Euler operator. 357 358 Field order is [rho E rho_u_x rho_u_y ...]. 359 """

360 - def __init__(self, dimensions, gamma, bc):

361 self.dimensions = dimensions 362 self.gamma = gamma 363 self.bc = bc

364

365 - def rho(self, q):

366 return q[0]

367

368 - def e(self, q):

369 return q[1]

370

371 - def rho_u(self, q):

372 return q[2:2+self.dimensions]

373

374 - def u(self, q):

375 from hedge.tools import make_obj_array 376 return make_obj_array([ 377 rho_u_i/self.rho(q) 378 for rho_u_i in self.rho_u(q)])

379

380 - def op_template(self):

381 from hedge.optemplate import make_vector_field, \ 382 make_common_subexpression as cse 383 384 def u(q): 385 return cse(self.u(q))

386 387 def p(q): 388 return cse((self.gamma-1)*(self.e(q) - 0.5*numpy.dot(self.rho_u(q), u(q))))

389 390 def flux(q): 391 from pytools import delta 392 from hedge.tools import make_obj_array, join_fields 393 return [ # one entry for each flux direction 394 cse(join_fields( 395 # flux rho 396 self.rho_u(q)[i], 397 398 # flux E 399 cse(self.e(q)+p(q))*u(q)[i], 400 401 # flux rho_u 402 make_obj_array([ 403 self.rho_u(q)[i]*self.u(q)[j] + delta(i,j) * p(q) 404 for j in range(self.dimensions) 405 ]) 406 )) 407 for i in range(self.dimensions)] 408 409 from hedge.optemplate import make_nabla, InverseMassOperator, \ 410 ElementwiseMaxOperator 411 412 from pymbolic import var 413 sqrt = var("sqrt") 414 415 state = make_vector_field("q", self.dimensions+2) 416 bc_state = make_vector_field("bc_q", self.dimensions+2) 417 418 c = cse(sqrt(self.gamma*p(state)/self.rho(state))) 419 420 speed = sqrt(numpy.dot(u(state), u(state))) + c 421 422 from hedge.tools import make_lax_friedrichs_flux, join_fields 423 from hedge.mesh import TAG_ALL 424 return join_fields( 425 (- numpy.dot(make_nabla(self.dimensions), flux(state)) 426 + InverseMassOperator()*make_lax_friedrichs_flux( 427 wave_speed=ElementwiseMaxOperator()*c, 428 state=state, flux_func=flux, 429 bdry_tags_and_states=[ 430 (TAG_ALL, bc_state) 431 ], 432 strong=True 433 )), 434 speed) 435

436 - def bind(self, discr):

437 from hedge.mesh import TAG_ALL 438 bound_op = discr.compile(self.op_template()) 439 440 def wrap(t, q): 441 opt_result = bound_op( 442 q=q, 443 bc_q=self.bc.boundary_interpolant(t, discr, TAG_ALL)) 444 max_speed = opt_result[-1] 445 ode_rhs = opt_result[:-1] 446 return ode_rhs, numpy.max(max_speed)

447 448 return wrap 449

Source Code for Module hedge.pde