Integral is missing an integration domain.
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. 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
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
Thank you for spotting my mistakes! I have implemented the suggested corrections and added the no-slip boundary conditions. Now I get the following exception: [scfw@comp-pc3237 CH_advection]$ python navier_stokes_CH.py 0.0005 pyop2:INFO Solving linear variational problem... Traceback (most recent call last): File "navier_stokes_CH.py", line 103, in <module> solve(a == L, u, bcs=[bc0, bc1, bc2, bc3]) # File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve _solve_varproblem(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 147, in _solve_varproblem solver.solve() File "<decorator-gen-295>", line 2, in solve File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper return f(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", line 190, in solve solving_utils.check_snes_convergence(self.snes) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", line 62, in check_snes_convergence %s""" % (snes.getIterationNumber(), msg)) RuntimeError: Nonlinear solve failed to converge after 0 nonlinear iterations. Reason: Inner linear solve failed to converge after 0 iterations with reason: DIVERGED_NANORINF ---- New code: from __future__ import division 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 gravity = False # switch for gravity 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 = UnitSquareMesh(200, 200) 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 # Dirichlet boundary condition for velocity noslip = Constant((0, 0)) bc0 = DirichletBC(W.sub(0), noslip, 1) bc1 = DirichletBC(W.sub(0), noslip, 2) bc2 = DirichletBC(W.sub(0), noslip, 3) bc3 = DirichletBC(W.sub(0), noslip, 4) # 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 = derivative(a1 + a2 + a3 + a4, u) L1 = (1/dt)*omega*phi0*dx L2 = - (beta/Cn)*phi0*psi*dx + (1/Cn)*dfdphi*psi*dx L3 = rho*(1/dt)*inner(chi, v0)*dx if gravity: L3 += -(1/Fr2)*rho_hat*inner(chi, ey)*dx L = L1 + L2 + L3 v0, p0, phi0, mu0 = u0.split() # Initial conditions for phi 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) # Initial conditions for velocity vel_ini = Expression(("0.0", "0.0")) v0.interpolate(vel_ini) # initialise density rho.interpolate(0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0))) 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: t+=dt ; print t # 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)) # Solve the system of equations solve(a == L, u, bcs=[bc0, bc1, bc2, bc3]) # # update rho rho -= dt*div(rho*v0) u0.assign(u) outfile.write(v, p, phi, mu, time=t) ----
Some considerations for the firedrake(s): - It may be worthwhile to use a direct solver at first to see whether the scheme works at all; this may be slow but that is of later consideration. Note that the system of equations to be solved is linear (the system at time level ^n+1 is linear), even though the system of equations is nonlinear. - For a different iterative solver one may need some adequate preconditioning? Some remarks just in case this was not clear. --Onno ________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Fryderyk Wilczynski <scfw@leeds.ac.uk> Sent: Monday, April 4, 2016 2:59 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Integral is missing an integration domain. Thank you for spotting my mistakes! I have implemented the suggested corrections and added the no-slip boundary conditions. Now I get the following exception: [scfw@comp-pc3237 CH_advection]$ python navier_stokes_CH.py 0.0005 pyop2:INFO Solving linear variational problem... Traceback (most recent call last): File "navier_stokes_CH.py", line 103, in <module> solve(a == L, u, bcs=[bc0, bc1, bc2, bc3]) # File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve _solve_varproblem(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 147, in _solve_varproblem solver.solve() File "<decorator-gen-295>", line 2, in solve File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper return f(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", line 190, in solve solving_utils.check_snes_convergence(self.snes) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", line 62, in check_snes_convergence %s""" % (snes.getIterationNumber(), msg)) RuntimeError: Nonlinear solve failed to converge after 0 nonlinear iterations. Reason: Inner linear solve failed to converge after 0 iterations with reason: DIVERGED_NANORINF ---- New code: from __future__ import division 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 gravity = False # switch for gravity 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 = UnitSquareMesh(200, 200) 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 # Dirichlet boundary condition for velocity noslip = Constant((0, 0)) bc0 = DirichletBC(W.sub(0), noslip, 1) bc1 = DirichletBC(W.sub(0), noslip, 2) bc2 = DirichletBC(W.sub(0), noslip, 3) bc3 = DirichletBC(W.sub(0), noslip, 4) # 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 = derivative(a1 + a2 + a3 + a4, u) L1 = (1/dt)*omega*phi0*dx L2 = - (beta/Cn)*phi0*psi*dx + (1/Cn)*dfdphi*psi*dx L3 = rho*(1/dt)*inner(chi, v0)*dx if gravity: L3 += -(1/Fr2)*rho_hat*inner(chi, ey)*dx L = L1 + L2 + L3 v0, p0, phi0, mu0 = u0.split() # Initial conditions for phi 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) # Initial conditions for velocity vel_ini = Expression(("0.0", "0.0")) v0.interpolate(vel_ini) # initialise density rho.interpolate(0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0))) 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: t+=dt ; print t # 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)) # Solve the system of equations solve(a == L, u, bcs=[bc0, bc1, bc2, bc3]) # # update rho rho -= dt*div(rho*v0) u0.assign(u) outfile.write(v, p, phi, mu, time=t) ---- _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 4 Apr 2016, at 15:04, Onno Bokhove <O.Bokhove@leeds.ac.uk> wrote:
Some considerations for the firedrake(s): - It may be worthwhile to use a direct solver at first to see whether the scheme works at all;
Definitely. Try with not to many elements. You'll need to do the following: parameters["matnest"] = False And run with: solve(..., solver_parameters={'pc_type': 'lu'}) Be aware that you may need to select a "better" direct solver (e.g. MUMPS) using: 'pc_factor_mat_solver_package': 'mumps'
this may be slow but that is of later consideration. Note that the system of equations to be solved is linear (the system at time level ^n+1 is linear), even though the system of equations is nonlinear.
Aha, I had wondered where the nonlinearities had disappeared to!
- For a different iterative solver one may need some adequate preconditioning?
Without question. There are good schemes for the separate navier-stokes and cahn-hilliard blocks, there appears to be less literature on preconditioning for the coupled system (although a good start would be to have a block-preconditioner that ignores the couplings). Jessica Bosch has been doing work recently on preconditioning the coupled system (see e.g. http://meetings.siam.org/sess/dsp_talk.cfm?p=71937), but I do not think a paper has appeared. Cheers, Lawrence
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
Dear all, I want to define a piecewise expression, the derivative of which will be used on the RHS in my weak for and needs to be updated every timestep. So where I previously had: phi = variable(phi) f = 0.25*(phi**2 - 1)**2 dfdphi = diff(f, phi) I want to define f as follows: f = (phi +1)^2 if phi < -1 0.25* ( phi^2 - 1 )^2 if -1 < phi < 1 (phi - 1)^2 if phi > 1 I have tried the following phi0 = variable(phi0) class FreeEnergyPotential(Expression): def eval(self, values, phi0): if phi0 > 1: values[:] = (phi0-1)**2 elif phi0 < -1: values[:] = (phi0+1)**2 else: values[:] = 0.25*(phi0**2 -1)**2 f= Function(P1).interpolate(FreeEnergyPotential()) dfdphi = diff(f, phi0) But that didn't work and I get the familiar exception: This integral is missing an integration domain. Traceback (most recent call last): File "navier_stokes_CH.py", line 96, in <module> L2 = - (beta/Cn)*phi0*psi*dx + (1/Cn)*dfdphi*psi*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. How can I implement it correctly? Full code attached. Thanks, Fryderyk ________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: 04 April 2016 15:37 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Integral is missing an integration domain.
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
Dear Fryderyk, You are trying to do a symbolic operation on numeric value, and, of course, it doesn't work as you expect it. f = Function(P1).interpolate(FreeEnergyPotential()) This construct f as the numerical values in the given function space. dfdphi = diff(f, phi0) This tries to symbolically differentiate f according to phi0, but f does not depend on phi0, so you get zero, hence the error message. You should define f symbolically: f = conditional(phi0 > 1, (phi0 - 1)**2, conditional(phi0 < -1, (phi0 + 1)**2, 0.25*(phi0**2 - 1)**2)) This way the symbolic differentiation will work correctly. On 06/04/16 13:47, Fryderyk Wilczynski wrote:
Dear all,
I want to define a piecewise expression, the derivative of which will be used on the RHS in my weak for and needs to be updated every timestep. So where I previously had:
phi = variable(phi) f = 0.25*(phi**2 - 1)**2 dfdphi = diff(f, phi)
I want to define f as follows: f = (phi +1)^2 if phi < -1 0.25* ( phi^2 - 1 )^2 if -1 < phi < 1 (phi - 1)^2 if phi > 1
I have tried the following
phi0 = variable(phi0) class FreeEnergyPotential(Expression): def eval(self, values, phi0): if phi0 > 1: values[:] = (phi0-1)**2 elif phi0 < -1: values[:] = (phi0+1)**2 else: values[:] = 0.25*(phi0**2 -1)**2 f= Function(P1).interpolate(FreeEnergyPotential()) dfdphi = diff(f, phi0)
But that didn't work and I get the familiar exception: This integral is missing an integration domain. Traceback (most recent call last): File "navier_stokes_CH.py", line 96, in <module> L2 = - (beta/Cn)*phi0*psi*dx + (1/Cn)*dfdphi*psi*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.
How can I implement it correctly? Full code attached.
Thanks, Fryderyk
________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: 04 April 2016 15:37 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Integral is missing an integration domain.
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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Thank you, that conditional worked. I have realized that I should in fact have trial functions, (equations are linear in (n+1) terms), I have been mistakenly following example scripts for nonlinear problems. I have changed v, p, phi, mu = split(u) rho = Function(P1) to v, p, phi, mu = TrialFunctions(W) rho = TrialFunction(P1) and I now I do have a bilinear form so I changed a = derivative(a1 + a2 + a3 + a4, u) back to a = a1 + a2 + a3 + a4 and b1 = (1/dt)*eta*rho*dx b = derivative(b1, rho) to b = (1/dt)*eta*rho*dx Everything else I kept unchanged. Now when I try to run the code I get the following: Traceback (most recent call last): File "navier_stokes_CH2.py", line 115, in <module> solve(a == L, u, bcs=[bc0, bc1, bc2, bc3], solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'}) # File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve _solve_varproblem(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 145, in _solve_varproblem options_prefix=options_prefix) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", line 252, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", line 126, in __init__ ctx = solving_utils._SNESContext(problem) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", line 106, in __init__ for problem in problems) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", line 106, in <genexpr> for problem in problems) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", line 66, in assemble inverse=inverse, nest=nest) File "<decorator-gen-296>", line 2, in _assemble File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", line 100, in _assemble inverse=inverse) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/tsfc_interface.py", line 243, in compile_form number_map).kernels File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line 203, in __new__ obj = make_obj() File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line 193, in make_obj obj.__init__(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/tsfc_interface.py", line 165, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/tsfc/driver.py", line 42, in compile_form do_estimate_degrees=True) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/algorithms/compute_form_data.py", line 362, in compute_form_data check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", line 143, in check_form_arity check_integrand_arity(itg.integrand(), arguments) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", line 137, in check_integrand_arity args = map_expr_dag(rules, expr, compress=False) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/corealg/map_dag.py", line 84, in map_expr_dags r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands]) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", line 55, in product raise ArityMismatch("Multiplying expressions with overlapping form argument number {0}, argument is {1}.".format(x.number(), x)) ufl.algorithms.check_arities.ArityMismatch: Multiplying expressions with overlapping form argument number 1, argument is v_1. Exception AttributeError: "'LinearVariationalSolver' object has no attribute '_parameters'" in <bound method LinearVariationalSolver.__del__ of <firedrake.variational_solver.LinearVariationalSolver object at 0x457c750>> ignored ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Miklós Homolya <m.homolya14@imperial.ac.uk> Sent: 06 April 2016 14:50 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Integral is missing an integration domain. Dear Fryderyk, You are trying to do a symbolic operation on numeric value, and, of course, it doesn't work as you expect it. f = Function(P1).interpolate(FreeEnergyPotential()) This construct f as the numerical values in the given function space. dfdphi = diff(f, phi0) This tries to symbolically differentiate f according to phi0, but f does not depend on phi0, so you get zero, hence the error message. You should define f symbolically: f = conditional(phi0 > 1, (phi0 - 1)**2, conditional(phi0 < -1, (phi0 + 1)**2, 0.25*(phi0**2 - 1)**2)) This way the symbolic differentiation will work correctly. On 06/04/16 13:47, Fryderyk Wilczynski wrote: Dear all, I want to define a piecewise expression, the derivative of which will be used on the RHS in my weak for and needs to be updated every timestep. So where I previously had: phi = variable(phi) f = 0.25*(phi**2 - 1)**2 dfdphi = diff(f, phi) I want to define f as follows: f = (phi +1)^2 if phi < -1 0.25* ( phi^2 - 1 )^2 if -1 < phi < 1 (phi - 1)^2 if phi > 1 I have tried the following phi0 = variable(phi0) class FreeEnergyPotential(Expression): def eval(self, values, phi0): if phi0 > 1: values[:] = (phi0-1)**2 elif phi0 < -1: values[:] = (phi0+1)**2 else: values[:] = 0.25*(phi0**2 -1)**2 f= Function(P1).interpolate(FreeEnergyPotential()) dfdphi = diff(f, phi0) But that didn't work and I get the familiar exception: This integral is missing an integration domain. Traceback (most recent call last): File "navier_stokes_CH.py", line 96, in <module> L2 = - (beta/Cn)*phi0*psi*dx + (1/Cn)*dfdphi*psi*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. How can I implement it correctly? Full code attached. Thanks, Fryderyk ________________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> <firedrake-bounces@imperial.ac.uk><mailto:firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk><mailto:lawrence.mitchell@imperial.ac.uk> Sent: 04 April 2016 15:37 To: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Subject: Re: [firedrake] Integral is missing an integration domain. On 4 Apr 2016, at 14:59, Fryderyk Wilczynski <scfw@leeds.ac.uk><mailto: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 _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 06/04/16 17:14, Fryderyk Wilczynski wrote:
Thank you, that conditional worked.
I have realized that I should in fact have trial functions, (equations are linear in (n+1) terms), I have been mistakenly following example scripts for nonlinear problems.
I have changed
v, p, phi, mu = split(u)
rho = Function(P1)
to
v, p, phi, mu = TrialFunctions(W)
rho = TrialFunction(P1)
You have a term in a3: inner(chi, phi*grad(mu))*dx but phi and mu are both trial functions. Did you mean phi0? Lawrence
Yes, that should be phi0. Thank you. ________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: 06 April 2016 17:30 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Integral is missing an integration domain. On 06/04/16 17:14, Fryderyk Wilczynski wrote:
Thank you, that conditional worked.
I have realized that I should in fact have trial functions, (equations are linear in (n+1) terms), I have been mistakenly following example scripts for nonlinear problems.
I have changed
v, p, phi, mu = split(u)
rho = Function(P1)
to
v, p, phi, mu = TrialFunctions(W)
rho = TrialFunction(P1)
You have a term in a3: inner(chi, phi*grad(mu))*dx but phi and mu are both trial functions. Did you mean phi0? Lawrence
How is it going? Harald wrote this about rho eqn: We hebben het algoritme in sterke vorm gepresenteerd voor de helderheid. Vgl (54) kan direct: rho zit dan in dezelfde FEM space als phi. Vgl (57) doen we via L2 projectie. rho is in same space as phi. Eqn (57) via L2-projection. Check whether phi goes over 1 or below -1: that can cause troubles but should go better on finer grids. ________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Fryderyk Wilczynski <scfw@leeds.ac.uk> Sent: Wednesday, April 6, 2016 5:43 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Integral is missing an integration domain. Yes, that should be phi0. Thank you. ________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: 06 April 2016 17:30 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Integral is missing an integration domain. On 06/04/16 17:14, Fryderyk Wilczynski wrote:
Thank you, that conditional worked.
I have realized that I should in fact have trial functions, (equations are linear in (n+1) terms), I have been mistakenly following example scripts for nonlinear problems.
I have changed
v, p, phi, mu = split(u)
rho = Function(P1)
to
v, p, phi, mu = TrialFunctions(W)
rho = TrialFunction(P1)
You have a term in a3: inner(chi, phi*grad(mu))*dx but phi and mu are both trial functions. Did you mean phi0? Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (4)
-
Fryderyk Wilczynski
-
Lawrence Mitchell
-
Miklós Homolya
-
Onno Bokhove