Dear All,
I have a fenics code whose results depend on the number of processes when ran in parallel. Not surprisingly I am deeply concerned about this. After my analysis I am not sure whether it is a bug. Anyhow, if somebody knew about a workaround that would be great.
I am solving a steady state advection diffusion equation in 3D with a 2D velocity u(x,y). I solve the velocity field in a separate program, load it from disk, and use Expression to define a 3D velocity g(x,y,z)=u(x,y) which I can use when I define the bilinear form.
The general idea is outlined below. I played with my code and found that if I replace the line which is marked by #(*) below (e. g., I do not call u_old(x,y)), the problem does not occur.
from dolfin import *
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
#Some global variables
Order = 2
def MakeUFunction2(FlowFieldMesh="TestFlowMesh",FlowField="TestFlowField"):
old_mesh=Mesh(FlowFieldMesh+".xml")
V_old = FunctionSpace(old_mesh,"Lagrange",Order)
u_old = Function(V_old,FlowField+".xml")
u_old.set_allow_extrapolation(True)
class MyExpression(Expression):
def eval(self,value,x):
value[0] = u_old(x[0],x[1]) #(*)
#value[0]=0 #Now it is independent of number of processes
def value_shape(self):
return (1,)
g = MyExpression(degree=3)
print "Example output: rank "+str(rank) + ".....u_old(0,0)" +str(u_old(0,0)) #(**)
##########Debug lines to make sure that g/u_old was loaded correctly.
uip = interpolate(g,V_old)
fd = File("FlowDump.pvd")
d<<uip
##########End Debug.
return g
def SolveAdvectionDiffusionEquation(...): #Example
#(...)
a = inner(grad(c),grad(v))*dx+U*u*c.dx(2)*v*dx
#(...)
L = f*v*dx() + g*v*ds()
#(...)
I suspected that it is an issue with extrapolating u_old, so I switched off the extrapolation - set_allow_extrapolation(False)). In this case, the code crashes complaining that I tried to evaluate u_old at a value outside of its domain and suggests to switch extrapolation on.
Even more suspicious now, I added the line marked by #(**). If extrapolation is not allowed, the code crashes if the number of processes is >1, if it is allowed and the number of processes is >1 most processes obtain values which do not make sense (I expect zero).
Any ideas?
I suspect the problem is the following:
The function u_old is partitioned between the processes in 2D, since the mesh is 2D. Each process can evaluate u_old (and therefore g) only for those x,y-values for which it has the partition. That's why it crashes at line #(**) when extrapolation is switched off and gives some pretty far off values when extrapolation is switched on.
I believe that (based on the vtu-files when saved for paraview) the 3D mesh for the advection diffusion equation is partitioned along the z-axes between the processes. Thus, when each process evaluates g(x,y,z) = u_old(x,y) the parts of u_old(x,y) which are not part of the processes' share have to extrapolated, which leads to inaccurate results.
Working example: the function MakeUFunction2() should work, mesh and flow field are given here
