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 ]
