[frames] | no frames]

# Source Code for Module 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
7
9  This program is free software: you can redistribute it and/or modify
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,
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
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
140          from hedge.optemplate import make_normal
142
145
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)
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,
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,
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
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
298          from hedge.optemplate import make_normal
300
304
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)
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,
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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Sat Aug 29 14:33:47 2009 http://epydoc.sourceforge.net