Dear all,

I'm experiencing some problems setting up weak forms for the system of equations I'm trying to solve.

Here's what I'm trying to do:

-----

from firedrake import *

# Model parameters

dt = 5.0e-04 # time step

Cn=0.01 # Cahn

We=0.45 # Weber

Pe=4.5 # Peclet

Re=100 # Reynolds

g_switch = 0 # switch for gravity 1 = ON, 0 = Off

Fr2=0.1 # Froude number squared

beta = 2 # stabilising term

rho1 = Constant(1)

rho2 = Constant(1)

drho_hat = Constant(0)

drho_hat.assign(0.5*(1 - (rho2/rho1))) # d rho_hat / d phi

alpha = (rho2-rho1)/(rho2+rho1)

# Create mesh and define function spaces

mesh = Mesh("square_mesh.msh")

P1 = FunctionSpace(mesh, "Lagrange", 1)

P2 = VectorFunctionSpace(mesh, "Lagrange", 2)

# Taylor-Hood elements for velocity and pressure, P1-P1 elements for phi and mu

W = P2*P1*P1*P1

# Define test functions

chi, theta, omega, psi = TestFunctions(W)

# Define functions

u = Function(W) # current solution

u0 = Function(W) # solution from previous converged step

# Split mixed functions

v, p, phi, mu = split(u)

v0, p0, phi0, mu0 = split(u0)

# Compute the chemical potential df/dphi

phi = variable(phi)

f = 0.25*(phi**2 - 1)**2

dfdphi = diff(f, phi)

# create a functions for the two definitions of density

rho_hat = Function(P1)

rho = Function(P1)

ex, ey = unit_vectors(2)

# Weak statement of the equations

a1 = (1/dt)*omega*phi*dx + omega*div(v*phi0)*dx + (1/Pe)*inner(grad(omega),grad(mu))*dx

a2 = psi*mu*dx - (beta/Cn)*phi*psi*dx - Cn*inner(grad(phi),grad(psi))*dx - (1/rho_hat)*drho_hat*psi*p*dx

a3 = rho*(1/dt)*inner(chi, v)*dx + rho*inner(dot(v0,nabla_grad(v)), chi)*dx +(1/We)*(-p*div(chi) + inner(chi, phi*grad(mu))-(1/rho_hat)*drho_hat*p*inner(chi, grad(phi0)) )*dx + (1/Re)*inner((nabla_grad(v)[i,j] + nabla_grad(v)[j,i]), nabla_grad(chi)[i,j] )*dx - (2/(3*Re))*div(v)*div(chi)*dx

a4 = theta*div(v)*dx + (alpha/Pe)*inner(grad(theta), grad(mu))*dx

a = a1 + a2 + a3 + a4

L1 = (1/dt)*omega*phi0*dx

L2 = - (beta/Cn)*phi0*psi*dx + (1/Cn)*dfdphi*psi*dx

L3 = - (g_switch/Fr2)*rho_hat*inner(chi,ey)*dx + rho*(1/dt)*inner(chi, v0)*dx

L = L1 + L2 + L3

v0, p0, phi0, mu0 = u0.split()

# Initial conditions for phi and velocity

phi_ini = Expression('''1 - tanh( ( sqrt( pow(x[0]-x0, 2)+pow(x[1]-y0, 2) ) - r0 )/C*sqrt(2) )- tanh( ( sqrt( pow(x[0]-x1, 2)+pow(x[1]-y1, 2) ) - r1 )/C*sqrt(2) )''', x0=0.4,y0=0.5,r0 = 0.25,x1=0.78,y1=0.5,r1 = 0.1,C=0.01 )

phi0.interpolate(phi_ini)

vel_ini = Expression(("0.0", "0.0"))

v0.interpolate(vel_ini)

# initialise density

rho.interpolate( Expression('''0.5*((1+phi_ini) + (rho2/rho1)*(1-phi_ini))''') )

outfile = File("./Results_NSCH/nsch_res.pvd")

t = 0

tend = 1.0

outfile.write(u0)

while t <= tend:

t+=dt ; print t

# compute density using algebraic expression rho_hat = 0.5*( (1+phi) + (rho2/rho1)*(1-phi))

rho_hat.interpolate(Expression(' 0.5*( (1+phi0) + (rho2/rho1)*(1-phi0))'))

# Solve the system of equations

solve(a == L, u)

# update rho

rho -= dt*div(rho*v0)

u0.assign(u)

outfile.write(u0)

-----

When I try to run the code I get an exception:

This integral is missing an integration domain.

Traceback (most recent call last):

File "navier_stokes_CH.py", line 54, in <module>

File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/measure.py", line 398, in __rmul__

error("This integral is missing an integration domain.")

File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/log.py", line 158, in error

raise self._exception_type(self._format_raw(*message))

ufl.log.UFLException: This integral is missing an integration domain.

Any help on fixing that would be greatly appreciated

I assume that this is not the only bug in my code, so I'd also appreciate other corrections that need to be made.

I have attached a pdf with a description of exactly what I want to do.

Thanks,

Fryderyk