Differences between revisions 17 and 18
Revision 17 as of 2012-10-19 16:21:27
Size: 0
Editor: 46
Comment:
Revision 18 as of 2013-01-16 23:15:42
Size: 14945
Comment:
Deletions are marked like this. Additions are marked like this.
Line 1: Line 1:
## page was renamed from SciComp/Hedge/Exercises/SampleSolutions
## page was renamed from SciComp/HedgeExercises/SampleSolutions
= Sample Solutions =
<<TableOfContents>>

== Dependency Mapper ==
{{{#!python
import pymbolic as p

expr = p.parse("x**r*y+z+5")

print expr
print repr(expr)

from pymbolic.mapper.dependency import DependencyMapper

print "Dependency Mapper:", DependencyMapper()(expr)

import pymbolic.primitives as primitives

def mapper(expression):
    variables = set()
    if isinstance(expression, primitives.Variable):
        variables.add(expression)
        return variables
    elif primitives.is_constant(expression):
        return variables
    elif isinstance(expression, primitives.Power):
        variables = variables.union(mapper(expression.base))
        variables = variables.union(mapper(expression.exponent))
        return variables
    for i in expression.children:
        variables = variables.union(mapper(i))
    return variables

print "--------------------"
print mapper(expr)
print "--------------------"
}}}

== Evaluator ==
=== visitor pattern approach ===
{{{#!python
list = dict({'x': 2, 'y': 4, 'z': 5})
print "Dict.:", list
expr = p.parse("x**(2*z)+y/2+z")
print "Expression:",expr

import pymbolic.primitives as primitives

class EvaluationMapper_new:
    def __init__(self, context={}):
        self.context = context
        self.common_subexp_cache = {}

    def __call__(self, expr, *args):
        try:
            method = expr.get_mapper_method(self)
        except AttributeError:
            if isinstance(expr, primitives.Expression):
                return self.handle_unsupported_expression(expr, *args)
            else:
                return self.map_foreign(expr, *args)
        else:
            return method(expr, *args)

    def rec(self, expr, *args):
        try:
            method = expr.get_mapper_method(self)
        except AttributeError:
            if isinstance(expr, primitives.Expression):
                return self.handle_unsupported_expression(expr, *args)
            else:
                return self.map_foreign(expr, *args)
        else:
            return method(expr, *args)
    
    # ======================================
    # Methods of Evaluation:
    # ======================================

    def map_constant(self, expr):
        return expr

    def map_variable(self, expr):
        try:
     # Schluesselstelle: Wenn eine Variable gefunden wurde, wird
     # deren Name in mittels dictionairy mit einer Zahl verknuepft.
     # Diese geht dann in den Summator - map_sum - oder
     # Produktbilder - map_product - zurueck.
            return self.context[expr.name]
        except KeyError:
            raise UnknownVariableError, expr.name

    def map_sum(self, expr):
        return sum(self.rec(child) for child in expr.children)

    def map_product(self, expr):
 import operator
        return reduce(operator.mul, (self.rec(child) for child in expr.children), 1)
        #product_intern(self.rec(child) for child in expr.children)

    def map_quotient(self, expr):
        return self.rec(expr.numerator) / self.rec(expr.denominator)

    def map_power(self, expr):
        return self.rec(expr.base) ** self.rec(expr.exponent)
    
    def map_foreign(self, expr, *args):
        if isinstance(expr, primitives.VALID_CONSTANT_CLASSES):
            return self.map_constant(expr, *args)
        elif isinstance(expr, list):
            return self.map_list(expr, *args)
        elif is_numpy_array(expr):
            return self.map_numpy_array(expr, *args)
        else:
            raise ValueError, "%s encountered invalid foreign object: %s" % (
                    self.__class__, repr(expr))

def evaluate_new(expression, context={}):
    return EvaluationMapper_new(context)(expression)
}}}

=== isinstance approach ===
{{{#!python
list = dict({'x': 2, 'y': 4, 'z': 5})
print "Dict.:", list
expr = p.parse("x**(2*z)+y/2+z")
print "Function:",expr

import pymbolic.primitives as primitives

class EvaluationMapper_new:
    def __init__(self, context={}):
        self.context = context
        self.common_subexp_cache = {}

    def __call__(self, expr, *args):
     if isinstance(expr, primitives.Variable):
                return self.map_variable(expr, *args)
            elif isinstance(expr, primitives.Sum):
         return self.map_sum(expr, *args)
            elif isinstance(expr, primitives.Product):
         return self.map_product(expr, *args)
            elif isinstance(expr, primitives.Quotient):
         return self.map_quotient(expr, *args)
            elif isinstance(expr, primitives.Power):
         return self.map_power(expr, *args)
            elif isinstance(expr, primitives.Expression):
                return self.handle_unsupported_expression(expr, *args)
            else:
                return self.map_foreign(expr, *args)

    def rec(self, expr, *args):
     if isinstance(expr, primitives.Variable):
                return self.map_variable(expr, *args)
            elif isinstance(expr, primitives.Sum):
         return self.map_sum(expr, *args)
            elif isinstance(expr, primitives.Product):
         return self.map_product(expr, *args)
            elif isinstance(expr, primitives.Quotient):
         return self.map_quotient(expr, *args)
            elif isinstance(expr, primitives.Power):
         return self.map_power(expr, *args)
            elif isinstance(expr, primitives.Expression):
                return self.handle_unsupported_expression(expr, *args)
            else:
                return self.map_foreign(expr, *args)


    def map_constant(self, expr):
        return expr

    def map_variable(self, expr):
        try:
            return self.context[expr.name]
        except KeyError:
            raise UnknownVariableError, expr.name

    def map_sum(self, expr):
        return sum(self.rec(child) for child in expr.children)

    def map_product(self, expr):
 import operator
        return reduce(operator.mul, (self.rec(child) for child in expr.children), 1)
        #product_intern(self.rec(child) for child in expr.children)

    def map_quotient(self, expr):
        return self.rec(expr.numerator) / self.rec(expr.denominator)

    def map_power(self, expr):
        return self.rec(expr.base) ** self.rec(expr.exponent)
    
    def map_foreign(self, expr, *args):
        if isinstance(expr, primitives.VALID_CONSTANT_CLASSES):
            return self.map_constant(expr, *args)
        elif isinstance(expr, list):
            return self.map_list(expr, *args)
        elif is_numpy_array(expr):
            return self.map_numpy_array(expr, *args)
        else:
            raise ValueError, "%s encountered invalid foreign object: %s" % (
                    self.__class__, repr(expr))

def evaluate_new(expression, context={}):
    return EvaluationMapper_new(context)(expression)

print evaluate_new(expr,list)
}}}

== Advection Equation ==

Final space-time dependent solution:


== Wave Equation ==
`/hedge/examples/wave/wave-min.py`:
{{{#!python
# Hedge - the Hybrid'n'Easy DG Environment
# Copyright (C) 2009 Andreas Stock
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.




from __future__ import division
import numpy
import numpy.linalg as la




def main() :
    from math import sin, exp, sqrt

    from hedge.mesh import make_disk_mesh
    mesh = make_disk_mesh(r=0.5,faces=100,max_area=0.008)
    
    from hedge.backends.jit import Discretization
    discr = Discretization(mesh, order=3)

    from hedge.visualization import VtkVisualizer
    vis = VtkVisualizer(discr, None, "fld")

    def source_u(x, el):
        return exp(-numpy.dot(x, x)*128)

    source_u_vec = discr.interpolate_volume_function(source_u)
    
    def source_vec_getter(t):
        from math import sin
        return source_u_vec*sin(10*t)
    
    def bound_u(pt, el):
        x, y = pt
 if x > 0:
     return 1
 else:
     return 0

    # Building the boundary-vector - a line-vector - which contains less -
    # but the requested - information than the volume-vector of the entire
    # discretization space.
    bound_u_vec = discr.interpolate_boundary_function(bound_u)

    def bound_vec_getter(t):
        if t > 0.5:
     return bound_u_vec * 0
 else:
     return bound_u_vec

    from hedge.pde import StrongWaveOperator
    from hedge.mesh import TAG_ALL, TAG_NONE

    op = StrongWaveOperator(1, discr.dimensions,
            source_vec_getter,
     flux_type="upwind",
            dirichlet_tag=TAG_ALL,
     dirichlet_bc_f=bound_vec_getter,
            neumann_tag=TAG_NONE,
            radiation_tag=TAG_NONE)

    from hedge.tools import join_fields
    fields = join_fields(discr.volume_zeros(),
            [discr.volume_zeros() for i in range(discr.dimensions)])

    dt = discr.dt_factor(op.max_eigenvalue()) / 10
    nsteps = int(1/dt)

    # timestep loop -----------------------------------------------------------
    from hedge.timestep import AdamsBashforthTimeStepper
    stepper = AdamsBashforthTimeStepper(3,None)
    
    rhs = op.bind(discr)
    for step in range(nsteps):
        t = step*dt

        if step % 10 == 0:
            print step, t
            visf = vis.make_file("fld-%04d" % step)

            vis.add_data(visf,
                    [ ("u", fields[0]), ("v", fields[1:]), ],
                    time=t, step=step)
            visf.close()

        fields = stepper(fields, t, dt, rhs)

    discr.close()
    vis.close()




if __name__ == "__main__":
    main()
}}}

`/hedge/hedge/pde.py`:
{{{#!python
class StrongWaveOperator:
    """This operator discretizes the Wave equation S{part}tt u = c^2 S{Delta} u.

    To be precise, we discretize the hyperbolic system

      * S{part}t u - c div v = 0
      * S{part}t v - c grad u = 0

    Note that this is not unique--we could also choose a different sign for M{v}.
    """

    def __init__(self, c, dimensions, source_f=None,
            flux_type="upwind",
            dirichlet_tag=hedge.mesh.TAG_ALL,
     dirichlet_bc_f=None,
            neumann_tag=hedge.mesh.TAG_NONE,
            radiation_tag=hedge.mesh.TAG_NONE):
        assert isinstance(dimensions, int)

        self.c = c
        self.dimensions = dimensions
        self.source_f = source_f

        if self.c > 0:
            self.sign = 1
        else:
            self.sign = -1

        self.dirichlet_tag = dirichlet_tag
        self.neumann_tag = neumann_tag
        self.radiation_tag = radiation_tag

 self.dirichlet_bc_f = dirichlet_bc_f

        self.flux_type = flux_type

    def flux(self):
        from hedge.flux import FluxVectorPlaceholder, make_normal

        dim = self.dimensions
        w = FluxVectorPlaceholder(1+dim)
        u = w[0]
        v = w[1:]
        normal = make_normal(dim)

        from hedge.tools import join_fields
        flux_weak = join_fields(
                numpy.dot(v.avg, normal),
                u.avg * normal)

        if self.flux_type == "central":
            pass
        elif self.flux_type == "upwind":
            # see doc/notes/hedge-notes.tm
            flux_weak -= self.sign*join_fields(
                    0.5*(u.int-u.ext),
                    0.5*(normal * numpy.dot(normal, v.int-v.ext)))
        else:
            raise ValueError, "invalid flux type '%s'" % self.flux_type

        flux_strong = join_fields(
                numpy.dot(v.int, normal),
                u.int * normal) - flux_weak

        return -self.c*flux_strong

    def op_template(self):
        from hedge.optemplate import \
  Field, \
                make_vector_field, \
                pair_with_boundary, \
                get_flux_operator, \
                make_nabla, \
                InverseMassOperator

        d = self.dimensions

        w = make_vector_field("w", d+1)
        u = w[0]
        v = w[1:]

        # boundary conditions -------------------------------------------------
        from hedge.flux import make_normal
        normal = make_normal(d)

        from hedge.tools import join_fields

        if self.dirichlet_bc_f is None:
            dir_bc = join_fields(-u, v)
        else:
            dir_bc = join_fields(-Field("dir_bc_u"), v)

        neu_bc = join_fields(u, -v)
        rad_bc = join_fields(
                0.5*(u - self.sign*numpy.dot(normal, v)),
                0.5*normal*(numpy.dot(normal, v) - self.sign*u)
                )

        # entire operator -----------------------------------------------------
        nabla = make_nabla(d)
        flux_op = get_flux_operator(self.flux())

        from hedge.tools import join_fields
        return (
                - join_fields(
                    -self.c*numpy.dot(nabla, v),
                    -self.c*(nabla*u)
                    )
                +
                InverseMassOperator() * (
                    flux_op*w
                    + flux_op * pair_with_boundary(w, dir_bc, self.dirichlet_tag)
                    + flux_op * pair_with_boundary(w, neu_bc, self.neumann_tag)
                    + flux_op * pair_with_boundary(w, rad_bc, self.radiation_tag)
                    ))

    
    def bind(self, discr):
        from hedge.mesh import check_bc_coverage
        check_bc_coverage(discr.mesh, [
            self.dirichlet_tag,
            self.neumann_tag,
            self.radiation_tag])

        compiled_op_template = discr.compile(self.op_template())

        def rhs(t, w):
            from hedge.tools import join_fields, ptwise_dot

            if self.dirichlet_bc_f is None:
                rhs = compiled_op_template(w=w)
            else:
                rhs = compiled_op_template(w=w, dir_bc_u=self.dirichlet_bc_f(t))

            if self.source_f is not None:
                rhs[0] += self.source_f(t)

            return rhs

        return rhs

    def max_eigenvalue(self):
        return abs(self.c)
}}}

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):
  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 "--------------------"

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)

Advection Equation

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
   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()

/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,
  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)

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