Introducing solver.System (#884)

.github/workflows/release.yml

@@ -11,7 +11,7 @@ jobs:
       - name: Checkout
         uses: actions/checkout@v4
       - name: Install dependencies
-        run: python3 -m pip install flit
+        run: python3 -m pip --break-system-packages install flit
       - name: Build package
         run: python3 -m flit build
       - name: Publish package to PyPI

.github/workflows/test.yaml

@@ -20,7 +20,7 @@ jobs:
       - name: Checkout
         uses: actions/checkout@v4
       - name: Install build dependencies
-        run: python3 -m pip install flit
+        run: python3 -m pip install --break-system-packages flit
       - name: Build package
         id: build
         run: |
@@ -183,8 +183,8 @@ jobs:
         uses: actions/checkout@v4
       - name: Install dependencies
         run: |
-          python3 -um pip install setuptools wheel
-          python3 -um pip install --upgrade --upgrade-strategy eager .[docs]
+          python3 -um pip install --break-system-packages setuptools wheel
+          python3 -um pip install --break-system-packages --upgrade --upgrade-strategy eager .[docs]
       - name: Build docs
         run: python3 -m sphinx -n -W --keep-going docs build/sphinx/html
   build-and-test-container-image:

examples/adaptivity.py

@@ -1,4 +1,5 @@
-from nutils import mesh, function, solver, export, testing
+from nutils import mesh, function, export, testing
+from nutils.solver import System
 from nutils.expression_v2 import Namespace
 import numpy
 import treelog
@@ -68,13 +69,13 @@ def main(etype: str = 'square',
         ns.du = 'u - uexact'
 
         sqr = domain.boundary['corner'].integral('u^2 dS' @ ns, degree=degree*2)
-        cons = solver.optimize('u,', sqr, droptol=1e-15)
+        cons = System(sqr, trial='u').optimize(droptol=1e-15)
 
         sqr = domain.boundary.integral('du^2 dS' @ ns, degree=7)
-        cons = solver.optimize('u,', sqr, droptol=1e-15, constrain=cons)
+        cons = System(sqr, trial='u').optimize(droptol=1e-15, constrain=cons)
 
         res = domain.integral('∇_k(v) ∇_k(u) dV' @ ns, degree=degree*2)
-        args = solver.solve_linear('u:v', res, constrain=cons)
+        args = System(res, trial='u', test='v').solve(constrain=cons)
 
         ndofs = len(args['u'])
         error = numpy.sqrt(domain.integral(['du^2 dV', '(du^2 + ∇_k(du) ∇_k(du)) dV'] @ ns, degree=7)).eval(**args)

examples/burgers.py

@@ -1,4 +1,5 @@
-from nutils import mesh, function, solver, export, testing
+from nutils import mesh, function, export, testing
+from nutils.solver import System
 from nutils.expression_v2 import Namespace
 import treelog as log
 import numpy
@@ -40,8 +41,8 @@ def main(nelems: int = 40,
     ns.x = geom
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
     ns.add_field(('u', 'u0', 'v'), domain.basis(btype, degree=degree))
-    ns.dt = timestep
-    ns.dudt = '(u - u0) / dt'
+    ns.add_field(('t', 't0'))
+    ns.dudt = '(u - u0) / (t - t0)'
     ns.f = '.5 u^2'
     ns.C = 1
     ns.uinit = 'exp(-25 x^2)'
@@ -50,18 +51,20 @@ def main(nelems: int = 40,
     res -= domain.interfaces.integral('[v] n ({f} - .5 C [u] n) dS' @ ns, degree=degree*2)
 
     sqr = domain.integral('(u - uinit)^2 dV' @ ns, degree=max(degree*2, 5))
-    args = solver.optimize('u,', sqr)
+    args = System(sqr, trial='u').solve()
+    args['t'] = 0.
+
+    system = System(res, trial='u', test='v')
 
     bezier = domain.sample('bezier', 7)
-    with log.iter.plain('timestep', itertools.count(step=timestep)) as times:
-        for t in times:
-            log.info('time:', round(t, 10))
+    with log.iter.plain('timestep', itertools.count()) as steps:
+        for _ in steps:
+            log.info('time:', round(args['t'], 10))
             x, u = bezier.eval(['x', 'u'] @ ns, **args)
             export.triplot('solution.png', x[:,numpy.newaxis], u, tri=bezier.tri, hull=bezier.hull, clim=(0, 1))
-            if t >= endtime:
+            if args['t'] >= endtime:
                 break
-            args['u0'] = args['u']
-            args = solver.newton('u:v', res, arguments=args).solve(newtontol)
+            args = system.step(timestep=timestep, arguments=args, timetarget='t', historysuffix='0', tol=newtontol)
 
     return args
 

examples/cahnhilliard.py

@@ -1,4 +1,5 @@
-from nutils import mesh, function, solver, numeric, export, testing
+from nutils import mesh, function, numeric, export, testing
+from nutils.solver import System
 from nutils.expression_v2 import Namespace
 from nutils.SI import Length, Time, Density, Tension, Energy, Pressure, Velocity, parse
 import numpy
@@ -178,6 +179,9 @@ def main(size: Length = parse('10cm'),
     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.ε = epsilon
     ns.σ = stens
     ns.σmean = (wtensp + wtensn) / 2
@@ -188,7 +192,7 @@ def main(size: Length = parse('10cm'),
     ns.δψ = stab.value
     ns.M = mobility
     ns.J_i = '-M ∇_i(η)'
-    ns.dt = timestep
+    ns.dt = 't - t0'
 
     nrg_mix = domain.integral('(ψ σ / ε) dV' @ ns, degree=degree*4)
     nrg_iface = domain.integral('.5 σ ε ∇_k(φ) ∇_k(φ) dV' @ ns, degree=degree*4)
@@ -198,14 +202,15 @@ def main(size: Length = parse('10cm'),
     numpy.random.seed(seed)
     args = dict(φ=numpy.random.normal(0, .5, basis.shape)) # initial condition
 
+    system = System(nrg / tol, trial='φ,η')
+
     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['φ0'] = args['φ']
-            args = solver.optimize(['φ', 'η'], nrg / tol, arguments=args, tol=1)
+            args = system.step(timestep=1., timetarget='t', historysuffix='0', arguments=args, tol=1, maxiter=5)
 
             with export.mplfigure('phase.png') as fig:
                 ax = fig.add_subplot(aspect='equal', xlabel='[mm]', ylabel='[mm]')

examples/coil.py

@@ -1,4 +1,5 @@
-from nutils import export, function, mesh, solver, testing
+from nutils import export, function, mesh, testing
+from nutils.solver import System
 from nutils.expression_v2 import Namespace
 import functools
 import numpy
@@ -115,7 +116,7 @@ def main(nelems: int = 50,
     res = REV.integral(RZ.integral('-∇_j(Atest_i) ∇_j(A_i) dV' @ ns, degree=2*degree), degree=0)
     res += REV.integral(RZ['coil'].integral('μ0 Atest_i J_i dV' @ ns, degree=2*degree), degree=0)
 
-    args = solver.solve_linear('A:Atest,', res)
+    args = System(res, trial='A', test='Atest').solve()
 
     # Since the coordinate transformation is singular at r=0 we can't evaluate
     # `B` (the curl of `A`) at r=0. We circumvent this problem by projecting `B`
@@ -124,7 +125,7 @@ def main(nelems: int = 50,
     ns.Borig = ns.B
     ns.add_field(('B', 'Btest'), RZ.basis('spline', degree=degree), ns.rot, dtype=complex)
     res = REV.integral(RZ.integral('Btest_i (B_i - Borig_i) dV' @ ns, degree=2*degree), degree=0)
-    args = solver.solve_linear('B:Btest,', res, arguments=args)
+    args = System(res, trial='B', test='Btest').solve(arguments=args)
 
     with export.mplfigure('magnetic-potential-1.png', dpi=300) as fig:
         ax = fig.add_subplot(111, aspect='equal', xlabel='$x_0$', ylabel='$x_2$', adjustable='datalim')

examples/cylinderflow.py

@@ -1,4 +1,5 @@
-from nutils import mesh, function, solver, export, testing, numeric
+from nutils import mesh, function, export, testing, numeric
+from nutils.solver import System
 from nutils.expression_v2 import Namespace
 import itertools
 import numpy
@@ -118,8 +119,8 @@ 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.dt = timestep
-    ns.DuDt_i = '(u_i - u0_i) / dt + ∇_j(u_i) u_j' # material derivative
+    ns.add_field(('t', 't0'))
+    ns.DuDt_i = '(u_i - u0_i) / (t - t0) + ∇_j(u_i) u_j' # material derivative
     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
@@ -128,10 +129,10 @@ def main(nelems: int = 99,
     ns.uwall_i = 'rotation ε_ij x_j' # clockwise positive rotation
 
     sqr = domain.boundary['inflow'].integral('Σ_i (u_i - uinf_i)^2 dS' @ ns, degree=degree*2)
-    cons = solver.optimize('u,', sqr, droptol=1e-15) # constrain inflow boundary to unit horizontal flow
+    cons = System(sqr, trial='u').optimize(droptol=1e-15) # constrain inflow boundary to unit horizontal flow
 
     sqr = domain.integral('(.5 Σ_i (u_i - uinf_i)^2 - ∇_k(u_k) p) dV' @ ns, degree=degree*2)
-    args = solver.optimize('u,p', sqr, constrain=cons) # set initial condition to potential flow
+    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)
@@ -142,10 +143,11 @@ def main(nelems: int = 99,
     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')
+
     for _ in steps:
         treelog.info(f'velocity divergence: {div.eval(**args):.0e}')
-        args['u0'] = args['u']
-        args = solver.newton('u:v,p:q', residual=res, arguments=args, constrain=cons).solve(1e-10)
+        args = system.step(timestep=timestep, timetarget='t', historysuffix='0', arguments=args, constrain=cons, tol=1e-10)
         postprocess(args)
 
     return args, numpy.sqrt(domain.integral('∇_k(u_k)^2 dV' @ ns, degree=2))
@@ -161,7 +163,6 @@ def __init__(self, topo, ns, region=4., aspect=16/9, figscale=7.2, spacing=.05,
         self.spacing = spacing
         self.pstep = pstep
         self.vortlim = vortlim
-        self.t = 0.
         self.initialize_xgrd()
         self.topo = topo
 
@@ -199,10 +200,8 @@ def __call__(self, args):
             ax.tricontour(*x.T, self.bezier.tri, p, numpy.arange(numpy.min(p)//1, numpy.max(p), self.pstep), colors='k', linestyles='solid', linewidths=.5, zorder=9)
             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, -self.t*self.ns.rotation.eval()*180/numpy.pi-90), markersize=20)
-        dt = self.ns.dt.eval()
-        self.t += dt
-        self.xgrd += ugrd * dt
+            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.regularize_xgrd()
 
 

examples/drivencavity.py

@@ -1,4 +1,5 @@
-from nutils import mesh, function, solver, export, testing
+from nutils import mesh, function, export, testing
+from nutils.solver import System, LinesearchNewton
 from nutils.expression_v2 import Namespace
 import treelog as log
 import numpy
@@ -117,12 +118,12 @@ def main(nelems: int = 32,
 
     # strong enforcement of non-penetrating boundary conditions
     sqr = domain.boundary.integral('(u_k n_k)^2 dS' @ ns, degree=degree*2)
-    cons = solver.optimize('u,', sqr, droptol=1e-15)
+    cons = System(sqr, trial='u').optimize(droptol=1e-15)
 
     if strongbc:
         # strong enforcement of tangential boundary conditions
         sqr = domain.boundary.integral('(ε_ij n_i (u_j - uwall_j))^2 dS' @ ns, degree=degree*2)
-        tcons = solver.optimize('u,', sqr, droptol=1e-15)
+        tcons = System(sqr, trial='u').optimize(droptol=1e-15)
         cons['u'] = numpy.choose(numpy.isnan(cons['u']), [cons['u'], tcons['u']])
     else:
         # weak enforcement of tangential boundary conditions via Nitsche's method
@@ -131,7 +132,7 @@ def main(nelems: int = 32,
         res += domain.boundary.integral('(nitsche_i (u_i - uwall_i) - v_i σ_ij n_j) dS' @ ns, degree=2*degree)
 
     with log.context('stokes'):
-        args0 = solver.solve_linear('u:v,p:q', res, constrain=cons)
+        args0 = System(res, trial='u,p', test='v,q').solve(constrain=cons)
         postprocess(domain, ns, **args0)
 
     # change to Navier-Stokes by adding convection
@@ -141,7 +142,7 @@ def main(nelems: int = 32,
         res += domain.integral('.5 u_i v_i ∇_j(u_j) dV' @ ns, degree=degree*3)
 
     with log.context('navier-stokes'):
-        args1 = solver.newton('u:v,p:q', res, arguments=args0, constrain=cons).solve(tol=1e-10)
+        args1 = System(res, trial='u,p', test='v,q').solve(arguments=args0, constrain=cons, tol=1e-10, method=LinesearchNewton())
         postprocess(domain, ns, **args1)
 
     return args0, args1
@@ -155,7 +156,7 @@ def postprocess(domain, ns, **arguments):
 
     # reconstruct velocity streamlines
     sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(ψ))^2 dV' @ ns, degree=4)
-    arguments = solver.optimize('ψ,', sqr, arguments=arguments)
+    arguments = System(sqr, trial='ψ').solve(arguments=arguments)
 
     bezier = domain.sample('bezier', 9)
     x, u, p, ψ = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'ψ'] @ ns, **arguments)

examples/elasticity.py

@@ -1,4 +1,5 @@
-from nutils import mesh, function, solver, export, testing
+from nutils import mesh, function, export, testing
+from nutils.solver import System
 from nutils.expression_v2 import Namespace
 import treelog as log
 import numpy
@@ -50,20 +51,18 @@ def main(nelems: int = 24,
     ns.q_i = '-δ_i1'
 
     sqr = domain.boundary['top'].integral('u_k u_k dS' @ ns, degree=degree*2)
-    cons = solver.optimize(('u',), sqr, droptol=1e-15)
+    cons = System(sqr, trial='u').optimize(droptol=1e-15)
 
     # solve for equilibrium configuration
-    internal = domain.integral('E dV' @ ns, degree=degree*2)
-    external = domain.integral('u_i q_i dV' @ ns, degree=degree*2)
-    args = solver.optimize('u,', internal - external, constrain=cons)
+    energy = domain.integral('(E - u_i q_i) dV' @ ns, degree=degree*2)
+    args = System(energy, trial='u').solve(constrain=cons)
 
     # evaluate tractions and net force
     if direct:
         ns.t_i = 'σ_ij n_j' # <-- this is an inadmissible boundary term
     else:
-        external += domain.boundary['top'].integral('u_i t_i dS' @ ns, degree=degree*2)
-        invcons = dict(t=numpy.choose(numpy.isnan(cons['u']), [numpy.nan, 0.]))
-        args = solver.solve_linear(('t',), [(internal - external).derivative('u')], constrain=invcons, arguments=args)
+        energy -= domain.boundary['top'].integral('u_i t_i dS' @ ns, degree=degree*2)
+        args = System(energy, trial='t', test='u').solve(constrain={'t': numpy.isnan(cons['u'])}, arguments=args)
     F = domain.boundary['top'].integrate('t_i dS' @ ns, degree=degree*2, arguments=args)
     log.user('total clamping force:', F)
 
@@ -132,8 +131,8 @@ def test_poisson(self):
                 eNpjaGBAhSBAZTEAEKAUAQ==''')
         with self.subTest('solution'):
             self.assertAlmostEqual64(args['u'], '''
-                eNqTNig6vcVwwekjRuJn5Iy1zzIAwQs999MdBmWn+w0Zz7QYpoPF/v+3Pv3/f/5pBgbWMwwM8WAxiYvu
-                pyvOl50uPMd4puYcRN3T80Wnfc4tOG1zVvzMozMQ8wAOOSng''')
+                eNqTNig6vcVwwekjRuJn5Iy1zzIAwQs999MdBmWn+w0Zz7QYpoPFGBisTzMw5AMx6xkGhniwmMRF99MV
+                58tOF55jPFNzDqLu6fmi0z7nFpy2OSt+5tEZiHkAKRAl5A==''')
         with self.subTest('traction'):
             self.assertAlmostEqual64(args['t'], '''
                 eNpjYEAF/Sc+maMJMdw0emzGgAFiMdSpn8VUV2j+yRwAoCAJFw==''')

examples/finitestrain.py

@@ -1,4 +1,5 @@
-from nutils import mesh, function, solver, export, testing
+from nutils import mesh, function, export, testing
+from nutils.solver import System, Minimize
 from nutils.expression_v2 import Namespace
 import numpy
 
@@ -57,18 +58,18 @@ def main(nelems: int = 20,
 
     sqr = domain.boundary['left'].integral('u_k u_k dS' @ ns, degree=degree*2)
     sqr += domain.boundary['right'].integral('((u_0 - X_1 sin(2 angle) - cos(angle) + 1)^2 + (u_1 - X_1 (cos(2 angle) - 1) + sin(angle))^2) dS' @ ns, degree=degree*2)
-    cons = solver.optimize('u,', sqr, droptol=1e-15)
+    cons = System(sqr, trial='u').optimize(droptol=1e-15)
 
     energy = domain.integral('energy dV' @ ns, degree=degree*2)
-    args0 = solver.optimize('u,', energy, constrain=cons)
+    args0 = System(energy, trial='u').solve(constrain=cons)
     x, energy = bezier.eval(['x_i', 'energy'] @ ns, **args0)
     export.triplot('linear.png', x, energy, tri=bezier.tri, hull=bezier.hull, cmap='jet')
 
     ns.ε_ij = '.5 (∇_j(u_i) + ∇_i(u_j) + ∇_i(u_k) ∇_j(u_k))'
     ns.energy = 'λ ε_ii ε_jj + 2 μ ε_ij ε_ij'
 
     energy = domain.integral('energy dV' @ ns, degree=degree*2)
-    args1 = solver.minimize('u,', energy, arguments=args0, constrain=cons).solve(restol)
+    args1 = System(energy, trial='u').solve(arguments=args0, constrain=cons, method=Minimize(), tol=restol)
     x, energy = bezier.eval(['x_i', 'energy'] @ ns, **args1)
     export.triplot('nonlinear.png', x, energy, tri=bezier.tri, hull=bezier.hull, cmap='jet')
 

examples/laplace.py

@@ -1,4 +1,5 @@
-from nutils import mesh, function, solver, export, testing
+from nutils import mesh, function, export, testing
+from nutils.solver import System
 from nutils.expression_v2 import Namespace
 import treelog
 
@@ -71,13 +72,13 @@ def main(nelems: int = 10,
 
     sqr = domain.boundary['left'].integral('u^2 dS' @ ns, degree=degree*2)
     sqr += domain.boundary['top'].integral('(u - cosh(1) sin(x_0))^2 dS' @ ns, degree=degree*2)
-    cons = solver.optimize('u,', sqr, droptol=1e-15)
+    cons = System(sqr, trial='u').optimize(droptol=1e-15)
 
     # The unconstrained entries of `u` are to be such that the residual
     # evaluates to zero for all possible values of `v`. The resulting array `u`
     # matches `cons` in the constrained entries.
 
-    args = solver.solve_linear('u:v', res, constrain=cons)
+    args = System(res, trial='u', test='v').solve(constrain=cons)
 
     # Once all arguments are establised, the corresponding solution can be
     # vizualised by sampling values of `ns.u` along with physical coordinates