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