Matrix-free iterative solver with a linear operator (matrix-vector product)

[adtzlr]

like the already implemented solid.assemble.matrix(field) but do not assemble the matrix as a whole and multiply it by the (system-) vector. Instead, try to assemble the matrix-vector-product (a linear operator) directly.

The vector could be treated as a FieldContainer (like the real FieldContainer) and is then evaluated at the quadrature points (gradient or interpolation). Then this evaluated quantity is used inside the constitutive material formulation. The assembly is only a vector and not a matrix.

Benefits

✅ That should result in great performance 🚀.

Drawbacks

❌ Boundary conditions with prescribed displacements require some manual work.

Example

import numpy as np
import tensortrax as tr
import tensortrax.math as tm
from scipy.sparse.linalg import LinearOperator
from scipy.sparse.linalg import minres as spsolve

import felupe as fem

mesh = fem.Cube(n=16)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries = fem.dof.symmetry(field[0])
dof0, dof1 = fem.dof.partition(field, boundaries)


def W(F, μ=1, λ=2):
    C = F.T @ F
    lnJ = tm.log(tm.linalg.det(F))
    return μ * (tm.trace(C) / 2 - lnJ) + λ * lnJ**2 / 2


P = tr.gradient(W, ntrax=2)
dP = tr.hessian_vector_product(W, ntrax=2)

dx = np.zeros(field[0].values.size)
values = np.zeros(mesh.ndof)


def hvp(v):
    "Hessian-vector-product dP(F, dF) = d²W(F)/(dFdF) : dF."
    global F

    values[dof1] = v
    δx = fem.Field(region, values=values.reshape(-1, 3)).grad()
    dr = fem.IntegralForm([dP(F, δx=δx)], field, region.dV).assemble()
    return dr.toarray()[dof1]


for factor in np.arange(1, 11, 3):
    print(f"Factor={factor}")
    for iteration in range(8):
        F = field.extract()[0]

        r = fem.IntegralForm([P(F)], field, region.dV).assemble().toarray()[dof1]
        f = np.ones_like(r) * factor / (2 * mesh.ndof)

        b = -r + f
        A = LinearOperator((dof1.size, dof1.size), hvp)

        dx[dof1], residuals = spsolve(A, b, tol=1e-3, x0=dx[dof1])
        field += dx.reshape(-1, 3)

        norm = np.linalg.norm(b)
        print(f"    Iteration {1 + iteration}, Norm {norm:.3e}")
        if norm < 1e-2:
            break

    field.imshow()

Factor=1
    Iteration 1, Norm 4.367e-03
Factor=4
    Iteration 1, Norm 1.389e-02
    Iteration 2, Norm 5.295e-03
Factor=7
    Iteration 1, Norm 1.105e-02
    Iteration 2, Norm 3.894e-03
Factor=10
    Iteration 1, Norm 1.142e-02
    Iteration 2, Norm 4.751e-03