Source Code for Module hedge.log

  1  """Timeseries data gathering sensors.""" 
  2   
  3  from __future__ import division 
  4   
  5  __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" 
  6   
  7  __license__ = """ 
  8  This program is free software: you can redistribute it and/or modify 
  9  it under the terms of the GNU General Public License as published by 
 10  the Free Software Foundation, either version 3 of the License, or 
 11  (at your option) any later version. 
 12   
 13  This program is distributed in the hope that it will be useful, 
 14  but WITHOUT ANY WARRANTY; without even the implied warranty of 
 15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 16  GNU General Public License for more details. 
 17   
 18  You should have received a copy of the GNU General Public License 
 19  along with this program.  If not, see U{http://www.gnu.org/licenses/}. 
 20  """ 
 21   
 22   
 23   
 24  from pytools.log import LogQuantity, MultiLogQuantity 
 25  import numpy 
 26   
 27   
 28   
 29   
 30 -def axis_name(axis): 
 31      if axis == 0: return "x" 
 32      elif axis == 1: return "y" 
 33      elif axis == 2: return "z" 
 34      else: raise RuntimeError, "invalid axis index" 
 35   
 36   
 37   
 38   
 39 -class Integral(LogQuantity): 
 40      """Log the volume integral of a variable in a scope.""" 
 41   
 42 -    def __init__(self, getter, discr, name=None,  
 43              unit="1", description=None): 
 44          """Construct the integral logger. 
 45   
 46          @arg getter: a callable that returns the value of which to  
 47            take the integral. 
 48          @arg discr: a L{discretization.Discretization} to which the variable belongs. 
 49          @arg name: the name reported to the C{LogManager}. 
 50          @arg unit: the unit of measure for the log quantity. 
 51          @arg description: A description fed to the C{LogManager}. 
 52          """ 
 53          self.getter = getter 
 54   
 55          if name is None: 
 56              try: 
 57                  name = "int_%s" % self.getter.name() 
 58              except AttributeError: 
 59                  raise ValueError, "must specify a name" 
 60   
 61          LogQuantity.__init__(self, name, unit, description) 
 62   
 63          self.discr = discr 
 64   
 65      @property 
 66 -    def default_aggregator(self):  
 67          return sum 
 68   
 69 -    def __call__(self): 
 70          var = self.getter() 
 71   
 72          from hedge.tools import log_shape 
 73   
 74          if len(log_shape(var)) == 1: 
 75              return sum( 
 76                      self.discr.integral(numpy.abs(v)) 
 77                      for v in var) 
 78          else: 
 79              return self.discr.integral(var) 
 80   
 81   
 82   
 83   
 84 -class LpNorm(LogQuantity): 
 85      """Log the Lp norm of a variable in a scope.""" 
 86   
 87 -    def __init__(self, getter, discr, p=2, name=None,  
 88              unit="1", description=None): 
 89          """Construct the Lp norm logger. 
 90   
 91          @arg getter: a callable that returns the value of which to  
 92            take the norm. 
 93          @arg discr: a L{discretization.Discretization} to which the variable belongs. 
 94          @arg p: the power of the norm. 
 95          @arg name: the name reported to the C{LogManager}. 
 96          @arg unit: the unit of measure for the log quantity. 
 97          @arg description: A description fed to the C{LogManager}. 
 98          """ 
 99          self.getter = getter 
100          self.discr = discr 
101          self.p = p 
102   
103          if name is None: 
104              try: 
105                  name = "l%d_%s" % (int(p), self.getter.name()) 
106              except AttributeError: 
107                  raise ValueError, "must specify a name" 
108   
109          LogQuantity.__init__(self, name, unit, description) 
110   
111      @property 
112 -    def default_aggregator(self):  
113          from pytools import norm_inf, Norm 
114   
115          if self.p == numpy.Inf: 
116              return norm_inf 
117          else: 
118              return Norm(self.p) 
119   
120 -    def __call__(self): 
121          var = self.getter() 
122          return self.discr.norm(var, self.p) 
123   
124   
125   
126   
127  # electromagnetic quantities -------------------------------------------------- 
128 -class EMFieldGetter(object): 
129      """Makes E and H field accessible as self.e and self.h from a variable lookup. 
130      To be used with the EM log quantities in this module.""" 
131 -    def __init__(self, discr, maxwell_op, fgetter): 
132          self.discr = discr 
133          self.maxwell_op = maxwell_op 
134          self.fgetter = fgetter 
135   
136      @property 
137 -    def e(self): 
138          fields = self.fgetter() 
139          e, h = self.maxwell_op.split_eh(fields) 
140          return e 
141   
142      @property 
143 -    def h(self): 
144          fields = self.fgetter() 
145          e, h = self.maxwell_op.split_eh(fields) 
146          return h 
147   
148   
149   
150 -class ElectricFieldEnergy(LogQuantity): 
151 -    def __init__(self, fields, name="W_el"): 
152          LogQuantity.__init__(self, name, "J", "Energy of the electric field") 
153          self.fields = fields 
154   
155      @property 
156 -    def default_aggregator(self):  
157          from pytools import norm_2 
158          return norm_2 
159   
160 -    def __call__(self): 
161          max_op = self.fields.maxwell_op 
162   
163          e = self.fields.e 
164          d = max_op.epsilon * e 
165   
166          from hedge.tools import ptwise_dot 
167          energy_density = 1/2*(ptwise_dot(1, 1, e, d)) 
168          return self.fields.discr.integral(energy_density) 
169   
170   
171   
172   
173 -class MagneticFieldEnergy(LogQuantity): 
174 -    def __init__(self, fields, name="W_mag"): 
175          LogQuantity.__init__(self, name, "J", "Energy of the magnetic field") 
176          self.fields = fields 
177   
178      @property 
179 -    def default_aggregator(self):  
180          from pytools import norm_2 
181          return norm_2 
182   
183 -    def __call__(self): 
184          max_op = self.fields.maxwell_op 
185   
186          h = self.fields.h 
187          b = max_op.mu * h 
188   
189          from hedge.tools import ptwise_dot 
190          energy_density = 1/2*(ptwise_dot(1, 1, h, b)) 
191          return self.fields.discr.integral(energy_density) 
192   
193   
194   
195 -class EMFieldMomentum(MultiLogQuantity): 
196 -    def __init__(self, fields, c0, names=None): 
197          if names is None: 
198              names = ["p%s_field" % axis_name(i)  
199                      for i in range(3)] 
200   
201          vdim = len(names) 
202   
203          MultiLogQuantity.__init__(self, names,  
204              units=["N*s"] * vdim,  
205              descriptions=["Field Momentum"] * vdim) 
206   
207          self.fields = fields 
208          self.c0 = c0 
209   
210          e_subset = fields.maxwell_op.get_eh_subset()[0:3] 
211          h_subset = fields.maxwell_op.get_eh_subset()[3:6] 
212   
213          from hedge.tools import SubsettableCrossProduct 
214          self.poynting_cross = SubsettableCrossProduct( 
215                  op1_subset=e_subset,  
216                  op2_subset=h_subset,  
217                  ) 
218   
219 -    def __call__(self): 
220          max_op = self.fields.maxwell_op 
221   
222          e = self.fields.e 
223          h = self.fields.h 
224   
225          poynting_s = self.poynting_cross(e, h) 
226   
227          momentum_density = poynting_s/self.c0**2 
228          return self.fields.discr.integral(momentum_density) 
229   
230   
231   
232   
233 -class EMFieldDivergenceD(LogQuantity): 
234 -    def __init__(self, maxwell_op, fields, name="divD"): 
235          LogQuantity.__init__(self, name, "C", "Integral over div D") 
236   
237          self.fields = fields 
238   
239          from hedge.models.nd_calculus import DivergenceOperator 
240          div_op = DivergenceOperator(maxwell_op.dimensions, 
241                  maxwell_op.get_eh_subset()[:3]) 
242          self.bound_div_op = div_op.bind(self.fields.discr) 
243   
244 -    def __call__(self): 
245          max_op = self.fields.maxwell_op 
246          d = max_op.epsilon * self.fields.e 
247          div_d = self.bound_div_op(d) 
248           
249          return self.fields.discr.integral(div_d) 
250   
251   
252   
253   
254 -class EMFieldDivergenceB(MultiLogQuantity): 
255 -    def __init__(self, maxwell_op, fields, names=None): 
256          self.fields = fields 
257   
258          from hedge.models.nd_calculus import DivergenceOperator 
259          self.div_op = DivergenceOperator(maxwell_op.dimensions, 
260                  maxwell_op.get_eh_subset()[3:]).bind(self.fields.discr) 
261   
262          if names is None: 
263              names = ["divB", "err_divB_l1"] 
264   
265          MultiLogQuantity.__init__(self,  
266                  names=names, 
267                  units=["T/m", "T/m"],  
268                  descriptions=["Integral over div B", "Integral over |div B|"]) 
269   
270 -    def __call__(self): 
271          max_op = self.fields.maxwell_op 
272          b = max_op.mu * self.fields.h 
273          div_b = self.div_op(b) 
274           
275          return [self.fields.discr.integral(div_b),  
276                  self.fields.discr.integral(numpy.absolute(div_b))] 
277   
278   
279   
280   
281 -def add_em_quantities(mgr, maxwell_op, fields): 
282      mgr.add_quantity(ElectricFieldEnergy(fields)) 
283      mgr.add_quantity(MagneticFieldEnergy(fields)) 
284      mgr.add_quantity(EMFieldMomentum(fields, maxwell_op.c)) 
285      mgr.add_quantity(EMFieldDivergenceD(maxwell_op, fields)) 
286      mgr.add_quantity(EMFieldDivergenceB(maxwell_op, fields)) 
287