From 3d9ac3e5f64bde1f7dd57e92b49958c04f61df14 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Tue, 15 Oct 2024 16:07:55 +0200 Subject: [PATCH 01/10] Reuse nutils hash for compiled function name This patch modifies the hash used for the name of the function produced by evaluable.compile, in order to reuse it as nutils hash. --- nutils/evaluable.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index 2eff1a2e5..02414b5a0 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -6282,8 +6282,9 @@ def collect(evaluable): print(' ' + line, file=script) print(' return ' + ret_fmt.format(*[v.py_expr for v in py_funcs]), file=script) script = script.getvalue() + script_hash = hashlib.sha1(script.encode('utf-8')).digest() - name = 'compiled_{}'.format(hashlib.sha256(script.encode('utf-8')).hexdigest()) + name = f'compiled_{script_hash.hex()}' if debug_flags.compile: print(script) @@ -6295,7 +6296,7 @@ def collect(evaluable): # Compile. eval(builtins.compile(script, name, 'exec'), globals) compiled = globals['compiled'] - compiled.__nutils_hash__ = types.nutils_hash(script) + compiled.__nutils_hash__ = script_hash return compiled From 5d30ed3c5e39f5e105d2a2e0e6a3fa10b44b0629 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 23 Sep 2024 13:16:43 +0200 Subject: [PATCH 02/10] Pull inflate through loopsum This patch adds a simplification for LoopSum, that promotes inflations for which the dofmaps are independent of index. --- nutils/evaluable.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index 02414b5a0..99b0ff4d5 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -4935,6 +4935,9 @@ def _simplified(self): return zeros_like(self) elif self.index not in self.func.arguments: return self.func * astype(self.index.length, self.func.dtype) + for axis, parts in self.func._inflations: + if not any(self.index in dofmap.arguments for dofmap in parts): + return util.sum(_inflate(loop_sum(func, self.index), dofmap, self.shape[axis], axis) for dofmap, func in parts.items()) return self.func._loopsum(self.index) def _takediag(self, axis1, axis2): From 0206a6cc9d3776155dbd0a2d6510b7679856b235 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Tue, 24 Sep 2024 17:22:16 +0200 Subject: [PATCH 03/10] Disable graphviz node background if ncalls == 0 This patch removes the background colour of graphviz nodes that are not visited once. For functions that are called more then once, this helps to distinguish nodes with cached values from those that evaluate in a short time span. --- nutils/evaluable.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index 99b0ff4d5..c407e6ecc 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -6365,7 +6365,7 @@ def _log_stats(func, stats): maxtime = builtins.max(n.metadata[1].time for n in node.walk(set())) tottime = builtins.sum(n.metadata[1].time for n in node.walk(set())) aggstats = tuple((key, builtins.sum(v.time for v in values), builtins.sum(v.ncalls for v in values)) for key, values in util.gather(n.metadata for n in node.walk(set()))) - fill_color = (lambda node: '0,{:.2f},1'.format(node.metadata[1].time/maxtime)) if maxtime else None + fill_color = (lambda node: '0,{:.2f},1'.format(node.metadata[1].time/maxtime) if node.metadata[1].ncalls else None) if maxtime else None if graphviz: node.export_graphviz(fill_color=fill_color, dot_path=graphviz) log.info('total time: {:.0f}ms\n'.format(tottime/1e6) + '\n'.join('{:4.0f} {} ({} calls, avg {:.3f} per call)'.format(t / 1e6, k, n, t / (1e6*n)) From cdeec92971768be9f9f3dc5936e989dfe5c66290 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 2 Oct 2024 16:57:57 +0200 Subject: [PATCH 04/10] Add Evaluable.argument_degree This patch adds the argument_degree method to evaluable arrays, which returns the highest power in a specific argument if the evaluable is polynomial in that argument; or raises a NotPolynomial error otherwise. --- nutils/evaluable.py | 73 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index c407e6ecc..3c65555d5 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -114,6 +114,11 @@ class ExpensiveEvaluationWarning(warnings.NutilsInefficiencyWarning): pass +class NotPolynomal(Exception): + def __init__(self, array, argument): + super().__init__(f'{array} is not polynomial in argument {argument.name!r}') + + class Evaluable(types.DataClass): 'Base class' @@ -217,6 +222,24 @@ def _compile_expression(self, py_self: _pyast.Variable, *args: _pyast.Expression return py_self.get_attr('evalf').call(*args) + def argument_degree(self, argument): + '''return the highest power of argument of self is polynomial, + or raise NotPolynomal otherwise.''' + # IMPORTANT: since we are tracking only the highest power, we cannot + # ever lower a power, e.g. via division. To see this, consider the sum + # of a 0th and 1st power, registered as degree 1. This would be lowered + # via division to 0, thereby suggesting that the evaluable is constant + # while in reality it is not polynomial. + if argument not in self.arguments: + return 0 + n = self._argument_degree(argument) + if n is None: + raise NotPolynomal(self, argument) + return n + + def _argument_degree(self, argument): + pass + class Tuple(Evaluable): @@ -1093,6 +1116,10 @@ def _intbounds_impl(self): def _const_uniform(self): return self.func._const_uniform + def _argument_degree(self, argument): + if argument not in self.length.arguments: + return self.func._argument_degree(argument) + class Transpose(Array): @@ -1295,6 +1322,9 @@ def _optimized_for_numpy(self): axes = [n for axis in self.axes for n in range(offsets[axis], offsets[axis+1])] return Assemble(transpose(self.func.func, axes), tuple(self.func.indices[i] for i in self.axes), self.shape) + def _argument_degree(self, argument): + return self.func._argument_degree(argument) + class Product(Array): @@ -1618,6 +1648,10 @@ def _intbounds_impl(self): extrema = [b1 and b2 and b1 * b2 for b1 in func1._intbounds for b2 in func2._intbounds] return min(extrema), max(extrema) + def _argument_degree(self, argument): + func1, func2 = self.funcs + return func1.argument_degree(argument) + func2.argument_degree(argument) + class Add(Array): @@ -1768,6 +1802,10 @@ def _intbounds_impl(self): lowers, uppers = zip(*[f._intbounds for f in self._terms]) return builtins.sum(lowers), builtins.sum(uppers) + def _argument_degree(self, argument): + func1, func2 = self.funcs + return max(func1.argument_degree(argument), func2.argument_degree(argument)) + class Einsum(Array): @@ -1907,6 +1945,9 @@ def _take(self, index, axis): def _takediag(self, axis1, axis2): return sum(_takediag(self.func, axis1, axis2), -2) + def _argument_degree(self, argument): + return self.func._argument_degree(argument) + class TakeDiag(Array): @@ -2026,6 +2067,10 @@ def _sum(self, axis): def _intbounds_impl(self): return self.func._intbounds + def _argument_degree(self, argument): + if argument not in self.indices.arguments: + return self.func._argument_degree(argument) + class Power(Array): @@ -2100,6 +2145,11 @@ def _take(self, index, axis): def _unravel(self, axis, shape): return Power(unravel(self.func, axis, shape), unravel(self.power, axis, shape)) + def _argument_degree(self, argument): + power, where = unalign(self.power.simplified) + if argument not in self.power.arguments and not power.ndim and isinstance(power, Constant) and power.value >= 0 and int(power.value) == power.value: + return self.func._argument_degree(argument) * int(power.value) + class Pointwise(Array): ''' @@ -3263,6 +3313,10 @@ def _intbounds_impl(self): lower, upper = self.func._intbounds return min(lower, 0), max(upper, 0) + def _argument_degree(self, argument): + if argument not in self.dofmap.arguments and argument not in self.length.arguments: + return self.func._argument_degree(argument) + class SwapInflateTake(Evaluable): @@ -3497,6 +3551,9 @@ def _loopsum(self, index): def _assparse(self): return tuple((*indices, indices[-1], values) for *indices, values in self.func._assparse) + def _argument_degree(self, argument): + return self.func._argument_degree(argument) + class Guard(Array): 'bar all simplifications' @@ -3695,6 +3752,10 @@ def _node(self, cache, subgraph, times, unique_loop_ids): def arguments(self): return frozenset({self}) + def _argument_degree(self, argument): + assert self == argument + return 1 + class IdentifierDerivativeTarget(DerivativeTargetBase): '''Virtual derivative target distinguished by an identifier. @@ -3901,6 +3962,10 @@ def _sum(self, axis): def _assparse(self): return tuple((*indices[:-1], *divmod(indices[-1], appendaxes(self.shape[-1], values.shape)), values) for *indices, values in self.func._assparse) + def _argument_degree(self, argument): + if argument not in self.sh1.arguments and argument not in self.sh2.arguments: + return self.func._argument_degree(argument) + class RavelIndex(Array): @@ -4986,6 +5051,10 @@ def _assparse(self): chunks.append(tuple(loop_concatenate(arr, self.index) for arr in (*elem_indices, elem_values))) return tuple(chunks) + def _argument_degree(self, argument): + if argument not in self.length.arguments: + return self.func._argument_degree(argument) + class _SizesToOffsets(Array): @@ -5144,6 +5213,10 @@ def _assparse(self): def _intbounds_impl(self): return self.func._intbounds + def _argument_degree(self, argument): + if argument not in self.start.arguments and argument not in self.stop.arguments and argument not in self.concat_length.arguments: + return self.func._argument_degree(argument) + class SearchSorted(Array): '''Find index of evaluable array into sorted numpy array.''' From e93c3b96a3e15c62a267549e81070c1c5725e1ef Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 23 Sep 2024 13:48:34 +0200 Subject: [PATCH 05/10] Add function.factor, evaluable.factor This patch adds function.factor and evaluable.factor for the creation of equivalent sparse polynomials. See the documentation of evaluable.factor for details. --- nutils/SI.py | 3 +- nutils/evaluable.py | 175 ++++++++++++++++++++++++++++++++++++++++ nutils/function.py | 8 ++ tests/test_evaluable.py | 62 ++++++++++++++ 4 files changed, 247 insertions(+), 1 deletion(-) diff --git a/nutils/SI.py b/nutils/SI.py index bfe4ec874..284f7bdf7 100644 --- a/nutils/SI.py +++ b/nutils/SI.py @@ -330,7 +330,8 @@ def register(*names, __table=__DISPATCH_TABLE): 'trace', 'ptp', 'amax', 'amin', 'max', 'min', 'mean', 'take', 'broadcast_to', 'transpose', 'getitem', 'opposite', 'jump', 'replace_arguments', 'linearize', 'derivative', 'integral', - 'sample', 'scatter', 'kronecker', 'real', 'imag', 'conjugate') + 'sample', 'scatter', 'kronecker', 'real', 'imag', 'conjugate', + 'factor') def __unary(op, *args, **kwargs): (dim0, arg0), = Quantity.__unpack(args[0]) return dim0.wrap(op(arg0, *args[1:], **kwargs)) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index 3c65555d5..8c2b0f2a1 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -5385,6 +5385,161 @@ def as_csr(array): return values, rowptr, colidx, ncols +@util.shallow_replace +def zero_all_arguments(value): + '''Replace all function arguments by zeros.''' + + if isinstance(value, Argument): + return zeros_like(value) + + +@log.withcontext +def factor(array): + '''Convert array to a sparse polynomial. + + This function forms the equivalent polynomial of an evaluable array, of + which the coefficients are already evaluated, with the aim to speed up + repeated evaluations in e.g. time or Newton loops. + + For example, if ``array`` is a quadratic function of ``arg``, then the + ``factor`` function returns the equivalent Taylor series:: + + array(arg) -> array(0) + darray_darg(0) arg + .5 arg d2array_darg2 arg + + The coefficients ``array(0)``, ``darray_darg(0)`` and ``d2array_darg2`` + are evaluated as sparse tensors. As such, the new representation retains + all sparsity of the original.''' + + array = array.as_evaluable_array.simplified + degree = {arg: array.argument_degree(arg) for arg in array.arguments} + log.info(f'analysing function of {", ".join(arg.name for arg in degree)} with highest power {max(degree.values())}') + + # PREPARATION. We construct the equivalent polynomial to the input array, + # parameterized by the m_coeffs (monomial coefficients) and m_args + # (monomial arguments) lists. For every index i, m_coeffs[i] is an an + # evaluable array of shape array.shape (the shape of the calling argument) + # plus all the shapes of m_args[i], in order, where every element of + # m_args[i] is an Argument object. The same argument may be repeated to + # form higher powers. The polynomial (not formed) would result from + # contracting the arguments in m_args[i] with the corresponding axes of + # m_coeffs[i], and adding the resulting monomials. + + m_coeffs = [] + m_args = [] + queue = [((), array)] + for args, func in queue: + func = func.simplified + m_args.append(args) + m_coeffs.append(zero_all_arguments(func).simplified) + for arg in func.arguments: # as m_args grows, fewer arguments will remain in func + # We keep only m_args that are alphabetically ordered to + # deduplicate permutations of the same argument set: + if not args or arg.name >= args[-1].name: + n = args.count(arg) + 1 # new monomial power of arg + assert n <= degree[arg] + queue.append(((*args, arg), derivative(func, arg) / float(n))) + log.info(f'constructing sparse polynomial of degree {max(len(args) for args in m_args)} with {len(m_args)} monomials:', + ', '.join('*'.join(f'{arg.name}^{n}' if n > 1 else arg.name for arg, n in collections.Counter(args).items()) or '1' for args in m_args)) + + # EVALUATION. We now form the polynomial, by accumulating evaluable + # monomials in which the coefficients are evaluated. To this end we + # evaluate the items of m_coeffs in sparse COO form, and contract each + # result with the relevant arguments from m_args. The sparse contraction is + # performed by first taking from the (flattened) argument the (raveled) + # sparse indices, and contracting the resulting vector with the relevant + # axis of the sparse values array. + + polynomial = zeros_like(array) + nvals = 0 + for args, (values, indices, shape) in log.iter.fraction('monomial', m_args, eval_coo(m_coeffs)): + indices, values = _sort_and_prune(shape, indices, values) + if not len(values): + continue + nvals += len(values) + factors = [constant(values)] + i = array.ndim + for arg, repeated in _iter_repeated(args): + if arg.ndim == 0: + factor = arg + else: + j = i + arg.ndim + factor = Take(_flat(arg), constant(numpy.ravel_multi_index(indices[i:j], shape[i:j]))) + i = j + if repeated: + # With reference to the example of the doc string, derivative + # will naively generate an evaluable of the form array'(arg) + # darg = darray_darg(0) darg + .5 arg d2array_darg2 darg + .5 + # darg d2array_darg2 arg. By Schwarz's theorem d2array_darg2 is + # symmetric, and the latter two terms are equal. However, since + # the row and column indices of d2array_darg2 differ, we cannot + # detect this equality but rather need to embed the information + # explicitly. To this end, factor produces an evaluable of the + # form array(arg) == array(0) + darray_darg(0) arg + .5 + # SymProd(arg, d2array_darg2 arg), which produces the desired + # derivative array'(arg) == darray_darg(0) + d2array_darg2 arg. + factor = SymProd(factors.pop(), factor) + factors.append(factor) + term = util.product(factors) + assert i == len(indices) + if array.ndim == 0: + monomial = Sum(term) + else: + index = constant(numpy.ravel_multi_index(indices[:array.ndim], shape[:array.ndim])) + size = constant(numpy.prod(shape[:array.ndim], dtype=int)) + monomial = Inflate(term, index, size) + while monomial.ndim < array.ndim: + monomial = Unravel(monomial, constant(shape[monomial.ndim-1]), constant(numpy.prod(shape[monomial.ndim:array.ndim], dtype=int))) + polynomial += monomial + log.info(f'factored function contains {nvals:,} doubles ({nvals>>17:,}MB)') + return polynomial + + +class SymProd(Array): + '''Symmetric product. + + The ``SymProd(a, b)`` operation behaves like ``Multiply(a, b)`` except that + its derivative is ``2 a db`` rather than ``a db + da b``, effectively + encoding the information that the two terms are equal. Care should be taken + that the ``SymProd`` operation is used only where this property applies.''' + + a: Array + b: Array + + def __post_init__(self): + assert not _any_certainly_different(self.a.shape, self.b.shape) + assert self.a.dtype == self.b.dtype + assert self.a.dtype != bool + assert not isinstance(self.b, SymProd) + + @property + def shape(self): + return self.a.shape + + @property + def dtype(self): + return self.a.dtype + + @property + def dependencies(self): + return self.a, self.b + + @property + def power(self): + return self.a.power + 1 if isinstance(self.a, SymProd) else 2 + + def evalf(self, a, b): + return a * b + + def _derivative(self, var, seen): + return self.dtype(self.power) * appendaxes(self.a, var.shape) * self.b._derivative(var, seen) + + def _optimized_for_numpy(self): + return self.a * self.b + + def _argument_degree(self, argument): + return self.a.argument_degree(argument) + self.b.argument_degree(argument) + + # AUXILIARY FUNCTIONS (FOR INTERNAL USE) @@ -5458,6 +5613,26 @@ def replace_inner(obj, root=None): return tuple(map(replace_inner, funcs)) +def _sort_and_prune(shape, indices, values): + assert len(shape) == len(indices) + if not shape: + v = values.sum() + return (), v[None] if v else numpy.zeros(0, int) + flat_index, inverse = numpy.unique(numpy.ravel_multi_index(indices, shape), return_inverse=True) + values = numeric.accumulate(values, [inverse], flat_index.shape) + nonzero = values != 0 + if nonzero.all(): + nonzero = slice(None) + return numpy.unravel_index(flat_index[nonzero], shape), values[nonzero] + + +def _iter_repeated(iterator): + last = object() + for item in iterator: + yield item, item == last + last = item + + class _Stats: def __init__(self, ncalls: int = 0, time: int = 0) -> None: diff --git a/nutils/function.py b/nutils/function.py index 5057bf404..eb8e76b78 100644 --- a/nutils/function.py +++ b/nutils/function.py @@ -2368,6 +2368,14 @@ def dotarg(__argname: str, *arrays: IntoArray, shape: Tuple[int, ...] = (), dtyp result = numpy.sum(_append_axes(result.transpose((*range(1, result.ndim), 0)), array.shape[1:]) * array, result.ndim-1) return result + +@nutils_dispatch +def factor(arg): + assert not arg.spaces + return _Unlower(evaluable.factor(arg), + spaces=set(), arguments=arg.arguments, lower_args=LowerArgs((), {}, {})) + + # BASES diff --git a/tests/test_evaluable.py b/tests/test_evaluable.py index b0b108b3f..879c9ec4f 100644 --- a/tests/test_evaluable.py +++ b/tests/test_evaluable.py @@ -1347,3 +1347,65 @@ def test_grad_variable_ncoeffs(self): evaluable.PolyGrad(eval_coeffs, 2).eval(ncoeffs=numpy.array(6)), poly.grad(const_coeffs, 2), ) + + +class factor(TestCase): + + def setUp(self): + index = evaluable.loop_index('i', 4) + values = evaluable.constant([1., 2., -2., -1.]) + dofs = index + evaluable.Range(values.shape[0]) + length = evaluable.constant(8) + basis = evaluable.loop_sum(evaluable._inflate(values, dofs, length, axis=0), index) + self.varg = evaluable.Argument('v', basis.shape, float) + self.v = (basis * self.varg).sum(0) + self.targ = evaluable.Argument('t', (), float) + self.t = 10. * self.targ + + def assertFactoredEqual(self, f, v=0, t=0): + self.assertEqual(f.argument_degree(self.varg), v) + self.assertEqual(f.argument_degree(self.targ), t) + g = evaluable.factor(f) + if not t: + tryargs = dict(v=numpy.zeros(8)), dict(v=numpy.ones(8)), dict(v=numpy.arange(8, dtype=float)) + elif not v: + tryargs = dict(t=0.), dict(t=1.), dict(t=-5.) + else: + tryargs = dict(v=numpy.arange(8, dtype=float), t=0.), dict(v=numpy.zeros(8, dtype=float), t=5.), dict(v=numpy.arange(8, dtype=float), t=5.) + for vderiv in range(v): + if vderiv: + f = evaluable.derivative(f, self.varg) + g = evaluable.derivative(g, self.varg) + f_ = f + g_ = f + for tderiv in range(t): + if tderiv: + f_ = evaluable.derivative(f_, self.targ) + g_ = evaluable.derivative(g_, self.targ) + with self.subTest(f'derivative v{vderiv}, t{tderiv}'): + F = evaluable.compile(f_) + G = evaluable.compile(g_) + for args in tryargs: + self.assertAllAlmostEqual(F(**args), G(**args)) + + def test_linear(self): + self.assertFactoredEqual(1. + self.v, v=1) + self.assertFactoredEqual(2. * self.v - 5. * self.t, v=1, t=1) + self.assertFactoredEqual(2. * self.v - self.t, v=1, t=1) + self.assertFactoredEqual(2. * self.v * self.t, v=1, t=1) + + def test_quadratic(self): + self.assertFactoredEqual(1. + self.v - self.v**2., v=2) + self.assertFactoredEqual(3. * self.v**2., v=2) + self.assertFactoredEqual(5. * self.t**2., t=2) + self.assertFactoredEqual(self.v * self.t**2. - 2. * self.v**2. * self.t, t=2, v=2) + + def test_cubic(self): + self.assertFactoredEqual(1. + self.v - self.v**2. + self.v**3., v=3) + self.assertFactoredEqual(-self.v**3., v=3) + self.assertFactoredEqual(3. * self.t**3., t=3) + self.assertFactoredEqual(self.t**2. * (self.v**2. - 2. * self.t), t=3, v=2) + + def test_not_polynomial(self): + with self.assertRaisesRegex(evaluable.NotPolynomal, "nutils.evaluable.Sign is not polynomial in argument 'v'"): + evaluable.factor(evaluable.Sign(self.v)) From eb6a3aeedd24a824886494075f5952fe724578de Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Fri, 25 Oct 2024 13:56:42 +0200 Subject: [PATCH 06/10] Replace System.trialshapes by argshapes This patch replaces the trialshapes dictionary in System, which contains only the shapes of the trial arguments, by the argshapes dictionary, which contains all of the evaluable's arguments. The reason for the change is that the broader set of arguments can be used for all sorts of argument checks prior to initiating a solve. --- nutils/solver.py | 12 ++++++------ tests/test_solver.py | 12 ++++++------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/nutils/solver.py b/nutils/solver.py index b7838b850..926a47c2d 100644 --- a/nutils/solver.py +++ b/nutils/solver.py @@ -234,8 +234,8 @@ def __init__(self, residual: Union[evaluable.AsEvaluableArray,Tuple[evaluable.As residuals = [evaluable.derivative(functional, argobjects[t]) for t in tests] self.is_symmetric = self.trials == tests - self.trialshapes = dict(zip(self.trials, evaluable.compile(tuple(argobjects[t].shape for t in self.trials))())) - self.__trial_offsets = numpy.cumsum([0] + [numpy.prod(self.trialshapes[t], dtype=int) for t in self.trials]) + self.argshapes = dict(zip(argobjects.keys(), evaluable.compile(tuple(arg.shape for arg in argobjects.values()))())) + self.__trial_offsets = numpy.cumsum([0] + [numpy.prod(self.argshapes[t], dtype=int) for t in self.trials]) value = functional if self.is_symmetric else () block_vector = [evaluable._flat(res) for res in residuals] @@ -285,7 +285,7 @@ def prepare_solution_vector(self, arguments: Dict[str, numpy.ndarray], constrain x = numpy.empty(self.__trial_offsets[-1], self.dtype) iscons = numpy.empty(self.__trial_offsets[-1], dtype=bool) for trial, i, j in zip(self.trials, self.__trial_offsets, self.__trial_offsets[1:]): - trialshape = self.trialshapes[trial] + trialshape = self.argshapes[trial] trialarg = x[i:j].reshape(trialshape) trialarg[...] = arguments.get(trial, 0) c = constrain.get(trial, False) @@ -301,7 +301,7 @@ def prepare_solution_vector(self, arguments: Dict[str, numpy.ndarray], constrain @property def _trial_info(self): - return ' and '.join(t + ' (' + ','.join(map(str, shape)) + ')' for t, shape in self.trialshapes.items()) + return ' and '.join(t + ' (' + ','.join(map(str, self.argshapes[t])) + ')' for t in self.trials) MethodIter = Iterator[Tuple[Dict[str, numpy.ndarray], float, float]] @@ -580,7 +580,7 @@ def __str__(self): def __call__(self, system, *, arguments: Dict[str, numpy.ndarray] = {}, constrain: Dict[str, numpy.ndarray] = {}, linargs: Dict[str, Any] = {}) -> System.MethodIter: x, iscons, arguments = system.prepare_solution_vector(arguments, constrain) linargs = _copy_with_defaults(linargs, rtol=1-3) - djac = self._assemble_inertia_matrix(system.trialshapes, arguments) + djac = self._assemble_inertia_matrix([(t, system.argshapes[t]) for t in system.trials], arguments) first = _First() while True: @@ -592,7 +592,7 @@ def __call__(self, system, *, arguments: Dict[str, numpy.ndarray] = {}, constrai x -= (jac + djac / timestep).solve_leniently(res, constrain=iscons, **linargs) def _assemble_inertia_matrix(self, trialshapes, arguments): - argobjs = [evaluable.Argument(t, tuple(map(evaluable.constant, shape)), float) for t, shape in trialshapes.items()] + argobjs = [evaluable.Argument(t, tuple(map(evaluable.constant, shape)), float) for t, shape in trialshapes] djacobians = [[evaluable._flat(evaluable.derivative(evaluable._flat(res), argobj).simplified, 2) for argobj in argobjs] for res in self.inertia] djac_blocks = evaluable.compile(tuple(tuple(map(evaluable.as_csr, row)) for row in djacobians)) return matrix.assemble_block_csr(djac_blocks(**arguments)) diff --git a/tests/test_solver.py b/tests/test_solver.py index 2841768ce..a97249fab 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -345,7 +345,7 @@ def test_constant(self): self.assertFalse(sys.is_symmetric) self.assertTrue(sys.is_constant) self.assertTrue(sys.is_constant_matrix) - self.assertEqual(sys.trialshapes, {'u': (10,)}) + self.assertEqual(sys.argshapes, {'u': (10,), 'v': (10,)}) args = {'u': numpy.arange(10, dtype=float)} mat, vec, val = sys.assemble(arguments=args) self.assertAllAlmostEqual(mat.export('dense'), function.eval(f.derivative('v').derivative('u'), **args)) @@ -365,7 +365,7 @@ def test_constant_symmetric(self): self.assertTrue(sys.is_symmetric) self.assertTrue(sys.is_constant) self.assertTrue(sys.is_constant_matrix) - self.assertEqual(sys.trialshapes, {'u': (10,)}) + self.assertEqual(sys.argshapes, {'u': (10,)}) args = {'u': numpy.arange(10, dtype=float)} mat, vec, val = sys.assemble(arguments=args) self.assertAllAlmostEqual(mat.export('dense'), function.eval(f.derivative('u').derivative('u'), **args)) @@ -388,7 +388,7 @@ def test_constant_matrix(self): self.assertFalse(sys.is_symmetric) self.assertFalse(sys.is_constant) self.assertTrue(sys.is_constant_matrix) - self.assertEqual(sys.trialshapes, {'u': (10,)}) + self.assertEqual(sys.argshapes, {'u': (10,), 'v': (10,), 't': ()}) args = {'u': numpy.arange(10, dtype=float), 't': 5} mat, vec, val = sys.assemble(arguments=args) self.assertAllAlmostEqual(mat.export('dense'), function.eval(f.derivative('v').derivative('u'), **args)) @@ -410,7 +410,7 @@ def test_linear(self): self.assertFalse(sys.is_symmetric) self.assertFalse(sys.is_constant) self.assertFalse(sys.is_constant_matrix) - self.assertEqual(sys.trialshapes, {'u': (10,)}) + self.assertEqual(sys.argshapes, {'u': (10,), 'v': (10,), 't': ()}) args = {'u': numpy.arange(10, dtype=float), 't': 5} mat, vec, val = sys.assemble(arguments=args) self.assertAllAlmostEqual(mat.export('dense'), function.eval(f.derivative('v').derivative('u'), **args)) @@ -428,7 +428,7 @@ def test_nonlinear(self): self.assertFalse(sys.is_symmetric) self.assertFalse(sys.is_constant) self.assertFalse(sys.is_constant_matrix) - self.assertEqual(sys.trialshapes, {'u': (10,)}) + self.assertEqual(sys.argshapes, {'u': (10,), 'v': (10,), 't': ()}) args = {'u': numpy.arange(10, dtype=float), 't': 5} mat, vec, val = sys.assemble(arguments=args) self.assertAllAlmostEqual(mat.export('dense'), function.eval(f.derivative('v').derivative('u'), **args)) @@ -444,7 +444,7 @@ def test_nonlinear_symmetric(self): self.assertTrue(sys.is_symmetric) self.assertFalse(sys.is_constant) self.assertFalse(sys.is_constant_matrix) - self.assertEqual(sys.trialshapes, {'u': (10,)}) + self.assertEqual(sys.argshapes, {'u': (10,)}) args = {'u': numpy.arange(10, dtype=float)} mat, vec, val = sys.assemble(arguments=args) self.assertAllAlmostEqual(mat.export('dense'), function.eval(f.derivative('u').derivative('u'), **args)) From 8a4ffdd81fe0ffeb5b3ae54a3f00ba6769503535 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 28 Oct 2024 16:16:37 +0100 Subject: [PATCH 07/10] Rename timetarget->timearg, historysuffix->suffix This patch renames the System.step parameters timetarget to timearg and historysuffix to suffix for brevity. --- examples/burgers.py | 2 +- examples/cahnhilliard.py | 2 +- examples/cylinderflow.py | 2 +- nutils/solver.py | 28 ++++++++++++++-------------- tests/test_solver.py | 2 +- 5 files changed, 18 insertions(+), 18 deletions(-) diff --git a/examples/burgers.py b/examples/burgers.py index a5d4586b4..7e072efea 100644 --- a/examples/burgers.py +++ b/examples/burgers.py @@ -64,7 +64,7 @@ def main(nelems: int = 40, export.triplot('solution.png', x[:,numpy.newaxis], u, tri=bezier.tri, hull=bezier.hull, clim=(0, 1)) if args['t'] >= endtime: break - args = system.step(timestep=timestep, arguments=args, timetarget='t', historysuffix='0', tol=newtontol) + args = system.step(timestep=timestep, arguments=args, timearg='t', suffix='0', tol=newtontol) return args diff --git a/examples/cahnhilliard.py b/examples/cahnhilliard.py index 7d55fdfde..8f45e3944 100644 --- a/examples/cahnhilliard.py +++ b/examples/cahnhilliard.py @@ -210,7 +210,7 @@ def main(size: Length = parse('10cm'), E = numpy.stack(function.eval([nrg_mix, nrg_iface, nrg_wall], **args)) log.user('energy: {0:,.0μJ/m} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(numpy.sum(E), 100*E/numpy.sum(E))) - args = system.step(timestep=1., timetarget='t', historysuffix='0', arguments=args, tol=1, maxiter=5) + args = system.step(timestep=1., timearg='t', suffix='0', arguments=args, tol=1, maxiter=5) with export.mplfigure('phase.png') as fig: ax = fig.add_subplot(aspect='equal', xlabel='[mm]', ylabel='[mm]') diff --git a/examples/cylinderflow.py b/examples/cylinderflow.py index 7e6ec1b01..856b6ebc6 100644 --- a/examples/cylinderflow.py +++ b/examples/cylinderflow.py @@ -147,7 +147,7 @@ def main(nelems: int = 99, for _ in steps: treelog.info(f'velocity divergence: {div.eval(**args):.0e}') - args = system.step(timestep=timestep, timetarget='t', historysuffix='0', arguments=args, constrain=cons, tol=1e-10) + args = system.step(timestep=timestep, timearg='t', suffix='0', arguments=args, constrain=cons, tol=1e-10) postprocess(args) return args, numpy.sqrt(domain.integral('∇_k(u_k)^2 dV' @ ns, degree=2)) diff --git a/nutils/solver.py b/nutils/solver.py index 926a47c2d..9f1afaa10 100644 --- a/nutils/solver.py +++ b/nutils/solver.py @@ -365,21 +365,21 @@ def solve(self, *, arguments: Dict[str, numpy.ndarray] = {}, constrain: Dict[str log.info(f'optimal value: {val:.1e}') return arguments - def step(self, *, timestep: float, timetarget: str, historysuffix: str, arguments: Dict[str, numpy.ndarray], **solveargs) -> Dict[str, numpy.ndarray]: + def step(self, *, timestep: float, timearg: str, suffix: str, arguments: Dict[str, numpy.ndarray], **solveargs) -> Dict[str, numpy.ndarray]: '''Advance a time step. - This method is best described by an example. Let ``timetarget`` equal - 't' and ``historysuffix`` '0', and let our system's trial arguments be - 'u' and 'v'. This method then creates argument copies 't0', 'u0', 'v0', - advances 't' by ``timestep``, and solves for the new 'u' and 'v'. + This method is best described by an example. Let ``timearg`` equal 't' + and ``suffix`` '0', and let our system's trial arguments be 'u' and 'v'. + This method then creates argument copies 't0', 'u0', 'v0', advances 't' + by ``timestep``, and solves for the new 'u' and 'v'. Parameters ---------- - timetarget + timearg Name of the scalar argument that tracks the time. timestep Size of the time increment. - historysuffix + suffix String suffix to add to argument names to denote their value at the beginning of the time step. arguments @@ -392,18 +392,18 @@ def step(self, *, timestep: float, timetarget: str, historysuffix: str, argument arguments = arguments.copy() for trial in self.trials: if trial in arguments: - arguments[trial + historysuffix] = arguments[trial] - time = arguments.get(timetarget, 0.) - arguments[timetarget + historysuffix] = time - arguments[timetarget] = time + timestep + arguments[trial + suffix] = arguments[trial] + time = arguments.get(timearg, 0.) + arguments[timearg + suffix] = time + arguments[timearg] = time + timestep try: return self.solve(arguments=arguments, **solveargs) except (SolverError, matrix.MatrixError) as e: log.error(f'error: {e}; retrying with timestep {timestep/2}') with log.context('tic'): - halfway_arguments = self.step(timestep=timestep/2, timetarget=timetarget, historysuffix=historysuffix, arguments=arguments, **solveargs) + halfway_arguments = self.step(timestep=timestep/2, timearg=timearg, suffix=suffix, arguments=arguments, **solveargs) with log.context('toc'): - return self.step(timestep=timestep/2, timetarget=timetarget, historysuffix=historysuffix, arguments=halfway_arguments, **solveargs) + return self.step(timestep=timestep/2, timearg=timearg, suffix=suffix, arguments=halfway_arguments, **solveargs) @cache.function @log.withcontext @@ -857,7 +857,7 @@ def thetamethod(target, residual, inertia, timestep: float, theta: float, *, lhs linesearch = newtonargs.pop('linesearch', NormBased()) method = None if system.is_linear else Newton() if linesearch is None else LinesearchNewton(strategy=linesearch, **newtonargs) return _thetamethod(system, arguments, - timestep=timestep, timetarget=timetarget, historysuffix=historysuffix, constrain=constrain or {}, tol=newtontol, method=method) + timestep=timestep, timearg=timetarget, suffix=historysuffix, constrain=constrain or {}, tol=newtontol, method=method) def _thetamethod(system, arguments, **stepargs): diff --git a/tests/test_solver.py b/tests/test_solver.py index a97249fab..13f7b5822 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -562,5 +562,5 @@ def setUp(self): self.system = solver.System(residual, trial='u', test='v') def test_step(self): - args = self.system.step(timestep=100, timetarget='t', historysuffix='0', arguments=self.arguments, method=solver.LinesearchNewton(), tol=1e-10) + args = self.system.step(timestep=100, timearg='t', suffix='0', arguments=self.arguments, method=solver.LinesearchNewton(), tol=1e-10) self.assertAlmostEqual64(args['u'], 'eNpzNBA1NjHuNHQ3FDsTfCbAuNz4nUGZgeyZiDOZxlONmQwU9W3OFJ/pNQAADZIOPA==') From 3793dcb914c827176a1fb688803a5a14f9bcfc11 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Tue, 15 Oct 2024 15:22:18 +0200 Subject: [PATCH 08/10] Add timesteptarget argument to System.step --- nutils/solver.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/nutils/solver.py b/nutils/solver.py index 9f1afaa10..77d285f50 100644 --- a/nutils/solver.py +++ b/nutils/solver.py @@ -365,7 +365,7 @@ def solve(self, *, arguments: Dict[str, numpy.ndarray] = {}, constrain: Dict[str log.info(f'optimal value: {val:.1e}') return arguments - def step(self, *, timestep: float, timearg: str, suffix: str, arguments: Dict[str, numpy.ndarray], **solveargs) -> Dict[str, numpy.ndarray]: + def step(self, *, timestep: float, suffix: str, arguments: Dict[str, numpy.ndarray], timearg: Optional[str] = None, timesteparg: Optional[str] = None, **solveargs) -> Dict[str, numpy.ndarray]: '''Advance a time step. This method is best described by an example. Let ``timearg`` equal 't' @@ -389,21 +389,28 @@ def step(self, *, timestep: float, timearg: str, suffix: str, arguments: Dict[st Remaining keyword arguments are passed on to the ``solver`` method. ''' + if not timearg and not timesteparg: + raise ValueError('either timearg or timesteparg should be specified') arguments = arguments.copy() for trial in self.trials: if trial in arguments: arguments[trial + suffix] = arguments[trial] - time = arguments.get(timearg, 0.) - arguments[timearg + suffix] = time - arguments[timearg] = time + timestep + if timesteparg: + arguments[timesteparg] = timestep + if timearg: + time = arguments.get(timearg, 0.) + arguments[timearg + suffix] = time + arguments[timearg] = time + timestep try: return self.solve(arguments=arguments, **solveargs) except (SolverError, matrix.MatrixError) as e: + if timearg not in self.argshapes and timesteparg not in self.argshapes: + raise log.error(f'error: {e}; retrying with timestep {timestep/2}') with log.context('tic'): - halfway_arguments = self.step(timestep=timestep/2, timearg=timearg, suffix=suffix, arguments=arguments, **solveargs) + halfway_arguments = self.step(timestep=timestep/2, timearg=timearg, timesteparg=timesteparg, suffix=suffix, arguments=arguments, **solveargs) with log.context('toc'): - return self.step(timestep=timestep/2, timearg=timearg, suffix=suffix, arguments=halfway_arguments, **solveargs) + return self.step(timestep=timestep/2, timearg=timearg, timesteparg=timesteparg, suffix=suffix, arguments=halfway_arguments, **solveargs) @cache.function @log.withcontext From bf63f37a758a658e1c8b1e600150f27594d1a2de Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Tue, 15 Oct 2024 15:24:11 +0200 Subject: [PATCH 09/10] Modify cahnhilliard example --- examples/cahnhilliard.py | 126 ++++++++++++++++----------------------- 1 file changed, 52 insertions(+), 74 deletions(-) diff --git a/examples/cahnhilliard.py b/examples/cahnhilliard.py index 8f45e3944..9485acb22 100644 --- a/examples/cahnhilliard.py +++ b/examples/cahnhilliard.py @@ -3,19 +3,8 @@ from nutils.expression_v2 import Namespace from nutils.SI import Length, Time, Density, Tension, Energy, Pressure, Velocity, parse import numpy -import itertools -import enum import treelog as log -# NOTE: This script uses dimensional quantities and requires the nutils.SI -# module, which is installed separate from the the core nutils. - - -class stab(enum.Enum): - nonlinear = '.25 dφ^2 (1 - φ^2 + φ dφ (2 / 3) - dφ^2 / 6)' - linear = '.25 dφ^2 (2 + 2 abs(φ0) - (φ + φ0)^2)' - none = '0' # not unconditionally stable - def main(size: Length = parse('10cm'), epsilon: Length = parse('1mm'), @@ -26,12 +15,12 @@ def main(size: Length = parse('10cm'), nelems: int = 0, etype: str = 'rectilinear', degree: int = 1, - timestep: Time = parse('.5s'), + timestep: Time = parse('.1s'), tol: Energy/Length = parse('1nJ/m'), endtime: Time = parse('1min'), seed: int = 0, circle: bool = True, - stab: stab = stab.linear, + stable: bool = False, showflux: bool = True): '''Unmixing of immiscible fluids @@ -78,7 +67,7 @@ def main(size: Length = parse('10cm'), For stability we wish for the perturbation `δψ` to be such that the time discrete system preserves the energy dissipation property `E(φ) ≤ E(φ0)` for - any timestep `dt`. To derive suitable perturbation terms to this effect, we + any timestep `dt`. To derive a suitable perturbation term to this effect, we define without loss of generality `δψ'(φ, φ0) = .5 (φ - φ0) f(φ, φ0)` and derive the following condition for unconditional stability: @@ -88,20 +77,8 @@ def main(size: Length = parse('10cm'), The inequality holds true if the perturbation `f` is bounded from below such that `f(φ, φ0) ≥ 1 - φ² - .5 (φ + φ0)²`. To keep the energy minima at the pure phases we additionally impose that `f(±1, φ0) = 0`, and select `1 - φ²` - as a suitable upper bound which we will call "nonlinear". - - We next observe that `η` is linear in `φ` if `f(φ, φ0) = g(φ0) - φ² - (φ + - φ0)²` for any function `g`, which dominates if `g(φ0) ≥ 1 + .5 (φ + φ0)²`. - While this cannot be made to hold true for all `φ`, we make it hold for `-√2 - ≤ φ, φ0 ≤ √2` by defining `g(φ0) = 2 + 2 |φ0| + φ0²`, which we will call - "linear". This scheme further satisfies a weak minima preservation `f(±1, - ±|φ0|) = 0`. - - We have thus arrived at the three stabilization schemes implemented here: - - - nonlinear: `f(φ, φ0) = 1 - φ²` - - linear: `f(φ, φ0) = 2 + 2 |φ0| - 2 φ (φ + φ0)` - - none: `f(φ, φ0) = 0` (not unconditionally stable) + as a suitable upper bound. For unconditional stability we thus obtained the + perturbation gradient `δψ'(φ, φ0) = .5 (φ - φ0) (1 - φ²)`. Finally, we observe that the weak formulation: @@ -112,8 +89,8 @@ def main(size: Length = parse('10cm'), F(φ, φ0, η) := E(φ) + ∫_Ω [ .5 dt J(η)·∇η + δψ(φ, φ0) σ / ε - η (φ - φ0) ] - For this reason, the `stab` enum in this script defines the stabilizing term - `δψ`, rather than `f`, allowing Nutils to construct the residuals through + For this reason, this script defines the stabilizing term `δψ`, rather than + its derivative `δψ'`, allowing Nutils to construct the residuals through symbolic differentiation. Parameters @@ -148,10 +125,10 @@ def main(size: Length = parse('10cm'), Random seed for the initial condition. circle Select circular domain as opposed to a unit square. - stab - Stabilization method (linear/nonlinear/none). + stable + Enable unconditional stability at the expense of dissipation. showflux - Overlay flux vectors on phase plot + Overlay flux vectors on phase plot. ''' nmin = round(size / epsilon) @@ -169,19 +146,13 @@ def main(size: Length = parse('10cm'), else: domain, geom = mesh.unitsquare(nelems, etype) - basis = domain.basis('std', degree=degree) - bezier = domain.sample('bezier', 5) # sample for surface plots - if showflux: - grid = domain.locate(geom, numeric.simplex_grid([1, 1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5) # sample for quivers - ns = Namespace() ns.x = size * geom ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) - ns.add_field(('φ', 'φ0'), basis) - ns.add_field('η', basis * stens / epsilon) # basis scaling to give η the required unit - ns.add_field(('t', 't0')) - ns.t *= timestep - ns.t0 *= timestep + ns.add_field(('φ', 'φ0'), domain.basis('std', degree=degree)) + ns.add_field('η', domain.basis('std', degree=degree) * stens / epsilon) # basis scaling to give η the required unit + ns.add_field('dt') + ns.dt *= timestep ns.ε = epsilon ns.σ = stens ns.σmean = (wtensp + wtensn) / 2 @@ -189,39 +160,46 @@ def main(size: Length = parse('10cm'), ns.σwall = 'σmean + φ σdiff' ns.dφ = 'φ - φ0' ns.ψ = '.25 (φ^2 - 1)^2' - ns.δψ = stab.value + ns.δψ = '.25 dφ^2 (1 - φ0^2 / 6 - φ φ0 / 3 - φ^2 / 2)' if stable else '0' ns.M = mobility ns.J_i = '-M ∇_i(η)' - ns.dt = 't - t0' + ns.v_i = 'φ J_i' - nrg_mix = domain.integral('(ψ σ / ε) dV' @ ns, degree=degree*4) - nrg_iface = domain.integral('.5 σ ε ∇_k(φ) ∇_k(φ) dV' @ ns, degree=degree*4) - nrg_wall = domain.boundary.integral('σwall dS' @ ns, degree=degree*2) - nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(δψ σ / ε - η dφ + .5 dt J_k ∇_k(η)) dV' @ ns, degree=degree*4) + nrg_mix = function.factor(domain.integral('(ψ σ / ε) dV' @ ns, degree=degree*4)) + nrg_iface = function.factor(domain.integral('.5 σ ε ∇_k(φ) ∇_k(φ) dV' @ ns, degree=degree*4)) + nrg_wall = function.factor(domain.boundary.integral('σwall dS' @ ns, degree=degree*2)) + nrg = nrg_mix + nrg_iface + nrg_wall + function.factor(domain.integral('(δψ σ / ε - η dφ + .5 dt J_k ∇_k(η)) dV' @ ns, degree=degree*4)) - numpy.random.seed(seed) - args = dict(φ=numpy.random.normal(0, .5, basis.shape)) # initial condition + bezier = domain.sample('bezier', 5) # sample for surface plots + bezier_x = bezier.eval(ns.x) + bezier_φ = function.factor(bezier(ns.φ)) + if showflux: + grid = domain.locate(geom, numeric.simplex_grid([1, 1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5) # sample for quivers + grid_x = grid.eval(ns.x) + grid_v = function.factor(grid(ns.v)) system = System(nrg / tol, trial='φ,η') + numpy.random.seed(seed) + args = dict(φ=numpy.random.normal(0, .5, system.argshapes['φ'])) # initial condition + with log.iter.fraction('timestep', range(round(endtime / timestep))) as steps: for istep in steps: E = numpy.stack(function.eval([nrg_mix, nrg_iface, nrg_wall], **args)) log.user('energy: {0:,.0μJ/m} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(numpy.sum(E), 100*E/numpy.sum(E))) - args = system.step(timestep=1., timearg='t', suffix='0', arguments=args, tol=1, maxiter=5) + args = system.step(timestep=1., timesteparg='dt', suffix='0', arguments=args, tol=1, maxiter=5) with export.mplfigure('phase.png') as fig: ax = fig.add_subplot(aspect='equal', xlabel='[mm]', ylabel='[mm]') - x, φ = bezier.eval(['x_i', 'φ'] @ ns, **args) - im = ax.tripcolor(*(x/'mm').T, bezier.tri, φ, shading='gouraud', cmap='coolwarm') + im = ax.tripcolor(*(bezier_x/'mm').T, bezier.tri, function.eval(bezier_φ, **args), shading='gouraud', cmap='coolwarm') im.set_clim(-1, 1) fig.colorbar(im) if showflux: - x, v = grid.eval(['x_i', 'φ J_i'] @ ns, **args) + v = function.eval(grid_v, **args) log.info('largest flux: {:.2mm/s}'.format(numpy.max(numpy.hypot(v[:, 0], v[:, 1])))) - ax.quiver(*(x/'mm').T, *(v/'m/s').T) + ax.quiver(*(grid_x/'mm').T, *(v/'m/s').T) ax.autoscale(enable=True, axis='both', tight=True) return args @@ -240,33 +218,33 @@ def test_square(self): args = main(epsilon=parse('5cm'), mobility=parse('1μL*s/kg'), nelems=3, degree=2, timestep=parse('1h'), endtime=parse('2h'), circle=False) with self.subTest('concentration'): self.assertAlmostEqual64(args['φ'], ''' - eNoBYgCd/zE1EjX1NaA2+TXiMxkz0TS9NL01ajaRNZoxYNElNRM1LDUlNZQw0cqgysI1nTWcNN4xLsuk - ybDJvDWaNTQ07s7nysnJ6ckPNQY1CzNozKjK58kOysQ0zTQKM73M3coVyjfKR9cuPg==''') + eNoBYgCd/y41EjX2NZ829DXcMxUz0jTANL41ajaNNZox/9EoNRY1LDUkNZAw1cqnysI1njWdNNkxMMuk + ybDJuTWXNTE0587oysjJ58kSNQM1ATNqzKjK58kNytA00DQJM8bM38oTyjfKbwku0w==''') with self.subTest('chemical-potential'): self.assertAlmostEqual64(args['η'], ''' - eNoBYgCd/3TIkccBNkQ6IDqIN4HMF8cSx9Y02DmdOVHLMcecxxLIEjQUOAHOa8a1xWw3izb1M9kzPMc0 - xmnGpzibODY359ETyJHHp8hbyWU2xzZSydfIOsrNyo3GjMjAyyXIm8hkzD3K1ggxvA==''') + eNoBYgCd/1PIicccNko6IzqRNwzMEccMx/M05TmfOTTLMceexwzI1TMiOMLNZsa3xXU3rjZdNE4zO8cr + xlrGoziaOEA3os8VyJLHk8hlyTw2sDZXydPISMoPy5zGe8i7yzfIncgAzGLKwgYwXw==''') def test_multipatchcircle(self): args = main(epsilon=parse('5cm'), mobility=parse('1μL*s/kg'), nelems=3, etype='multipatch', degree=2, timestep=parse('1h'), endtime=parse('2h')) with self.subTest('concentration'): self.assertAlmostEqual64(args['φ'], ''' - eNoNwk1I02EcB/AGnYQCQcvDBqYYVGz7P8/v93QIyhdYK8lyDktKx8ouQ0hz6EEdjOpWdIqCXrAONrAg - AhGGBSHhxe/v+f//8yVpWIGXgk5DSAjf+HzyZsKkTZv5xdf5PN01A8aYCndwgrKq2hw0bzjEtfTdmfIK - /JBX6Z2eiN6zz+Uf+bSt76g635WoFHVO59SfyFE3KxsIqU2nyZmMOF6rjclbeSTDcmNfv3RLFEEoXEIe - s/iKT3iMJ3iJaczgPUL2s4wISQVFvEBKG92oLkSSblGCEuezvEgnKKOueWXp46tczz/pMI2qW94DkzPt - 5pDJ8yx16Sqzw/M8xt+orCvObb7IB7hAazql5laDtK4zuqSGnd3lTrcuHCid9pLumh20W7JQivkfXbFL - ckq+oDV6LHzTv+8esc0iOKOeOqlIjT9oe6SMBBE906/V35N52y9JjvMPukxLaqh03DasxN0tGZD/2JaU - /MYCpvABWZzDFaQxjl60YAgFdKMT7XiFgOwBQRG+vg==''') + eNoNz01IlFEUBmByEcVsWkiBoKHYoh9nvnvPOa5GcCE1gqNjDZOBBUM1iSYYEf2JEGZE0SoIokWMCCYk + 0qZaJKE5Qvae+33XkUhwFi4iWhSKGEELc/vsnhG5IxekVSp8lk/SkPQLyQZ3cJbumwNSJUU+zLW019T5 + Io/xV3pvvyc+uln9RSX6a++aCb+lKZ2yA3bQ/IgH4UP9g2rzMzgUPIhT1O4yOqlP9Lqe1169qFll1KMZ + GYziHUqYx1M8x0tM4y1m0Ojm9baKbuPDrp2zCVtjjsZPh7Nar60svEANlDc90brmuItreJX20zVzJRqV + YUlJTIb5DXXZmOzwHN9kT5FdD/o4zXu4SF9si8mv1FLFFuxnkwselZPh+ImSPxKlQ+8KblMXltv863DR + eW3SRXQmkk0t/lIYc0n1aDdTQUe8HOVdVis4Q0LP7Atz9diIK2iG23iN0lQ2Y8vH3Y1vneG29us/HHQD + +htLeIVPuIdTyOEyHqMPKdza/ebRg25MYJ/+BxBNvrM=''') with self.subTest('chemical-potential'): self.assertAlmostEqual64(args['η'], ''' - eNoN0F8onWEcB/AjQ1pNoZb8u8Ki0HvO+zy/F4uJ2eTG7I9W1q602toOdS4pIXZSs+zUpJN2ywXJboTE - /H7f53lf08nalZqkJldGmXYzbj+Xn2X6StM0RKf6v0rtvKS3dI8yqFOd2W5k0i3K0MfhLPNeiuQ+NVC6 - O4NK6ec33ESKYu4/PJdJjnGzJm0jIRvIT45zcnc19SVo9PPMX3FkClOYxBg+YxuMqmt5KsOyJ9UoQ1Jm - RORCbuI2DuQAh/iNfbxCCPOSr+v0kKp1e/yQKUSH99g7olZKKt8HHE956d4G/dLFasXGKUHvqIXiulnN - +7m6TEf1pm7TnaokcMMPw2nK1eeqVtUEp86WMxuei9REup1Rm/P9xk6j32XJ1tmPZs2kTML0mj6sSYVs - mxPzw6pgwNQjKgmzZKr8wVSBncaieOq1unSb3PLgkSnFChlqp096ws3y15GNdowgiif4w3fkhYzLkjzg - AV7kM74rz/gDf+M8iV0fLvAxt8mWXAGfucYq''') + eNoNzzFLW1EUAOAtDiLFB0ohggUVTG2i9HnvPS9DpkjAIg7WBG0GWxwEoRUqtkJQEETR4CKCSDsIYjEY + nIriYojlnHPvTSqablVBOgSpCCFQF8HK9wu+AziErzAHt+pOLhTH4TP0wIOKyr8myU+gEWpVxXX0KrVS + FMJQL76xoBR+xJcQgDHh06O0jtPoqoDa764zv+gCV3DgbDdULgzZRv2PumiZ07zE87zOyDlup04apHkq + UQc38xfaoB9UpRp2+JzO+fLRKY/wPWWoRcVUWoKI23v2c7+X8K6hD7blsUUWHng+D6GsgjJvVh8HkxCD + BdUhP9iAiqgZlVdCPZPb9rd75zbJoCrJG7FhX7iO+9Pd7a641a7XZrG4UwjZpAmbkJnVZ7qs1/SATnKG + muiP9pmsGbVh7ed39F2XdIPdKp7oT7xJw3JROvKNeF6I6aecgxOIw47aFwGb4xqO8Ht+y6+4im0Upzna + o15MYRYrGKEEpvEIHZqiCczgFUYpT/8Bk47KLA==''') if __name__ == '__main__': From ebd55be255c91d233dba1b60a28231ab703566b5 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Tue, 15 Oct 2024 15:24:31 +0200 Subject: [PATCH 10/10] Modify cylinderflow example --- examples/cylinderflow.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/cylinderflow.py b/examples/cylinderflow.py index 856b6ebc6..3a1f7ace1 100644 --- a/examples/cylinderflow.py +++ b/examples/cylinderflow.py @@ -119,8 +119,7 @@ def main(nelems: int = 99, domain.basis('spline', degree=(degree, degree-1), removedofs=((0,), None)), domain.basis('spline', degree=(degree-1, degree))]) @ J.T / detJ) ns.add_field(('p', 'q'), domain.basis('spline', degree=degree-1) / detJ) - ns.add_field(('t', 't0')) - ns.DuDt_i = '(u_i - u0_i) / (t - t0) + ∇_j(u_i) u_j' # material derivative + ns.add_field(('t', 'dt')) ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij' ns.ω = 'ε_ij ∇_j(u_i)' ns.N = 10 * degree / elemangle # Nitsche constant based on element size = elemangle/2 @@ -134,20 +133,21 @@ def main(nelems: int = 99, sqr = domain.integral('(.5 Σ_i (u_i - uinf_i)^2 - ∇_k(u_k) p) dV' @ ns, degree=degree*2) args = System(sqr, trial='u,p').solve(constrain=cons) # set initial condition to potential flow - res = domain.integral('(v_i DuDt_i + ∇_j(v_i) σ_ij + q ∇_k(u_k)) dV' @ ns, degree=degree*3) - res += domain.boundary['inner'].integral('(nitsche_i (u_i - uwall_i) - v_i σ_ij n_j) dS' @ ns, degree=degree*2) - div = numpy.sqrt(domain.integral('∇_k(u_k)^2 dV' @ ns, degree=2)) # L2 norm of velocity divergence + res = domain.integral('v_i (u_i - u0_i) dV' @ ns, degree=degree*3) + res += domain.integral('(v_i ∇_j(u_i) u_j + ∇_j(v_i) σ_ij + q ∇_k(u_k)) dt dV' @ ns, degree=degree*3) + res += domain.boundary['inner'].integral('(nitsche_i (u_i - uwall_i) - v_i σ_ij n_j) dt dS' @ ns, degree=degree*2) + div = numpy.sqrt(abs(function.factor(domain.integral('∇_k(u_k)^2 dV' @ ns, degree=2)))) # L2 norm of velocity divergence postprocess = PostProcessor(domain, ns) steps = treelog.iter.fraction('timestep', range(round(endtime / timestep))) if endtime < float('inf') \ else treelog.iter.plain('timestep', itertools.count()) - system = System(res, trial='u,p', test='v,q') + system = System(function.factor(res), trial='u,p', test='v,q') for _ in steps: - treelog.info(f'velocity divergence: {div.eval(**args):.0e}') - args = system.step(timestep=timestep, timearg='t', suffix='0', arguments=args, constrain=cons, tol=1e-10) + treelog.info(f'velocity divergence: {function.eval(div, **args):.0e}') + args = system.step(timestep=timestep, timearg='t', timesteparg='dt', suffix='0', arguments=args, constrain=cons, tol=1e-10) postprocess(args) return args, numpy.sqrt(domain.integral('∇_k(u_k)^2 dV' @ ns, degree=2)) @@ -201,7 +201,7 @@ def __call__(self, args): export.plotlines_(ax, x.T, self.bezier.hull, colors='k', linewidths=.1, alpha=.5) ax.quiver(*self.xgrd.T, *ugrd.T, angles='xy', width=1e-3, headwidth=3e3, headlength=5e3, headaxislength=2e3, zorder=9, alpha=.5, pivot='tip') ax.plot(0, 0, 'k', marker=(3, 2, -args['t']*self.ns.rotation.eval()*180/numpy.pi-90), markersize=20) - self.xgrd += ugrd * (args['t'] - args['t0']) + self.xgrd += ugrd * args['dt'] self.regularize_xgrd()