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):
  17         variables.add(expression)
  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 "--------------------"

   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)

   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)

   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
   5 # it under the terms of the GNU General Public License as published by
   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,
  75             radiation_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 -----------------------------------------------------------
  85     from hedge.timestep import AdamsBashforthTimeStepper
  86     stepper = AdamsBashforthTimeStepper(3,None)
  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 
  96             vis.add_data(visf,
  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()

   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,
  17             radiation_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
  31         self.radiation_tag = radiation_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)
  94         rad_bc = join_fields(
  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)
 114                     + flux_op * pair_with_boundary(w, rad_bc, self.radiation_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,
 123             self.radiation_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)