Package sfc :: Package representation :: Module integralrepresentation
[hide private]
[frames] | no frames]

Source Code for Module sfc.representation.integralrepresentation

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  This module contains representation classes for integrals. 
  4  """ 
  5   
  6  # Copyright (C) 2008-2009 Martin Sandve Alnes and Simula Resarch Laboratory 
  7  # 
  8  # This file is part of SyFi. 
  9  # 
 10  # SyFi 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 2 of the License, or 
 13  # (at your option) any later version. 
 14  # 
 15  # SyFi 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 SyFi. If not, see <http://www.gnu.org/licenses/>. 
 22  # 
 23  # Modified by Kent-Andre Mardal, 2010. 
 24  # 
 25  # First added:  2008-08-13 
 26  # Last changed: 2009-03-12 
 27   
 28  from itertools import izip 
 29   
 30  #import SyFi 
 31  import swiginac 
 32   
 33  #import ufl 
 34  from ufl.permutation import compute_indices 
 35  from ufl.algorithms import Graph, partition, expand_indices2, expand_indices, strip_variables, tree_format 
 36  from ufl.classes import Expr, Terminal, UtilityType, SpatialDerivative, Indexed 
 37   
 38  from sfc.common import sfc_assert, sfc_warning, sfc_debug, sfc_error 
 39  from sfc.common.utilities import indices_subset 
 40  from sfc.codegeneration.codeformatting import indent, row_major_index_string 
 41  from sfc.symbolic_utils import symbols, symbol 
 42  from sfc.representation.swiginac_eval import SwiginacEvaluator, EvaluateAsSwiginac 
 43   
44 -def find_locals(Vin, Vout, P):
45 # Count the number of vertices depending on each vertex 46 Vcount = [len(vi) for vi in Vin] 47 # Count down all dependencies internally in the partition, O(|E| log|P|) 48 Pset = set(P) 49 for i in P: # v_i is in our partition 50 for j in Vout[i]: # v_i depends on v_j 51 if j in Pset: # is v_j in our partition? 52 Vcount[j] -= 1 # don't count with internal dependencies in partition 53 # Return mapping (v_i -> is_local_to_P) for all v_i in P 54 return dict((i, (Vcount[i] == 0)) for i in P)
55
56 -def is_simple(e):
57 return (e.nops() == 0) 58 # TODO: Try something like this: 59 #return (e.nops() == 0) or \ 60 # (isinstance(e, (swiginac.add, swiginac.mul)) \ 61 # and all(o.nops() for o in e.ops())) 62
63 -class IntegralData:
64 - def __init__(self, itgrep, integral):
65 self.itgrep = itgrep 66 self.integral = integral 67 68 # TODO: Override these and other options with metadata 69 metadata = integral.measure().metadata() 70 self.integration_method = itgrep.options.integration_method 71 self.integration_order = itgrep.options.integration_order 72 73 self.quad_rule = None 74 self.facet_quad_rules = None 75 76 self.G = None 77 self.Vdeps = None 78 self.Vin_count = None 79 self.partitions = {} 80 self.evaluator = None 81 self.evaluate = None 82 83 self.V_symbols = {} 84 self.vertex_data_set = {}
85
86 - def store(self, value, i, basis_functions):
87 #print "STORING AT", i, basis_functions, str(value) 88 vd = self.vertex_data_set.get(basis_functions) 89 if vd is None: 90 vd = {} 91 self.vertex_data_set[basis_functions] = vd 92 vd[i] = value
93
94 - def get_vd(self, j, basis_functions):
95 jkeep = [(("v%d" % k) in self.Vdeps[j]) for k in range(self.itgrep.formrep.rank)] 96 iota, = indices_subset([basis_functions], jkeep) 97 if all(i is None for i in iota): 98 iota = () 99 vd = self.vertex_data_set[iota] 100 return vd
101
102 - def fetch_storage(self, j, basis_functions):
103 vd = self.get_vd(j, basis_functions) 104 return vd[j]
105
106 - def free_storage(self, j, basis_functions):
107 "Delete stored data for j." 108 vd = self.get_vd(j, basis_functions)
109 #del vd[j] # FIXME: Enable when counts are corrected 110
111 -class IntegralRepresentation(object):
112 - def __init__(self, integrals, formrep, on_facet):
113 sfc_debug("Entering IntegralRepresentation.__init__") 114 115 self.integrals = integrals 116 self.formrep = formrep 117 self._on_facet = on_facet 118 self.classname = formrep.itg_names[integrals[0]] 119 120 self.options = self.formrep.options.code.integral 121 122 # Shape of element tensor 123 A_shape = [] 124 for i in range(self.formrep.rank): 125 element = self.formrep.formdata.elements[i] 126 rep = self.formrep.element_reps[element] 127 A_shape.append(rep.local_dimension) 128 self.A_shape = tuple(A_shape) 129 130 # All permutations of element tensor indices TODO: Not same for interior_facet_integrals, move to subclasses 131 self.indices = compute_indices(self.A_shape) 132 133 # Symbols for output tensor entry A[iota] 134 Asym = symbols("A[%s]" % row_major_index_string(i, self.A_shape) for i in self.indices) 135 self.A_sym = dict(izip(self.indices, Asym)) 136 137 # Data structures for token symbols 138 self._free_symbols = [] 139 self._symbol_counter = 0 140 141 # Data structures to hold integrals for each measure 142 self.symbolic_integral = None # Only a single symbolic integral is allowed 143 self.quadrature_integrals = [] # Multiple quadrature integrals are possible, each with different quadrature rules 144 self.integral_data = {} # Convenient mapping from measure to IntegralData object 145 146 # Compute data about integrals for each measure 147 fd = self.formrep.formdata 148 for integral in self.integrals: 149 data = IntegralData(self, integral) 150 151 if data.integration_method == "symbolic": 152 self.symbolic_integral = integral 153 data.evaluator = SwiginacEvaluator(self.formrep, use_symbols=False, on_facet=self._on_facet) 154 155 elif data.integration_method == "quadrature": 156 157 # FIXME: Build a separate quadrature rule for each integral, if necessary 158 data.quad_rule = self.formrep.quad_rule 159 data.facet_quad_rules = self.formrep.facet_quad_rules 160 161 self.quadrature_integrals.append(integral) 162 if self.options.safemode: 163 # Produce naive code, good for debugging and verification 164 data.evaluator = SwiginacEvaluator(self.formrep, use_symbols=True, on_facet=self._on_facet) 165 166 else: 167 # Create an evaluation multifunction 168 data.evaluate = EvaluateAsSwiginac(self.formrep, self, data, on_facet=self._on_facet) 169 170 # Expand all indices before building graph 171 integrand = integral.integrand() 172 integrand = strip_variables(integrand) 173 if self.options.use_expand_indices2: 174 integrand = expand_indices2(integrand) 175 else: 176 integrand = expand_indices(integrand) 177 178 # Build linearized computational graph 179 data.G = Graph(integrand) 180 181 # Count dependencies for each node 182 V = data.G.V() 183 n = len(V) 184 data.Vin_count = [0]*n 185 for i, vs in enumerate(data.G.Vin()): 186 for j in vs: 187 vj = V[j] 188 # ... FIXME: Need to store one counter for each iota relevant for i, and to add "size(iota_space_j)" instead of 1 189 data.Vin_count[i] += 1 190 #data.Vin_count[-1] += ... # FIXME: Increment further for integrand root vertex 191 192 # Partition graph based on dependencies 193 #basis_function_deps = [set(("v%d" % i, "x")) \ 194 # for (i, bf) in enumerate(fd.basis_functions)] 195 #function_deps = [set(("w", "c", "x",)) \ 196 # for (i, bf) in enumerate(fd.functions)] 197 # TODO: Input actual dependencies of functions and 198 # basis functions and their derivatives here 199 data.partitions, data.Vdeps = partition(data.G) 200 201 self.integral_data[integral.measure()] = data 202 203 sfc_debug("Leaving IntegralRepresentation.__init__")
204 205 # --- Printing --- 206
207 - def __str__(self):
208 s = "" 209 s += "IntegralRepresentation:\n" 210 s += " classname: %s\n" % self.classname 211 s += " A_shape: %s\n" % self.A_shape 212 s += " indices: %s\n" % self.indices 213 s += " A_sym: %s\n" % self.A_sym 214 215 # Not computing whole A at any point because of memory usage 216 #s += " A:\n" 217 #for ii in self.indices: 218 # s += " %s = %s\n" % (self.A_sym[ii], self.A[ii]) 219 220 s += " UFL Integrals:\n" 221 for integral in self.integrals: 222 s += indent(str(integral)) + "\n\n" 223 return s.strip("\n")
224 225 # --- Utilities for allocation of symbols and association of data with graph vertices --- 226
227 - def allocate(self, data, j):
228 """Allocate symbol(s) for vertex j. Either gets symbol from the free 229 symbol set or creates a new symbol and increases the symbol counter.""" 230 if self._free_symbols: 231 return self._free_symbols.pop() 232 s = symbol("s[%d]" % self._symbol_counter) 233 self._symbol_counter += 1 234 data.V_symbols[j] = s 235 return s
236
237 - def free_symbols(self, data, j):
238 "Delete stored data for j and make its symbols available again." 239 s = data.V_symbols.get(j) 240 if s is None: 241 sfc_debug("Trying to deallocate symbols that are not allocated!") 242 return 243 self._free_symbols.append(s) 244 del data.V_symbols[j]
245 246 # --- Generators for tokens of various kinds --- 247
248 - def iter_partition(self, data, deps, basis_functions=()):
249 sfc_debug("Entering IntegralRepresentation.iter_partition") 250 251 deps = frozenset(deps) 252 253 P = data.partitions.get(deps) 254 if not P: 255 sfc_debug("Leaving IntegralRepresentation.iter_partition, empty") 256 return 257 258 data.evaluate.current_basis_function = basis_functions 259 260 # Get graph components 261 G = data.G 262 V = G.V() 263 E = G.E() 264 Vin = G.Vin() 265 Vout = G.Vout() 266 267 # Figure out which variables are only needed within the partition 268 is_local = find_locals(Vin, Vout, P) 269 270 # For all vertices in partition P 271 for i in P: 272 v = V[i] 273 274 if isinstance(v, UtilityType): 275 # Skip labels and indices 276 continue 277 278 if v.shape(): 279 # Skip expressions with shape: they will be evaluated when indexed. 280 if not isinstance(v, (Terminal, SpatialDerivative)): 281 print "="*30 282 print "type:", type(v) 283 print "str:", str(v) 284 print "child types:", [type(V[j]) for j in Vout[i]] 285 print "child str:" 286 print "\n".join( " vertex %d: %s" % (j, str(V[j])) for j in Vout[i] ) 287 print "number of parents:", len(Vin[i]) 288 if len(Vin[i]) < 5: 289 print "parent types:", [type(V[j]) for j in Vin[i]] 290 #print "parent str:" 291 #print "\n".join( " vertex %d: %s" % (j, str(V[j])) for j in Vin[i] ) 292 sfc_error("Expecting all indexing to have been propagated to terminals?") 293 continue 294 295 if Vin[i] and all(isinstance(V[j], SpatialDerivative) for j in Vin[i]): 296 # Skipping expressions that only occur later in 297 # differentiated form: we don't need their non-differentiated 298 # expression so they will be evaluated when indexed. 299 if not isinstance(v, (Terminal, SpatialDerivative, Indexed)): 300 print "="*30 301 print type(v) 302 print str(v) 303 sfc_error("Expecting all indexing to have been propagated to terminals?") 304 continue 305 306 if isinstance(v, (Indexed, SpatialDerivative)): 307 ops = v.operands() 308 309 # Evaluate vertex v with given mapped ops 310 if not all(isinstance(o, (Expr, swiginac.basic)) for o in ops): 311 print ";"*80 312 print tree_format(v) 313 print str(v) 314 print type(ops) 315 print str(ops) 316 print repr(ops) 317 print "types:" 318 print "\n".join(str(type(o)) for o in ops) 319 print ";"*80 320 321 e = data.evaluate(v, *ops) 322 323 else: 324 # Get already computed operands 325 # (if a vertex isn't already computed 326 # at this point, that is a bug) 327 ops = [] 328 for j in Vout[i]: 329 try: 330 # Fetch expression or symbol for vertex j 331 e = data.fetch_storage(j, basis_functions) 332 except: 333 print "Failed to fetch expression for vertex %d," % j 334 print " V[%d] = %s" % (j, repr(V[j])) 335 print " parent V[%d] = %s" % (i, repr(V[i])) 336 raise RuntimeError 337 ops.append(e) 338 ops = tuple(ops) 339 e = data.evaluate(v, *ops) 340 341 # TODO: Make sure that dependencies of skipped expressions 342 # are counted down when evaluated later! 343 # Or do we need that? We don't make symbols for them anyway. 344 345 # Count down number of uses of v's dependencies 346 for j in Vout[i]: # For each vertex j that depend on vertex i 347 data.Vin_count[j] -= 1 348 if False: # data.Vin_count[j] == 0: # FIXME: This doesn't work yet, counters are wrong 349 # Make the symbols for j available again 350 self.free_symbols(data, j) 351 data.free_storage(j, basis_functions) 352 353 # Store some result for vertex i and yield token based on some heuristic 354 if is_simple(e): 355 # Store simple expression associated with 356 # vertex j, no code generation necessary 357 data.store(e, i, basis_functions) 358 else: 359 if is_local[i]: 360 pass # TODO: Use this information to make a local set of symbols "sl[%d]" for each partition 361 362 # Allocate a symbol and remember it, just throw away the expression 363 s = self.allocate(data, i) 364 data.store(s, i, basis_functions) 365 # Only yield when we have a variable to generate code for! 366 yield (s, e) 367 368 sfc_debug("Leaving IntegralRepresentation.iter_partition")
369 370 # Some stuff from the old code: 371 #bf = self._current_basis_function # FIXME: Filter depending on deps! 372 #key = tuple(chain((count,), component, index_values, bf)) 373 #compstr = "_".join("%d" % k for k in key) 374 #vname = "s_%s" % compstr 375 376 # Register token (s, e) with variable v in evaluator, 377 # such that it can return it from evaluator.variable 378 #self.evaluator._variable2symbol[key] = vsym 379
380 - def iter_member_quad_tokens(self, data): # FIXME: Make this into precomputation of a static array of constants
381 "Return an iterator over member tokens dependent of spatial variables. Overload in subclasses!" 382 383 assert data.integration_method == "quadrature" 384 385 # TODO: Precompute basis function values symbolically here instead of generating a loop in the constructor 386 387 # TODO: Skip what's not needed! 388 389 # Precompute all basis functions 390 fr = self.formrep 391 fd = fr.formdata 392 generated = set() 393 for iarg in range(fr.rank + fr.num_coefficients): 394 elm = fd.elements[iarg] 395 rep = fr.element_reps[elm] 396 for i in range(rep.local_dimension): 397 for component in rep.value_components: 398 # Yield basis function itself 399 s = fr.v_sym(iarg, i, component, self._on_facet) 400 if not (s == 0 or s in generated): 401 e = fr.v_expr(iarg, i, component) 402 t = (s, e) 403 yield t 404 generated.add(s) 405 # Yield its derivatives w.r.t. local coordinates 406 for d in range(fr.cell.nsd): 407 der = (d,) 408 s = fr.dv_sym(iarg, i, component, der, self._on_facet) 409 if not (s == 0 or s in generated): 410 e = fr.dv_expr(iarg, i, component, der) 411 t = (s, e) 412 yield t 413 generated.add(s)
414
415 - def iter_geometry_tokens(self):
416 "Return an iterator over geometry tokens independent of spatial variables. Overload in subclasses!" 417 fr = self.formrep 418 419 # TODO: Skip what's not needed! 420 421 # vx 422 for (ss,ee) in zip(fr.vx_sym, fr.vx_expr): 423 for i in range(ss.nops()): 424 yield (ss.op(i), ee.op(i)) 425 426 # G 427 (ss,ee) = (fr.G_sym, fr.G_expr) 428 for i in range(ss.nops()): 429 yield (ss.op(i), ee.op(i)) 430 431 # detG 432 yield (fr.detGtmp_sym, fr.detGtmp_expr) 433 yield (fr.detG_sym, fr.detG_expr) 434 435 # Ginv 436 (ss,ee) = (fr.Ginv_sym, fr.Ginv_expr) 437 for i in range(ss.nops()): 438 yield (ss.op(i), ee.op(i)) 439 440 if self._on_facet: 441 # Needed for normal vector TODO: Skip if not needed 442 yield (fr.detG_sign_sym, fr.detG_sign_expr) 443 else: 444 if self.symbolic_integral is not None: 445 # Scaling by cell volume factor, determinant of coordinate mapping 446 yield (fr.D_sym, fr.detG_sym)
447
448 - def iter_runtime_quad_tokens(self, data):
449 "Return an iterator over runtime tokens dependent of spatial variables. Overload in subclasses!" 450 assert data.integration_method == "quadrature" 451 452 # TODO: yield geometry tokens like G etc here instead if they depend on x,y,z 453 454 # TODO: Skip what's not needed! 455 456 # Generate all basis function derivatives 457 fr = self.formrep 458 fd = fr.formdata 459 generated = set() 460 for iarg in range(fr.rank + fr.num_coefficients): 461 elm = fd.elements[iarg] 462 rep = fr.element_reps[elm] 463 for i in range(rep.local_dimension): 464 for component in rep.value_components: 465 # Yield first order derivatives w.r.t. global coordinates 466 for d in range(fr.cell.nsd): 467 der = (d,) 468 s = fr.Dv_sym(iarg, i, component, der, self._on_facet) 469 if not (s == 0 or s in generated): 470 e = fr.Dv_expr(iarg, i, component, der, True, self._on_facet) 471 t = (s, e) 472 yield t 473 generated.add(s) 474 475 # Currently placing all w in here: 476 generated = set() 477 for iarg in range(fr.num_coefficients): 478 elm = fd.elements[fr.rank+iarg] 479 rep = fr.element_reps[elm] 480 for component in rep.value_components: 481 # Yield coefficient function itself 482 s = fr.w_sym(iarg, component) 483 if not s in generated: 484 e = fr.w_expr(iarg, component, True, self._on_facet) 485 t = (s, e) 486 yield t 487 generated.add(s) 488 # Yield first order derivatives w.r.t. global coordinates 489 for d in range(fr.cell.nsd): 490 der = (d,) 491 s = fr.Dw_sym(iarg, component, der) 492 if not (s == 0 or s in generated): 493 e = fr.Dw_expr(iarg, component, der, True, self._on_facet) 494 t = (s, e) 495 yield t 496 generated.add(s) 497 498 # Scaling factor 499 if self._on_facet: 500 # TODO: yield tokens that depend on both x and facet here 501 502 # Scaling by facet area factor 503 D_expr = fr.quad_weight_sym*fr.facet_D_sym 504 else: 505 # Scaling by cell volume factor, determinant of coordinate mapping 506 D_expr = fr.quad_weight_sym*fr.detG_sym 507 yield (fr.D_sym, D_expr)
508 509 # --- Element tensor producers 510
511 - def iter_A_tokens(self, data, facet=None):
512 "Iterate over all A[iota] tokens." 513 for iota in self.indices: 514 A_sym = self.A_sym[iota] 515 A_expr = self.compute_A(data, iota, facet) 516 yield (A_sym, A_expr)
517
518 - def compute_A(self, data, iota, facet=None):
519 "Compute expression for A[iota]. Overload in subclasses!" 520 raise NotImplementedError
521