Dear all, I would like some advice to decrease the running time in my code 3D_NL.py . Details are given in the attached pdf and below is a quick summary: In the attached files, I run a 3D nonlinear code for potential flow water wave equations. The vertical and horizontal space discretisation are separated, so that Firedrake solves the equations in the horizontal plane only, and the vertical dependency is involved through matrices computed by evaluating the z-dependent integrals with the functions defined in the file 'discr_sigma.py'. Basically, the potential flow phi(x,y,z,t) is defined as phi(x,y,t)phi(z), where phi(z) is discretised as a n^th order Lagrange expansion in one element: phi(z_i) = product_{k=0}^n (z-z_k)/(z_i-z_k) for k!= i There are therefore 3 unknowns, which are : - h (the total depth of water) defined in the horizontal plane - phi_s (the potential flow at the surface) defined in the horizontal plane - hat_psi a vector function containing all the phi_i for i = 1:n , which are the values of phi in the vertical nodes (excluding the surface ones), where n is the order of the Lagrange expansion. Each phi_i is defined in the horizontal plane. The time scheme is defined like this : step 1 : update phi_s and hat_psi simultaneously (that is, sum of n+1 weak formulations) step 2 : update h and hat_psi simultaneously (that is, sum of n+1 weak formulations) step 3 : update phi_s (one weak formulation). (those steps are written in the pdf file as well) At the moment, the equations are solved in a small domain with resolution 0.1, but it's running very slowly. Running in parallel does not improve it since the number of nodes is too small. Is there another way to increase the simulation speed? Maybe by changing the iteration criteria in the solver ? I hope this is clear enough, otherwise I can send you more details. Best regards, Floriane