PETSc error code 63: Argument out of range / alternatives for project()?

Hello,

im trying to calculate dynamic hyperelasticity using the alpha-method and am running into PETSc error 63 when using my avg()-function.

In my extendedd code this function looks like this (in the minimal example attached it is simplified):

def avg(x_old, x_new, alpha):
    return alpha*x_old + (1-alpha)*x_new

x_old and x_new are Functions of the VectorFunctionSpace V, alpha is a scalar
The shown avg()-function returns a Sum-object and therefore i used project() to return a function again, but this leads to the PETSc error mentioned.

def avg(x_old, x_new, alpha):
    return project(alpha*x_old + (1-alpha)*x_new, V)

Using the docs i found that there is an additional argument for project() called ‘function’ and when i pass ‘function=u’ the error does not appear, but the results are not what would be expected. Also i found no further information what this argument means.

Would you guys have an idea of how to resolve this error?
Is there an alternative for project()-function?

from fenics import *
from mshr import *
import numpy as np

E = Constant(330e6)                             #Young's Modulus in Pa
nu = Constant(0.25)                             #Poisson's ratio
mu = Constant(E/(2.0*(1.0+nu)))
lmbda = Constant(E*nu/((1.0+nu)*(1.0-2.0*nu)))

mesh = UnitCubeMesh(20, 20, 20)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 0.0) and on_boundary

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 1.0) and on_boundary
    
V = VectorFunctionSpace(mesh, "CG", 1)
du = TrialFunction(V) 
u_ = TestFunction(V)
u = Function(V, name="Displacement")

boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_subdomains.set_all(0) #mark all with '0'
top = Top()
top.mark(boundary_subdomains, 2)
bottom = Bottom()
bottom.mark(boundary_subdomains, 3)
dss = ds(subdomain_data=boundary_subdomains)
zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, top)
bcs = [bc]

def GreenLagrangeStrain(u):
    d = len(u)
    I = Identity(d)
    F = I+grad(u)
    C = F.T * F
    E = 0.5*(C - I)
    return E
    
def PK2_GLS(u):
    E = variable(GreenLagrangeStrain(u))
    psi = lmbda/2*(tr(E)**2) + mu*tr(E*E) #St. Venant Kirchhoff material
    S = variable(diff(psi, E))
    return S

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
p = Constant((0, 0, 200000))

def k(u, u_):
    dE = derivative(GreenLagrangeStrain(u), u, u_)
    return inner(PK2_GLS(u), dE)*dx 

def Wext(u_):
    return dot(u_, p)*dss(3)

def avg(x_old, alpha):
    #return alpha*x_old +(1-alpha)*x_old
    return project(x_old*alpha+(1-alpha)*x_old, V)

alp = Constant(0.5)

#res = k(u, u_) - Wext(u_)
res = k(avg(u, alp), u_) - Wext(u_)

J = derivative(res, u)

problem = NonlinearVariationalProblem(res, u, bcs=bcs, J=J)
NLsolver = NonlinearVariationalSolver(problem)
NLsolver.parameters['newton_solver']['relative_tolerance'] = 1e-6
NLsolver.parameters['newton_solver']['linear_solver'] = 'mumps'
NLsolver.solve()

Output:

                        Solving nonlinear variational problem.
                          Newton iteration 0: r (abs) = 9.750e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_743/2749866848.py in <module>
     73 NLsolver.parameters['newton_solver']['relative_tolerance'] = 1e-6
     74 NLsolver.parameters['newton_solver']['linear_solver'] = 'mumps'
---> 75 NLsolver.solve()
     76 

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason:  PETSc error code is: 63 (Argument out of range).
*** Where:   This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1633628080759/work/dolfin/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  7e99df3564aea9dc05729ddb64498ae47b6bc15d
*** -------------------------------------------------------------------------



Thanks a lot!
Josef