On 4 Apr 2016, at 14:59, Fryderyk Wilczynski <scfw@leeds.ac.uk> wrote:
Thank you for spotting my mistakes!
I have implemented the suggested corrections and added the no-slip boundary conditions.
...
New code:
rho_hat = Function(P1)
Here you make a function rho_hat which is initialised to zero.
a2 = psi*mu*dx - (beta/Cn)*phi*psi*dx - Cn*inner(grad(phi),grad(psi))*dx - (1/rho_hat)*drho_hat*psi*p*dx
This term has a 1/rho_hat -> 1/0.0 -> NAN
# compute density using algebraic expression rho_hat = 0.5*( (1+phi) + (rho2/rho1)*(1-phi)) rho_hat = 0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0))
... This just assigns the name rho_hat the expression on the RHS, it does not magically go back and substitute it in to the expression for a2 above. My suggestion was to define the expression before using it when defining the variational forms: rho_hat = 0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0)) a1 = ... a2 = ... a3 = ... Then no need to do anything in the while loop.
# Solve the system of equations solve(a == L, u, bcs=[bc0, bc1, bc2, bc3]) # # update rho rho -= dt*div(rho*v0)
This line doesn't do what you think it does. You cannot compute div(rho*v0) point wise. This is your update for equation (10), right? rho^{n+1} = rho^n - dt div (rho^n v^n) You will need to solve this equation weakly in the normal way. Else wise (although I am not sure if this will destroy the conservation properties of the scheme) you could interpolate the expression: tmp_rho = Function(P1) while ...: ... tmp_rho.interpolate(rho - dt*div(rho*v0)) rho.assign(tmp_rho) Thanks, Lawrence