Projections in Viscoelasticity extremely slow (Quadrature/DG elements)

I am working with a viscoelastic model which requires updating the internal variable (C^v) at each time step inside the domain. Moreover the viscosity is a complex nonlinear function of the applied deformation C = F^T F where F is the deformation gradient tensor. Keeping it short, I have something like this

# Defining the weak form a_uv and the problem
Jac = derivative(a_uv, u, utrial)
F = NonlinearVariationalProblem(a_uv, u, bcs, J=Jac)
solver = NonlinearVariationalSolver(F)

# Loop over time steps: 
for i, t in enumerate(timeArray):
    # Assign the loading
    # Update the internal variable
     Cv = Cv + some_nonlinear_function(u, Cv)

    # Using a local projection after the update
    Cv = local_project(Cv, V)    # this process is extremely slow

I am using @bleyerj 's idea of LocalSolver as mentioned here and @kamensky 's suggestions here . I have tried using Quadrature as well as DG elements for my internal variable but there is no significant performance gain (even after using LocalSolver). Even a simple UnitCubeMesh(2,2,2) takes > 30 minutes in a serial computation.


  • If I use constant viscosity (a simpler update for the internal variable) then the code runs pretty fast (as expected).
  • I cannot use a monolithic scheme as suggested here as I need a different update rule than implicit euler.

Any pointers as to how the projection could be made faster/more efficient?

Thanks in advance.

1 Like

Hello Bhavesh,
I am surprised by the projection being extremely slow, I used this approach for 3D elastoplastic problems in large strains and it was OK.
What is exactly your some_nonlinear_function(u, Cv) is it a UFL non-linear expression ? In that case did you enforce the number of quadrature points used in the form computation ? maybe your expression is highly non-linear and FEniCS automatically chooses a large number of quadrature points.

Hi Jeremy,
Thanks for your reply. Indeed it is a complex function

C^v_{n+1} = C^v_n + \Delta t /90 * \sum\limits_{i=1, i \neq 2}^6 \omega_i k_i(u_{n+1}, u_n, C^v_n) as in this model, pages 108-109. Basically it is using an “explicit” higher order time integration. I must add that the I had already specified

parameters['form_compiler']['quadrature_degree'] = 3

and used the same degree for Quadrature elements, as well as switched UFL representation to quadrature (tried tsfc too but not much changed) inside the local_project function and back to uflacs outside it, like the suggestion posted here

Individually local_project-ing the sum \sum\limits_{i=1, i\neq 2}^6 \omega_ik_i(u_{n+1}, u_n, C^v_n) before doing the update C^v_{n+1} = \dots sped up things. This was the step I was missing perhaps. Thanks!