Differences between revisions 17 and 18

# Sample Solutions

## Dependency Mapper

```   1 import pymbolic as p
2
3 expr = p.parse("x**r*y+z+5")
4
5 print expr
6 print repr(expr)
7
8 from pymbolic.mapper.dependency import DependencyMapper
9
10 print "Dependency Mapper:", DependencyMapper()(expr)
11
12 import pymbolic.primitives as primitives
13
14 def mapper(expression):
15     variables = set()
16     if isinstance(expression, primitives.Variable):
18         return variables
19     elif primitives.is_constant(expression):
20         return variables
21     elif isinstance(expression, primitives.Power):
22         variables = variables.union(mapper(expression.base))
23         variables = variables.union(mapper(expression.exponent))
24         return variables
25     for i in expression.children:
26         variables = variables.union(mapper(i))
27     return variables
28
29 print "--------------------"
30 print mapper(expr)
31 print "--------------------"
```

## Evaluator

### visitor pattern approach

```   1 list = dict({'x': 2, 'y': 4, 'z': 5})
2 print "Dict.:", list
3 expr = p.parse("x**(2*z)+y/2+z")
4 print "Expression:",expr
5
6 import pymbolic.primitives as primitives
7
8 class EvaluationMapper_new:
9     def __init__(self, context={}):
10         self.context = context
11         self.common_subexp_cache = {}
12
13     def __call__(self, expr, *args):
14         try:
15             method = expr.get_mapper_method(self)
16         except AttributeError:
17             if isinstance(expr, primitives.Expression):
18                 return self.handle_unsupported_expression(expr, *args)
19             else:
20                 return self.map_foreign(expr, *args)
21         else:
22             return method(expr, *args)
23
24     def rec(self, expr, *args):
25         try:
26             method = expr.get_mapper_method(self)
27         except AttributeError:
28             if isinstance(expr, primitives.Expression):
29                 return self.handle_unsupported_expression(expr, *args)
30             else:
31                 return self.map_foreign(expr, *args)
32         else:
33             return method(expr, *args)
34
35     # ======================================
36     # Methods of Evaluation:
37     # ======================================
38
39     def map_constant(self, expr):
40         return expr
41
42     def map_variable(self, expr):
43         try:
44             # Schluesselstelle: Wenn eine Variable gefunden wurde, wird
45             # deren Name in mittels dictionairy mit einer Zahl verknuepft.
46             # Diese geht dann in den Summator - map_sum - oder
47             # Produktbilder - map_product - zurueck.
48             return self.context[expr.name]
49         except KeyError:
50             raise UnknownVariableError, expr.name
51
52     def map_sum(self, expr):
53         return sum(self.rec(child) for child in expr.children)
54
55     def map_product(self, expr):
56         import operator
57         return reduce(operator.mul, (self.rec(child) for child in expr.children), 1)
58         #product_intern(self.rec(child) for child in expr.children)
59
60     def map_quotient(self, expr):
61         return self.rec(expr.numerator) / self.rec(expr.denominator)
62
63     def map_power(self, expr):
64         return self.rec(expr.base) ** self.rec(expr.exponent)
65
66     def map_foreign(self, expr, *args):
67         if isinstance(expr, primitives.VALID_CONSTANT_CLASSES):
68             return self.map_constant(expr, *args)
69         elif isinstance(expr, list):
70             return self.map_list(expr, *args)
71         elif is_numpy_array(expr):
72             return self.map_numpy_array(expr, *args)
73         else:
74             raise ValueError, "%s encountered invalid foreign object: %s" % (
75                     self.__class__, repr(expr))
76
77 def evaluate_new(expression, context={}):
78     return EvaluationMapper_new(context)(expression)
```

### isinstance approach

```   1 list = dict({'x': 2, 'y': 4, 'z': 5})
2 print "Dict.:", list
3 expr = p.parse("x**(2*z)+y/2+z")
4 print "Function:",expr
5
6 import pymbolic.primitives as primitives
7
8 class EvaluationMapper_new:
9     def __init__(self, context={}):
10         self.context = context
11         self.common_subexp_cache = {}
12
13     def __call__(self, expr, *args):
14             if isinstance(expr, primitives.Variable):
15                 return self.map_variable(expr, *args)
16             elif isinstance(expr, primitives.Sum):
17                 return self.map_sum(expr, *args)
18             elif isinstance(expr, primitives.Product):
19                 return self.map_product(expr, *args)
20             elif isinstance(expr, primitives.Quotient):
21                 return self.map_quotient(expr, *args)
22             elif isinstance(expr, primitives.Power):
23                 return self.map_power(expr, *args)
24             elif isinstance(expr, primitives.Expression):
25                 return self.handle_unsupported_expression(expr, *args)
26             else:
27                 return self.map_foreign(expr, *args)
28
29     def rec(self, expr, *args):
30             if isinstance(expr, primitives.Variable):
31                 return self.map_variable(expr, *args)
32             elif isinstance(expr, primitives.Sum):
33                 return self.map_sum(expr, *args)
34             elif isinstance(expr, primitives.Product):
35                 return self.map_product(expr, *args)
36             elif isinstance(expr, primitives.Quotient):
37                 return self.map_quotient(expr, *args)
38             elif isinstance(expr, primitives.Power):
39                 return self.map_power(expr, *args)
40             elif isinstance(expr, primitives.Expression):
41                 return self.handle_unsupported_expression(expr, *args)
42             else:
43                 return self.map_foreign(expr, *args)
44
45
46     def map_constant(self, expr):
47         return expr
48
49     def map_variable(self, expr):
50         try:
51             return self.context[expr.name]
52         except KeyError:
53             raise UnknownVariableError, expr.name
54
55     def map_sum(self, expr):
56         return sum(self.rec(child) for child in expr.children)
57
58     def map_product(self, expr):
59         import operator
60         return reduce(operator.mul, (self.rec(child) for child in expr.children), 1)
61         #product_intern(self.rec(child) for child in expr.children)
62
63     def map_quotient(self, expr):
64         return self.rec(expr.numerator) / self.rec(expr.denominator)
65
66     def map_power(self, expr):
67         return self.rec(expr.base) ** self.rec(expr.exponent)
68
69     def map_foreign(self, expr, *args):
70         if isinstance(expr, primitives.VALID_CONSTANT_CLASSES):
71             return self.map_constant(expr, *args)
72         elif isinstance(expr, list):
73             return self.map_list(expr, *args)
74         elif is_numpy_array(expr):
75             return self.map_numpy_array(expr, *args)
76         else:
77             raise ValueError, "%s encountered invalid foreign object: %s" % (
78                     self.__class__, repr(expr))
79
80 def evaluate_new(expression, context={}):
81     return EvaluationMapper_new(context)(expression)
82
83 print evaluate_new(expr,list)
```

Final space-time dependent solution:

## Wave Equation

/hedge/examples/wave/wave-min.py:

```   1 # Hedge - the Hybrid'n'Easy DG Environment
2 # Copyright (C) 2009 Andreas Stock
3 #
4 # This program is free software: you can redistribute it and/or modify
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
16
17
18
19
20 from __future__ import division
21 import numpy
22 import numpy.linalg as la
23
24
25
26
27 def main() :
28     from math import sin, exp, sqrt
29
30     from hedge.mesh import make_disk_mesh
31     mesh = make_disk_mesh(r=0.5,faces=100,max_area=0.008)
32
33     from hedge.backends.jit import Discretization
34     discr = Discretization(mesh, order=3)
35
36     from hedge.visualization import VtkVisualizer
37     vis = VtkVisualizer(discr, None, "fld")
38
39     def source_u(x, el):
40         return exp(-numpy.dot(x, x)*128)
41
42     source_u_vec = discr.interpolate_volume_function(source_u)
43
44     def source_vec_getter(t):
45         from math import sin
46         return source_u_vec*sin(10*t)
47
48     def bound_u(pt, el):
49         x, y = pt
50         if x > 0:
51             return 1
52         else:
53             return 0
54
55     # Building the boundary-vector - a line-vector - which contains less -
56     # but the requested - information than the volume-vector of the entire
57     # discretization space.
58     bound_u_vec = discr.interpolate_boundary_function(bound_u)
59
60     def bound_vec_getter(t):
61         if t > 0.5:
62             return bound_u_vec * 0
63         else:
64             return bound_u_vec
65
66     from hedge.pde import StrongWaveOperator
67     from hedge.mesh import TAG_ALL, TAG_NONE
68
69     op = StrongWaveOperator(1, discr.dimensions,
70             source_vec_getter,
71             flux_type="upwind",
72             dirichlet_tag=TAG_ALL,
73             dirichlet_bc_f=bound_vec_getter,
74             neumann_tag=TAG_NONE,
76
77     from hedge.tools import join_fields
78     fields = join_fields(discr.volume_zeros(),
79             [discr.volume_zeros() for i in range(discr.dimensions)])
80
81     dt = discr.dt_factor(op.max_eigenvalue()) / 10
82     nsteps = int(1/dt)
83
84     # timestep loop -----------------------------------------------------------
87
88     rhs = op.bind(discr)
89     for step in range(nsteps):
90         t = step*dt
91
92         if step % 10 == 0:
93             print step, t
94             visf = vis.make_file("fld-%04d" % step)
95
97                     [ ("u", fields[0]), ("v", fields[1:]), ],
98                     time=t, step=step)
99             visf.close()
100
101         fields = stepper(fields, t, dt, rhs)
102
103     discr.close()
104     vis.close()
105
106
107
108
109 if __name__ == "__main__":
110     main()
```

/hedge/hedge/pde.py:

```   1 class StrongWaveOperator:
2     """This operator discretizes the Wave equation S{part}tt u = c^2 S{Delta} u.
3
4     To be precise, we discretize the hyperbolic system
5
6       * S{part}t u - c div v = 0
7       * S{part}t v - c grad u = 0
8
9     Note that this is not unique--we could also choose a different sign for M{v}.
10     """
11
12     def __init__(self, c, dimensions, source_f=None,
13             flux_type="upwind",
14             dirichlet_tag=hedge.mesh.TAG_ALL,
15             dirichlet_bc_f=None,
16             neumann_tag=hedge.mesh.TAG_NONE,
18         assert isinstance(dimensions, int)
19
20         self.c = c
21         self.dimensions = dimensions
22         self.source_f = source_f
23
24         if self.c > 0:
25             self.sign = 1
26         else:
27             self.sign = -1
28
29         self.dirichlet_tag = dirichlet_tag
30         self.neumann_tag = neumann_tag
32
33         self.dirichlet_bc_f = dirichlet_bc_f
34
35         self.flux_type = flux_type
36
37     def flux(self):
38         from hedge.flux import FluxVectorPlaceholder, make_normal
39
40         dim = self.dimensions
41         w = FluxVectorPlaceholder(1+dim)
42         u = w[0]
43         v = w[1:]
44         normal = make_normal(dim)
45
46         from hedge.tools import join_fields
47         flux_weak = join_fields(
48                 numpy.dot(v.avg, normal),
49                 u.avg * normal)
50
51         if self.flux_type == "central":
52             pass
53         elif self.flux_type == "upwind":
54             # see doc/notes/hedge-notes.tm
55             flux_weak -= self.sign*join_fields(
56                     0.5*(u.int-u.ext),
57                     0.5*(normal * numpy.dot(normal, v.int-v.ext)))
58         else:
59             raise ValueError, "invalid flux type '%s'" % self.flux_type
60
61         flux_strong = join_fields(
62                 numpy.dot(v.int, normal),
63                 u.int * normal) - flux_weak
64
65         return -self.c*flux_strong
66
67     def op_template(self):
68         from hedge.optemplate import \
69                 Field, \
70                 make_vector_field, \
71                 pair_with_boundary, \
72                 get_flux_operator, \
73                 make_nabla, \
74                 InverseMassOperator
75
76         d = self.dimensions
77
78         w = make_vector_field("w", d+1)
79         u = w[0]
80         v = w[1:]
81
82         # boundary conditions -------------------------------------------------
83         from hedge.flux import make_normal
84         normal = make_normal(d)
85
86         from hedge.tools import join_fields
87
88         if self.dirichlet_bc_f is None:
89             dir_bc = join_fields(-u, v)
90         else:
91             dir_bc = join_fields(-Field("dir_bc_u"), v)
92
93         neu_bc = join_fields(u, -v)
95                 0.5*(u - self.sign*numpy.dot(normal, v)),
96                 0.5*normal*(numpy.dot(normal, v) - self.sign*u)
97                 )
98
99         # entire operator -----------------------------------------------------
100         nabla = make_nabla(d)
101         flux_op = get_flux_operator(self.flux())
102
103         from hedge.tools import join_fields
104         return (
105                 - join_fields(
106                     -self.c*numpy.dot(nabla, v),
107                     -self.c*(nabla*u)
108                     )
109                 +
110                 InverseMassOperator() * (
111                     flux_op*w
112                     + flux_op * pair_with_boundary(w, dir_bc, self.dirichlet_tag)
113                     + flux_op * pair_with_boundary(w, neu_bc, self.neumann_tag)
115                     ))
116
117
118     def bind(self, discr):
119         from hedge.mesh import check_bc_coverage
120         check_bc_coverage(discr.mesh, [
121             self.dirichlet_tag,
122             self.neumann_tag,
124
125         compiled_op_template = discr.compile(self.op_template())
126
127         def rhs(t, w):
128             from hedge.tools import join_fields, ptwise_dot
129
130             if self.dirichlet_bc_f is None:
131                 rhs = compiled_op_template(w=w)
132             else:
133                 rhs = compiled_op_template(w=w, dir_bc_u=self.dirichlet_bc_f(t))
134
135             if self.source_f is not None:
136                 rhs[0] += self.source_f(t)
137
138             return rhs
139
140         return rhs
141
142     def max_eigenvalue(self):
143         return abs(self.c)
```

Hedge/Exercises/SampleSolutions (last edited 2013-01-16 23:15:42 by AndreasKloeckner)