hedge.flux

1 """Building blocks for flux computation. Flux compilation.""" 2 3 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" 4 5 __license__ = """ 6 This program is free software: you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation, either version 3 of the License, or 9 (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with this program. If not, see U{http://www.gnu.org/licenses/}. 18 """ 19 20 21 22 23 import numpy 24 import hedge._internal as _internal 25 import pymbolic.primitives 26 import pymbolic.mapper.collector 27 import pymbolic.mapper.expander 28 import pymbolic.mapper.flattener 29 import pymbolic.mapper.substitutor 30 import pymbolic.mapper.constant_folder 31 import pymbolic.mapper.flop_counter 32 33 34 35 36 FluxFace = _internal.FluxFace

37 38 39 40 41 # python fluxes --------------------------------------------------------------- 42 -class Flux(pymbolic.primitives.AlgebraicLeaf):

43 - def stringifier(self):

44 return FluxStringifyMapper

45

46 47 48 49 -class FluxScalarParameter(pymbolic.primitives.Variable):

50 - def __init__(self, name, is_complex=False):

51 pymbolic.primitives.Variable.__init__(self, name) 52 self.is_complex = is_complex

53

54 - def get_mapper_method(self, mapper):

55 return mapper.map_scalar_parameter

56

57 58 59 60 -class FieldComponent(Flux):

61 - def __init__(self, index, is_local):

62 self.index = index 63 self.is_local = is_local

64

65 - def is_equal(self, other):

66 return (isinstance(other, FieldComponent) 67 and self.index == other.index 68 and self.is_local == other.is_local 69 )

70

71 - def __getinitargs__(self):

72 return self.index, self.is_local

73

74 - def get_hash(self):

75 return hash(( 76 self.__class__, 77 self.index, 78 self.is_local))

79

80 - def get_mapper_method(self, mapper):

81 return mapper.map_field_component

82

83 84 85 86 -class Normal(Flux):

87 - def __init__(self, axis):

88 self.axis = axis

89

90 - def __getinitargs__(self):

91 return self.axis,

92

93 - def is_equal(self, other):

94 return isinstance(other, Normal) and self.axis == other.axis

95

96 - def get_hash(self):

97 return hash(( 98 self.__class__, 99 self.axis))

100

101 - def get_mapper_method(self, mapper):

102 return mapper.map_normal

103

104 105 106 107 -class PenaltyTerm(Flux):

108 - def __init__(self, power=1):

109 self.power = power

110

111 - def __getinitargs__(self):

112 return self.power,

113

114 - def is_equal(self, other):

115 return isinstance(other, PenaltyTerm) and self.power == other.power

116

117 - def get_hash(self):

118 return hash(( 119 self.__class__, 120 self.power))

121

122 - def get_mapper_method(self, mapper):

123 return mapper.map_penalty_term

124

125 126 127 128 -class IfPositive(Flux):

129 - def __init__(self, criterion, then, else_):

130 self.criterion = criterion 131 self.then = then 132 self.else_ = else_

133

134 - def __getinitargs__(self):

135 return self.criterion, self.then, self.else_

136

137 - def is_equal(self, other):

138 return (isinstance(other, IfPositive) 139 and self.criterion == other.criterion 140 and self.then == other.then 141 and self.else_ == other.else_)

142

143 - def get_hash(self):

144 return hash(( 145 self.__class__, 146 self.criterion, 147 self.then, 148 self.else_))

149

150 - def get_mapper_method(self, mapper):

151 return mapper.map_if_positive

152

153 154 155 156 -class FluxFunctionSymbol(pymbolic.primitives.FunctionSymbol):

157 pass

158

159 -class Abs(FluxFunctionSymbol):

160 arg_count = 1

161

162 -class Max(FluxFunctionSymbol):

163 arg_count = 2

164

165 -class Min(FluxFunctionSymbol):

166 arg_count = 2

167 168 flux_abs = Abs() 169 flux_max = Max() 170 flux_min = Min()

171 172 173 174 -def norm(v):

175 return numpy.dot(v, v)**0.5

176

177 178 179 180 181 -def make_normal(dimensions):

182 return numpy.array([Normal(i) for i in range(dimensions)], dtype=object)

183

184 185 186 187 188 -class FluxZeroPlaceholder(object):

189 @property

190 - def int(self):

191 return 0

192 193 @property

194 - def ext(self):

195 return 0

196 197 @property

198 - def avg(self):

199 return 0

200

201 202 203 204 -class FluxScalarPlaceholder(object):

205 - def __init__(self, component=0):

206 self.component = component

207 208 @property

209 - def int(self):

210 return FieldComponent(self.component, True)

211 212 @property

213 - def ext(self):

214 return FieldComponent(self.component, False)

215 216 @property

217 - def avg(self):

218 return 0.5*(self.int+self.ext)

219

220 221 222 223 -class FluxVectorPlaceholder(object):

224 - def __init__(self, components=None, scalars=None):

225 if not (components is not None or scalars is not None): 226 raise ValueError, "either components or scalars must be specified" 227 if components is not None and scalars is not None: 228 raise ValueError, "only one of components and scalars may be specified" 229 230 # make them arrays for the better indexing 231 if components: 232 self.scalars = numpy.array([ 233 FluxScalarPlaceholder(i) 234 for i in range(components)]) 235 else: 236 self.scalars = numpy.array(scalars)

237

238 - def __len__(self):

239 return len(self.scalars)

240

241 - def __getitem__(self, idx):

242 if isinstance(idx, int): 243 return self.scalars[idx] 244 else: 245 return FluxVectorPlaceholder(scalars=self.scalars.__getitem__(idx))

246 247 @property

248 - def int(self):

249 return numpy.array([scalar.int for scalar in self.scalars])

250 251 @property

252 - def ext(self):

253 return numpy.array([scalar.ext for scalar in self.scalars])

254 255 @property

256 - def avg(self):

257 return numpy.array([scalar.avg for scalar in self.scalars])

258

259 260 261 262 # internal flux wrangling ----------------------------------------------------- 263 -class FluxIdentityMapperMixin(object):

264 - def map_field_component(self, expr):

265 return expr

266

267 - def map_normal(self, expr):

268 return expr

269

270 - def map_penalty_term(self, expr):

271 return expr.__class__(self.rec(expr.power))

272

273 - def map_if_positive(self, expr):

274 return expr.__class__( 275 self.rec(expr.criterion), 276 self.rec(expr.then), 277 self.rec(expr.else_), 278 )

279

280 - def map_scalar_parameter(self, expr):

281 return expr

282

283 284 285 286 -class FluxIdentityMapper( 287 pymbolic.mapper.IdentityMapper, 288 FluxIdentityMapperMixin):

289 pass

290

291 -class FluxSubstitutionMapper(pymbolic.mapper.substitutor.SubstitutionMapper, 292 FluxIdentityMapperMixin):

293 - def map_field_component(self, expr):

294 result = self.subst_func(expr) 295 if result is not None: 296 return result 297 else: 298 return expr

299 300 map_normal = map_field_component

301

302 303 304 305 -class FluxStringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):

306 - def map_field_component(self, expr, enclosing_prec):

307 if expr.is_local: 308 return "Int[%d]" % expr.index 309 else: 310 return "Ext[%d]" % expr.index

311

312 - def map_normal(self, expr, enclosing_prec):

313 return "Normal(%d)" % expr.axis

314

315 - def map_penalty_term(self, expr, enclosing_prec):

316 return "Penalty(%s)" % (expr.power)

317

318 - def map_if_positive(self, expr, enclosing_prec):

319 return "IfPositive(%s, %s, %s)" % (expr.criterion, expr.then, expr.else_)

320

321 322 323 324 -class FluxFlattenMapper(pymbolic.mapper.flattener.FlattenMapper, 325 FluxIdentityMapperMixin):

326 pass

327

328 329 330 331 332 -class FluxDependencyMapper(pymbolic.mapper.dependency.DependencyMapper):

333 - def map_field_component(self, expr):

334 return set([expr])

335

336 - def map_normal(self, expr):

337 return set()

338

339 - def map_penalty_term(self, expr):

340 return set()

341

342 - def map_if_positive(self, expr):

343 return self.rec(expr.criterion) | self.rec(expr.then) | self.rec(expr.else_)

344

345 - def map_scalar_parameter(self, expr):

346 return set([expr])

347

348 349 350 351 352 -class FluxTermCollector(pymbolic.mapper.collector.TermCollector, 353 FluxIdentityMapperMixin):

354 pass

355

356 357 358 359 -class FluxAllDependencyMapper(FluxDependencyMapper):

360 - def map_normal(self, expr):

361 return set([expr])

362

363 - def map_penalty_term(self, expr):

364 return set([expr])

365

366 367 368 369 -class FluxNormalizationMapper(pymbolic.mapper.collector.TermCollector, 370 FluxIdentityMapperMixin):

371 - def get_dependencies(self, expr):

372 return FluxAllDependencyMapper()(expr)

373

374 - def map_constant_flux(self, expr):

375 if expr.local_c == expr.neighbor_c: 376 return expr.local_c 377 else: 378 return expr

379

380 381 382 383 -class FluxCCFMapper(pymbolic.mapper.constant_folder.CommutativeConstantFoldingMapper, 384 FluxIdentityMapperMixin):

385 - def is_constant(self, expr):

386 return not bool(FluxAllDependencyMapper()(expr))

387

388 389 390 391 -class FluxExpandMapper(pymbolic.mapper.expander.ExpandMapper, 392 FluxIdentityMapperMixin):

393 - def __init__(self):

394 pymbolic.mapper.expander.ExpandMapper.__init__(self, 395 FluxNormalizationMapper())

396

397 398 399 400 -class FluxFlipper(FluxIdentityMapper):

401 - def map_normal(self, expr):

402 return -expr

403

404 - def map_field_component(self, expr):

405 return expr.__class__(expr.index, not expr.is_local)

406

407 408 409 410 411 -class FluxFlopCounter(pymbolic.mapper.flop_counter.FlopCounter):

412 - def map_penalty_term(self, expr):

413 return 0

414

415 - def map_normal(self, expr):

416 return 0

417

418 - def map_field_component(self, expr):

419 return 0

420

421 - def map_if_positive(self, expr):

422 return self.rec(expr.criterion) + max( 423 self.rec(expr.then), 424 self.rec(expr.else_))

425

426 - def map_function_symbol(self, expr):

427 return 1

428

429 - def map_scalar_parameter(self, expr):

430 return 0

431

432 433 434 435 436 -def normalize_flux(flux):

437 return FluxCCFMapper()(FluxNormalizationMapper()( 438 FluxFlattenMapper()(FluxExpandMapper()(flux))))

439

Source Code for Module hedge.flux