Source Code for Module hedge.timestep.ab

  1  # -*- coding: utf8 -*- 
  2   
  3  """Adams-Bashforth ODE solvers.""" 
  4   
  5  from __future__ import division 
  6   
  7  __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" 
  8   
  9  __license__ = """ 
 10  This program is free software: you can redistribute it and/or modify 
 11  it under the terms of the GNU General Public License as published by 
 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 --------------------------------------------------------------- 
123 -class AdamsBashforthTimeStepper(TimeStepper): 
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: 
148              from operator import add 
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