1 """Just-in-time compiling backend: Flux code generation."""
2
3 from __future__ import division
4
5 __copyright__ = "Copyright (C) 2008 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
25 from pymbolic.mapper.c_code import CCodeMapper
26 from hedge.flux import FluxIdentityMapper
27
28
29
30
31
34 self.flux_idx = flux_idx
35 self.flux_var_info = fvi
36
38 if expr.is_local:
39 where = "loc"
40 else:
41 where = "opp"
42
43 arg_name = self.flux_var_info.flux_idx_and_dep_to_arg_name[
44 self.flux_idx, expr]
45
46 if not arg_name:
47 return 0
48 else:
49 from pymbolic import var
50 return var(arg_name+"_it")[var(where+"_idx")]
51
53 from pymbolic import var
54 return var("args._scalar_arg_%d"
55 % self.flux_var_info.scalar_parameters.index(expr))
62 return "uncomplex_type(fp.loc.normal[%d])" % (expr.axis)
63
65 return ("uncomplex_type(pow(fp.loc.order*fp.loc.order/fp.loc.h, %(pwr)r))"
66 % {"pwr": expr.power})
67
79
81 import numpy
82 if isinstance(x, complex):
83 return "std::complex<uncomplex_type>(%s, %s)" % (
84 repr(x.real), repr(x.imag))
85 else:
86 return "uncomplex_type(%s)" % repr(x)
92 -def flux_to_code(f2c, is_flipped, flux_idx, fvi, flux, prec):
98
99
100
101
102
104 from pytools import Record
105 class FluxVariableInfo(Record):
106 pass
107
108 scalar_parameters = set()
109
110 fvi = FluxVariableInfo(
111 scalar_parameters=None,
112 arg_specs=[],
113 arg_names=[],
114 flux_idx_and_dep_to_arg_name={},
115 )
116
117 field_expr_to_arg_name = {}
118
119 from hedge.flux import \
120 FieldComponent, FluxDependencyMapper, \
121 FluxScalarParameter
122
123 from hedge.optemplate import BoundaryPair
124
125 for flux_idx, flux_binding in enumerate(fluxes):
126 for dep in FluxDependencyMapper(include_calls=False)(flux_binding.op.flux):
127 if isinstance(dep, FluxScalarParameter):
128 scalar_parameters.add(dep)
129 elif isinstance(dep, FieldComponent):
130 is_bdry = isinstance(flux_binding.field, BoundaryPair)
131 if is_bdry:
132 if dep.is_local:
133 this_field_expr = flux_binding.field.field
134 else:
135 this_field_expr = flux_binding.field.bfield
136 else:
137 this_field_expr = flux_binding.field
138
139 from hedge.tools import is_obj_array
140 if is_obj_array(this_field_expr):
141 fc_field_expr = this_field_expr[dep.index]
142 else:
143 assert dep.index == 0
144 fc_field_expr = this_field_expr
145
146 def set_or_check(dict_instance, key, value):
147 try:
148 existing_value = dict_instance[key]
149 except KeyError:
150 dict_instance[key] = value
151 else:
152 assert existing_value == value
153
154 from pymbolic.primitives import is_zero
155 if is_zero(fc_field_expr):
156 fvi.flux_idx_and_dep_to_arg_name[flux_idx, dep] = 0
157 else:
158 if fc_field_expr not in field_expr_to_arg_name:
159 arg_name = "arg%d" % len(fvi.arg_specs)
160 field_expr_to_arg_name[fc_field_expr] = arg_name
161
162 fvi.arg_names.append(arg_name)
163 fvi.arg_specs.append((fc_field_expr, dep.is_local))
164 else:
165 arg_name = field_expr_to_arg_name[fc_field_expr]
166
167 set_or_check(
168 fvi.flux_idx_and_dep_to_arg_name,
169 (flux_idx, dep),
170 arg_name)
171
172 if not is_bdry:
173
174
175
176 set_or_check(
177 fvi.flux_idx_and_dep_to_arg_name,
178 (flux_idx,
179 FieldComponent(dep.index, not dep.is_local)),
180 arg_name)
181 else:
182 raise ValueError("unknown flux dependency type: %s" % dep)
183
184 fvi.scalar_parameters = list(scalar_parameters)
185
186 return fvi
187
188
189
190
203
204
205
206
208 from codepy.cgen import \
209 FunctionDeclaration, FunctionBody, \
210 Const, Reference, Value, MaybeUnused, Typedef, POD, \
211 Statement, Include, Line, Block, Initializer, Assign, \
212 CustomLoop, For, Struct
213
214 from codepy.bpl import BoostPythonModule
215 mod = BoostPythonModule()
216
217 from pytools import to_uncomplex_dtype, flatten
218
219 S = Statement
220 mod.add_to_preamble([
221 Include("cstdlib"),
222 Include("algorithm"),
223 Line(),
224 Include("boost/foreach.hpp"),
225 Line(),
226 Include("hedge/face_operators.hpp"),
227 ])
228
229 mod.add_to_module([
230 S("using namespace hedge"),
231 S("using namespace pyublas"),
232 Line(),
233 Typedef(POD(dtype, "value_type")),
234 Typedef(POD(to_uncomplex_dtype(dtype), "uncomplex_type")),
235 Line(),
236 ])
237
238 arg_struct = Struct("arg_struct", [
239 Value("numpy_array<value_type>", "flux%d_on_faces" % i)
240 for i in range(len(fluxes))
241 ]+[
242 Value("numpy_array<value_type>", arg_name)
243 for arg_name in fvi.arg_names
244 ]+[
245 Value("value_type" if scalar_par.is_complex else "uncomplex_type",
246 "_scalar_arg_%d" % i)
247 for i, scalar_par in enumerate(fvi.scalar_parameters)
248 ])
249
250 mod.add_struct(arg_struct, "ArgStruct")
251 mod.add_to_module([Line()])
252
253 fdecl = FunctionDeclaration(
254 Value("void", "gather_flux"),
255 [
256 Const(Reference(Value("face_group", "fg"))),
257 Reference(Value("arg_struct", "args"))
258 ])
259
260 from pymbolic.mapper.stringifier import PREC_PRODUCT
261
262 def gen_flux_code():
263 f2cm = FluxToCodeMapper(repr, reverse=False)
264
265 result = [
266 Assign("fof%d_it[%s_fof_base+%s]" % (flux_idx, where, tgt_idx),
267 "uncomplex_type(fp.loc.face_jacobian) * " +
268 flux_to_code(f2cm, is_flipped, flux_idx, fvi, flux.op.flux, PREC_PRODUCT))
269 for flux_idx, flux in enumerate(fluxes)
270 for where, is_flipped, tgt_idx in [
271 ("loc", False, "i"),
272 ("opp", True, "opp_write_map[i]")
273 ]]
274
275 return [
276 Initializer(Value("value_type", f2cm.cse_prefix+str(i)), cse)
277 for i, cse in f2cm.cses] + result
278
279 fbody = Block([
280 Initializer(
281 Const(Value("numpy_array<value_type>::iterator", "fof%d_it" % i)),
282 "args.flux%d_on_faces.begin()" % i)
283 for i in range(len(fluxes))
284 ]+[
285 Initializer(
286 Const(Value("numpy_array<value_type>::const_iterator", "%s_it" % arg_name)),
287 "args.%s.begin()" % arg_name)
288 for arg_name in fvi.arg_names
289 ]+[
290 Line(),
291 CustomLoop("BOOST_FOREACH(const face_pair &fp, fg.face_pairs)", Block(
292 list(flatten([
293 Initializer(Value("node_number_t", "%s_ebi" % where),
294 "fp.%s.el_base_index" % where),
295 Initializer(Value("index_lists_t::const_iterator", "%s_idx_list" % where),
296 "fg.index_list(fp.%s.face_index_list_number)" % where),
297 Initializer(Value("node_number_t", "%s_fof_base" % where),
298 "fg.face_length()*(fp.%(where)s.local_el_number*fg.face_count"
299 " + fp.%(where)s.face_id)" % {"where": where}),
300 Line(),
301 ]
302 for where in ["loc", "opp"]
303 ))+[
304 Initializer(Value("index_lists_t::const_iterator", "opp_write_map"),
305 "fg.index_list(fp.opp_native_write_map)"),
306 Line(),
307 For(
308 "unsigned i = 0",
309 "i < fg.face_length()",
310 "++i",
311 Block(
312 [
313 Initializer(MaybeUnused(Value("node_number_t", "%s_idx" % where)),
314 "%(where)s_ebi + %(where)s_idx_list[i]"
315 % {"where": where})
316 for where in ["loc", "opp"]
317 ]+gen_flux_code()
318 )
319 )
320 ]))
321 ])
322 mod.add_function(FunctionBody(fdecl, fbody))
323
324
325
326
327
328 return mod.compile(get_flux_toolchain(discr, fluxes),
329 wait_on_error="jit_wait_on_compile_error" in discr.debug)
330
331
332
333
335 from codepy.cgen import \
336 FunctionDeclaration, FunctionBody, Typedef, Struct, \
337 Const, Reference, Value, POD, MaybeUnused, \
338 Statement, Include, Line, Block, Initializer, Assign, \
339 CustomLoop, For
340
341 from pytools import to_uncomplex_dtype, flatten
342
343 from codepy.bpl import BoostPythonModule
344 mod = BoostPythonModule()
345
346 mod.add_to_preamble([
347 Include("cstdlib"),
348 Include("algorithm"),
349 Line(),
350 Include("boost/foreach.hpp"),
351 Line(),
352 Include("hedge/face_operators.hpp"),
353 ])
354
355 S = Statement
356 mod.add_to_module([
357 S("using namespace hedge"),
358 S("using namespace pyublas"),
359 Line(),
360 Typedef(POD(dtype, "value_type")),
361 Typedef(POD(to_uncomplex_dtype(dtype), "uncomplex_type")),
362 ])
363
364 arg_struct = Struct("arg_struct", [
365 Value("numpy_array<value_type>", "flux%d_on_faces" % i)
366 for i in range(len(fluxes))
367 ]+[
368 Value("numpy_array<value_type>", arg_name)
369 for arg_name in fvi.arg_names
370 ])
371
372 mod.add_struct(arg_struct, "ArgStruct")
373 mod.add_to_module([Line()])
374
375 fdecl = FunctionDeclaration(
376 Value("void", "gather_flux"),
377 [
378 Const(Reference(Value("face_group", "fg"))),
379 Reference(Value("arg_struct", "args"))
380 ])
381
382 from pymbolic.mapper.stringifier import PREC_PRODUCT
383
384 def gen_flux_code():
385 f2cm = FluxToCodeMapper(repr, reverse=False)
386
387 result = [
388 Assign("fof%d_it[loc_fof_base+i]" % flux_idx,
389 "uncomplex_type(fp.loc.face_jacobian) * " +
390 flux_to_code(f2cm, False, flux_idx, fvi, flux.op.flux, PREC_PRODUCT))
391 for flux_idx, flux in enumerate(fluxes)
392 ]
393
394 return [
395 Initializer(Value("value_type", f2cm.cse_prefix+str(i)), cse)
396 for i, cse in enumerate(f2cm.cses)] + result
397
398 fbody = Block([
399 Initializer(
400 Const(Value("numpy_array<value_type>::iterator", "fof%d_it" % i)),
401 "args.flux%d_on_faces.begin()" % i)
402 for i in range(len(fluxes))
403 ]+[
404 Initializer(
405 Const(Value("numpy_array<value_type>::const_iterator",
406 "%s_it" % arg_name)),
407 "args.%s.begin()" % arg_name)
408 for arg_name in fvi.arg_names
409 ]+[
410 Line(),
411 CustomLoop("BOOST_FOREACH(const face_pair &fp, fg.face_pairs)", Block(
412 list(flatten([
413 Initializer(Value("node_number_t", "%s_ebi" % where),
414 "fp.%s.el_base_index" % where),
415 Initializer(Value("index_lists_t::const_iterator", "%s_idx_list" % where),
416 "fg.index_list(fp.%s.face_index_list_number)" % where),
417 Line(),
418 ]
419 for where in ["loc", "opp"]
420 ))+[
421 Line(),
422 Initializer(Value("node_number_t", "loc_fof_base"),
423 "fg.face_length()*(fp.%(where)s.local_el_number*fg.face_count"
424 " + fp.%(where)s.face_id)" % {"where": "loc"}),
425 Line(),
426 For(
427 "unsigned i = 0",
428 "i < fg.face_length()",
429 "++i",
430 Block(
431 [
432 Initializer(MaybeUnused(
433 Value("node_number_t", "%s_idx" % where)),
434 "%(where)s_ebi + %(where)s_idx_list[i]"
435 % {"where": where})
436 for where in ["loc", "opp"]
437 ]+gen_flux_code()
438 )
439 )
440 ]))
441 ])
442
443 mod.add_function(FunctionBody(fdecl, fbody))
444
445
446
447
448
449 return mod.compile(get_flux_toolchain(discr, fluxes),
450 wait_on_error="jit_wait_on_compile_error" in discr.debug)
451