[frames] | no frames]

# Source Code for Module hedge.timestep.ab

```  1  # -*- coding: utf8 -*-
2
4
5  from __future__ import division
6
8
10  This program is free software: you can redistribute it and/or modify
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14
15  This program is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18  GNU General Public License for more details.
19
20  You should have received a copy of the GNU General Public License
21  along with this program.  If not, see U{http://www.gnu.org/licenses/}.
22  """
23
24
25
26  import numpy
27  import numpy.linalg as la
28  from pytools import memoize
29  from hedge.timestep.base import TimeStepper
30
31
32
33
34  # coefficient generators ------------------------------------------------------
35 -def make_generic_ab_coefficients(levels, int_start, tap):
36      """Find coefficients (αᵢ) such that
37         ∑ᵢ αᵢ F(tᵢ) = ∫[int_start..tap] f(t) dt."""
38
39      # explanations --------------------------------------------------------------
40      # To calculate the AB coefficients this method makes use of the interpolation
41      # connection of the Vandermonde matrix:
42      #
43      #  Vᵀ * α = fe(t₀₊₁),                                    (1)
44      #
45      # with Vᵀ as the transposed Vandermonde matrix (with monomial base: xⁿ),
46      #
47      #  α = (..., α₋₂, α₋₁,α₀)ᵀ                               (2)
48      #
49      # a vector of interpolation coefficients and
50      #
51      #  fe(t₀₊₁) = (t₀₊₁⁰, t₀₊₁¹, t₀₊₁²,...,t₀₊₁ⁿ)ᵀ           (3)
52      #
53      # a vector of the evaluated interpolation polynomial f(t) at t₀₊₁ = t₀ ∓ h
54      # (h being any arbitrary stepsize).
55      #
56      # Solving the system (1) by knowing Vᵀ and fe(t₀₊₁) receiving α makes it
57      # possible for any function F(t) - the function which gets interpolated
58      # by the interpolation polynomial f(t) - to calculate f(t₀₊₁) by:
59      #
60      # f(t₀₊₁) =  ∑ᵢ αᵢ F(tᵢ)                                 (5)
61      #
62      # with F(tᵢ) being the values of F(t) at the sampling points tᵢ.
63      # --------------------------------------------------------------------------
64      # The Adams-Bashforth method is defined by:
65      #
66      #  y(t₀₊₁) = y(t₀) + Δt * ∫₀⁰⁺¹ f(t) dt                  (6)
67      #
68      # with:
69      #
70      #  ∫₀⁰⁺¹ f(t) dt = ∑ᵢ ABcᵢ F(tᵢ),                        (8)
71      #
72      # with ABcᵢ = [AB coefficients], f(t) being the interpolation polynomial,
73      # and F(tᵢ) being the values of F (= RHS) at the sampling points tᵢ.
74      # --------------------------------------------------------------------------
75      # For the AB method (1) becomes:
76      #
77      #  Vᵀ * ABc = ∫₀⁰⁺¹ fe(t₀₊₁)                             (7)
78      #
79      # with ∫₀⁰⁺¹ fe(t₀₊₁) being a vector evalueting the integral of the
80      # interpolation polynomial in the form oft
81      #
82      #  1/(n+1)*(t₀₊₁⁽ⁿ⁾-t₀⁽ⁿ⁾)                               (8)
83      #
84      #  for n = 0,1,...,N sampling points, and
85      #
86      # ABc = [c₀,c₁, ... , cn]ᵀ                               (9)
87      #
88      # being the AB coefficients.
89      #
90      # For example ∫₀⁰⁺¹ f(t₀₊₁) evaluated for the timestep [t₀,t₀₊₁] = [0,1]
91      # is:
92      #
93      #  point_eval_vec = [1, 0.5, 0.333, 0.25, ... ,1/n]ᵀ.
94      #
95      # For substep levels the bounds of the integral has to be adapted to the
96      # size and position of the substep interval:
97      #
98      #  [t₀,t₀₊₁] = [substep_int_start, substep_int_end]
99      #
100      # which is equal to the implemented [int_start, tap].
101      #
102      # Since Vᵀ and ∫₀⁰⁺¹ f(t₀₊₁) is known the AB coefficients c can be
103      # predicted by solving system (7) and calculating:
104      #
105      #  ∫₀⁰⁺¹ f(t) dt = ∑ᵢ ABcᵢ F(tᵢ),
106
107      from hedge.polynomial import monomial_vdm
108      point_eval_vec = numpy.array([
109          1/(n+1)*(tap**(n+1)-int_start**(n+1)) for n in range(len(levels))])
110      return la.solve(monomial_vdm(levels).T, point_eval_vec)
111
112
113
114
115  @memoize
116 -def make_ab_coefficients(order):
117      return make_generic_ab_coefficients(numpy.arange(0, -order, -1), 0, 1)
118
119
120
121
122  # time steppers ---------------------------------------------------------------
124      dt_fudge_factor = 0.95
125
126 -    def __init__(self, order, startup_stepper=None):
127          self.coefficients = make_ab_coefficients(order)
128          self.f_history = []
129
130          if startup_stepper is not None:
131              self.startup_stepper = startup_stepper
132          else:
133              from hedge.timestep.rk4 import RK4TimeStepper
134              self.startup_stepper = RK4TimeStepper()
135
136 -    def __call__(self, y, t, dt, rhs):
137          if len(self.f_history) == 0:
138              # insert IC
139              self.f_history.append(rhs(t, y))
140
141          if len(self.f_history) < len(self.coefficients):
142              ynew = self.startup_stepper(y, t, dt, rhs)
143              if len(self.f_history) == len(self.coefficients) - 1:
144                  # here's some memory we won't need any more
145                  del self.startup_stepper
146
147          else:
149
150              assert len(self.coefficients) == len(self.f_history)
151              ynew = y + dt * reduce(add,
152                      (coeff * f
153                          for coeff, f in
154                          zip(self.coefficients, self.f_history)))
155
156              self.f_history.pop()
157
158          self.f_history.insert(0, rhs(t+dt, ynew))
159          return ynew
160
<!--
expandto(location.href);
// -->

```

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