I'm glad to hear that worked for you.

On Mar 21, 2015, at 7:36 PM, [email protected] wrote:

> I made the changes that you suggested. The program works perfectly now.
> Thank you for the help. This should shave hours off of my simulations.
> 
> Thanks again,
> Tony
> 
> 
>> It doesn't have anything to do with the solver.
>> 
>> The issue is that multiprocessing pickles the arguments to the
>> multiprocessed function. FiPy has known issues with pickling (see
>> https://github.com/usnistgov/fipy/issues/358 and
>> https://github.com/usnistgov/fipy/wiki/MeshPickling#sour-pickles (which
>> didn't survive migration from
>> http://matforge.org/fipy/wiki/MeshPickling#sourpickles)).
>> 
>> In this case, the problem is that our pickling code is unaware that you
>> have changed the mesh origin.
>> This situation is easily resolved by replacing
>> 
>>> mesh = fipy.Grid2D(dx=dx, dy=dp, nx=2*nx, ny=2*np, communicator=comm)
>>> mesh.origin = numpy.array([[-xMax],[-pMax]])
>> 
>> with
>> 
>>> mesh = fipy.Grid2D(dx=dx, dy=dp, nx=2*nx, ny=2*np, communicator=comm) +
>>> [[-xMax],[-pMax]]
>> 
>> 
>> As a final remark, I would encourage trimming a *lot* more to get
>> diagnostic examples. A single iteration is sufficient to demonstrate the
>> problem and is much easier to follow. I set n1 = n2 = 1 and looped 'for k
>> in range(1)'. Nothing more is needed. You can probably throw a lot more
>> away than that.
>> 
>> 
>> On Mar 13, 2015, at 4:16 PM, [email protected] wrote:
>> 
>>> 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 ]
>> 
>> 
>> _______________________________________________
>> 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 ]

Reply via email to