From efbb7df5872d22087588c9704e0fbba058aa85b8 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 21 Oct 2026 11:00:04 +0100 Subject: [PATCH 01/19] compiler: Preserve CSE scheduling independently of temporary names --- devito/passes/clusters/cse.py | 12 +++++------- tests/test_cse.py | 19 ++++++++++++++++++- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/devito/passes/clusters/cse.py b/devito/passes/clusters/cse.py index 9c899eb23e..75fe68d636 100644 --- a/devito/passes/clusters/cse.py +++ b/devito/passes/clusters/cse.py @@ -293,15 +293,13 @@ def _toposort(exprs): if e0.lhs in e1.rhs.free_symbols: dag.add_edge(e0, e1, force_add=True) + ordering = {e: n for n, e in enumerate(exprs)} + def choose_element(queue, scheduled): + # Preserve expression order among ready temporaries, independent of names tmps = [i for i in queue if isinstance(i.lhs, CTemp)] - if tmps: - # Try to honor temporary names as much as possible - first = sorted(tmps, key=lambda i: i.lhs.name).pop(0) - queue.remove(first) - else: - first = sorted(queue, key=lambda i: exprs.index(i)).pop(0) - queue.remove(first) + first = min(tmps or queue, key=ordering.__getitem__) + queue.remove(first) return first processed = dag.topological_sort(choose_element) diff --git a/tests/test_cse.py b/tests/test_cse.py index 6ab68f05a0..4fc439c49b 100644 --- a/tests/test_cse.py +++ b/tests/test_cse.py @@ -8,7 +8,7 @@ ConditionalDimension, Eq, Function, Grid, Operator, TimeFunction, cos, sin ) from devito.finite_differences.differentiable import diffify -from devito.ir import Conditional, DummyEq, FindNodes, FindSymbols +from devito.ir import Conditional, DummyEq, Expression, FindNodes, FindSymbols from devito.ir.support import generator from devito.passes.clusters.cse import CTemp, _cse from devito.symbolics import DefFunction, FieldFromPointer, IndexedPointer, indexify @@ -131,6 +131,23 @@ def test_temp_order(): assert type(args[2]) is CTemp +@pytest.mark.parametrize('names', [('r2', 'r3'), ('r9', 'r10'), ('later', 'earlier')]) +@pytest.mark.parametrize('mode', ['smartsort', 'advanced']) +def test_schedule_independent_of_temp_names(names, mode): + grid = Grid(shape=(4, 4)) + f = Function(name='f', grid=grid) + g = Function(name='g', grid=grid) + u = Function(name='u', grid=grid) + v = Function(name='v', grid=grid) + a, b = [CTemp(name=name, dtype=np.float32) for name in names] + + op = Operator([Eq(a, cos(f)), Eq(b, sin(g)), Eq(u, a + b), Eq(v, a*b)], + opt=('advanced', {'cse-algo': mode})) + + exprs = [e.expr for e in FindNodes(Expression).visit(op)] + assert [e.lhs for e in exprs if isinstance(e.lhs, CTemp)][:2] == [a, b] + + def test_deffunction_not_captured(): """ Opaque C-level calls (DefFunction) must not be CSE-captured: their return From b0b82c37833d0ffc29ac40cea3ab92f3f5b4330d Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 21 Oct 2026 11:14:38 +0100 Subject: [PATCH 02/19] compiler: Preserve lowered arithmetic during substitutions --- devito/ir/cgen/printer.py | 2 +- devito/symbolics/manipulation.py | 14 ++++++++++++-- tests/test_symbolics.py | 13 ++++++++++++- 3 files changed, 25 insertions(+), 4 deletions(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index 27ee0e3112..adfc7651aa 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -289,7 +289,7 @@ def _print_Mul(self, expr): term = self.parenthesize(args[0], precedence(expr)) if neg: - return f'-{term}' + return term[1:] if term.startswith('-') else f'-{term}' else: return term diff --git a/devito/symbolics/manipulation.py b/devito/symbolics/manipulation.py index b545290a15..c1e2ec4169 100644 --- a/devito/symbolics/manipulation.py +++ b/devito/symbolics/manipulation.py @@ -53,8 +53,9 @@ def uxreplace(expr, rule): be avoided at all costs (e.g., to prevent SymPy from unpicking Devito transformations, such as factorization). - The canonical ordering of the arguments is however guaranteed; where this is - not possible, a re-evaluation will be enforced. + Evaluated expressions retain canonical argument ordering. Lowered, + unevaluated arithmetic retains its operand order and grouping unless + substitution makes it completely numeric. By avoiding re-evaluations, this function is typically much quicker than SymPy's xreplace. @@ -191,6 +192,15 @@ def _(expr, args, kwargs): return expr.func(*args, evaluate=evaluate) +@_uxreplace_handle.register(UnevalAdd) +@_uxreplace_handle.register(UnevalMul) +def _(expr, args, kwargs): + # Lowered arithmetic must survive substitutions without reassociation + if all(i.is_Number for i in args): + return expr.func.__base__(*args) + return expr.func(*args) + + @_uxreplace_handle.register(Add) def _(expr, args, kwargs): args = [i for i in args if i != 0] diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index 51a759a705..db96c6a842 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -21,9 +21,10 @@ retrieve_derivatives, retrieve_functions, retrieve_indexed, subs_if_composite, uxreplace, xreplace_indices ) +from devito.symbolics.unevaluation import Add as UnevalAdd, Mul as UnevalMul from devito.tools import CustomDtype, as_tuple, dtypes_vector_mapper from devito.types import ( - Array, Bundle, ComponentAccess, FIndexed, LocalObject, Object, StencilDimension + Array, Bundle, ComponentAccess, FIndexed, LocalObject, Object, StencilDimension, Temp ) from devito.types import Symbol as dSymbol from devito.types.basic import AbstractSymbol @@ -1123,6 +1124,16 @@ def test_halo_touch(self): assert ht1.args == (g[x, y],) + @pytest.mark.parametrize('op', [UnevalAdd, UnevalMul]) + def test_lowered_arithmetic_order(self, op): + a, b, c = [Temp(name=f'r{i}') for i in (9, 10, 11)] + x, y, z = [Temp(name=f'r{i}') for i in (0, 1, 2)] + expr = op(c, op(b, a)) + + rebuilt = uxreplace(expr, {a: x, b: y, c: z}) + + assert rebuilt.args == (z, op(y, x)) + def test_minmax(): grid = Grid(shape=(5, 5)) From c2c8f7b28243318563465b7fd4499312d3df18dd Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 21 Oct 2026 11:14:38 +0100 Subject: [PATCH 03/19] compiler: Substitute dynamic shared-memory dereference offsets --- devito/ir/iet/visitors.py | 3 ++- tests/test_iet.py | 16 +++++++++++++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/devito/ir/iet/visitors.py b/devito/ir/iet/visitors.py index 6186ff4e76..8a825c7619 100644 --- a/devito/ir/iet/visitors.py +++ b/devito/ir/iet/visitors.py @@ -1506,7 +1506,8 @@ def visit_PointerCast(self, o): def visit_Dereference(self, o): pointee = self.mapper.get(o.pointee, o.pointee) pointer = self.mapper.get(o.pointer, o.pointer) - return reuse_if_unchanged(o, pointee=pointee, pointer=pointer) + offset = uxreplace(o.offset, self.mapper) + return reuse_if_unchanged(o, pointee=pointee, pointer=pointer, offset=offset) def visit_Pragma(self, o): arguments = [uxreplace(i, self.mapper) for i in o.arguments] diff --git a/tests/test_iet.py b/tests/test_iet.py index 99a94ed24c..d1e97009ed 100644 --- a/tests/test_iet.py +++ b/tests/test_iet.py @@ -12,7 +12,7 @@ from devito.ir.iet import ( Call, Callable, CGen, Conditional, Definition, Dereference, DeviceCall, DummyExpr, ElementalFunction, FindNodes, FindSymbols, Iteration, KernelLaunch, Lambda, List, - Switch, Transformer, filter_iterations, make_callable, make_efunc, + Switch, Transformer, Uxreplace, filter_iterations, make_callable, make_efunc, retrieve_iteration_tree ) from devito.ir.iet.visitors import sorted_efuncs @@ -120,6 +120,20 @@ def test_nested_calls_cgen(): assert str(code) == 'foo(bar());' +def test_uxreplace_dereference_offset(): + d = CustomDimension(name='d', symbolic_size=4) + array = Array(name='a0', dimensions=(d,), shape=(4,), scope='shared') + pointer = Pointer(name='sm7') + old, new = Symbol(name='off7'), Symbol(name='off0') + node = Dereference(array, pointer, offset=old + 32) + + replaced = Uxreplace({old: new}).visit(node) + + assert replaced.offset == new + 32 + assert replaced.pointee is array + assert replaced.pointer is pointer + + @pytest.mark.parametrize('mode,expected', [ ('basics', '["x"]'), ('symbolics', '["f"]') From 6113dc541ee27c347207c656a7248bde5935cdce Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 21 Oct 2026 11:14:39 +0100 Subject: [PATCH 04/19] compiler: Finalize IET symbols before extracting macros --- devito/operator/operator.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/devito/operator/operator.py b/devito/operator/operator.py index 8f3c02f983..3ec62fbdae 100644 --- a/devito/operator/operator.py +++ b/devito/operator/operator.py @@ -472,6 +472,11 @@ def _specialize_iet(cls, graph, **kwargs): """ return graph + @classmethod + def _finalize_iet(cls, graph, **kwargs): + """Backend hook after lowering, before macro and signature extraction.""" + return graph + @classmethod @timed_pass(name='lowering.uIET') def _lower_uiet(cls, stree, profiler=None, **kwargs): @@ -515,15 +520,17 @@ def _lower_iet(cls, uiet, **kwargs): # specialization further Sections may be introduced cls._Target.instrument(graph, **kwargs) - # Extract the necessary macros from the symbolic objects - generate_macros(graph, **kwargs) - # Target-specific lowering lower_dtypes(graph, **kwargs) # Target-independent optimizations minimize_symbols(graph) + graph = cls._finalize_iet(graph, **kwargs) + + # Extract macros only after the symbolic objects have their final names + generate_macros(graph, **kwargs) + # Finalize helper signatures after all IET transformations have settled. finalize_args(graph) From 91ee6a00151335f514f4f72a930050e08c269e45 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 21 Oct 2026 16:39:45 +0100 Subject: [PATCH 05/19] compiler: Avoid generated-name collisions with application symbols --- devito/ir/support/symregistry.py | 10 ++++++++-- devito/operator/operator.py | 9 ++++++++- tests/test_operator.py | 13 +++++++++++++ 3 files changed, 29 insertions(+), 3 deletions(-) diff --git a/devito/ir/support/symregistry.py b/devito/ir/support/symregistry.py index d7fca20f7c..552034cd61 100644 --- a/devito/ir/support/symregistry.py +++ b/devito/ir/support/symregistry.py @@ -15,6 +15,7 @@ class SymbolRegistry: def __init__(self): # {name -> generator()} -- to create unique names for symbols, functions, ... self.counters = {} + self.reserved = set() # Special symbols self.nthreads = NThreads() @@ -42,10 +43,15 @@ def make_name(self, prefix=None, increment_first=True): counter = self.counters.setdefault(prefix, generator()) # Only increment symbol names after the first encountered - if not increment_first: + if not increment_first and prefix not in self.reserved: + self.reserved.add(prefix) return prefix - return f'{prefix}{counter()}' + name = f'{prefix}{counter()}' + while name in self.reserved: + name = f'{prefix}{counter()}' + self.reserved.add(name) + return name def make_npthreads(self, size): name = self.make_name(prefix='npthreads') diff --git a/devito/operator/operator.py b/devito/operator/operator.py index 3ec62fbdae..c4a7bc0394 100644 --- a/devito/operator/operator.py +++ b/devito/operator/operator.py @@ -36,7 +36,9 @@ lower_dtypes, lower_index_derivatives, lower_local_sums, minimize_symbols, optimize_pows, unevaluate ) -from devito.symbolics import estimate_cost, subs_op_args +from devito.symbolics import ( + estimate_cost, retrieve_functions, retrieve_symbols, subs_op_args +) from devito.tools import ( DAG, CacheInstances, MemoryEstimate, OrderedSet, ReducerMap, Signer, as_mapper, as_tuple, contains_val, filter_sorted, flatten, frozendict, is_integer, memoized_func, @@ -293,6 +295,11 @@ def _lower(cls, expressions, **kwargs): expressions = as_tuple(expressions) + # Protect application identifiers before compiler temporaries are introduced + sregistry = kwargs['sregistry'] + sregistry.reserved.update(i.name for i in retrieve_functions(expressions)) + sregistry.reserved.update(i.name for i in retrieve_symbols(expressions)) + # Enable recursive lowering # This may be used by a compilation pass that constructs a new # expression for which a partial or complete lowering is desired diff --git a/tests/test_operator.py b/tests/test_operator.py index e6fc86463a..8b34a93bc9 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -198,6 +198,19 @@ def test_jit_compiler_flags(self, envvar, flags, monkeypatch): assert op0._soname != op1._soname assert op1._soname == op2._soname + @pytest.mark.parametrize('name', ['r0', 'r1']) + def test_generated_names_avoid_inputs(self, name): + grid = Grid(shape=(8, 8), extent=(7., 7.)) + f = Function(name=name, grid=grid, space_order=4) + u = Function(name='u', grid=grid) + x, y = np.indices(f.data_with_halo.shape, dtype=np.float32) + f.data_with_halo[:] = x + 2*y + + op = Operator(Eq(u, f.dx + f.dy), name='names_inputs') + op.apply() + + np.testing.assert_allclose(u.data, 3., rtol=1.e-5, atol=1.e-6) + class TestCodeGen: From f1c7ddc20ae5424624272ef1f35e17a09566bbab Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 21 Oct 2026 16:50:17 +0100 Subject: [PATCH 06/19] compiler: Preserve scoped bindings during symbolic substitutions --- devito/passes/iet/misc.py | 5 +++-- devito/symbolics/manipulation.py | 24 +++++++++++++++++++++-- tests/test_symbolics.py | 33 ++++++++++++++++++++++++++++++-- 3 files changed, 56 insertions(+), 6 deletions(-) diff --git a/devito/passes/iet/misc.py b/devito/passes/iet/misc.py index 1b404daca1..a94a754124 100644 --- a/devito/passes/iet/misc.py +++ b/devito/passes/iet/misc.py @@ -188,8 +188,9 @@ def _generate_macros_findexeds(iet, sregistry=None, tracker=None, **kwargs): subs = {} for i in findexeds: + key = (i.base, i.strides_map) try: - v = tracker[i.base].v + v = tracker[key].v subs[i] = v.func(v.base, *i.indices) continue except KeyError: @@ -199,7 +200,7 @@ def _generate_macros_findexeds(iet, sregistry=None, tracker=None, **kwargs): header, v = i.bind(pname) subs[i] = v - tracker[i.base] = Bunch(header=header, v=v) + tracker[key] = Bunch(header=header, v=v) iet = Uxreplace(subs).visit(iet) diff --git a/devito/symbolics/manipulation.py b/devito/symbolics/manipulation.py index c1e2ec4169..82fffcfb9a 100644 --- a/devito/symbolics/manipulation.py +++ b/devito/symbolics/manipulation.py @@ -12,7 +12,7 @@ EvalDerivative, IndexDerivative, LocalSum ) from devito.symbolics.extended_dtypes import LONG -from devito.symbolics.extended_sympy import DefFunction, rfunc +from devito.symbolics.extended_sympy import CallFromPointer, DefFunction, rfunc from devito.symbolics.queries import q_leaf from devito.symbolics.search import ( retrieve_functions, retrieve_indexed, retrieve_terminals @@ -22,11 +22,12 @@ from devito.symbolics.unevaluation import Pow as UnevalPow from devito.symbolics.unevaluation import UnevaluableMixin from devito.tools import ( - EnrichedTuple, as_list, as_tuple, flatten, split, transitive_closure + EnrichedTuple, as_list, as_tuple, flatten, frozendict, split, transitive_closure ) from devito.types.array import ComponentAccess from devito.types.basic import Basic, Indexed from devito.types.equation import Eq +from devito.types.misc import FIndexed from devito.types.relational import Ge, Gt, Le, Lt __all__ = [ @@ -118,6 +119,17 @@ def _uxreplace(expr, rule): return expr, False +@_uxreplace.register(CallFromPointer) +def _(expr, rule): + if expr in rule or expr.__class__ in rule: + return _uxreplace.dispatch(object)(expr, rule) + + args, changed = _uxreplace_dispatch((expr.pointer, expr.params), rule) + if changed: + return expr.func(expr.call, *args), True + return expr, False + + @_uxreplace.register(LocalSum) def _(expr, rule): # Only the sum owns these guards; ordinary ConditionalDimensions stay atomic @@ -166,6 +178,7 @@ def _(iterable, rule): @_uxreplace_dispatch.register(dict) +@_uxreplace_dispatch.register(frozendict) def _(mapper, rule): ret = {} changed = False @@ -184,6 +197,12 @@ def _uxreplace_handle(expr, args, kwargs): return expr.func(*args) +@_uxreplace_handle.register(CallFromPointer) +def _(expr, args, kwargs): + # C member names belong to their compound type, rather than the local scope + return expr.func(expr.call, *args[1:]) + + @_uxreplace_handle.register(Min) @_uxreplace_handle.register(Max) @_uxreplace_handle.register(Pow) @@ -262,6 +281,7 @@ def dispatchable(self, obj): _uxreplace_registry.register(Eq) _uxreplace_registry.register(DefFunction) _uxreplace_registry.register(ComponentAccess) +_uxreplace_registry.register(FIndexed) _uxreplace_registry.register(LocalSum) diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index db96c6a842..25ef9b2cac 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -15,8 +15,9 @@ from devito.ir.support.guards import GuardExpr, pairwise_or, simplify_and from devito.mpi.halo_scheme import HaloTouch from devito.symbolics import ( # noqa - INT, BaseCast, CallFromPointer, Cast, DefFunction, FieldFromComposite, - FieldFromPointer, IntDiv, ListInitializer, Namespace, ReservedWord, RoundUp, Rvalue, + INT, BaseCast, CallFromComposite, CallFromPointer, Cast, DefFunction, + FieldFromComposite, FieldFromPointer, IntDiv, ListInitializer, Namespace, + ReservedWord, RoundUp, Rvalue, SizeOf, VectorAccess, evalrel, pow_to_mul, q_negative, q_positive, retrieve_derivatives, retrieve_functions, retrieve_indexed, subs_if_composite, uxreplace, xreplace_indices @@ -1134,6 +1135,34 @@ def test_lowered_arithmetic_order(self, op): assert rebuilt.args == (z, op(y, x)) + @pytest.mark.parametrize('cls,converted', [(CallFromPointer, CallFromComposite), + (CallFromComposite, CallFromPointer), + (FieldFromPointer, FieldFromComposite), + (FieldFromComposite, FieldFromPointer)]) + def test_bound_member_names(self, cls, converted): + a, b = [dSymbol(name=f'r{i}') for i in (9, 0)] + p, q = [dSymbol(name=f'p{i}') for i in (9, 0)] + expr = cls(a, p, [a]) + + assert uxreplace(expr, {a: b, p: q}) == cls(a, q, [b]) + assert uxreplace(expr, {expr: cls(b, q, [b])}) == cls(b, q, [b]) + assert uxreplace(expr, {cls: converted, a: b, p: q}) == converted(a, q, [b]) + + bound = cls(a, p) + assert uxreplace(bound, {a: b}) is bound + + def test_findexed_strides(self): + grid = Grid(shape=(4, 4)) + x, y = grid.dimensions + f = Function(name='f', grid=grid) + a, b = [dSymbol(name=f'y_stride{i}') for i in (9, 0)] + expr = FIndexed(f.indexed, x, y, strides_map={y: a}) + + rebuilt = uxreplace(expr, {a: b}) + + assert rebuilt.strides == (b,) + assert expr.strides == (a,) + def test_minmax(): grid = Grid(shape=(5, 5)) From 2390a7e2e2a430e1e08b5cf51add67710c94488f Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 21 Oct 2026 16:50:17 +0100 Subject: [PATCH 07/19] compiler: Abstract equivalent symbol instances once --- devito/passes/iet/engine.py | 2 ++ tests/test_iet.py | 11 +++++++---- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/devito/passes/iet/engine.py b/devito/passes/iet/engine.py index fe9a78b65b..8615582da5 100644 --- a/devito/passes/iet/engine.py +++ b/devito/passes/iet/engine.py @@ -564,6 +564,8 @@ def abstract_objects(objects0): mapper = {} sregistry = SymbolRegistry() for i in objects: + if i in mapper: + continue abstract_object(i, mapper, sregistry) return mapper diff --git a/tests/test_iet.py b/tests/test_iet.py index d1e97009ed..a1af790f77 100644 --- a/tests/test_iet.py +++ b/tests/test_iet.py @@ -556,14 +556,17 @@ def test_codegen_quality0(): assert foo1.parameters[0] is a -def test_reuse_efuncs_natural_numbering(): +@pytest.mark.parametrize('copy_symbols', [False, True]) +def test_reuse_efuncs_natural_numbering(copy_symbols): - def make_foo(name, start): + def make_foo(name, start, copies=False): pointers = [Pointer(name=f'p{i}') for i in range(start, start + 3)] - return make_callable(name, Call('bar', pointers)) + aliases = [p._rebuild() for p in pointers] if copies else pointers + return make_callable(name, List(body=[Call('bar', pointers), + Call('bar', aliases)])) foo0 = make_foo('foo0', 0) - foo1 = make_foo('foo1', 9) + foo1 = make_foo('foo1', 9, copies=copy_symbols) root = make_callable('root', List(body=[ Call(foo0.name, foo0.parameters), From 423877879d8b6d96e01365fcaee72efbc787c08c Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 22 Oct 2026 11:05:55 +0100 Subject: [PATCH 08/19] compiler: Encapsulate input identifier reservation --- devito/ir/support/symregistry.py | 11 ++++++++++- devito/operator/operator.py | 9 ++------- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/devito/ir/support/symregistry.py b/devito/ir/support/symregistry.py index 552034cd61..01ceaf3bdf 100644 --- a/devito/ir/support/symregistry.py +++ b/devito/ir/support/symregistry.py @@ -1,4 +1,5 @@ -from devito.tools import generator +from devito.symbolics import retrieve_functions +from devito.tools import as_tuple, generator from devito.types import ( DeviceID, NPThreads, NThreads, NThreadsNested, NThreadsNonaffine, ThreadID ) @@ -32,6 +33,14 @@ def __init__(self): # passes, to maximize symbol (especially Dimension) reuse self.caches = {} + def reserve(self, exprs): + """Reserve identifiers already used by the input expressions.""" + exprs = as_tuple(exprs) + self.reserved.update(i.name for e in exprs for i in e.free_symbols) + + # Before indexification, a Function's free symbols omit its own name + self.reserved.update(i.name for i in retrieve_functions(exprs)) + def make_name(self, prefix=None, increment_first=True): # By default we're creating a new symbol if prefix is None: diff --git a/devito/operator/operator.py b/devito/operator/operator.py index c4a7bc0394..e36c48a351 100644 --- a/devito/operator/operator.py +++ b/devito/operator/operator.py @@ -36,9 +36,7 @@ lower_dtypes, lower_index_derivatives, lower_local_sums, minimize_symbols, optimize_pows, unevaluate ) -from devito.symbolics import ( - estimate_cost, retrieve_functions, retrieve_symbols, subs_op_args -) +from devito.symbolics import estimate_cost, subs_op_args from devito.tools import ( DAG, CacheInstances, MemoryEstimate, OrderedSet, ReducerMap, Signer, as_mapper, as_tuple, contains_val, filter_sorted, flatten, frozendict, is_integer, memoized_func, @@ -295,10 +293,7 @@ def _lower(cls, expressions, **kwargs): expressions = as_tuple(expressions) - # Protect application identifiers before compiler temporaries are introduced - sregistry = kwargs['sregistry'] - sregistry.reserved.update(i.name for i in retrieve_functions(expressions)) - sregistry.reserved.update(i.name for i in retrieve_symbols(expressions)) + kwargs['sregistry'].reserve(expressions) # Enable recursive lowering # This may be used by a compilation pass that constructs a new From 19a7c03b57b8e09cb941b13adb320f584dc88d6d Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 22 Oct 2026 14:14:01 +0100 Subject: [PATCH 09/19] tests: Fix import sorting and continuation indentation --- tests/test_iet.py | 2 +- tests/test_symbolics.py | 19 ++++++++++--------- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/tests/test_iet.py b/tests/test_iet.py index a1af790f77..2834ac2712 100644 --- a/tests/test_iet.py +++ b/tests/test_iet.py @@ -563,7 +563,7 @@ def make_foo(name, start, copies=False): pointers = [Pointer(name=f'p{i}') for i in range(start, start + 3)] aliases = [p._rebuild() for p in pointers] if copies else pointers return make_callable(name, List(body=[Call('bar', pointers), - Call('bar', aliases)])) + Call('bar', aliases)])) foo0 = make_foo('foo0', 0) foo1 = make_foo('foo1', 9, copies=copy_symbols) diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index 25ef9b2cac..23c74b5525 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -17,17 +17,18 @@ from devito.symbolics import ( # noqa INT, BaseCast, CallFromComposite, CallFromPointer, Cast, DefFunction, FieldFromComposite, FieldFromPointer, IntDiv, ListInitializer, Namespace, - ReservedWord, RoundUp, Rvalue, - SizeOf, VectorAccess, evalrel, pow_to_mul, q_negative, q_positive, - retrieve_derivatives, retrieve_functions, retrieve_indexed, subs_if_composite, - uxreplace, xreplace_indices + ReservedWord, RoundUp, Rvalue, SizeOf, VectorAccess, evalrel, pow_to_mul, q_negative, + q_positive, retrieve_derivatives, retrieve_functions, retrieve_indexed, + subs_if_composite, uxreplace, xreplace_indices ) -from devito.symbolics.unevaluation import Add as UnevalAdd, Mul as UnevalMul +from devito.symbolics.unevaluation import Add as UnevalAdd +from devito.symbolics.unevaluation import Mul as UnevalMul from devito.tools import CustomDtype, as_tuple, dtypes_vector_mapper from devito.types import ( - Array, Bundle, ComponentAccess, FIndexed, LocalObject, Object, StencilDimension, Temp + Array, Bundle, ComponentAccess, FIndexed, LocalObject, Object, StencilDimension ) from devito.types import Symbol as dSymbol +from devito.types import Temp from devito.types.basic import AbstractSymbol @@ -1136,9 +1137,9 @@ def test_lowered_arithmetic_order(self, op): assert rebuilt.args == (z, op(y, x)) @pytest.mark.parametrize('cls,converted', [(CallFromPointer, CallFromComposite), - (CallFromComposite, CallFromPointer), - (FieldFromPointer, FieldFromComposite), - (FieldFromComposite, FieldFromPointer)]) + (CallFromComposite, CallFromPointer), + (FieldFromPointer, FieldFromComposite), + (FieldFromComposite, FieldFromPointer)]) def test_bound_member_names(self, cls, converted): a, b = [dSymbol(name=f'r{i}') for i in (9, 0)] p, q = [dSymbol(name=f'p{i}') for i in (9, 0)] From 6ca953e3de67b8938610f0b8768524ee8e2f6c84 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 22 Oct 2026 14:52:11 +0100 Subject: [PATCH 10/19] docs: Refresh generated code in GPU tutorial --- examples/performance/01_gpu.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/performance/01_gpu.ipynb b/examples/performance/01_gpu.ipynb index dc096aa0eb..d37e956cd6 100644 --- a/examples/performance/01_gpu.ipynb +++ b/examples/performance/01_gpu.ipynb @@ -309,7 +309,7 @@ " {\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " uL0(time + 1, x + 2, y + 2) = dt*(c*(r1*uL0(time, x + 1, y + 2) + r1*uL0(time, x + 3, y + 2) + r2*uL0(time, x + 2, y + 1) + r2*uL0(time, x + 2, y + 3) - 2.0F*(r1*uL0(time, x + 2, y + 2) + r2*uL0(time, x + 2, y + 2))) + r0*uL0(time, x + 2, y + 2));\n", + " uL0(time + 1, x + 2, y + 2) = dt*(c*(r1*uL0(time, x + 1, y + 2) + r1*uL0(time, x + 3, y + 2) + r2*uL0(time, x + 2, y + 1) + r2*uL0(time, x + 2, y + 3) + (-2.0F)*(r1*uL0(time, x + 2, y + 2) + r2*uL0(time, x + 2, y + 2))) + r0*uL0(time, x + 2, y + 2));\n", " }\n", " }\n", " STOP(section0,timers)\n", From b2c4e178749a263a82c99488a658dd7e2e55c3f4 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 22 Oct 2026 17:44:13 +0100 Subject: [PATCH 11/19] compiler: Preserve lowered rounding functions during reconstruction --- devito/symbolics/manipulation.py | 9 +++++++-- devito/symbolics/unevaluation.py | 10 +++++++++- tests/test_pickle.py | 7 ++++--- 3 files changed, 20 insertions(+), 6 deletions(-) diff --git a/devito/symbolics/manipulation.py b/devito/symbolics/manipulation.py index 82fffcfb9a..f1a7524ffb 100644 --- a/devito/symbolics/manipulation.py +++ b/devito/symbolics/manipulation.py @@ -3,7 +3,7 @@ from functools import singledispatch import numpy as np -from sympy import Add, Max, Min, Mul, Pow, S, SympifyError, Tuple, sympify +from sympy import Add, Max, Min, Mul, Pow, S, SympifyError, Tuple, ceiling, floor, sympify from sympy.core import Basic as SympyBasic from sympy.core.add import _addsort from sympy.core.mul import _mulsort @@ -21,6 +21,8 @@ from devito.symbolics.unevaluation import Mul as UnevalMul from devito.symbolics.unevaluation import Pow as UnevalPow from devito.symbolics.unevaluation import UnevaluableMixin +from devito.symbolics.unevaluation import ceiling as UnevalCeil +from devito.symbolics.unevaluation import floor as UnevalFloor from devito.tools import ( EnrichedTuple, as_list, as_tuple, flatten, frozendict, split, transitive_closure ) @@ -615,7 +617,10 @@ def evalrel(func=min, input=None, assumptions=None): return rfunc(func, *input) -uneval_mapper = {Add: UnevalAdd, Mul: UnevalMul, Pow: UnevalPow} +# Rounding functions also need freezing: reconstruction otherwise reevaluates +# their lowered arguments, including during Operator unpickling. +uneval_mapper = {Add: UnevalAdd, Mul: UnevalMul, Pow: UnevalPow, + ceiling: UnevalCeil, floor: UnevalFloor} def unevaluate(expr): diff --git a/devito/symbolics/unevaluation.py b/devito/symbolics/unevaluation.py index 6a860c70b5..7e348ec518 100644 --- a/devito/symbolics/unevaluation.py +++ b/devito/symbolics/unevaluation.py @@ -1,7 +1,7 @@ import sympy from sympy import sympify -__all__ = ['Add', 'Mod', 'Mul', 'Pow'] +__all__ = ['Add', 'Mod', 'Mul', 'Pow', 'ceiling', 'floor'] class UnevaluableMixin: @@ -38,3 +38,11 @@ def __new__(cls, base, exp, evaluate=None, **kwargs): return sympy.S.One else: return cls.__base__.__new__(cls, base, exp, evaluate=False, **kwargs) + + +class ceiling(sympy.ceiling, UnevaluableMixin): + __new__ = UnevaluableMixin.__new__ + + +class floor(sympy.floor, UnevaluableMixin): + __new__ = UnevaluableMixin.__new__ diff --git a/tests/test_pickle.py b/tests/test_pickle.py index 8bd989d31f..838dd2ea80 100644 --- a/tests/test_pickle.py +++ b/tests/test_pickle.py @@ -9,7 +9,7 @@ from devito import ( MPI, ConditionalDimension, Constant, Dimension, Eq, Function, Grid, IncrDimension, Min, Operator, PrecomputedSparseTimeFunction, SparseFunction, SteppingDimension, - SubDimension, SubDomain, TimeDimension, TimeFunction, floor, solve + SubDimension, SubDomain, TimeDimension, TimeFunction, ceiling, floor, solve ) from devito.data import LEFT, OWNED from devito.finite_differences.tools import centered, direct, left, right, transpose @@ -905,14 +905,15 @@ def test_collected_coeffs(self, pickle): assert str(op) == str(new_op) - def test_elemental(self, pickle): + @pytest.mark.parametrize('func', [ceiling, floor]) + def test_elemental(self, pickle, func): """ Tests that elementary functions don't get reconstructed differently. """ grid = Grid(shape=(4, 4)) f = Function(name='f', grid=grid) - op = Operator(Eq(f, floor(f + 0.5))) + op = Operator(Eq(f, func(f + 0.5))) pkl_op = pickle.dumps(op) new_op = pickle.loads(pkl_op) From 37960a8c078f0d0687f2d1fd4fd5676806e33cfa Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 22 Oct 2026 18:03:00 +0100 Subject: [PATCH 12/19] compiler: Preserve lowered operand order when printing code --- devito/ir/cgen/printer.py | 16 ++++++++++++++++ examples/cfd/01_convection.ipynb | 2 +- examples/cfd/01_convection_revisited.ipynb | 2 +- examples/performance/00_overview.ipynb | 20 ++++++++++---------- examples/performance/01_gpu.ipynb | 2 +- tests/test_symbolics.py | 21 +++++++++++++++++++++ 6 files changed, 50 insertions(+), 13 deletions(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index adfc7651aa..e7910b4745 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -19,6 +19,7 @@ from devito.symbolics.extended_sympy import BitwiseAnd from devito.symbolics.inspection import has_integer_args, sympy_dtype from devito.symbolics.queries import q_leaf +from devito.symbolics.unevaluation import UnevaluableMixin from devito.tools import ( ctypes_to_cstr, ctypes_vector_mapper, dtype_to_ctype, memoized_func ) @@ -59,6 +60,21 @@ class BasePrinter(CodePrinter): _namespaces = [] _headers = [('_POSIX_C_SOURCE', '200809L')] + def __init__(self, settings=None): + super().__init__(settings) + self._order = self.order + + def _print(self, expr, **kwargs): + # Lowered arithmetic already owns its operand order. Reset the setting + # for each child so ordinary index expressions keep their usual form. + order = self._settings['order'] + self._settings['order'] = ('none' if isinstance(expr, UnevaluableMixin) + else self._order) + try: + return super()._print(expr, **kwargs) + finally: + self._settings['order'] = order + @property def dtype(self): try: diff --git a/examples/cfd/01_convection.ipynb b/examples/cfd/01_convection.ipynb index 13ee3ac72e..91ff571c90 100644 --- a/examples/cfd/01_convection.ipynb +++ b/examples/cfd/01_convection.ipynb @@ -451,7 +451,7 @@ " {\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " u[t1][x + 1][y + 1] = dt*(-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt);\n", + " u[t1][x + 1][y + 1] = (-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt)*dt;\n", " }\n", " }\n", " STOP(section0,timers)\n", diff --git a/examples/cfd/01_convection_revisited.ipynb b/examples/cfd/01_convection_revisited.ipynb index ae762a6d16..a5c9f0b80f 100644 --- a/examples/cfd/01_convection_revisited.ipynb +++ b/examples/cfd/01_convection_revisited.ipynb @@ -427,7 +427,7 @@ " {\n", " for (int y = y_m + y_ltkn0; y <= y_M - y_rtkn0; y += 1)\n", " {\n", - " u[t1][x + 1][y + 1] = dt*(-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt);\n", + " u[t1][x + 1][y + 1] = (-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt)*dt;\n", " }\n", " }\n", " STOP(section0,timers)\n", diff --git a/examples/performance/00_overview.ipynb b/examples/performance/00_overview.ipynb index 917b75df65..35e5e3d472 100644 --- a/examples/performance/00_overview.ipynb +++ b/examples/performance/00_overview.ipynb @@ -605,7 +605,7 @@ " {\n", " float r0 = 1.0F/h_y;\n", "- u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n", - "+ u[t1][x + 4][y + 4][z + 4] = (r0*r0)*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(6.94444444e-3F*(u[t0][x + 4][y][z + 4] + u[t0][x + 4][y + 8][z + 4]) + 4.44444444e-1F*(u[t0][x + 4][y + 2][z + 4] + u[t0][x + 4][y + 6][z + 4]) + 1.11111111e-1F*(-u[t0][x + 4][y + 1][z + 4] + u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4] - u[t0][x + 4][y + 7][z + 4]) - 9.02777778e-1F*u[t0][x + 4][y + 4][z + 4])*sinf(f[x + 1][y + 1][z + 1]);\n", + "+ u[t1][x + 4][y + 4][z + 4] = (r0*r0)*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(6.94444444e-3F*(u[t0][x + 4][y + 8][z + 4] + u[t0][x + 4][y][z + 4]) + 1.11111111e-1F*(-u[t0][x + 4][y + 1][z + 4] - u[t0][x + 4][y + 7][z + 4] + u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]) + 4.44444444e-1F*(u[t0][x + 4][y + 2][z + 4] + u[t0][x + 4][y + 6][z + 4]) - 9.02777778e-1F*u[t0][x + 4][y + 4][z + 4])*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -1223,7 +1223,7 @@ " #pragma omp simd aligned(u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " r2[ys][z] = r1*(8.33333333e-2F*(u[t0][x + 4][y + 2][z + 4] - u[t0][x + 4][y + 6][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", + " r2[ys][z] = r1*(8.33333333e-2F*(-u[t0][x + 4][y + 6][z + 4] + u[t0][x + 4][y + 2][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", " }\n", " }\n", " for (int y = y0_blk0, ys = 0; y <= MIN(y_M, y0_blk0 + y0_blk0_size - 1); y += 1, ys += 1)\n", @@ -1231,7 +1231,7 @@ " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(r2[ys][z] - r2[ys + 4][z]) + 6.66666667e-1F*(-r2[ys + 1][z] + r2[ys + 3][z]))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(-r2[ys + 4][z] + r2[ys][z]) + 6.66666667e-1F*(-r2[ys + 1][z] + r2[ys + 3][z]))*r0[x][y][z];\n", " }\n", " }\n", " }\n", @@ -1350,13 +1350,13 @@ " #pragma omp simd aligned(u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " r2[yi][z] = r1*(8.33333333e-2F*(u[t0][x + 4][y + yc + 2][z + 4] - u[t0][x + 4][y + yc + 6][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + yc + 3][z + 4] + u[t0][x + 4][y + yc + 5][z + 4]));\n", + " r2[yi][z] = r1*(8.33333333e-2F*(-u[t0][x + 4][y + yc + 6][z + 4] + u[t0][x + 4][y + yc + 2][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + yc + 3][z + 4] + u[t0][x + 4][y + yc + 5][z + 4]));\n", " }\n", " }\n", " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(r2[yr0][z] - r2[yr2][z]) + 6.66666667e-1F*(r2[yr1][z] - r2[yr3][z]))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(-r2[yr2][z] + r2[yr0][z]) + 6.66666667e-1F*(-r2[yr3][z] + r2[yr1][z]))*r0[x][y][z];\n", " }\n", " }\n", " }\n", @@ -1518,7 +1518,7 @@ " #pragma omp simd aligned(u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " r2[y + 2][z] = r1*(8.33333333e-2F*(u[t0][x + 4][y + 2][z + 4] - u[t0][x + 4][y + 6][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", + " r2[y + 2][z] = r1*(8.33333333e-2F*(-u[t0][x + 4][y + 6][z + 4] + u[t0][x + 4][y + 2][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", " }\n", " }\n", " for (int y = y_m; y <= y_M; y += 1)\n", @@ -1526,7 +1526,7 @@ " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(r2[y][z] - r2[y + 4][z]) + 6.66666667e-1F*(-r2[y + 1][z] + r2[y + 3][z]))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(-r2[y + 4][z] + r2[y][z]) + 6.66666667e-1F*(-r2[y + 1][z] + r2[y + 3][z]))*r0[x][y][z];\n", " }\n", " }\n", " }\n", @@ -1665,8 +1665,8 @@ " #pragma omp simd aligned(u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " r3[x + 2][y + 2][z] = r1*(8.33333333e-2F*(u[t0][x + 2][y + 4][z + 4] - u[t0][x + 6][y + 4][z + 4]) + 6.66666667e-1F*(-u[t0][x + 3][y + 4][z + 4] + u[t0][x + 5][y + 4][z + 4]));\n", - " r4[x + 2][y + 2][z] = r2*(8.33333333e-2F*(u[t0][x + 4][y + 2][z + 4] - u[t0][x + 4][y + 6][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", + " r3[x + 2][y + 2][z] = r1*(8.33333333e-2F*(-u[t0][x + 6][y + 4][z + 4] + u[t0][x + 2][y + 4][z + 4]) + 6.66666667e-1F*(-u[t0][x + 3][y + 4][z + 4] + u[t0][x + 5][y + 4][z + 4]));\n", + " r4[x + 2][y + 2][z] = r2*(8.33333333e-2F*(-u[t0][x + 4][y + 6][z + 4] + u[t0][x + 4][y + 2][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", " }\n", " }\n", " }\n", @@ -1687,7 +1687,7 @@ " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(r1*(8.33333333e-2F*(r3[x][y + 2][z] - r3[x + 4][y + 2][z]) + 6.66666667e-1F*(-r3[x + 1][y + 2][z] + r3[x + 3][y + 2][z])) + r2*(8.33333333e-2F*(r4[x + 2][y][z] - r4[x + 2][y + 4][z]) + 6.66666667e-1F*(-r4[x + 2][y + 1][z] + r4[x + 2][y + 3][z])))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(r1*(8.33333333e-2F*(-r3[x + 4][y + 2][z] + r3[x][y + 2][z]) + 6.66666667e-1F*(-r3[x + 1][y + 2][z] + r3[x + 3][y + 2][z])) + r2*(8.33333333e-2F*(-r4[x + 2][y + 4][z] + r4[x + 2][y][z]) + 6.66666667e-1F*(-r4[x + 2][y + 1][z] + r4[x + 2][y + 3][z])))*r0[x][y][z];\n", " }\n", " }\n", " }\n", diff --git a/examples/performance/01_gpu.ipynb b/examples/performance/01_gpu.ipynb index d37e956cd6..6f64552c6a 100644 --- a/examples/performance/01_gpu.ipynb +++ b/examples/performance/01_gpu.ipynb @@ -309,7 +309,7 @@ " {\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " uL0(time + 1, x + 2, y + 2) = dt*(c*(r1*uL0(time, x + 1, y + 2) + r1*uL0(time, x + 3, y + 2) + r2*uL0(time, x + 2, y + 1) + r2*uL0(time, x + 2, y + 3) + (-2.0F)*(r1*uL0(time, x + 2, y + 2) + r2*uL0(time, x + 2, y + 2))) + r0*uL0(time, x + 2, y + 2));\n", + " uL0(time + 1, x + 2, y + 2) = (r0*uL0(time, x + 2, y + 2) + ((r1*uL0(time, x + 2, y + 2) + r2*uL0(time, x + 2, y + 2))*(-2.0F) + r1*uL0(time, x + 1, y + 2) + r1*uL0(time, x + 3, y + 2) + r2*uL0(time, x + 2, y + 1) + r2*uL0(time, x + 2, y + 3))*c)*dt;\n", " }\n", " }\n", " STOP(section0,timers)\n", diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index 23c74b5525..b3aae10efb 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -339,6 +339,27 @@ def test_extended_sympy_arithmetic(): assert ccode(-1 + bar) == '-1 + o->bar' +@pytest.mark.parametrize('cls, separator', [(UnevalAdd, ' + '), (UnevalMul, '*')]) +def test_lowered_arithmetic_printing(cls, separator): + a, b, c = [dSymbol(i, dtype=np.float32) for i in ('a', 'b', 'c')] + + # Lowered operands carry an order that must survive identifier changes + assert ccode(cls(c, a, b)) == separator.join(('c', 'a', 'b')) + renamed = uxreplace(cls(c, a, b), {c: a, a: c}) + assert ccode(renamed) == separator.join(('a', 'c', 'b')) + + # Ordinary index arithmetic retains its usual printing inside lowered math + grid = Grid(shape=(8,)) + x, = grid.dimensions + f = Function(name='f', grid=grid) + assert ccode(cls(c, f.indexed[x + 1], a)) == \ + separator.join(('c', 'f[x + 1]', 'a')) + assert ccode(x + 1) == 'x + 1' + + # Each nested lowered operation retains its own operands and precedence + assert ccode(UnevalMul(b, UnevalAdd(c, a))) == 'b*(c + a)' + + def test_integer_abs(): i1 = Dimension(name="i1") assert ccode(Abs(i1 - 1)) == "abs(i1 - 1)" From e3b183ba867e7b8b0bb06c76537b919462146a4d Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 26 Oct 2026 10:14:34 +0000 Subject: [PATCH 13/19] compiler: Address kernel naming review feedback Retrieve Functions uniquely when reserving names, share the default uxreplace traversal through an explicit helper, and remove the unneeded Operator name from the input-name collision test. Keep the CSE position map to avoid repeated linear scans of exprs. In a comparison of the actual topological sort with 128 independent CTemps, including map construction, cached lookup took 39 ms versus 1255 ms for list.index. Add a comment explaining the cached lookup. Validation: focused C/OpenMP regressions (64 passed, 1 xfailed) and C++/OpenMP regressions (62 passed, 1 xfailed). The full CPU audit found separate numerical failures already present before the rebase. --- devito/ir/support/symregistry.py | 2 +- devito/passes/clusters/cse.py | 1 + devito/symbolics/manipulation.py | 8 ++++++-- tests/test_operator.py | 2 +- 4 files changed, 9 insertions(+), 4 deletions(-) diff --git a/devito/ir/support/symregistry.py b/devito/ir/support/symregistry.py index 01ceaf3bdf..e32e37acd6 100644 --- a/devito/ir/support/symregistry.py +++ b/devito/ir/support/symregistry.py @@ -39,7 +39,7 @@ def reserve(self, exprs): self.reserved.update(i.name for e in exprs for i in e.free_symbols) # Before indexification, a Function's free symbols omit its own name - self.reserved.update(i.name for i in retrieve_functions(exprs)) + self.reserved.update(i.name for i in retrieve_functions(exprs, mode='unique')) def make_name(self, prefix=None, increment_first=True): # By default we're creating a new symbol diff --git a/devito/passes/clusters/cse.py b/devito/passes/clusters/cse.py index 75fe68d636..f00bf05e04 100644 --- a/devito/passes/clusters/cse.py +++ b/devito/passes/clusters/cse.py @@ -293,6 +293,7 @@ def _toposort(exprs): if e0.lhs in e1.rhs.free_symbols: dag.add_edge(e0, e1, force_add=True) + # Cache positions to avoid repeated linear scans while choosing ready nodes ordering = {e: n for n, e in enumerate(exprs)} def choose_element(queue, scheduled): diff --git a/devito/symbolics/manipulation.py b/devito/symbolics/manipulation.py index f1a7524ffb..0f068077e1 100644 --- a/devito/symbolics/manipulation.py +++ b/devito/symbolics/manipulation.py @@ -81,6 +81,10 @@ def uxreplace(expr, rule): @singledispatch def _uxreplace(expr, rule): + return _uxreplace_default(expr, rule) + + +def _uxreplace_default(expr, rule): if expr in rule: v = rule[expr] if not isinstance(v, dict): @@ -124,7 +128,7 @@ def _uxreplace(expr, rule): @_uxreplace.register(CallFromPointer) def _(expr, rule): if expr in rule or expr.__class__ in rule: - return _uxreplace.dispatch(object)(expr, rule) + return _uxreplace_default(expr, rule) args, changed = _uxreplace_dispatch((expr.pointer, expr.params), rule) if changed: @@ -143,7 +147,7 @@ def _(expr, rule): mapper[d] = d._rebuild(**kwargs) # Apply the same dimension substitutions to the summand and its bindings - return _uxreplace.dispatch(object)(expr, {**mapper, **rule}) + return _uxreplace_default(expr, {**mapper, **rule}) @singledispatch diff --git a/tests/test_operator.py b/tests/test_operator.py index 8b34a93bc9..a2cf17861a 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -206,7 +206,7 @@ def test_generated_names_avoid_inputs(self, name): x, y = np.indices(f.data_with_halo.shape, dtype=np.float32) f.data_with_halo[:] = x + 2*y - op = Operator(Eq(u, f.dx + f.dy), name='names_inputs') + op = Operator(Eq(u, f.dx + f.dy)) op.apply() np.testing.assert_allclose(u.data, 3., rtol=1.e-5, atol=1.e-6) From 0325dd2d74b42f734e6b217a1759494d7c1ed149 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 26 Oct 2026 10:14:34 +0000 Subject: [PATCH 14/19] tests: Match lowered operand order in conditional dimension output --- tests/test_dimension.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_dimension.py b/tests/test_dimension.py index 3e28210fe5..72badad3cf 100644 --- a/tests/test_dimension.py +++ b/tests/test_dimension.py @@ -1535,11 +1535,11 @@ def test_stepping_dim_in_condition_lowering(self): # ``` # if (g[t0][x + 1][y + 1] <= 10) # { - # g[t1][x + 1][y + 1] = g[t0][x + 1][y + 1] + 1 + # g[t1][x + 1][y + 1] = 1 + g[t0][x + 1][y + 1] # ``` # is in the generated code, but indentation etc. seems to vary. part1 = 'if (g[t0][x + 1][y + 1] <= 10)\n' - part2 = 'g[t1][x + 1][y + 1] = g[t0][x + 1][y + 1] + 1' + part2 = 'g[t1][x + 1][y + 1] = 1 + g[t0][x + 1][y + 1]' whole_code = str(op.ccode) try: @@ -1549,7 +1549,7 @@ def test_stepping_dim_in_condition_lowering(self): except AssertionError: # Try the alternative string part1 = 'if (gL0(t0, x + 1, y + 1) <= 10)\n' - part2 = 'gL0(t1, x + 1, y + 1) = gL0(t0, x + 1, y + 1) + 1' + part2 = 'gL0(t1, x + 1, y + 1) = 1 + gL0(t0, x + 1, y + 1)' loc = whole_code.find(part1) assert loc != -1 assert whole_code.find(part2, loc + len(part1)) != -1 From 5f09d1f70d2c0ce30877d69e2ed52a8a2b963ca2 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 26 Oct 2026 11:18:31 +0000 Subject: [PATCH 15/19] compiler: Print integer offsets after symbolic operands Stably place numeric offsets last in commutative integer sums. This emits bx_size - 1 without sorting symbolic operands by name or changing floating-point summation order. Preserve opaque noncommutative expressions. Cover positive and negative offsets with 32- and 64-bit integers and mixed floating-point addends. Refresh two stale floating-point code-display lines in the MPI overview, keeping notebook source and comparison rules unchanged. Validation: both MPI notebooks pass all 35 cell checks on a private four-rank cluster; affected C/OpenMP compiler tests pass 554 cases, and C++/OpenMP symbolic/IET tests pass 192 cases. --- devito/ir/cgen/printer.py | 3 +++ examples/mpi/overview.ipynb | 4 ++-- tests/test_symbolics.py | 14 ++++++++++++++ 3 files changed, 19 insertions(+), 2 deletions(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index e7910b4745..a339f980c9 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -365,6 +365,9 @@ def _print_Add(self, expr, order=None): Print an addition. """ terms = self._as_ordered_terms(expr, order=order) + if expr.is_commutative and terms[0].is_Number and has_integer_args(*expr.args): + # Keep integer offsets last without sorting symbolic operands + terms = sorted(terms, key=lambda i: i.is_Number) PREC = precedence(expr) l = [] diff --git a/examples/mpi/overview.ipynb b/examples/mpi/overview.ipynb index 8c209d3bc6..194e3746aa 100644 --- a/examples/mpi/overview.ipynb +++ b/examples/mpi/overview.ipynb @@ -362,7 +362,7 @@ " #pragma omp simd aligned(u:64)\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " u[t1][x + 2][y + 2] = u[t0][x + 2][y + 2] + 1;\n", + " u[t1][x + 2][y + 2] = 1 + u[t0][x + 2][y + 2];\n", " }\n", " }\n", " STOP(section0,timers)\n", @@ -471,7 +471,7 @@ " #pragma omp simd aligned(u:64)\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " u[t1][x + 2][y + 2] = -r0*u[t0][x + 2][y + 2] + r0*u[t0][x + 3][y + 2] + 1;\n", + " u[t1][x + 2][y + 2] = 1 - r0*u[t0][x + 2][y + 2] + r0*u[t0][x + 3][y + 2];\n", " }\n", " }\n", " STOP(section0,timers)\n", diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index b3aae10efb..73d7498887 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -360,6 +360,20 @@ def test_lowered_arithmetic_printing(cls, separator): assert ccode(UnevalMul(b, UnevalAdd(c, a))) == 'b*(c + a)' +@pytest.mark.parametrize('dtype', [np.int32, np.int64]) +@pytest.mark.parametrize('offset, expected', [(-1, ' - 1'), (1, ' + 1')]) +def test_lowered_integer_offset(dtype, offset, expected): + size = dSymbol('bx_size', dtype=dtype) + other = dSymbol('a_size', dtype=dtype) + + assert ccode(UnevalAdd(offset, size)) == 'bx_size' + expected + assert ccode(UnevalAdd(offset, size, other)) == 'bx_size + a_size' + expected + + # A floating-point addend requires preserving the original summation order + value = dSymbol('value', dtype=np.float32) + assert ccode(UnevalAdd(offset, size, value)) == f'{offset} + bx_size + value' + + def test_integer_abs(): i1 = Dimension(name="i1") assert ccode(Abs(i1 - 1)) == "abs(i1 - 1)" From 5e98cfb7fafcb9c9eed228a8cc0d9d5c91e237aa Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 26 Oct 2026 11:26:03 +0000 Subject: [PATCH 16/19] compiler: Require integer addends when formatting offsets --- devito/ir/cgen/printer.py | 4 +++- tests/test_symbolics.py | 2 ++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index a339f980c9..b6daffd4ce 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -365,7 +365,9 @@ def _print_Add(self, expr, order=None): Print an addition. """ terms = self._as_ordered_terms(expr, order=order) - if expr.is_commutative and terms[0].is_Number and has_integer_args(*expr.args): + # Integer inputs alone do not prove an integer result (e.g. sin(n)) + if (expr.is_commutative and terms[0].is_Number and + all(has_integer_args(i) for i in expr.args)): # Keep integer offsets last without sorting symbolic operands terms = sorted(terms, key=lambda i: i.is_Number) diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index 73d7498887..815d502ed1 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -372,6 +372,8 @@ def test_lowered_integer_offset(dtype, offset, expected): # A floating-point addend requires preserving the original summation order value = dSymbol('value', dtype=np.float32) assert ccode(UnevalAdd(offset, size, value)) == f'{offset} + bx_size + value' + assert ccode(UnevalAdd(offset, size, sympy.sin(size))) == \ + f'{offset} + bx_size + sin(bx_size)' def test_integer_abs(): From cac600d7316791be62960330bce85d5d609afeee Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 26 Oct 2026 14:39:34 +0000 Subject: [PATCH 17/19] compiler: Freeze arithmetic order before local renaming Choose emitted Add/Mul order at the existing freezing boundary, using capture-avoiding definition-order comparison symbols for local Temps. Preserve original operands and groups through replacement, serialization and printing. Remove the per-node printer setting mutation and integer result-type formatting workaround. Keep frozen sums grouped and preserve factor order during negative coefficient extraction. Retain ordinary symbolic display, refresh generated notebook output, and avoid fixed operand positions in the visitor regression. Preserve immutable mapping types during symbolic reconstruction, and explain why the concrete frozen Add/Mul dispatch registrations are required by MRO. --- devito/ir/cgen/printer.py | 29 ++------ devito/passes/clusters/unevaluate.py | 22 ++++-- devito/symbolics/manipulation.py | 23 ++++++- devito/symbolics/unevaluation.py | 7 ++ examples/cfd/01_convection.ipynb | 2 +- examples/cfd/01_convection_revisited.ipynb | 2 +- examples/mpi/overview.ipynb | 4 +- examples/performance/00_overview.ipynb | 44 ++++++------ examples/performance/01_gpu.ipynb | 2 +- tests/test_dimension.py | 6 +- tests/test_symbolics.py | 79 ++++++++++++++++++---- tests/test_visitors.py | 5 +- 12 files changed, 151 insertions(+), 74 deletions(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index b6daffd4ce..c770e2a25b 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -19,7 +19,7 @@ from devito.symbolics.extended_sympy import BitwiseAnd from devito.symbolics.inspection import has_integer_args, sympy_dtype from devito.symbolics.queries import q_leaf -from devito.symbolics.unevaluation import UnevaluableMixin +from devito.symbolics.unevaluation import Add as UnevalAdd from devito.tools import ( ctypes_to_cstr, ctypes_vector_mapper, dtype_to_ctype, memoized_func ) @@ -60,21 +60,6 @@ class BasePrinter(CodePrinter): _namespaces = [] _headers = [('_POSIX_C_SOURCE', '200809L')] - def __init__(self, settings=None): - super().__init__(settings) - self._order = self.order - - def _print(self, expr, **kwargs): - # Lowered arithmetic already owns its operand order. Reset the setting - # for each child so ordinary index expressions keep their usual form. - order = self._settings['order'] - self._settings['order'] = ('none' if isinstance(expr, UnevaluableMixin) - else self._order) - try: - return super()._print(expr, **kwargs) - finally: - self._settings['order'] = order - @property def dtype(self): try: @@ -299,6 +284,11 @@ def _print_Mul(self, expr): args = [a for a in expr.args if a != -1] neg = (len(expr.args) - len(args)) % 2 + # SymPy's sign extraction rebuilds a plain Mul, losing frozen factor order + if args and args[0].is_Number and args[0].is_negative: + args[0] = -args[0] + neg ^= 1 + if len(args) > 1: term = super()._print_Mul(expr.func(*args, evaluate=False)) else: @@ -365,17 +355,12 @@ def _print_Add(self, expr, order=None): Print an addition. """ terms = self._as_ordered_terms(expr, order=order) - # Integer inputs alone do not prove an integer result (e.g. sin(n)) - if (expr.is_commutative and terms[0].is_Number and - all(has_integer_args(i) for i in expr.args)): - # Keep integer offsets last without sorting symbolic operands - terms = sorted(terms, key=lambda i: i.is_Number) PREC = precedence(expr) l = [] for term in terms: t = self._print(term) - if precedence(term) < PREC: + if precedence(term) < PREC or (isinstance(expr, UnevalAdd) and term.is_Add): l.extend(["+", f"({t})"]) elif t.startswith('-'): l.extend(["-", t[1:]]) diff --git a/devito/passes/clusters/unevaluate.py b/devito/passes/clusters/unevaluate.py index b895e1eafa..26e65553e9 100644 --- a/devito/passes/clusters/unevaluate.py +++ b/devito/passes/clusters/unevaluate.py @@ -1,11 +1,23 @@ -from devito.ir import cluster_pass +from devito.ir import SymbolRegistry from devito.symbolics import unevaluate as _unevaluate +from devito.types import Temp __all__ = ['unevaluate'] -@cluster_pass(mode='all') -def unevaluate(cluster): - exprs = [_unevaluate(e) for e in cluster.exprs] +def unevaluate(clusters): + # Abstract local temporaries in definition order before choosing operand order + exprs = [e for c in clusters for e in c.exprs] + sregistry = SymbolRegistry() + sregistry.reserve(exprs) + # Lexical comparison must respect definition order across digit boundaries + width = len(str(len(exprs))) + mapper = {e.lhs: e.lhs._rebuild(name=sregistry.make_name( + prefix=f'_order{n:0{width}}_')) + for n, e in enumerate(exprs) if isinstance(e.lhs, Temp)} - return cluster.rebuild(exprs=exprs) + processed = [] + for c in clusters: + exprs = [_unevaluate(e, mapper) for e in c.exprs] + processed.append(c.rebuild(exprs=exprs)) + return processed diff --git a/devito/symbolics/manipulation.py b/devito/symbolics/manipulation.py index 0f068077e1..0df752e197 100644 --- a/devito/symbolics/manipulation.py +++ b/devito/symbolics/manipulation.py @@ -192,7 +192,7 @@ def _(mapper, rule): vx, flag = _uxreplace_dispatch(v, rule) ret[k] = vx changed |= flag - return (ret, True) if changed else (mapper, False) + return (mapper.__class__(ret), True) if changed else (mapper, False) @singledispatch @@ -217,6 +217,7 @@ def _(expr, args, kwargs): return expr.func(*args, evaluate=evaluate) +# The SymPy bases precede UnevaluableMixin in these classes' dispatch MRO @_uxreplace_handle.register(UnevalAdd) @_uxreplace_handle.register(UnevalMul) def _(expr, args, kwargs): @@ -627,11 +628,27 @@ def evalrel(func=min, input=None, assumptions=None): ceiling: UnevalCeil, floor: UnevalFloor} -def unevaluate(expr): +def unevaluate(expr, mapper=None): + """ + Freeze code-generation operand order without reevaluating arithmetic. + + `mapper` supplies comparison symbols; original operands and grouping survive. + """ if q_leaf(expr): return expr - args = [unevaluate(a) for a in expr.args] + args = [unevaluate(a, mapper) for a in expr.args] + if expr.func in (Add, Mul): + # Order comparison copies; preserve the original operands and grouping + keys = [uxreplace(a, mapper) for a in args] + originals = dict(zip(keys, args, strict=True)) + handle = expr.func(*keys, evaluate=False) + if expr.func is Add: + keys = handle.as_ordered_terms() + else: + keys = handle.as_ordered_factors() + # Factor ordering can split a negative coefficient into -1 and its magnitude + args = [originals.get(a, a) for a in keys] try: return uneval_mapper[expr.func](*args) diff --git a/devito/symbolics/unevaluation.py b/devito/symbolics/unevaluation.py index 7e348ec518..cb49840fcc 100644 --- a/devito/symbolics/unevaluation.py +++ b/devito/symbolics/unevaluation.py @@ -13,6 +13,10 @@ def __new__(cls, *args, evaluate=None, **kwargs): class Add(sympy.Add, UnevaluableMixin): __new__ = UnevaluableMixin.__new__ + def as_ordered_terms(self, order=None, data=False): + # Lowered arithmetic already carries its code-generation order + return self.as_terms() if data else list(self.args) + class Mod(sympy.Mod, UnevaluableMixin): __new__ = UnevaluableMixin.__new__ @@ -26,6 +30,9 @@ def __new__(cls, *args, evaluate=None, **kwargs): return super().__new__(cls, *args, evaluate=False, **kwargs) + def as_ordered_factors(self, order=None): + return list(self.args) + class Pow(sympy.Pow, UnevaluableMixin): diff --git a/examples/cfd/01_convection.ipynb b/examples/cfd/01_convection.ipynb index 91ff571c90..13ee3ac72e 100644 --- a/examples/cfd/01_convection.ipynb +++ b/examples/cfd/01_convection.ipynb @@ -451,7 +451,7 @@ " {\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " u[t1][x + 1][y + 1] = (-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt)*dt;\n", + " u[t1][x + 1][y + 1] = dt*(-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt);\n", " }\n", " }\n", " STOP(section0,timers)\n", diff --git a/examples/cfd/01_convection_revisited.ipynb b/examples/cfd/01_convection_revisited.ipynb index a5c9f0b80f..ae762a6d16 100644 --- a/examples/cfd/01_convection_revisited.ipynb +++ b/examples/cfd/01_convection_revisited.ipynb @@ -427,7 +427,7 @@ " {\n", " for (int y = y_m + y_ltkn0; y <= y_M - y_rtkn0; y += 1)\n", " {\n", - " u[t1][x + 1][y + 1] = (-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt)*dt;\n", + " u[t1][x + 1][y + 1] = dt*(-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt);\n", " }\n", " }\n", " STOP(section0,timers)\n", diff --git a/examples/mpi/overview.ipynb b/examples/mpi/overview.ipynb index 194e3746aa..5376881b0d 100644 --- a/examples/mpi/overview.ipynb +++ b/examples/mpi/overview.ipynb @@ -362,7 +362,7 @@ " #pragma omp simd aligned(u:64)\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " u[t1][x + 2][y + 2] = 1 + u[t0][x + 2][y + 2];\n", + " u[t1][x + 2][y + 2] = u[t0][x + 2][y + 2] + 1;\n", " }\n", " }\n", " STOP(section0,timers)\n", @@ -471,7 +471,7 @@ " #pragma omp simd aligned(u:64)\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " u[t1][x + 2][y + 2] = 1 - r0*u[t0][x + 2][y + 2] + r0*u[t0][x + 3][y + 2];\n", + " u[t1][x + 2][y + 2] = (-r0*u[t0][x + 2][y + 2] + r0*u[t0][x + 3][y + 2]) + 1;\n", " }\n", " }\n", " STOP(section0,timers)\n", diff --git a/examples/performance/00_overview.ipynb b/examples/performance/00_overview.ipynb index 35e5e3d472..230fe642bf 100644 --- a/examples/performance/00_overview.ipynb +++ b/examples/performance/00_overview.ipynb @@ -208,7 +208,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -283,7 +283,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -341,7 +341,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -431,7 +431,7 @@ " {\n", " for (int z = z0_blk1; z <= MIN(z_M, z0_blk1 + z0_blk1_size - 1); z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -522,7 +522,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 5][z + 4]) + (-8.33333333e-2F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 8][z + 4]) + (8.33333333e-2F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 7][z + 4]))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 7][z + 4]) + (-8.33333333e-2F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 8][z + 4]) + (-6.66666667e-1F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 5][z + 4]))*r0[x][y][z];\n", " }\n", " }\n", " }\n", @@ -567,9 +567,9 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - "- u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + "- u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", "+ float r0 = 1.0F/h_y;\n", - "+ u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n", + "+ u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -604,8 +604,8 @@ " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", " float r0 = 1.0F/h_y;\n", - "- u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n", - "+ u[t1][x + 4][y + 4][z + 4] = (r0*r0)*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(6.94444444e-3F*(u[t0][x + 4][y + 8][z + 4] + u[t0][x + 4][y][z + 4]) + 1.11111111e-1F*(-u[t0][x + 4][y + 1][z + 4] - u[t0][x + 4][y + 7][z + 4] + u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]) + 4.44444444e-1F*(u[t0][x + 4][y + 2][z + 4] + u[t0][x + 4][y + 6][z + 4]) - 9.02777778e-1F*u[t0][x + 4][y + 4][z + 4])*sinf(f[x + 1][y + 1][z + 1]);\n", + "- u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n", + "+ u[t1][x + 4][y + 4][z + 4] = (r0*r0)*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(6.94444444e-3F*(u[t0][x + 4][y][z + 4] + u[t0][x + 4][y + 8][z + 4]) + 4.44444444e-1F*(u[t0][x + 4][y + 2][z + 4] + u[t0][x + 4][y + 6][z + 4]) + 1.11111111e-1F*(-u[t0][x + 4][y + 1][z + 4] + u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4] - u[t0][x + 4][y + 7][z + 4]) - 9.02777778e-1F*u[t0][x + 4][y + 4][z + 4])*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -747,7 +747,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*r0[x][y + 2][z]/h_x - 6.66666667e-1F*r0[x + 1][y + 2][z]/h_x + 6.66666667e-1F*r0[x + 3][y + 2][z]/h_x - 8.33333333e-2F*r0[x + 4][y + 2][z]/h_x + 8.33333333e-2F*r1[x + 2][y][z]/h_y - 6.66666667e-1F*r1[x + 2][y + 1][z]/h_y + 6.66666667e-1F*r1[x + 2][y + 3][z]/h_y - 8.33333333e-2F*r1[x + 2][y + 4][z]/h_y)*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r0[x][y + 2][z]/h_x - 6.66666667e-1F*r0[x + 1][y + 2][z]/h_x + 6.66666667e-1F*r0[x + 3][y + 2][z]/h_x - 8.33333333e-2F*r0[x + 4][y + 2][z]/h_x) + (8.33333333e-2F*r1[x + 2][y][z]/h_y - 6.66666667e-1F*r1[x + 2][y + 1][z]/h_y + 6.66666667e-1F*r1[x + 2][y + 3][z]/h_y - 8.33333333e-2F*r1[x + 2][y + 4][z]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -818,7 +818,7 @@ " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*r0[x][y][z]/h_x - 6.66666667e-1F*r0[x + 1][y][z]/h_x + 6.66666667e-1F*r0[x + 3][y][z]/h_x - 8.33333333e-2F*r0[x + 4][y][z]/h_x + 8.33333333e-2F*r1[y][z]/h_y - 6.66666667e-1F*r1[y + 1][z]/h_y + 6.66666667e-1F*r1[y + 3][z]/h_y - 8.33333333e-2F*r1[y + 4][z]/h_y)*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r0[x][y][z]/h_x - 6.66666667e-1F*r0[x + 1][y][z]/h_x + 6.66666667e-1F*r0[x + 3][y][z]/h_x - 8.33333333e-2F*r0[x + 4][y][z]/h_x) + (8.33333333e-2F*r1[y][z]/h_y - 6.66666667e-1F*r1[y + 1][z]/h_y + 6.66666667e-1F*r1[y + 3][z]/h_y - 8.33333333e-2F*r1[y + 4][z]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -882,7 +882,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -1020,7 +1020,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -1223,7 +1223,7 @@ " #pragma omp simd aligned(u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " r2[ys][z] = r1*(8.33333333e-2F*(-u[t0][x + 4][y + 6][z + 4] + u[t0][x + 4][y + 2][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", + " r2[ys][z] = r1*(8.33333333e-2F*(u[t0][x + 4][y + 2][z + 4] - u[t0][x + 4][y + 6][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", " }\n", " }\n", " for (int y = y0_blk0, ys = 0; y <= MIN(y_M, y0_blk0 + y0_blk0_size - 1); y += 1, ys += 1)\n", @@ -1231,7 +1231,7 @@ " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(-r2[ys + 4][z] + r2[ys][z]) + 6.66666667e-1F*(-r2[ys + 1][z] + r2[ys + 3][z]))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(r2[ys][z] - r2[ys + 4][z]) + 6.66666667e-1F*(-r2[ys + 1][z] + r2[ys + 3][z]))*r0[x][y][z];\n", " }\n", " }\n", " }\n", @@ -1350,13 +1350,13 @@ " #pragma omp simd aligned(u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " r2[yi][z] = r1*(8.33333333e-2F*(-u[t0][x + 4][y + yc + 6][z + 4] + u[t0][x + 4][y + yc + 2][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + yc + 3][z + 4] + u[t0][x + 4][y + yc + 5][z + 4]));\n", + " r2[yi][z] = r1*(8.33333333e-2F*(u[t0][x + 4][y + yc + 2][z + 4] - u[t0][x + 4][y + yc + 6][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + yc + 3][z + 4] + u[t0][x + 4][y + yc + 5][z + 4]));\n", " }\n", " }\n", " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(-r2[yr2][z] + r2[yr0][z]) + 6.66666667e-1F*(-r2[yr3][z] + r2[yr1][z]))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(r2[yr0][z] - r2[yr2][z]) + 6.66666667e-1F*(r2[yr1][z] - r2[yr3][z]))*r0[x][y][z];\n", " }\n", " }\n", " }\n", @@ -1518,7 +1518,7 @@ " #pragma omp simd aligned(u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " r2[y + 2][z] = r1*(8.33333333e-2F*(-u[t0][x + 4][y + 6][z + 4] + u[t0][x + 4][y + 2][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", + " r2[y + 2][z] = r1*(8.33333333e-2F*(u[t0][x + 4][y + 2][z + 4] - u[t0][x + 4][y + 6][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", " }\n", " }\n", " for (int y = y_m; y <= y_M; y += 1)\n", @@ -1526,7 +1526,7 @@ " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(-r2[y + 4][z] + r2[y][z]) + 6.66666667e-1F*(-r2[y + 1][z] + r2[y + 3][z]))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = r1*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*(r2[y][z] - r2[y + 4][z]) + 6.66666667e-1F*(-r2[y + 1][z] + r2[y + 3][z]))*r0[x][y][z];\n", " }\n", " }\n", " }\n", @@ -1665,8 +1665,8 @@ " #pragma omp simd aligned(u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " r3[x + 2][y + 2][z] = r1*(8.33333333e-2F*(-u[t0][x + 6][y + 4][z + 4] + u[t0][x + 2][y + 4][z + 4]) + 6.66666667e-1F*(-u[t0][x + 3][y + 4][z + 4] + u[t0][x + 5][y + 4][z + 4]));\n", - " r4[x + 2][y + 2][z] = r2*(8.33333333e-2F*(-u[t0][x + 4][y + 6][z + 4] + u[t0][x + 4][y + 2][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", + " r3[x + 2][y + 2][z] = r1*(8.33333333e-2F*(u[t0][x + 2][y + 4][z + 4] - u[t0][x + 6][y + 4][z + 4]) + 6.66666667e-1F*(-u[t0][x + 3][y + 4][z + 4] + u[t0][x + 5][y + 4][z + 4]));\n", + " r4[x + 2][y + 2][z] = r2*(8.33333333e-2F*(u[t0][x + 4][y + 2][z + 4] - u[t0][x + 4][y + 6][z + 4]) + 6.66666667e-1F*(-u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4]));\n", " }\n", " }\n", " }\n", @@ -1687,7 +1687,7 @@ " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(r1*(8.33333333e-2F*(-r3[x + 4][y + 2][z] + r3[x][y + 2][z]) + 6.66666667e-1F*(-r3[x + 1][y + 2][z] + r3[x + 3][y + 2][z])) + r2*(8.33333333e-2F*(-r4[x + 2][y + 4][z] + r4[x + 2][y][z]) + 6.66666667e-1F*(-r4[x + 2][y + 1][z] + r4[x + 2][y + 3][z])))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(r1*(8.33333333e-2F*(r3[x][y + 2][z] - r3[x + 4][y + 2][z]) + 6.66666667e-1F*(-r3[x + 1][y + 2][z] + r3[x + 3][y + 2][z])) + r2*(8.33333333e-2F*(r4[x + 2][y][z] - r4[x + 2][y + 4][z]) + 6.66666667e-1F*(-r4[x + 2][y + 1][z] + r4[x + 2][y + 3][z])))*r0[x][y][z];\n", " }\n", " }\n", " }\n", diff --git a/examples/performance/01_gpu.ipynb b/examples/performance/01_gpu.ipynb index 6f64552c6a..dc096aa0eb 100644 --- a/examples/performance/01_gpu.ipynb +++ b/examples/performance/01_gpu.ipynb @@ -309,7 +309,7 @@ " {\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " uL0(time + 1, x + 2, y + 2) = (r0*uL0(time, x + 2, y + 2) + ((r1*uL0(time, x + 2, y + 2) + r2*uL0(time, x + 2, y + 2))*(-2.0F) + r1*uL0(time, x + 1, y + 2) + r1*uL0(time, x + 3, y + 2) + r2*uL0(time, x + 2, y + 1) + r2*uL0(time, x + 2, y + 3))*c)*dt;\n", + " uL0(time + 1, x + 2, y + 2) = dt*(c*(r1*uL0(time, x + 1, y + 2) + r1*uL0(time, x + 3, y + 2) + r2*uL0(time, x + 2, y + 1) + r2*uL0(time, x + 2, y + 3) - 2.0F*(r1*uL0(time, x + 2, y + 2) + r2*uL0(time, x + 2, y + 2))) + r0*uL0(time, x + 2, y + 2));\n", " }\n", " }\n", " STOP(section0,timers)\n", diff --git a/tests/test_dimension.py b/tests/test_dimension.py index 72badad3cf..3e28210fe5 100644 --- a/tests/test_dimension.py +++ b/tests/test_dimension.py @@ -1535,11 +1535,11 @@ def test_stepping_dim_in_condition_lowering(self): # ``` # if (g[t0][x + 1][y + 1] <= 10) # { - # g[t1][x + 1][y + 1] = 1 + g[t0][x + 1][y + 1] + # g[t1][x + 1][y + 1] = g[t0][x + 1][y + 1] + 1 # ``` # is in the generated code, but indentation etc. seems to vary. part1 = 'if (g[t0][x + 1][y + 1] <= 10)\n' - part2 = 'g[t1][x + 1][y + 1] = 1 + g[t0][x + 1][y + 1]' + part2 = 'g[t1][x + 1][y + 1] = g[t0][x + 1][y + 1] + 1' whole_code = str(op.ccode) try: @@ -1549,7 +1549,7 @@ def test_stepping_dim_in_condition_lowering(self): except AssertionError: # Try the alternative string part1 = 'if (gL0(t0, x + 1, y + 1) <= 10)\n' - part2 = 'gL0(t1, x + 1, y + 1) = 1 + gL0(t0, x + 1, y + 1)' + part2 = 'gL0(t1, x + 1, y + 1) = gL0(t0, x + 1, y + 1) + 1' loc = whole_code.find(part1) assert loc != -1 assert whole_code.find(part2, loc + len(part1)) != -1 diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index 815d502ed1..11cb229971 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -19,11 +19,11 @@ FieldFromComposite, FieldFromPointer, IntDiv, ListInitializer, Namespace, ReservedWord, RoundUp, Rvalue, SizeOf, VectorAccess, evalrel, pow_to_mul, q_negative, q_positive, retrieve_derivatives, retrieve_functions, retrieve_indexed, - subs_if_composite, uxreplace, xreplace_indices + subs_if_composite, unevaluate, uxreplace, xreplace_indices ) from devito.symbolics.unevaluation import Add as UnevalAdd from devito.symbolics.unevaluation import Mul as UnevalMul -from devito.tools import CustomDtype, as_tuple, dtypes_vector_mapper +from devito.tools import CustomDtype, as_tuple, dtypes_vector_mapper, frozendict from devito.types import ( Array, Bundle, ComponentAccess, FIndexed, LocalObject, Object, StencilDimension ) @@ -358,6 +358,33 @@ def test_lowered_arithmetic_printing(cls, separator): # Each nested lowered operation retains its own operands and precedence assert ccode(UnevalMul(b, UnevalAdd(c, a))) == 'b*(c + a)' + assert ccode(cls(b, cls(c, a))) == f'b{separator}(c{separator}a)' + assert ccode(UnevalAdd(b, sympy.Add(c, a, evaluate=False))) == 'b + (a + c)' + + +def test_lowered_negative_product(): + a, b, c = [dSymbol(i, dtype=np.float32) for i in ('a', 'b', 'c')] + + # Extracting a negative coefficient must retain the frozen factors + assert ccode(UnevalMul(-2, c, a, b)) == '-2*c*a*b' + expr = -2*c*a*b + assert ccode(unevaluate(expr)) == ccode(expr) + + +@pytest.mark.parametrize('cls, lowered, separator', [ + (sympy.Add, UnevalAdd, ' + '), + (sympy.Mul, UnevalMul, '*') +]) +def test_freeze_grouping(cls, lowered, separator): + a, b, c = [Temp(name=n, dtype=np.float32) for n in ('r9', 'r10', 'r11')] + x, y, z = [Temp(name=n, dtype=np.float32) for n in ('r0', 'r1', 'r2')] + expr = cls(c, cls(b, a, evaluate=False), evaluate=False) + + frozen = unevaluate(expr, {a: x, b: y, c: z}) + + assert frozen.args == (c, lowered(a, b)) + assert ccode(frozen) == f'r11{separator}(r9{separator}r10)' + assert unevaluate(frozen) is frozen @pytest.mark.parametrize('dtype', [np.int32, np.int64]) @@ -366,14 +393,35 @@ def test_lowered_integer_offset(dtype, offset, expected): size = dSymbol('bx_size', dtype=dtype) other = dSymbol('a_size', dtype=dtype) - assert ccode(UnevalAdd(offset, size)) == 'bx_size' + expected - assert ccode(UnevalAdd(offset, size, other)) == 'bx_size + a_size' + expected + lowered = unevaluate(size + offset) + assert lowered.args == (size, offset) + assert ccode(lowered) == 'bx_size' + expected + + # Renaming retains the order chosen before freezing + lowered = unevaluate(size + other + offset) + renamed = uxreplace(lowered, {size: other, other: size}) + assert ccode(renamed) == 'bx_size + a_size' + expected - # A floating-point addend requires preserving the original summation order + # Result types do not affect how the arithmetic order is frozen value = dSymbol('value', dtype=np.float32) - assert ccode(UnevalAdd(offset, size, value)) == f'{offset} + bx_size + value' - assert ccode(UnevalAdd(offset, size, sympy.sin(size))) == \ - f'{offset} + bx_size + sin(bx_size)' + for addend in (value, sympy.sin(size)): + expr = size + addend + offset + assert ccode(unevaluate(expr)) == ccode(expr) + + +@pytest.mark.parametrize('cls, separator', [(sympy.Add, ' + '), (sympy.Mul, '*')]) +@pytest.mark.parametrize('names', [('r9', 'r10', 'r11'), ('r21', 'r20', 'r19'), + ('later', 'middle', 'earlier')]) +def test_lowered_temporary_order(cls, separator, names): + grid = Grid(shape=(4,)) + f, g, h, out = [Function(name=n, grid=grid) for n in ('f', 'g', 'h', 'out')] + a, b, c = [Temp(name=n, dtype=np.float32) for n in names] + op = Operator([Eq(a, cos(f)), Eq(b, sin(g)), Eq(c, cos(h)), + Eq(out, cls(a, b, c))], opt='noop') + expr, = [i for i in FindNodes(Expression).visit(op) if i.write == out] + + assert expr.expr.rhs.args == (a, b, c) + assert str(expr).endswith(f' = {separator.join(names)};') def test_integer_abs(): @@ -1077,17 +1125,20 @@ def test_xreplace_indices_reuses_untouched_sequence(self): assert xreplace_indices(exprs, {z: z + 1}) is exprs assert xreplace_indices(exprs, {x: x + 1}) is not exprs - def test_custom_reconstructable(self): + @pytest.mark.parametrize('mapping_type', [dict, frozendict]) + def test_custom_reconstructable(self, mapping_type): class MyDefFunction(DefFunction): __rargs__ = ('name', 'arguments') - __rkwargs__ = ('p0', 'p1', 'p2') + __rkwargs__ = ('p0', 'p1', 'p2', 'p3') - def __new__(cls, name=None, arguments=None, p0=None, p1=None, p2=None): + def __new__(cls, name=None, arguments=None, p0=None, p1=None, p2=None, + p3=None): obj = super().__new__(cls, name=name, arguments=arguments) obj.p0 = p0 obj.p1 = as_tuple(p1) obj.p2 = p2 + obj.p3 = p3 return obj grid = Grid(shape=(4, 4)) @@ -1096,7 +1147,7 @@ def __new__(cls, name=None, arguments=None, p0=None, p1=None, p2=None): g = Function(name='g', grid=grid) func = MyDefFunction(name='foo', arguments=f.indexify(), - p0=f, p1=f, p2='bar') + p0=f, p1=f, p2='bar', p3=mapping_type({'field': f})) mapper = {f: g, f.indexify(): g.indexify()} func1 = uxreplace(func, mapper) @@ -1105,6 +1156,10 @@ def __new__(cls, name=None, arguments=None, p0=None, p1=None, p2=None): assert func1.p0 is g assert func1.p1 == (g,) assert func1.p2 == 'bar' + assert type(func1.p3) is mapping_type + assert func1.p3['field'] is g + assert func.p3['field'] is f + assert uxreplace(func1, mapper) is func1 def test_custom_def_function_reconstruction_no_aliasing(self): diff --git a/tests/test_visitors.py b/tests/test_visitors.py index bdbd5dd6da..e59d8798fb 100644 --- a/tests/test_visitors.py +++ b/tests/test_visitors.py @@ -11,7 +11,7 @@ FindNodes, FindSections, FindSymbols, FindWithin, IsPerfectIteration, Iteration, List, MapNodes, Node, Transformer, Uxreplace, printAST ) -from devito.symbolics import ListInitializer +from devito.symbolics import ListInitializer, retrieve_indexed from devito.types import Array, LocalObject, SpaceDimension, Symbol @@ -422,7 +422,8 @@ def test_find_symbols_with_duplicates(): # objects because they are constructed in two different places during # the CIRE pass r0a = exprs[0].output - r0b = exprs[1].expr.rhs.args[0].args[0] + r0b, = [i for i in retrieve_indexed(exprs[1].expr.rhs) + if i.name == r0a.name] assert r0a.function is r0b.function assert r0a.name == r0b.name == "r0" assert hash(r0a) == hash(r0b) From 928f4d7ef55daccc74f41582957676fb75d7c645 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 26 Oct 2026 17:10:16 +0000 Subject: [PATCH 18/19] compiler: Keep expression freezing separate from local naming Restore recursive unevaluation without comparison symbols or operand sorting. Keep ordinary Add/Mul printing and avoid new parentheses around nested sums. Preserve the earlier scoped substitution and reconstruction fixes. Extract numeric factor signs locally in the Devito printer, including signs outside the leading factor and products containing only minus-one factors. Cover pure freezing, equal operand instances and canonical emitted code; bring tutorial snapshots back to ordinary printing. Validate CPU/C++ symbolic and numerical cases, 91 strict notebook checks, CUDA naming/pickling and the reported CUDA/OpenACC compilation failures. The existing min-storage and bitcomp precision failures remain unchanged in scope; no numerical tolerances are relaxed. --- devito/ir/cgen/printer.py | 20 +++--- devito/passes/clusters/unevaluate.py | 22 ++---- devito/symbolics/manipulation.py | 21 +----- devito/symbolics/unevaluation.py | 7 -- examples/mpi/overview.ipynb | 2 +- examples/performance/00_overview.ipynb | 24 +++---- tests/test_symbolics.py | 94 +++++++++++++++++--------- 7 files changed, 95 insertions(+), 95 deletions(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index c770e2a25b..d83a0cb421 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -19,7 +19,6 @@ from devito.symbolics.extended_sympy import BitwiseAnd from devito.symbolics.inspection import has_integer_args, sympy_dtype from devito.symbolics.queries import q_leaf -from devito.symbolics.unevaluation import Add as UnevalAdd from devito.tools import ( ctypes_to_cstr, ctypes_vector_mapper, dtype_to_ctype, memoized_func ) @@ -281,13 +280,16 @@ def _print_Mod(self, expr): return '%'.join(args) def _print_Mul(self, expr): - args = [a for a in expr.args if a != -1] - neg = (len(expr.args) - len(args)) % 2 - - # SymPy's sign extraction rebuilds a plain Mul, losing frozen factor order - if args and args[0].is_Number and args[0].is_negative: - args[0] = -args[0] - neg ^= 1 + args, neg = [], False + # Extract numeric signs before SymPy can rebuild a plain Mul + for a in expr.args: + if a.is_Number and a.is_negative: + a = -a + neg = not neg + if a == 1: + continue + args.append(a) + args = args or [S.One] if len(args) > 1: term = super()._print_Mul(expr.func(*args, evaluate=False)) @@ -360,7 +362,7 @@ def _print_Add(self, expr, order=None): l = [] for term in terms: t = self._print(term) - if precedence(term) < PREC or (isinstance(expr, UnevalAdd) and term.is_Add): + if precedence(term) < PREC: l.extend(["+", f"({t})"]) elif t.startswith('-'): l.extend(["-", t[1:]]) diff --git a/devito/passes/clusters/unevaluate.py b/devito/passes/clusters/unevaluate.py index 26e65553e9..b895e1eafa 100644 --- a/devito/passes/clusters/unevaluate.py +++ b/devito/passes/clusters/unevaluate.py @@ -1,23 +1,11 @@ -from devito.ir import SymbolRegistry +from devito.ir import cluster_pass from devito.symbolics import unevaluate as _unevaluate -from devito.types import Temp __all__ = ['unevaluate'] -def unevaluate(clusters): - # Abstract local temporaries in definition order before choosing operand order - exprs = [e for c in clusters for e in c.exprs] - sregistry = SymbolRegistry() - sregistry.reserve(exprs) - # Lexical comparison must respect definition order across digit boundaries - width = len(str(len(exprs))) - mapper = {e.lhs: e.lhs._rebuild(name=sregistry.make_name( - prefix=f'_order{n:0{width}}_')) - for n, e in enumerate(exprs) if isinstance(e.lhs, Temp)} +@cluster_pass(mode='all') +def unevaluate(cluster): + exprs = [_unevaluate(e) for e in cluster.exprs] - processed = [] - for c in clusters: - exprs = [_unevaluate(e, mapper) for e in c.exprs] - processed.append(c.rebuild(exprs=exprs)) - return processed + return cluster.rebuild(exprs=exprs) diff --git a/devito/symbolics/manipulation.py b/devito/symbolics/manipulation.py index 0df752e197..11e085722a 100644 --- a/devito/symbolics/manipulation.py +++ b/devito/symbolics/manipulation.py @@ -628,27 +628,12 @@ def evalrel(func=min, input=None, assumptions=None): ceiling: UnevalCeil, floor: UnevalFloor} -def unevaluate(expr, mapper=None): - """ - Freeze code-generation operand order without reevaluating arithmetic. - - `mapper` supplies comparison symbols; original operands and grouping survive. - """ +def unevaluate(expr): + """Protect lowered arithmetic from reevaluation during reconstruction.""" if q_leaf(expr): return expr - args = [unevaluate(a, mapper) for a in expr.args] - if expr.func in (Add, Mul): - # Order comparison copies; preserve the original operands and grouping - keys = [uxreplace(a, mapper) for a in args] - originals = dict(zip(keys, args, strict=True)) - handle = expr.func(*keys, evaluate=False) - if expr.func is Add: - keys = handle.as_ordered_terms() - else: - keys = handle.as_ordered_factors() - # Factor ordering can split a negative coefficient into -1 and its magnitude - args = [originals.get(a, a) for a in keys] + args = [unevaluate(a) for a in expr.args] try: return uneval_mapper[expr.func](*args) diff --git a/devito/symbolics/unevaluation.py b/devito/symbolics/unevaluation.py index cb49840fcc..7e348ec518 100644 --- a/devito/symbolics/unevaluation.py +++ b/devito/symbolics/unevaluation.py @@ -13,10 +13,6 @@ def __new__(cls, *args, evaluate=None, **kwargs): class Add(sympy.Add, UnevaluableMixin): __new__ = UnevaluableMixin.__new__ - def as_ordered_terms(self, order=None, data=False): - # Lowered arithmetic already carries its code-generation order - return self.as_terms() if data else list(self.args) - class Mod(sympy.Mod, UnevaluableMixin): __new__ = UnevaluableMixin.__new__ @@ -30,9 +26,6 @@ def __new__(cls, *args, evaluate=None, **kwargs): return super().__new__(cls, *args, evaluate=False, **kwargs) - def as_ordered_factors(self, order=None): - return list(self.args) - class Pow(sympy.Pow, UnevaluableMixin): diff --git a/examples/mpi/overview.ipynb b/examples/mpi/overview.ipynb index 5376881b0d..8c209d3bc6 100644 --- a/examples/mpi/overview.ipynb +++ b/examples/mpi/overview.ipynb @@ -471,7 +471,7 @@ " #pragma omp simd aligned(u:64)\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " u[t1][x + 2][y + 2] = (-r0*u[t0][x + 2][y + 2] + r0*u[t0][x + 3][y + 2]) + 1;\n", + " u[t1][x + 2][y + 2] = -r0*u[t0][x + 2][y + 2] + r0*u[t0][x + 3][y + 2] + 1;\n", " }\n", " }\n", " STOP(section0,timers)\n", diff --git a/examples/performance/00_overview.ipynb b/examples/performance/00_overview.ipynb index 230fe642bf..3345a1ac4e 100644 --- a/examples/performance/00_overview.ipynb +++ b/examples/performance/00_overview.ipynb @@ -208,7 +208,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -283,7 +283,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -341,7 +341,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -431,7 +431,7 @@ " {\n", " for (int z = z0_blk1; z <= MIN(z_M, z0_blk1 + z0_blk1_size - 1); z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -522,7 +522,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 7][z + 4]) + (-8.33333333e-2F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 8][z + 4]) + (-6.66666667e-1F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 5][z + 4]))*r0[x][y][z];\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 5][z + 4]) + (-8.33333333e-2F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 8][z + 4]) + (8.33333333e-2F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r1)*(8.33333333e-2F*r1*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r1*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r1*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r1*u[t0][x + 4][y + 7][z + 4]))*r0[x][y][z];\n", " }\n", " }\n", " }\n", @@ -567,9 +567,9 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - "- u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + "- u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", "+ float r0 = 1.0F/h_y;\n", - "+ u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n", + "+ u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -604,7 +604,7 @@ " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", " float r0 = 1.0F/h_y;\n", - "- u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n", + "- u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n", "+ u[t1][x + 4][y + 4][z + 4] = (r0*r0)*(f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(6.94444444e-3F*(u[t0][x + 4][y][z + 4] + u[t0][x + 4][y + 8][z + 4]) + 4.44444444e-1F*(u[t0][x + 4][y + 2][z + 4] + u[t0][x + 4][y + 6][z + 4]) + 1.11111111e-1F*(-u[t0][x + 4][y + 1][z + 4] + u[t0][x + 4][y + 3][z + 4] + u[t0][x + 4][y + 5][z + 4] - u[t0][x + 4][y + 7][z + 4]) - 9.02777778e-1F*u[t0][x + 4][y + 4][z + 4])*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", @@ -747,7 +747,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r0[x][y + 2][z]/h_x - 6.66666667e-1F*r0[x + 1][y + 2][z]/h_x + 6.66666667e-1F*r0[x + 3][y + 2][z]/h_x - 8.33333333e-2F*r0[x + 4][y + 2][z]/h_x) + (8.33333333e-2F*r1[x + 2][y][z]/h_y - 6.66666667e-1F*r1[x + 2][y + 1][z]/h_y + 6.66666667e-1F*r1[x + 2][y + 3][z]/h_y - 8.33333333e-2F*r1[x + 2][y + 4][z]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*r0[x][y + 2][z]/h_x - 6.66666667e-1F*r0[x + 1][y + 2][z]/h_x + 6.66666667e-1F*r0[x + 3][y + 2][z]/h_x - 8.33333333e-2F*r0[x + 4][y + 2][z]/h_x + 8.33333333e-2F*r1[x + 2][y][z]/h_y - 6.66666667e-1F*r1[x + 2][y + 1][z]/h_y + 6.66666667e-1F*r1[x + 2][y + 3][z]/h_y - 8.33333333e-2F*r1[x + 2][y + 4][z]/h_y)*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -818,7 +818,7 @@ " #pragma omp simd aligned(f,u:32)\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F*r0[x][y][z]/h_x - 6.66666667e-1F*r0[x + 1][y][z]/h_x + 6.66666667e-1F*r0[x + 3][y][z]/h_x - 8.33333333e-2F*r0[x + 4][y][z]/h_x) + (8.33333333e-2F*r1[y][z]/h_y - 6.66666667e-1F*r1[y + 1][z]/h_y + 6.66666667e-1F*r1[y + 3][z]/h_y - 8.33333333e-2F*r1[y + 4][z]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*(8.33333333e-2F*r0[x][y][z]/h_x - 6.66666667e-1F*r0[x + 1][y][z]/h_x + 6.66666667e-1F*r0[x + 3][y][z]/h_x - 8.33333333e-2F*r0[x + 4][y][z]/h_x + 8.33333333e-2F*r1[y][z]/h_y - 6.66666667e-1F*r1[y + 1][z]/h_y + 6.66666667e-1F*r1[y + 3][z]/h_y - 8.33333333e-2F*r1[y + 4][z]/h_y)*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -882,7 +882,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", @@ -1020,7 +1020,7 @@ " {\n", " for (int z = z_m; z <= z_M; z += 1)\n", " {\n", - " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y) + (8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", + " u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n", " }\n", " }\n", " }\n", diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index 11cb229971..7681fef792 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -339,68 +339,99 @@ def test_extended_sympy_arithmetic(): assert ccode(-1 + bar) == '-1 + o->bar' -@pytest.mark.parametrize('cls, separator', [(UnevalAdd, ' + '), (UnevalMul, '*')]) -def test_lowered_arithmetic_printing(cls, separator): +@pytest.mark.parametrize('cls, regular, separator', [ + (UnevalAdd, sympy.Add, ' + '), + (UnevalMul, sympy.Mul, '*') +]) +def test_lowered_arithmetic_printing(cls, regular, separator): a, b, c = [dSymbol(i, dtype=np.float32) for i in ('a', 'b', 'c')] - # Lowered operands carry an order that must survive identifier changes - assert ccode(cls(c, a, b)) == separator.join(('c', 'a', 'b')) + # Lowering retains ordinary canonical printing + assert ccode(cls(c, a, b)) == separator.join(('a', 'b', 'c')) renamed = uxreplace(cls(c, a, b), {c: a, a: c}) - assert ccode(renamed) == separator.join(('a', 'c', 'b')) + assert ccode(renamed) == separator.join(('a', 'b', 'c')) # Ordinary index arithmetic retains its usual printing inside lowered math grid = Grid(shape=(8,)) x, = grid.dimensions f = Function(name='f', grid=grid) assert ccode(cls(c, f.indexed[x + 1], a)) == \ - separator.join(('c', 'f[x + 1]', 'a')) + ccode(regular(c, f.indexed[x + 1], a, evaluate=False)) assert ccode(x + 1) == 'x + 1' - # Each nested lowered operation retains its own operands and precedence - assert ccode(UnevalMul(b, UnevalAdd(c, a))) == 'b*(c + a)' - assert ccode(cls(b, cls(c, a))) == f'b{separator}(c{separator}a)' - assert ccode(UnevalAdd(b, sympy.Add(c, a, evaluate=False))) == 'b + (a + c)' + # Nested sums retain main's printing policy + assert ccode(UnevalMul(b, UnevalAdd(c, a))) == 'b*(a + c)' + assert ccode(cls(b, cls(c, a))) == \ + ccode(regular(b, regular(c, a, evaluate=False), evaluate=False)) + assert ccode(UnevalAdd(b, sympy.Add(c, a, evaluate=False))) == 'b + a + c' -def test_lowered_negative_product(): +@pytest.mark.parametrize('coefficient', [-2, sympy.Float(-2), sympy.Rational(-1, 2)]) +@pytest.mark.parametrize('position', [0, 1, 3]) +def test_lowered_negative_product(coefficient, position): a, b, c = [dSymbol(i, dtype=np.float32) for i in ('a', 'b', 'c')] + args = [c, a, b] + args.insert(position, coefficient) - # Extracting a negative coefficient must retain the frozen factors - assert ccode(UnevalMul(-2, c, a, b)) == '-2*c*a*b' - expr = -2*c*a*b + assert ccode(UnevalMul(*args)) == ccode(coefficient*a*b*c) + expr = coefficient*c*a*b assert ccode(unevaluate(expr)) == ccode(expr) -@pytest.mark.parametrize('cls, lowered, separator', [ - (sympy.Add, UnevalAdd, ' + '), - (sympy.Mul, UnevalMul, '*') +@pytest.mark.parametrize('coefficient', [-1, -2, sympy.Float(-2), sympy.Rational(-1, 2)]) +def test_lowered_product_sign(coefficient): + a, b, c = [dSymbol(i, dtype=np.float32) for i in ('a', 'b', 'c')] + expr = UnevalMul(-1, coefficient, c, a, b) + + assert ccode(expr) == ccode(-coefficient*a*b*c) + + +def test_lowered_sign_identity(): + assert ccode(UnevalMul(-1, -1)) == '1' + + +@pytest.mark.parametrize('cls, lowered', [ + (sympy.Add, UnevalAdd), + (sympy.Mul, UnevalMul) ]) -def test_freeze_grouping(cls, lowered, separator): +def test_freeze_grouping(cls, lowered): a, b, c = [Temp(name=n, dtype=np.float32) for n in ('r9', 'r10', 'r11')] - x, y, z = [Temp(name=n, dtype=np.float32) for n in ('r0', 'r1', 'r2')] expr = cls(c, cls(b, a, evaluate=False), evaluate=False) - frozen = unevaluate(expr, {a: x, b: y, c: z}) + frozen = unevaluate(expr) - assert frozen.args == (c, lowered(a, b)) - assert ccode(frozen) == f'r11{separator}(r9{separator}r10)' + assert frozen.args == (c, lowered(b, a)) + assert ccode(frozen) == ccode(expr) assert unevaluate(frozen) is frozen +def test_freeze_equal_instances(): + a = Symbol.__xnew__(Symbol, 'equal_frozen_operand') + b = Symbol.__xnew__(Symbol, 'equal_frozen_operand') + expr = sympy.Add(a, b, evaluate=False) + + frozen = unevaluate(expr) + + assert a == b and a is not b + assert frozen.args[0] is a + assert frozen.args[1] is b + + @pytest.mark.parametrize('dtype', [np.int32, np.int64]) @pytest.mark.parametrize('offset, expected', [(-1, ' - 1'), (1, ' + 1')]) def test_lowered_integer_offset(dtype, offset, expected): size = dSymbol('bx_size', dtype=dtype) other = dSymbol('a_size', dtype=dtype) - lowered = unevaluate(size + offset) - assert lowered.args == (size, offset) + expr = size + offset + lowered = unevaluate(expr) + assert lowered.args == expr.args assert ccode(lowered) == 'bx_size' + expected - # Renaming retains the order chosen before freezing + # Renaming preserves the tree, while printing remains canonical lowered = unevaluate(size + other + offset) renamed = uxreplace(lowered, {size: other, other: size}) - assert ccode(renamed) == 'bx_size + a_size' + expected + assert ccode(renamed) == 'a_size + bx_size' + expected # Result types do not affect how the arithmetic order is frozen value = dSymbol('value', dtype=np.float32) @@ -409,19 +440,20 @@ def test_lowered_integer_offset(dtype, offset, expected): assert ccode(unevaluate(expr)) == ccode(expr) -@pytest.mark.parametrize('cls, separator', [(sympy.Add, ' + '), (sympy.Mul, '*')]) +@pytest.mark.parametrize('cls', [sympy.Add, sympy.Mul]) @pytest.mark.parametrize('names', [('r9', 'r10', 'r11'), ('r21', 'r20', 'r19'), ('later', 'middle', 'earlier')]) -def test_lowered_temporary_order(cls, separator, names): +def test_noop_temporary_order(cls, names): grid = Grid(shape=(4,)) f, g, h, out = [Function(name=n, grid=grid) for n in ('f', 'g', 'h', 'out')] a, b, c = [Temp(name=n, dtype=np.float32) for n in names] + rhs = cls(a, b, c) op = Operator([Eq(a, cos(f)), Eq(b, sin(g)), Eq(c, cos(h)), - Eq(out, cls(a, b, c))], opt='noop') + Eq(out, rhs)], opt='noop') expr, = [i for i in FindNodes(Expression).visit(op) if i.write == out] - assert expr.expr.rhs.args == (a, b, c) - assert str(expr).endswith(f' = {separator.join(names)};') + assert tuple(i.name for i in expr.expr.rhs.args) == tuple(i.name for i in rhs.args) + assert ccode(expr.expr.rhs) == ccode(rhs) def test_integer_abs(): From 89a123f2ae8e265ede9bcaaff1e4cc7591444cbd Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Tue, 27 Oct 2026 07:46:36 +0000 Subject: [PATCH 19/19] compiler: Explain the local multiplication sign workaround Document that SymPy calls the standalone _keep_coeff helper, so UnevalMul cannot override its reconstruction. Explain why Devito extracts signs before delegating instead of patching SymPy or copying its multiplication printer. Comment only; Python syntax trees are identical. Ruff and diff checks pass. --- devito/ir/cgen/printer.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py index d83a0cb421..f7fcae920b 100644 --- a/devito/ir/cgen/printer.py +++ b/devito/ir/cgen/printer.py @@ -281,7 +281,11 @@ def _print_Mod(self, expr): def _print_Mul(self, expr): args, neg = [], False - # Extract numeric signs before SymPy can rebuild a plain Mul + # SymPy's CodePrinter extracts negative coefficients with the standalone + # _keep_coeff helper, which can rebuild a plain SymPy Mul. Since this + # is not a method, UnevalMul cannot override it. Extract signs before + # delegating to preserve the product's class without patching SymPy or + # copying its multiplication printer. for a in expr.args: if a.is_Number and a.is_negative: a = -a