There are three possible ways to execute the code. One could execute
everything in serial, or use multiprocessing with processes=1, or
multiprocessing with processes=N. With my code, I'm finding that the two
multiprocessing results are the same but do not agree with the serial
result.
I'm attaching a trimmed down version of the program which still displays
this behavior. It is currently set up to test the impulse responses both
in serial and parallel with processes=1, forcing both methods to use
LinearPCGSolver. It also prints out lists of the impulse response after a
small time evolution.
Thank you,
Tony
# Import Modules
import fipy
import numpy
import multiprocessing
############################## Functions ############################
def evolution(phi, tStart, tEnd, mesh, xMax, pMax):
tEvo = tStart
X, P = mesh.cellCenters()
xFace, pFace = mesh.faceCenters()
dx = mesh.dx
dp = mesh.dy
nHalf = mesh.shape[1]/2
while(tEvo<tEnd):
timeStepDuration=((1./(nHalf))*(1./numpy.sqrt(2)))
if(tEvo+timeStepDuration>tEnd):
timeStepDuration = tEnd - tEvo
# Create Diffusion Term
gxx = 0.*X
gxp = 0.*X
gpx = 0.*X
gpp = 0.*X + 2.
dTerm = fipy.DiffusionTerm(fipy.CellVariable(mesh=mesh,
value=[[gxx,gxp],[gpx,gpp]]))
del gxx, gxp, gpx, gpp
# Create convection term
uX = 1.*pFace
uP = -1.*xFace - 2.*pFace
cCoeff = fipy.FaceVariable(mesh=mesh, value=[uX,uP])
cTerm = fipy.ExponentialConvectionTerm(cCoeff)
del uX, uP, cCoeff
# Create evolution equation
eq = fipy.TransientTerm()+cTerm==dTerm
# Evolve system
solver = fipy.solvers.LinearPCGSolver(tolerance=10**-20, \
iterations=1000, precon=None)
eq.solve(var=phi, dt=timeStepDuration, solver=solver)
norm = fipy.numerix.sum(phi.value, axis=0)*dx*dp
tEvo = tEvo + timeStepDuration
if(abs(norm-1)>10**(-10)):
s1 = 'Error: Distribution is no longer normalized.\n'
s2 = 'Abs(1-Norm) = '
s3 = repr(norm-1)
print(s1 + s2 + s3)
del s1, s2, s3
break
del timeStepDuration, cTerm, dTerm, eq, norm, solver
del tEvo, X, P, xFace, pFace, dx, dp, nHalf
del tStart, tEnd, mesh, xMax, pMax
return phi
def impulseResponse(argArray):
tStart = argArray[0]
tFinal = argArray[1]
mesh = argArray[2]
xMax = argArray[3]
pMax = argArray[4]
k = argArray[5]
dx = mesh.dx
dp = mesh.dy
nSide = mesh.nx
phi = fipy.CellVariable(mesh=mesh, name=r"$\rho(x,p)$")
temp = numpy.zeros(((nSide)**2,))
temp[k]=1.
norm = dx*dp
phi.setValue(temp/norm)
del dx, dp, nSide, temp, norm, argArray
response = evolution(phi, tStart, tFinal, mesh, xMax, pMax)
del phi, tStart, tFinal, mesh, xMax, pMax
print(repr(k)+' , '+repr(response.value[k]))
del k
return response
#####################################################################
#####################################################################
#################### User Specified Parameters ######################
# Initial time for evolution
tStart = 0.
# Final time for evolution
tFinal = 0.001
# Number of divisions per thermal length in phase space
n1 = 2
# Number of thermal lengths in phase space domain
n2 = 3
#####################################################################
#####################################################################
# Initialize phase-space lattice
nx = n1*n2
np = nx
dx = 1./n1
dp = 1./n1
xMax = dx*nx
pMax = dp*np
comm = fipy.tools.serial
mesh = fipy.Grid2D(dx=dx, dy=dp, nx=2*nx, ny=2*np, communicator=comm)
mesh.origin = numpy.array([[-xMax],[-pMax]])
del np, dx, dp
# Find all impulse response functions in parallel
pool = multiprocessing.Pool(processes=1)
responseStore = pool.map(impulseResponse, \
[[tStart, tFinal, mesh, xMax, pMax, k]\
for k in range((2*nx)**2)])
pool.close()
pool.join()
del pool
# Find all impulse response functions in serial
responseStoreSer = numpy.zeros(((2*nx)**2,(2*nx)**2))
for k in range((2*nx)**2):
responseStoreSer[k,:] = impulseResponse([tStart, tFinal, mesh, \
xMax, pMax, k])
> I'm not aware that we've encountered this, but I also don't think we've
> ever tried.
>
> Just to clarify, the issue is that the answer changes between processes=1
> or processes=2+?
>
> I don't think that any of FiPy's parallel code (e.g., using mpirun to
> invoke PyTrilinos in parallel) should be triggered by this use of
> multiprocessing, but we've never tried.
>
> I'll try to do some experiments. If you'd like to send me a complete
> (non)working example, I'll be happy to take a look at it.
>
>
> On Mar 13, 2015, at 2:46 PM, Tony Bartolotta <[email protected]> wrote:
>
>> FiPy Developers,
>>
>> I am currently using FiPy 3.1 with Python 2.7 on a laptop running Ubuntu
>> 14.04. I am attempting to find the impulse responses of a Fokker-Planck
>> equation. In the below snippets of code, funcCall(args) calculate an
>> impulse
>> response with FiPy using a non-parallelized solver.
>>
>>
>> Serial Calls
>>
>> for k in range(N):
>> responses[k,:] = funcCall(args)
>>
>>
>> Parallel Calls
>>
>> pool = multiprocessing.Pool(processes=1)
>> responses = pool.map(funcCall, [args for k in range(N)])
>> pool.close()
>> pool.join()
>>
>>
>> I am finding that FiPy returns different numerical results when I call
>> this
>> method in serial or in parallel. Strangely, even if the parallel calls
>> are
>> carried out by only one worker using the same solver, the results are
>> still
>> different.
>>
>> Using a debugging tool, all arguments appear to be correctly passed to
>> FiPy.
>> After considering a few test cases, the serial calls appear to arrive at
>> the
>> correct numerical solutions while the parallel version does not. I am
>> unsure
>> what is causing this behavior. Are you aware of any potential cause of
>> this
>> behavior?
>>
>> If necessary, I can post a more complete version of the code.
>>
>> Thank you,
>> Tony
>>
>> _______________________________________________
>> fipy mailing list
>> [email protected]
>> http://www.ctcms.nist.gov/fipy
>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
