On 3 Apr 2016, at 12:06, Fryderyk Wilczynski <scfw@leeds.ac.uk> wrote:
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> 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 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.
This has occurred because you've been bitten by python division of integers. You have: Re = 100 # an integer and then write: ... (1/Re)*...*dx But 1/100 = 0 (not 0.01!) And so you end up with 0*dx and UFL complains, because it doesn't know what mesh to integrate over (because the 0 provides no information). There are two ways you can fix this: 1. Write Re = 100.0 2. At the *very top* of your python file include the line: from __future__ import division Now the normal "/" division acts "mathematically": In [1]: 1/100 Out[1]: 0 In [2]: from __future__ import division In [3]: 1/100 Out[3]: 0.01 If, for some reason, you need integer division, you can get that using "//": In [4]: 1//100 Out[4]: 0 There is also a similar problem when you use g_switch to select the gravitational term. Since again you end up with 0*...*dx I would redo this to be: gravity = False ... L3 = rho*(1/dt)*inner(chi, v0)*dx if gravity: L3 += -(1/Fr2)*rho_hat*inner(chi, ey)*dx With this done, there are a few other problems. You can't output the mixed function u0 directly, you must "split" it: v, p, phi, mu = u0.split() and output the components With modern firedrake, all of these can be output to the same file: outfile = File("./Results_NSCH/nsch_res.pvd") t = 0 tend = 1.0 v, p, phi, mu = u0.split() v.rename("Velocity") p.rename("Pressure") phi.rename("Order parameter") mu.rename("Chemical potential") outfile.write(v, p, phi, mu, time=t) while t < tend: ... outfile.write(v, p, phi, mu, time=t) Note how I have "renamed" the split functions to give them descriptive names that will appear in the output file. Also, your use of Expressions does not do what you want it to. You don't have access to python variables inside the C code strings. If you really want to do this, please instead use interpolation of UFL expressions as documented here: http://firedrakeproject.org/interpolation.html So to compute rho I would do: rho.interpolate(0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0))) For rho_hat, I note that it is always just a function of phi0, so I would just substitute in symbolically: rho_hat = 0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0)) And then use rho_hat in all your forms. Finally, your definition of "a" does not produce a bilinear form, since it contains test functions, but no trial functions. If I read the formulation correctly, I think you can compute the bilinear form by writing: a = derivative(a1 + a2 + a3 + a4, u) Additionally, I don't see where you're applying the strong (dirichlet) boundary condition on v. The natural conditions on the Cahn-Hilliard equations result in the surface integrals vanishing (as you have) but I think you need to apply v=0 strongly. Cheers, Lawrence