Source code for yade.mpy

# 2018 © Bruno Chareyre <bruno.chareyre@grenoble-inp.fr> 

'''
This module defines mpirun(), a parallel implementation of run() using a distributed memory approach. Message passing is done with mpi4py mainly, however some messages are also handled in c++ (with openmpi).

.. note:: Many internals of the mpy module listed on this page are not helpful to the user. Instead, please find :ref:`introductory material on mpy module<mpy>` in user manual.


Logic:
------

The logic for an initially centralized scene is as follows:

1. Instanciate a complete, ordinary, yade scene
2. Insert subdomains as special yade bodies. This is somehow similar to adding a clump body on the top of clump members
3. Broadcast this scene to all workers. In the initialization phase the workers will:

	- define the bounding box of their assigned bodies and return it to other workers
	- detect which assigned bodies are virtually in interaction with other domains (based on their bounding boxes) and communicate the lists to the relevant workers
	- erase the bodies which are neither assigned nor virtually interacting with the subdomain

4. Run a number of 'regular' iterations without re-running collision detection (verlet dist mechanism). In each regular iteration the workers will:

	- calculate internal and cross-domains interactions
	- execute Newton on assigned bodies (modified Newton skips other domains)
	- send updated positions to other workers and partial force on floor to master
	
5. When one worker triggers collision detection all workers will follow. It will result in updating the intersections between subdomains.

6. If enabled, bodies may be re-allocated to different domains just after a collision detection, based on a filter. Custom filters are possible. One is predidefined here (medianFilter)

Rules:
------

	#- intersections[0] has 0-bodies (to which we need to send force)
	#- intersections[thisDomain] has ids of the other domains overlapping the current ones
	#- intersections[otherDomain] has ids of bodies in _current_ domain which are overlapping with other domain (for which we need to send updated pos/vel)

Hints:
------

	#- handle subD.intersections with care (same for mirrorIntersections). subD.intersections.append() will not reach the c++ object. subD.intersections can only be assigned (a list of list of int)

'''

import sys,os,inspect
import time
import numpy as np
import yade.runtime
this = sys.modules[__name__]


#from yade import *
from yade.utils import *
from yade.wrapper import *
#import yade.runtime
from yade import timing; timing.mpi={} #prepare a dictionnary for mpi-related stats

# for coloring bodies
import colorsys

## Config flags

ACCUMULATE_FORCES=True #control force summation on master's body. FIXME: if false master goes out of sync since nothing is blocking rank=0 thread
VERBOSE_OUTPUT=False
NO_OUTPUT=False
COLOR_OUTPUT=True # color mprint output per rank / turn off here for sphinx build mainly
MAX_RANK_OUTPUT=5 #: larger ranks will be skipped in mprint
TIMEOUT = 60 # after this delay self-kill, so we don't leave the OS with nasty orphan processes
SEND_SHAPES=False #if false only bodies' states are communicated between threads, else shapes as well (to be implemented)
ERASE_REMOTE = True # True is MANDATORY. Erase bodies not interacting wit a given subdomain? else keep dead clones of all bodies in each scene
ERASE_REMOTE_MASTER = True # erase remotes on master or keep them for fast merge (updating just b.state)
OPTIMIZE_COM=True
USE_CPP_MPI=True and OPTIMIZE_COM
YADE_TIMING=False #report timing.stats()?
MERGE_SPLIT = False
MERGE_W_INTERACTIONS = True
COPY_MIRROR_BODIES_WHEN_COLLIDE = True  # True is MANDATORY 
RESET_SUBDOMAINS_WHEN_COLLIDE = False
DOMAIN_DECOMPOSITION = False
NUM_MERGES = 0
SEND_BYTEARRAYS = True
ENABLE_PFACETS = False    #PFacets need special (and expensive) tricks, if PFacets are not used skip the tricks
DISTRIBUTED_INSERT = False  #if True each worker is supposed to "O.bodies.insertAtId" its own bodies
REALLOCATE_FREQUENCY = 0  # if >0 checkAndCollide() will automatically reallocate bodies to subdomains, if =1 realloc. happens each time collider is triggered, if >1 it happens every N trigger
REALLOCATE_FILTER = None # pointer to filtering function, will be set to 'medianFilter' hereafter, could point to other ones if implemented
FAIR_SHARE = True # try to keep equal no. bodies per subdomain when reallocating  
AUTO_COLOR = True
MINIMAL_INTERSECTIONS = False # Reduces the size of position/velocity comms (at the end of the colliding phase, we can exclude those bodies with no interactions besides body<->subdomain from intersections). 
REALLOCATE_MINIMAL = False # if true, intersections are minimized before reallocations, hence minimizing the number of reallocated bodies
MASTER_UPDATE_STATES = False # does master thread need to provide updated pos/vel? else they will be fixed or moved by the workers (e.g. MotionEngines)
fibreList = []
FLUID_COUPLING = False
fluidBodies = [] 
USE_CPP_REALLOC = True
USE_CPP_INTERS = False #sending intersections using mpi4py sometimes fails (dependent on mpi4py version, needs confirmation) (ERR : MPI_ERR_TRUNCATE)


### Internals

comm = None
comm_slave = None # will be MPI.Comm.Get_parent() at runtime, until we disconnect 
rank = None
numThreads = None
MPI = None # will be mpi4py.MPI after configure()

waitingCommands=False		# are workers currently interactive?
userScriptInCheckList=""	# detect if mpy is executed by checkList.py
caller_name=""			# name of the executed script (typically user script, alternatively yadeSphinx.py or checkList.py)
_REALLOC_COUNT=0
colorScale=None
# for coloring processes outputs differently
bcolors=['\x1b[95m','\x1b[94m','\x1b[93m','\x1b[92m','\x1b[91m','\x1b[90m','\x1b[95m','\x1b[93m','\x1b[91m','\x1b[1m','\x1b[4m','\x1b[0m']

#tags for mpi messages
_SCENE_=11
_SUBDOMAINSIZE_=12
_INTERSECTION_=13
_ID_STATE_SHAPE_=14
_FORCES_=15
_MIRROR_INTERSECTIONS_ = 16
_POS_VEL_ = 17
_BOUNDS_ = 18
_MASTER_COMMAND_ = 19
_RETURN_VALUE_ = 20
_ASSIGNED_IDS_ = 21
_GET_CONNEXION_= 22



	
[docs]def makeColorScale(n=numThreads): scale= [(0.3+random.random())*Vector3(colorsys.hsv_to_rgb(value*1.0/n, 1, 1)) for value in range(0, n)] from random import shuffle random.seed(1) shuffle(scale) return scale
[docs]def mprint(*args,force=False): #this one will print regardless of VERBOSE_OUTPUT """ Print with rank-reflecting color regardless of mpy.VERBOSE_OUTPUT, still limited to rank<=mpy.MAX_RANK_OUTPUT """ if (NO_OUTPUT or rank>MAX_RANK_OUTPUT) and not force: return m=bcolors[min(rank,len(bcolors)-2)] if COLOR_OUTPUT else '' resetFont='\x1b[0m' if COLOR_OUTPUT else '' if rank==0: m+='Master: ' else: m+='Worker'+str(rank)+": " for a in args: m+=str(a)+" " print (m+resetFont)
# sys.stdout.flush() # maybe?
[docs]def wprint(*args): """ Print with rank-reflecting color, *only if* mpy.VERBOSE_OUTPUT=True (else see :yref:`yade.mpy.mprint`), limited to rank<=mpy.MAX_RANK_OUTPUT """ if not VERBOSE_OUTPUT: return mprint(*args)
### Here we detect current MPI context *only if needed*, i.e. importing mpy without actually using it will not call configure() (then no MPI related warnings) #if yade.runtime.opts.mpi_mode: #configure()
[docs]def colorDomains(): ''' Apply color to body to reflect their subdomain idx ''' global colorScale if len(colorScale)<numThreads: colorScale = makeColorScale(numThreads) for b in O.bodies: b.shape.color=colorScale[b.subdomain]
## Initialization
[docs]def configure(): # calling this function will import mpi4py.MPI, ''' Import MPI and define context, configure will no spawn workers by itself, that is done by initialize() openmpi environment variables needs to be set before calling configure() ''' global comm, comm_slave, rank, numThreads, colorScale, MPI os.environ["OMPI_MCA_rmaps_base_oversubscribe"]="1" # needed here, after importing MPI is too late (or there is a way to update flags before the spawn?) os.environ["OMPI_MCA_pmix_server_max_wait"]="3600" os.environ["OMPI_MCA_pmix_base_exchange_timeout"]="3600" from mpi4py import MPI worldComm = MPI.COMM_WORLD color = 3; key =0; comm = worldComm.Split(color, key) # if OFOAM coupled, split communicator comm_slave = comm.Get_parent() if comm_slave!=MPI.COMM_NULL: # if executor is a spawned worker merge comm with master comm=comm_slave.Merge() rank = comm.Get_rank() # set rank and numThreads if rank>0: sys.stderr.write=sys.stdout.write # so we see error messages from workers in terminal numThreads = comm.Get_size() colorScale= makeColorScale(numThreads)
configure() # needed at import time, else we can't check rank or numThreads right after import
[docs]def disconnect(): ''' Kill all mpi processes, leaving python interpreter to rank 0 as in single-threaded execution. The scenes in workers are lost since further reconnexion to mpi will just spawn new processes. The scene in master thread is left unchanged. ''' global comm,comm_slave,numThreads if comm==None: return #not connected if rank==0: # exit the interactive mode on master _after_ telling workers to exit wprint("sending exit command") sendCommand(executors="slaves",command="exit",wait=False) if comm!=MPI.COMM_NULL and comm_slave!=MPI.COMM_NULL: # only true after spawn/merge, false if executed with 'mpirun -np yade' wprint("disconnecting") # (1) should be a Disconnect(), no a Free(), but we have an issue with openmpi it seems # https://bitbucket.org/mpi4py/mpi4py/issues/176/disconnect-hangs-with-openmp31-python-38 # (2) for some reason even comm.Free() causes trouble: multiple parallel scripts can be # chained in 'yade --check' after free(), but not all of them. Ultimately it crashes like this on checkMPISilo.py (unbuntu20.04): # Signal: Segmentation fault (11) # Signal code: Address not mapped (1) # Failing at address: 0x49 # [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x46210)[0x7ff3982d3210] # [ 1] /usr/lib/x86_64-linux-gnu/openmpi/lib/openmpi3/mca_pml_ob1.so(mca_pml_ob1_add_comm+0x169)[0x7ff378060789] # commenting out "Free()" seems to workaround, so let it be. Unclear if it is a bug here or an issue with openmpi comm.barrier() comm.Free() comm_slave.Disconnect() if rank>0: # kill workers exit else: mprint("mpy already disconnected, size=", ) comm=None comm_slave=None numThreads=None yade.runtime.opts.mpi_mode=False O.splittedOnce=False O.splitted=False
[docs]def makeMpiArgv(): stack = inspect.stack() interactive=False yadeArgv = sys.yade_argv.copy() # the arguments to be passed to workers if '--check' in yadeArgv: yadeArgv.remove("--check") # escape this since the check scripts will be called from master thread already if '--checkall' in yadeArgv: yadeArgv.remove("--checkall") # 'mpy.py' will be the caller in first N lines, followed by a py script or some interactive command # we search backward the first line where caller is not mpy.py to identify the script calling initialize # in a way to handle execfiled scripts (e.g. yade --check) lastCalledScript = "" searchDepth=6 for line in range(searchDepth): if stack[line][3]=='run_code': # we are in interactive shell, no script to run interactive = True break for kArg in range(1,len(stack[line])): if str(stack[line][kArg])[-3:]=='.py': if str(stack[line][kArg])[-6:]=='mpy.py': break # jump to next line in the stack else: lastCalledScript = stack[line][kArg] if lastCalledScript != "": break if lastCalledScript != "": kArg=1 while kArg<len(yadeArgv): if yadeArgv[kArg][-3:]=='.py': yadeArgv[kArg]=lastCalledScript break kArg+=1 if kArg==len(yadeArgv): # yadeArgv empty of any script => execfiled from python shell yadeArgv.append(lastCalledScript) return yadeArgv,interactive
[docs]def initialize(np): global comm,comm_slave,rank,numThreads,userScriptInCheckList,colorScale,waitingCommands if (comm!=None and yade.runtime.opts.mpi_mode==True): if rank==0: disconnect() # reset to virgin context else: exit # kill myself, I'm part of deprecated pool if comm==None: configure() # should only happen after a despawn process_count = comm.Get_size() if rank==0: # MASTER yadeArgv,waitingCommands = makeMpiArgv() numThreads=np colorScale= makeColorScale(numThreads) mprint("will spawn ",numThreads-process_count," workers running:",yadeArgv[0], yadeArgv[1:]) comm_slave = MPI.COMM_WORLD.Spawn(yadeArgv[0], args=yadeArgv[1:],maxprocs=numThreads-process_count) comm = comm_slave.Merge() yade.runtime.opts.mpi_mode=True if waitingCommands: declareMasterInteractive() #in interactive session tell the workers to ignore timeout else: #WORKERS wprint("spawned") #initialize subdomains. For Master it will be used storage and comm only, for workers it will be over-written in the split operation O.subD=Subdomain() O.subD.comm = comm return rank,numThreads
[docs]def spawnedProcessWaitCommand(): global waitingCommands if waitingCommands: return if comm==None: configure() wprint("I'm now waiting") waitingCommands = True s=MPI.Status() while 1: maxTime = time.time()+TIMEOUT waitTime=0 while (TIMEOUT==0 or waitTime<maxTime) and not comm.Iprobe(source=MPI.ANY_SOURCE, tag=_MASTER_COMMAND_, status=s): time.sleep(0.001) waitTime=time.time() if TIMEOUT!=0 and waitTime>=maxTime: mprint("orphan yade worker disconnecting") break command = comm.recv(source=s.source,tag=_MASTER_COMMAND_) if command=="exit": #this is to terminate the waiting loop remotely comm.send(0,dest=s.source,tag=_RETURN_VALUE_) break wprint("will now execute ",command) try: exec(command) except: comm.send(None,dest=s.source,tag=_RETURN_VALUE_) mprint(sys.exc_info(),force=True) raise waitingCommands = False disconnect() # this will kill the workers
[docs]def declareMasterInteractive(): ''' This is to signal that we are in interactive session, so TIMEOUT will be reset to 0 (ignored) ''' global TIMEOUT TIMEOUT = 0 if rank==0: waitingCommands=True sendCommand("slaves","sys.modules['yade.mpy'].TIMEOUT=0",False)
[docs]def sendCommand(executors,command,wait=True,workerToWorker=False): ''' Send a command to a worker (or list of) from master or from another worker. Accepted executors are "i", "[i,j,k]", "slaves", "all" (then even master will execute the command). ''' start=time.time() if (rank>0 and not workerToWorker): wprint("sendCommand ignored by worker", rank, ", pass workerToWorker=True to force it"); return if (executors=="all"): executors=list(range(numThreads)) if (executors=="slaves"): executors=list(range(1,numThreads)) argIsList=isinstance(executors,list) toMaster = (argIsList and 0 in executors) or executors==0 if (toMaster and rank>0): mprint("workers cannot sendCommand to master (only master to itself)") if not argIsList: executors = [executors] #if 0 in executors: mprint("master does not accept mpi commands"); return if len(executors)>numThreads: mprint("executors > numThreads"); return if wait and not command=="exit":#trick command to make it return a result by mpi commandSent="resCommand="+command+";comm.send(resCommand,dest="+str(rank)+",tag=_RETURN_VALUE_)" else: commandSent=command reqs=[] for w in executors: #note: if the return from this isend() is not appended to a list we have random deadlock if (w>0): reqs.append(comm.isend(commandSent,dest=w,tag=_MASTER_COMMAND_) ) resCommand=[] if toMaster:#eval command on master since it wasn't done yet try: resCommand = [eval(command)] except: resCommand = [None] mprint(sys.exc_info(),force=True) raise if wait: resCommand=resCommand+ [comm.recv(source=w,tag=_RETURN_VALUE_) for w in executors if w>0] wprint("sendCommand returned in "+str(time.time()-start)+" s") return (resCommand if argIsList else resCommand[0]) else: for r in reqs: r.wait() # make sure everything has been sent return None
[docs]def probeRecvMessage(source, tag): msgStat = MPI.Status() comm.Probe(source=source, tag=tag, status=msgStat) if msgStat.tag == tag : print("message size recvd") data = bytearray(msgStat.Get_count(MPI.BYTE)) comm.Recv([data, MPI.BYTE], source=source, tag=tag) return data
### TIMING ###
[docs]def recordMpiTiming(name,val): ''' append val to a list of values defined by 'name' in the dictionnary timing.mpi ''' if not name in yade.timing.mpi.keys(): yade.timing.mpi[name]=[] yade.timing.mpi[name].append(val)
[docs]def mpiStats(): if (waitingCommands and rank==0): sendCommand(list(range(1,numThreads)),"mpiStats()") dat = yade.timing.mpi dat = comm.gather(dat,root=0) return dat
yade.timing.mpiStats=mpiStats
[docs]class Timing_comm(): def __init__(self): self.timings={}
[docs] def clear(self): self.timings={}
[docs] def print_all(self): time.sleep((numThreads-rank)*0.001) message = "COMMUNICATION TIMINGS:\n" max_string_len = len(max(self.timings.keys(),key=len)) for k,v in sorted(self.timings.items(), key=lambda x: x[1][1], reverse=True): message += ("{:<"+str(max_string_len)+"}").format(k) + " " + str(v) + "\n" mprint(message)
[docs] def enable_timing(comm_function): def wrapper(self,timing_name, *args, **kwargs): #pre-exec ti=time.time() #exec rvalue=comm_function(self, *args, **kwargs) #post-exec if(not timing_name in self.timings.keys()): self.timings[timing_name]=[0,0] self.timings[timing_name][0]+=1 self.timings[timing_name][1]+=time.time()-ti return rvalue return wrapper
[docs] @enable_timing def send(self, *args, **kwargs): return comm.send(*args,**kwargs)
[docs] @enable_timing def recv(self, *args, **kwargs): return comm.recv(*args,**kwargs)
[docs] @enable_timing def bcast(self, *args, **kwargs): return comm.bcast(*args,**kwargs)
[docs] @enable_timing def allreduce(self, *args, **kwargs): return comm.allreduce(*args,**kwargs)
[docs] @enable_timing def Gather(self, *args, **kwargs): return comm.Gather(*args, **kwargs)
[docs] @enable_timing def Gatherv(self, *args, **kwargs): return comm.Gatherv(*args, **kwargs)
[docs] @enable_timing def Allgather(self, *args, **kwargs): return comm.Allgather(*args, **kwargs)
#this is to time the cpp messages
[docs] @enable_timing def mpiWaitReceived(self,*args, **kwargs): return O.subD.mpiWaitReceived(*args, **kwargs)
[docs] @enable_timing def mpiWait(self,req): return comm.wait(req)
[docs] @enable_timing def mpiSendStates(self,*args, **kwargs): return O.subD.mpiSendStates(*args, **kwargs)
timing_comm = Timing_comm()
[docs]def receiveForces(subdomains): ''' Accumulate forces from subdomains (only executed by master process), should happen after ForceResetter but before Newton and before any other force-dependent engine (e.g. StressController), could be inserted via yade's pyRunner. ''' if 0: #non-blocking: this should be replaced by isend+recv since irecv may have too small default buffer size, use blocking for now reqForces=[] #for sd in subdomains: for sd in O.subD.intersections[0]: #wprint( "master getting forces from "+str(b.subdomain)+"(id="+str(b.id)+")") reqForces.append(comm.irecv(None,sd, tag=_FORCES_)) #wprint( "master got forces from "+str(b.subdomain)+": "+str(forces)) for r in reqForces: forces=r.wait() for ft in forces: #wprint( "adding force "+str(ft[1])+" to body "+str(ft[0])) O.forces.addF(ft[0],ft[1]) O.forces.addT(ft[0],ft[2]) else: for sd in subdomains: forces=timing_comm.recv("isendRecvForces",source=sd, tag=_FORCES_) #wprint( "master got forces from "+str(sd)+": "+str(forces)+" iter="+str(O.iter)+" dt="+str(O.dt)) for ft in forces: #wprint( "adding force "+str(ft[1])+" to body "+str(ft[0])) O.forces.addF(ft[0],ft[1]) O.forces.addT(ft[0],ft[2])
[docs]def shrinkIntersections(): ''' Reduce intersections and mirrorIntersections to bodies effectively interacting with another statefull body form current subdomain This will reduce the number of updates in sendRecvStates Initial lists are backed-up and need to be restored (and all states updated) before collision detection (see checkAndCollide()) ''' if hasattr(O.subD,"fullIntersections") and O.subD.fullIntersections!=None: mprint("Problem HERE!!!!!!!!!!!") O.subD.fullIntersections = O.subD.intersections O.subD.fullMirrorIntersections = O.subD.mirrorIntersections if (rank==0): return 0,0 oriLen=sum([len(c) for c in O.subD.intersections]) res=O.subD.filterIntersections() ##ints = O.subD.mirrorIntersections ## What follows did not show any anomalous result, so we better skip communication ## of something every worker agree on splontaneously #reqs=[] #for other in O.subD.intersections[rank]: #if other==0: continue #reqs.append([other,comm.irecv(None, other, tag=_MIRROR_INTERSECTIONS_)]) #comm.send(O.subD.intersections[other],dest=other,tag=_MIRROR_INTERSECTIONS_) #for r in reqs: #ints[r[0]]=r[1].wait() #if ints[r[0]]!=O.subD.mirrorIntersections[r[0]]: #mprint("inconsistency in the filtering of intersections[",r[0],"]:",len(ints[r[0]]),"received vs.",len(O.subD.mirrorIntersections[r[0]])) #O.subD.mirrorIntersections = ints #that's because python wrapping only enable assignment return res,oriLen
[docs]def checkAndCollide(): ''' return true if collision detection needs activation in at least one SD, else false. If COPY_MIRROR_BODIES_WHEN_COLLIDE run collider when needed, and in that case return False. ''' global _REALLOC_COUNT needsCollide = int(typedEngine("InsertionSortCollider").isActivated()) if(needsCollide!=0): wprint("triggers collider at iter "+str(O.iter)) needsCollide = timing_comm.allreduce("checkcollider",needsCollide,op=MPI.SUM) if needsCollide: if(COPY_MIRROR_BODIES_WHEN_COLLIDE): if MINIMAL_INTERSECTIONS: if hasattr(O.subD,"fullIntersections") and O.subD.fullIntersections!=None: # if we have tricked intersections in previous steps we set them backto full content in order to update all positions before colliding O.subD.intersections = O.subD.fullIntersections O.subD.mirrorIntersections = O.subD.fullMirrorIntersections O.subD.fullMirrorIntersections=O.subD.fullIntersections=None #clear, so we know we can re-use #else: mprint("fullIntersections not initialized (first iteration or rank=0)") sendRecvStates() # triggers comm # parallel collision detection (incl. insertion of newly intersecting bodies) parallelCollide() # reallocate and/or minimize intersections if REALLOCATE_FREQUENCY>0: _REALLOC_COUNT+=1 if _REALLOC_COUNT>=REALLOCATE_FREQUENCY: #comm.barrier() #we will modify intersections while they can still be accessed by calls to mpi in parallelCollide() if (REALLOCATE_MINIMAL): # shrink mprint("don't use REALLOCATE_MINIMAL. It seems broken for the moment") r=shrinkIntersections() #if we filter before reallocation we minimize the reallocations #mprint("filtered out (1)",r[0],"of",r[1]) reallocateBodiesToSubdomains(REALLOCATE_FILTER,blocking=True) _REALLOC_COUNT=0 if (MINIMAL_INTERSECTIONS): #filter here, even if already done before, since realloc updated intersections #if rank>0: r=shrinkIntersections() recordMpiTiming("filteredInts",r[0]); recordMpiTiming("totInts",r[1]); recordMpiTiming("interactionsInts",len(O.interactions)); recordMpiTiming("iterInts",O.iter) #mprint("filtered out (2)",r[0],"of",r[1]) return False else: return True return False
[docs]def unboundRemoteBodies(): ''' Turn bounding boxes on/off depending on rank ''' for b in O.bodies:# unbound the bodies assigned to workers (not interacting directly with other bodies in master scene) if not b.isSubdomain and b.subdomain!=rank: b.bounded=False
[docs]def reboundRemoteBodies(ids): ''' update states of bodies handled by other workers, argument 'states' is a list of [id,state] (or [id,state,shape] conditionnaly) ''' if isinstance(ids,list): for id in ids: b = O.bodies[id] if b and not isinstance(b.shape,GridNode): b.bounded=True else: #when passing numpy array we need to convert 'np.int32' to 'int' for id in ids: b = O.bodies[id.item()] if b and not isinstance(b.shape,GridNode): b.bounded=True
[docs]def updateDomainBounds(subdomains): #subdomains is the list of subdomains by body ids ''' Update bounds of current subdomain, broadcast, and receive updated bounds from other subdomains Precondition: collider.boundDispatcher.__call__() ''' wprint( "Updating bounds: "+str(subdomains)) if(rank==0): send_buff=np.zeros(6)*np.nan else: subD=O.bodies[subdomains[rank-1]].shape #shorthand to shape of current subdomain send_buff=np.append(subD.boundsMin,subD.boundsMax) recv_buff = np.empty(6*numThreads) timing_comm.Allgather("updateDomainBounds",send_buff,recv_buff) for r in range(1,numThreads): O.bodies[subdomains[r-1]].shape.boundsMin = recv_buff[6*r:6*r+3] O.bodies[subdomains[r-1]].shape.boundsMax = recv_buff[3+6*r:6+6*r]
#if(VERBOSE_OUTPUT):#condition here to avoid concatenation overhead #mprint("Updated ", O.bodies[subdomains[r-1]].subdomain, " with min=", O.bodies[subdomains[r-1]].shape.boundsMin," and max=", O.bodies[subdomains[r-1]].shape.boundsMax) def maskedPFacet(pf, boolArray): ''' List bodies within a facet selectively, the ones marked 'True' in boolArray (i.e. already selected from another facet) are discarded ''' l=[] for id in [pf.node1.id, pf.node2.id, pf.node3.id, pf.conn1.id, pf.conn2.id, pf.conn3.id]: if not boolArray[id]: l.append(id) boolArray[id]=True
[docs]def maskedPFacet(b, boolArray): ''' List bodies within a facet selectively, the ones marked 'True' in boolArray (i.e. already selected from another facet) are discarded ''' l=[] pf=b.shape for id in [b.id,pf.node1.id, pf.node2.id, pf.node3.id, pf.conn1.id, pf.conn2.id, pf.conn3.id]: if not boolArray[id]: l.append(id) boolArray[id]=True return l
[docs]def maskedConnection(b, boolArray): ''' List bodies within a facet selectively, the ones marked 'True' in boolArray (i.e. already selected from another facet) are discarded ''' l=[] pf=b.shape for id in [b.id,pf.node1.id, pf.node2.id]: if not boolArray[id]: l.append(id) boolArray[id]=True return l
[docs]def genLocalIntersections(subdomains): ''' Defines sets of bodies within current domain overlapping with other domains. The structure of the data for domain 'k' is: [[id1, id2, ...], <----------- intersections[0] = ids of bodies in domain k interacting with master domain (subdomain k itself excluded) [id3, id4, ...], <----------- intersections[1] = ids of bodies in domain k interacting with domain rank=1 (subdomain k itself excluded) ... [domain1, domain2, domain3, ...], <---------- intersections[k] = ranks (not ids!) of external domains interacting with domain k ... ] ''' intersections=[[] for n in range(numThreads)] for sdId in subdomains: #grid nodes or grid connections could be appended twice or more, as they can participate in multiple pfacets and connexions #this bool list is used to append only once subdIdx=O.bodies[sdId].subdomain intrs=O.interactions.withBodyAll(sdId) if (ENABLE_PFACETS): appended = np.repeat([False],len(O.bodies)) #special case when we get interactions with current domain, only used to define interactions with master, otherwise some intersections would appear twice if subdIdx==rank: if not ENABLE_PFACETS: # no PFacets, same block duplicated below for the opposite case - only way to move ENABLE_PFACETS condition out of the loop for i in intrs: otherId=i.id1 if i.id2==sdId else i.id2 b=O.bodies[otherId] if not b:continue #in case the body was deleted if b.subdomain==0: intersections[0].append(otherId) else: for i in intrs: otherId=i.id1 if i.id2==sdId else i.id2 b=O.bodies[otherId] if not b:continue #in case the body was deleted if b.subdomain==0: if isinstance(b.shape,PFacet): intersections[0]+= maskedPFacet(b, appended); continue if isinstance(b.shape,GridConnection): intersections[0]+=maskedConnection(b, appended); continue #else (standalone body, normal case) intersections[0].append(otherId) if len(intersections[0])>0: intersections[subdIdx].append(0) continue # normal case if not ENABLE_PFACETS: # same as above for i in intrs: otherId=i.id1 if i.id2==sdId else i.id2 b=O.bodies[otherId] if not b:continue #in case the body was deleted if b.subdomain!=rank: continue if b.isSubdomain : intersections[rank].append(subdIdx) #intersecting subdomain (will need to receive updated positions from there) else: intersections[subdIdx].append(otherId) else: for i in intrs: otherId=i.id1 if i.id2==sdId else i.id2 b=O.bodies[otherId] if not b:continue #in case the body was deleted if b.subdomain!=rank: continue if b.isSubdomain : intersections[rank].append(subdIdx) #intersecting subdomain (will need to receive updated positions from there) else: if isinstance(b.shape,PFacet): intersections[subdIdx]+= maskedPFacet(b, appended); continue if isinstance(b.shape,GridConnection): intersections[subdIdx]+=maskedConnection(b, appended); continue #else (standalone body, normal case) intersections[subdIdx].append(otherId) #for master domain set list of interacting subdomains (could be handled above but for the sake of clarity complex if-else-if are avoided for now) if rank==0 and len(intersections[subdIdx])>0: intersections[0].append(subdIdx) #wprint( "found "+str(len(intrs))+" intersections"+str(intersections)) return intersections
[docs]def updateRemoteStates(states, setBounded=False): ''' update states of bodies handled by other workers, argument 'states' is a list of [id,state] (or [id,state,shape] conditionnaly) ''' ids=[] for bst in states: #print bst[0],O.bodies[bst[0]] ids.append(bst[0]) b=O.bodies[bst[0]] b.state=bst[1] #if SEND_SHAPES: b.shape=bst[2] if setBounded and not isinstance(b.shape,GridNode): b.bounded=True return ids
[docs]def genUpdatedStates(b_ids): ''' return list of [id,state] (or [id,state,shape] conditionnaly) to be sent to other workers ''' return [[id,O.bodies[id].state] for id in b_ids] if not SEND_SHAPES else [[id,O.bodies[id].state,O.bodies[id].shape] for id in b_ids]
############# COMMUNICATIONS ################" statesCommTime=0
[docs]def sendRecvStates(): global statesCommTime if rank==0 and not MASTER_UPDATE_STATES: comm.barrier() return # master has just nothing to do if workers don't need updated pos/vel start=time.time() #____1. get ready to receive positions from other subdomains pstates = [] buf = [] #heuristic guess, assuming number of intersecting is ~linear in the number of rows, needs if rank!=0: #the master process never receive updated states (except when gathering) for otherDomain in O.subD.intersections[rank]: if len(O.subD.mirrorIntersections[otherDomain])==0: continue #can happen if MINIMAL_INTERSECTIONS if otherDomain==0 and not MASTER_UPDATE_STATES: continue # skip, we don't need news from master if not USE_CPP_MPI: buf.append(bytearray(1<<22)) #FIXME: smarter size? this is for a few thousands states max (empirical); bytearray(1<<24) = 128 MB pstates.append( comm.irecv(buf[-1],otherDomain, tag=_ID_STATE_SHAPE_)) #warning leaving buffer size undefined crash for large subdomains (MPI_ERR_TRUNCATE: message truncated) else: O.subD.mpiIrecvStates(otherDomain) #use yade's messages (coded in cpp) #____2. broadcast new positions (should be non-blocking if n>2, else lock) - this includes subdomain bodies intersecting the current one reqs=[] for k in O.subD.intersections[rank]: if k==rank or k==0: continue #don't broadcast to itself... OTOH this list intersections[rank] will be used to receive if len(O.subD.intersections[k])==0: continue #can happen if MINIMAL_INTERSECTIONS #if len(b_ids)>0:#skip empty intersections, it means even the bounding boxes of the corresponding subdomains do not overlap wprint("sending "+str(len(O.subD.intersections[k]))+" states to "+str(k)) if not OPTIMIZE_COM: timing_comm.send("sendRecvStates",genUpdatedStates(O.subD.intersections[k]), dest=k, tag=_ID_STATE_SHAPE_) #should be non-blocking if n>2, else lock? else: if not USE_CPP_MPI: reqs.append(comm.isend(O.subD.getStateValues(k), dest=k, tag=_ID_STATE_SHAPE_)) #should be non-blocking if n>2, else lock? else: timing_comm.mpiSendStates("mpiSendStates",k) for r in reqs: r.wait() #empty if USE_CPP_MPI #____3. receive positions and update bodies if rank==0: return #positions sent from master, done. Will receive forces instead of states if not USE_CPP_MPI: nn=0 for ss in pstates: states=ss.wait() if not OPTIMIZE_COM: updateRemoteStates(states) else: O.subD.setStateValuesFromIds(O.subD.mirrorIntersections[O.subD.intersections[rank][nn]],states) nn+=1 else: for otherDomain in O.subD.intersections[rank]: if len(O.subD.mirrorIntersections[otherDomain])==0: continue #can happen if MINIMAL_INTERSECTIONS if otherDomain==0 and not MASTER_UPDATE_STATES: continue timing_comm.mpiWaitReceived("mpiWaitReceived(States)",otherDomain) O.subD.setStateValuesFromBuffer(otherDomain) comm.barrier() statesCommTime+=(time.time()-start)
[docs]def isendRecvForces(): ''' Communicate forces from subdomain to master Warning: the sending sides (everyone but master) must wait() the returned list of requests ''' O.freqs=[] #keep that one defined even if empty, it is accessed in other functions #TDOD: FORCES FROM FLUID DOMAIN BOXES!!!! if ACCUMULATE_FORCES: if rank!=0: if not 0 in O.subD.intersections[rank]: return if FLUID_COUPLING: forces0 = [] for id in O.subD.mirrorIntersections[0]: if not isinstance(O.bodies[id].shape, FluidDomainBbox): forces0.append([id, O.forces.f(id), O.forces.t(id)]) else: forces0=[[id,O.forces.f(id),O.forces.t(id)] for id in O.subD.mirrorIntersections[0]] #wprint ("worker "+str(rank)+": sending "+str(len(forces0))+" "+str("forces to 0 ")) #O.freqs.append(comm.isend(forces0, dest=0, tag=_FORCES_)) timing_comm.send("isendRecvForces",forces0, dest=0, tag=_FORCES_) else: #master receiveForces(O.subD.intersections[0])
[docs]def waitForces(): ''' wait until all forces are sent to master. O.freqs is empty for master, and for all threads if not ACCUMULATE_FORCES ''' for r in O.freqs: r.wait()
##### INITIALIZE MPI ######### # Flag used after import of this module, turned True after scene is distributed O.splitted=False O.splittedOnce=False #after the first split we have additional bodies (Subdomains) and engines in the merged scene, use this flag to know
[docs]def mergeScene(): if(rank==0 and waitingCommands): sendCommand("slaves","mergeScene()",False) if O.splitted: if MERGE_W_INTERACTIONS or ERASE_REMOTE_MASTER or DISTRIBUTED_INSERT: O.subD.mergeOp() sendRecvStatesRunner.dead = isendRecvForcesRunner.dead = waitForcesRunner.dead = collisionChecker.dead = True O.splitted=False #collider.doSort = True if (AUTO_COLOR): colorDomains() global NUM_MERGES; NUM_MERGES +=1; else: if rank>0: # Workers send_buff=np.asarray(O.subD.getStateBoundsValuesFromIds([b.id for b in O.bodies if b.subdomain==rank])) size=np.array(len(send_buff),dtype=int) else: #Master send_buff=np.array([0]) size=np.array(0,dtype=int) sizes=np.empty(numThreads,dtype=int) # Master get sizes from all workers timing_comm.Gather("mergeScene_sizes",size,sizes,root=0) if(rank==0): # MASTER # Alloc sizes for workers dat=np.ones(sizes.sum(),dtype=np.float64) # Displacement indexes where data should be stored/received in targeted array # dspl should be visible by everyone dspl=np.empty(numThreads, dtype=int) dspl[0] = 0 for i in range(1, len(sizes)): dspl[i] = dspl[i-1] + sizes[i-1]; else: dspl=None dat=None # data sent = [data, size of data] (for each worker) # data recv = [allocated target_array, array of different sizes, displacement, data type] timing_comm.Gatherv("mergeScene_data",[send_buff, size], [dat, sizes, dspl, MPI.DOUBLE], root=0) if(rank==0): #master for worker_id in range(1, numThreads): # generate corresponding ids (order is the same for both master and worker) #ids = [b.id for b in O.bodies if b.subdomain==worker_id] ids = O.bodies[O.subD.subdomains[worker_id-1]].shape.ids+[O.subD.subdomains[worker_id-1]] #faster than looping on all bodies #if (O.bodies[O.subD.subdomains[worker_id-1]].shape.ids+[O.subD.subdomains[worker_id-1]] != ids): #print("______________INCONSISTENCY!______________",ids," vs. ",O.bodies[O.subD.subdomains[worker_id-1]].shape.ids+[O.subD.subdomains[worker_id-1]]) shift = dspl[worker_id]; if (worker_id != numThreads-1): shift_plus_one = dspl[worker_id+1]; else: shift_plus_one = len(dat); O.subD.setStateBoundsValuesFromIds(ids,dat[shift: shift_plus_one]); reboundRemoteBodies(ids) # turn mpi engines off sendRecvStatesRunner.dead = isendRecvForcesRunner.dead = waitForcesRunner.dead = collisionChecker.dead = True O.splitted=False #collider.doSort = True if (AUTO_COLOR): colorDomains() if rank==0: O.engines = O.initialEngines
[docs]def splitScene(): ''' Split a monolithic scene into distributed scenes on threads. Precondition: the bodies have subdomain no. set in user script ''' if not COPY_MIRROR_BODIES_WHEN_COLLIDE: mprint("COPY_MIRROR_BODIES_WHEN_COLLIDE=False is not supported") #if not ERASE_REMOTE: mprint("ERASE_REMOTE=False is not supported") if not O.splittedOnce: O.initialEngines = O.engines if DOMAIN_DECOMPOSITION: #if not already partitionned by the user we partition here if rank == 0: import yade.bisectionDecomposition as dd decomposition = dd.decompBodiesSerial(comm) decomposition.partitionDomain(fibreList) maxid = len(O.bodies)-1 if DISTRIBUTED_INSERT: #find max id before inserting subdomains maxid = timing_comm.allreduce("splitScene",maxid,op=MPI.MAX) wprint("Splitting with maxId=",maxid) if rank == 0 or DISTRIBUTED_INSERT: subdomains=[] #list subdomains by body ids #insert "meta"-bodies for k in range(1,numThreads): domainBody=Body(shape=Subdomain(ids=[b.id for b in O.bodies if b.subdomain==k]),subdomain=k) #note: not clear yet how shape.subDomainIndex and body.subdomain should interact, currently equal values domainBody.isSubdomain=True if rank==k: O._sceneObj.subD=domainBody.shape subdomains.append(O.bodies.insertAtId(domainBody,maxid+k)) if rank==0: O._sceneObj.subD = Subdomain() # make sure it's initialized here O.subD = O._sceneObj.subD O.subD.subdomains = subdomains subD= O.subD #alias subD.comm=comm #make sure the c++ uses the merged intracommunicator masterBodies = [b.id for b in O.bodies if b.subdomain==0] #for VTKRecorderParallel, easier to loop through the owned bodies. subD.setIDstoSubdomain(masterBodies) #tell the collider how to handle this new thing collider = typedEngine("InsertionSortCollider") if FLUID_COUPLING: collider.boundDispatcher.functors = collider.boundDispatcher.functors+[Bo1_FluidDomainBbox_Aabb()] collider.boundDispatcher.functors=collider.boundDispatcher.functors+[Bo1_Subdomain_Aabb()] collider.targetInterv=0 collider.keepListsShort=True # probably not needed, O.bodies.insertAtId should turn it on automaticaly O.bodies.useRedirection=True # idem O.bodies.allowRedirection=False #BEGIN Garbage (should go to some init(), usually done in collider.__call__() but in the mpi case we want to collider.boundDispatcher.__call__() before collider.__call__() collider.boundDispatcher.sweepDist=collider.verletDist; collider.boundDispatcher.minSweepDistFactor=collider.minSweepDistFactor; collider.boundDispatcher.targetInterv=collider.targetInterv; collider.boundDispatcher.updatingDispFactor=collider.updatingDispFactor; #END Garbage if not DISTRIBUTED_INSERT: #we send scene from master to all workers sceneAsString= O.sceneToString() if rank==0 else None sceneAsString=timing_comm.bcast("splitScene",sceneAsString,root=0) if rank > 0: O.stringToScene(sceneAsString) #receive a scene pre-processed by master (i.e. with appropriate body.subdomain's) # as long as subD.subdomains isn't serialized we need to rebuild it here since it's lost domainBody=None subdomains=[] #list of subdomains by body id for b in O.bodies: if b.isSubdomain: subdomains.append(b.id) if b.subdomain==rank: domainBody=b if domainBody==None: wprint("SUBDOMAIN NOT FOUND FOR RANK=",rank) O._sceneObj.subD = domainBody.shape O.subD = O._sceneObj.subD O.subD.subdomains = subdomains subD = O.subD if FLUID_COUPLING: fluidCoupling = typedEngine("FoamCoupling") fluidCoupling.comm = comm fluidCoupling.setIdList(fluidBodies) fluidCoupling.couplingModeParallel = True O._sceneObj.subdomain = rank O.subD.comm=comm #make sure the c++ uses the merged intracommunicator O.subD.init() wprint("to parallel collide") parallelCollide() wprint("end parallel collide") # insert states communicator after newton idx = O.engines.index(typedEngine("NewtonIntegrator")) O.engines=O.engines[:idx+1]+[PyRunner(iterPeriod=1,initRun=True,command="sys.modules['yade.mpy'].sendRecvStates(); ",label="sendRecvStatesRunner")]+O.engines[idx+1:] # insert force communicator before Newton O.engines=O.engines[:idx]+[PyRunner(iterPeriod=1,initRun=True,command="sys.modules['yade.mpy'].isendRecvForces()",label="isendRecvForcesRunner")]+O.engines[idx:] # append engine waiting until forces are effectively sent to master O.engines=O.engines+[PyRunner(iterPeriod=1,initRun=True,command="pass",label="waitForcesRunner")] O.engines=O.engines+[PyRunner(iterPeriod=1,initRun=True,command="if sys.modules['yade.mpy'].checkAndCollide(): O.pause();",label="collisionChecker")] O.splittedOnce = True O.splittedEngines = O.engines else: if (DOMAIN_DECOMPOSITION and RESET_SUBDOMAINS_WHEN_COLLIDE): if rank == 0: import yade.bisectionDecomposition as dd decomposition = dd.decompBodiesSerial(comm) decomposition.partitionDomain() O.subD.splitBodiesToWorkers(RESET_SUBDOMAINS_WHEN_COLLIDE) parallelCollide() if rank == 0 : O.engines = O.splittedEngines O.interactions.clear() unboundRemoteBodies() if (ERASE_REMOTE and ERASE_REMOTE_MASTER): eraseRemote() sendRecvStatesRunner.dead = isendRecvForcesRunner.dead = waitForcesRunner.dead = collisionChecker.dead = False if MINIMAL_INTERSECTIONS: shrinkIntersections() O.splitted = True
[docs]def updateMirrorOwners(): mirrorInts=O.subD.mirrorIntersections for kk in range(len(mirrorInts)): for id in mirrorInts[kk]: if not O.bodies[id]: continue elif O.bodies[id].subdomain!=kk: O.bodies[id].subdomain = kk
[docs]def updateAllIntersections(): subD=O.subD subD.intersections=genLocalIntersections(subD.subdomains) if USE_CPP_INTERS: #to be used in case of possible issue with older mpi4py versions, this is slightly faster. subD.getMirrorIntrs() else: #update mirror intersections so we know message sizes in advance subD.mirrorIntersections=[[] for n in range(numThreads)] if rank==0:#master domain for worker in range(1,numThreads):#FIXME: we actually don't need so much data since at this stage the states are unchanged and the list is used to re-bound intersecting bodies, this is only done in the initialization phase, though #wprint("sending mirror intersections to "+str(worker)+" ("+str(len(subD.intersections[worker]))+" bodies), "+str(subD.intersections[worker])) m = O.intrsctToBytes(subD,worker,False) if SEND_BYTEARRAYS else subD.intersections[worker]; timing_comm.send("sendIntersections",m, dest=worker, tag=_MIRROR_INTERSECTIONS_) else: # from master b_ids=comm.recv(source=0, tag=_MIRROR_INTERSECTIONS_) wprint("Received mirrors from master: ",len(b_ids)) #FIXME: we are assuming that Body::id_t is 4 bytes here, not that portable... numInts0= int(len(b_ids)/4) if SEND_BYTEARRAYS else len(b_ids) #ints = 4 bytes if numInts0>0: if SEND_BYTEARRAYS: O.bufferFromIntrsct(subD,0,numInts0,True)[:]=b_ids b_ids=np.frombuffer(b_ids,dtype=np.int32) else: subD.mirrorIntersections= [b_ids]+subD.mirrorIntersections[1:] #reboundRemoteBodies(b_ids) # since interaction with 0-bodies couldn't be detected before, mirror intersections from master will # tell if we need to wait messages from master (and this is declared via intersections) if not 0 in subD.intersections[rank]: temp=subD.intersections[rank] temp+=[0] subD.intersections=subD.intersections[:rank]+[temp]+subD.intersections[rank+1:] else: if not O.splittedOnce: mprint("0 already in intersections (should not happen)") reqs=[] for worker in subD.intersections[rank]: if worker==0: continue #we do not send positions to master, only forces #wprint("sending "+str(len(subD.intersections[worker]))+" states to "+str(worker)) m = O.intrsctToBytes(subD,worker,False) if SEND_BYTEARRAYS else subD.intersections[worker]; reqs.append(comm.isend(m, dest=worker, tag=_MIRROR_INTERSECTIONS_)) for worker in subD.intersections[rank]: if worker==0: continue intrs=timing_comm.recv("recvAllIntersections",source=worker, tag=_MIRROR_INTERSECTIONS_) if SEND_BYTEARRAYS: #wprint("Received mirrors from: ", worker, " : ",int(len(intrs)/4)) O.bufferFromIntrsct(subD,worker,int(len(intrs)/4),True)[:]=intrs intrs=np.frombuffer(intrs,dtype=np.int32) else: subD.mirrorIntersections = subD.mirrorIntersections[0:worker]+[intrs]+subD.mirrorIntersections[worker+1:] for req in reqs: req.wait() updateMirrorOwners()
bodiesToImport=[]
[docs]def parallelCollide(): global bodiesToImport subD=O.subD start = time.time() if (not O.splitted): unboundRemoteBodies() eraseRemote() collider.boundDispatcher.__call__() updateDomainBounds(subD.subdomains) #triggers communications collider.__call__() #see [1] unboundRemoteBodies() #in splitted stage we exploit bounds to detect bodies which are no longer part of intersections (they will be left with no bounds after what follows) updateAllIntersections() #triggers communications if rank!=0: for l in subD.mirrorIntersections: if len(l)>0: reboundRemoteBodies(l) if(ERASE_REMOTE): eraseRemote() # erase the bodies which still have no bounds """ " NOTE: FK, what to do here: " 1- all threads loop on reqs, i.e the intersecting subdomains of the current subdomain. " 2- during this loop, check whether the current subdomain needs bodies from the intersecting ones " 3- in all cases, isend the ids needed, (if no ids needed send empty array) " 3.1- build a ranks list "requestedSomethingFrom" to loop later on to receive data " 4- loop again on reqs to get the ids needed by other subdomains (with blocking recv as we used isend) " 5- if the data recved is empty (nothing requested), do nothing. Else isend the bodies (c++) " 6- loop on "requestedSomethingFrom" ranks and recv the bodies (blocking, c++, using MPI_Probe to know the message size) " 7- comm.barrier(), just in case """ if(COPY_MIRROR_BODIES_WHEN_COLLIDE or MERGE_W_INTERACTIONS): requestedSomethingFrom=[] bodiesToImport=[[] for worker in range(numThreads)] sent=[] if rank>0: #master doesn't need bodies #reqs=subD.intersections[rank] for worker in subD.intersections[rank]: #worker=req[0] #if(worker==0):continue for mirrorBodyId in subD.mirrorIntersections[worker]: if O.bodies[mirrorBodyId]==None: bodiesToImport[worker]+=[mirrorBodyId] if(len(bodiesToImport[worker])>0): requestedSomethingFrom.append(worker) wprint("I request ids: ",len(bodiesToImport[worker]), " from ",worker) sent.append(comm.isend(bodiesToImport[worker], worker, tag=_MIRROR_INTERSECTIONS_)) wprint("will wait requests from ",subD.intersections[rank]) for worker in subD.intersections[rank]: if worker!=0: wprint("waiting requests from ",worker) requestedIds=timing_comm.recv("parallelCollide",source=worker, tag=_MIRROR_INTERSECTIONS_) wprint("requested: ",len(requestedIds),"from ",worker) if(len(requestedIds)>0): wprint("will now send ",len(requestedIds)," to ",worker) subD.sendBodies(worker,requestedIds) for worker in requestedSomethingFrom: subD.receiveBodies(worker) for s in sent: s.wait() subD.completeSendBodies(); if not collider.keepListsShort: collider.doSort = True collider.__call__() collider.execTime+=int((time.time()-start)*1e9) collider.execCount+=1 try: collisionChecker.execTime-=int((time.time()-start)*1e9) except: pass #maxVelocitySq is normally reset in NewtonIntegrator in the same iteration as bound dispatching, since Newton will not run before next iter in our case we force that value to avoid another collision detection at next step typedEngine("NewtonIntegrator").maxVelocitySq=0.5
[docs]def eraseRemote(): if rank>0 or ERASE_REMOTE_MASTER: # suppress external bodies from scene #numBodies = len(O.bodies) #for id in range(numBodies): #mprint("will erase ",[b.id for b in O.bodies if (not b.bounded and b.subdomain!=rank)]) for b in O.bodies: if not b.bounded and b.subdomain != rank: connected = False #a gridNode could be needed as part of interacting facet/connection even if not overlaping a specific subdomain. Assume connections are always bounded for now, we thus only check nodes. if isinstance(b.shape,GridNode): for f in b.shape.getPFacets(): if f.bounded: connected = True for c in b.shape.getConnections(): if c.bounded: connected = True if not connected: O.bodies.erase(b.id)
##### RUN MPI #########
[docs]def mpirun(nSteps,np=None,withMerge=False): ''' Parallel version of O.run() using MPI domain decomposition. Parameters nSteps : The numer of steps to compute np : number of mpi workers (master+subdomains), if=1 the function fallback to O.run() withMerge : wether subdomains should be merged into master at the end of the run (default False). If True the scene in the master process is exactly in the same state as after O.run(nSteps,True). The merge can be time consumming, it is recommended to activate only if post-processing or other similar tasks require it. ''' if comm==None: configure() if np==None: np=numThreads if(np==1): mprint("single-core, fall back to O.run()") O.run(nSteps,True) return if (np!=numThreads): if numThreads >1: mprint("warning: it is unsafe to change numThreads in consecutive executions of mpy.initialize/mpirun. In general it needs explicit disconnect() and rebuilding a scene in between.") initialize(np) #this will set numThreads if(rank==0 and waitingCommands): for w in range(1,numThreads): comm.send("yade.mpy.mpirun(nSteps="+str(nSteps)+",withMerge="+str(withMerge)+")",dest=w,tag=_MASTER_COMMAND_) wprint("mpirun command sent to ",w) if FLUID_COUPLING: fluidCoupling = typedEngine("FoamCoupling") fluidCoupling.comm = comm fluidCoupling.getFluidDomainBbox() #triggers the communication between yade procs and Yales2/openfoam procs, get's fluid domain bounding boxes from all fluid procs. # split if needed initStep = O.iter if not O.splitted: wprint("splitting") splitScene() wprint("splitted") O.timingEnabled=YADE_TIMING #turn it ON/OFF # run iterations if not (MERGE_SPLIT): O.run(nSteps,True) if withMerge: mergeScene() #will be useful to see evolution in QGLViewer, for instance else: #merge/split or body_copy for each collider update collisionChecker.dead=True while (O.iter-initStep)<nSteps: O.step() if checkAndCollide(): mergeScene() splitScene() mergeScene() # report performance if YADE_TIMING and (rank<=MAX_RANK_OUTPUT or rank>=(numThreads-MAX_RANK_OUTPUT)): timing_comm.print_all() from yade import timing time.sleep((numThreads-rank)*0.002) #avoid mixing the final output, timing.stats() is independent of the sleep mprint( "##### Worker "+str(rank)+" ######") timing.stats() #specific numbers for -n4 and gabion.py
####################################### #######  Bodies re-allocation ######## #######################################
[docs]def bodyErase(ids): ''' The parallel version of O.bodies.erase(id), should be called collectively else the distributed scenes become inconsistent with each other (even the subdomains which don't have 'id' can call safely). For performance, better call on a list: bodyErase([i,j,k]). ''' if not isinstance(ids,list): ids = [ids] erased=0 if not numThreads or numThreads<=1:# non-parallel case, fall back but handle lists for i in ids: erased += O.bodies.erase(i) return erased if(rank==0 and waitingCommands): # make it collective erased = sendCommand("slaves","bodyErase("+str(ids)+")",True) if hasattr(O.subD,"fullIntersections") and O.subD.fullIntersections!=None: #un-shrink O.subD.intersections = O.subD.fullIntersections O.subD.mirrorIntersections = O.subD.fullMirrorIntersections O.subD.fullMirrorIntersections=O.subD.fullIntersections=None intersections_ = O.subD.intersections mirrorIntersections_ = O.subD.mirrorIntersections ori = len(O.subD.ids) localIds=O.subD.ids for id in ids: if O.bodies[id]==None: #mprint("not erasing None"); continue b=O.bodies[id] owner = b.subdomain if rank==owner: if not id in localIds: mprint("a nasty bug somewhere") localIds.remove(id) removedFrom=[] for ii in O.interactions.withBodyAll(id): otherID = ii.id1 if ii.id2==id else ii.id2 otherSD = O.bodies[otherID].subdomain if otherSD != rank and not otherSD in removedFrom: #mprint("removing",id," from intersections with",O.bodies[otherID].subdomain) if id in intersections_[otherSD]: intersections_[otherSD].remove(id) removedFrom.append(otherSD) else: if id in mirrorIntersections_[owner]: mirrorIntersections_[owner].remove(id) #if hasattr(O.subD,"fullMirrorIntersections") and id in O.subD.fullMirrorIntersections[owner]: #O.subD.fullMirrorIntersections[owner].remove(id) # FIXME: same as above in the symmetric case #if len(O.subD.intersections[O.bodies[otherID].subdomain])==0: # if it was last interacting body then the two subdomains don't interact any more #O.subD.intersections[rank].remove(O.bodies[otherID].subdomain) erased+=O.bodies.erase(id) # write result back into the subdomain # other subdomain's ids are locally deprecated (should be safe since they are not used) O.subD.ids = localIds O.subD.intersections = intersections_ O.subD.mirrorIntersections = mirrorIntersections_ return erased
[docs]def runOnSynchronouslPairs(workers,command): ''' Locally (from one worker POV), this function runs interactive mpi tasks defined by 'command' on a list of other workers (typically the list of interacting subdomains). Overall, peer-to-peer connexions are established so so that 'command' is executed symmetrically and simultaneously on both sides of each worker pair. I.e. if worker "i" executes "command" with argument "j" (index of another worker), then by design "j" will execute the same thing with argument "i" *simultaneously*. In many cases a similar series of data exchanges can be obtained more simply (and fastly) with asynchronous irecv+send like below. for w in workers: m=comm.irecv(w) comm.send(data,dest=w) The above only works if the messages are all known in advance locally, before any communication. If the interaction with workers[1] depends on the result of a previous interaction with workers[0] OTOH, it needs synchronous execution, hence this function. Synchronicity is also required if more than one blocking call is present in 'command', else an obvious deadlock as if 'irecv' was replaced by 'recv' in that naive loop. Both cases occur with the 'medianFilter' algorithm, hence why we need this synchronous method. In this function pair connexions are established by the workers in a non-supervized and non-deterministic manner. Each time an interactive communication (i,j) is established 'command' is executed simultaneously by i and j. It is guaranted that all possible pairs are visited. The function can be used for all-to-all operations (N^2 pairs), but more interestingly it works with workers=intersections[rank] (O(N) pairs). It can be tested with the dummy funtion 'pairOp': runOnSynchronouslPairs(range(numThreads),pairOp) command: a function taking index of another worker as argument, can include blocking communications with the other worker since runOnSynchronouslPairs guarantee that the other worker will be running the command symmetrically. ''' global t1 t1=time.time() workersTemp = list(workers) #we will remove elements from this copy as we proceed if rank in workersTemp: workersTemp.remove(rank) #don't talk to yourself #skip master (master is useless typically if we are reallocating bodies, we skip it here although the function would work with master to) if True: if 0 in workersTemp: workersTemp.remove(0) #don't talk to master if rank==0: return #don't talk if you are master # now call the workers one by one until one of them reacts s=MPI.Status(); sentTo=-1 #last worker to which connexion request was sent while len(workersTemp)>0: other=workersTemp[-1] connected=False if rank>other and other!=sentTo: rs = comm.issend(None,dest=other,tag=_GET_CONNEXION_) sentTo = other while not (connected): if comm.Iprobe(source=MPI.ANY_SOURCE, status=s,tag=_GET_CONNEXION_): data = timing_comm.recv("runOnSynchronouslPairs",None,s.source,tag=_GET_CONNEXION_) talkTo = s.source connected = True elif rank>other: if rs.test()[0]: talkTo = other rs.wait() connected=True command(talkTo) workersTemp.remove(talkTo)
#mprint("TOTAL TIME(",rank,"):",time.time()-t1) # a dummy test function for runOnSynchronouslPairs(range(numThreads,pairOp)
[docs]def pairOp(talkTo): global t1 # send/recv data message="haha" comm.isend(message,talkTo,tag=_PAIR_OP_) feedback = timming_comm.recv("pair_op",None,talkTo,tag=_PAIR_OP_) # crunch feedback and numbers... time.sleep(0.01) # send/recv result comm.isend(message,talkTo,tag=_PAIR_OP_) feedback = timing_comm.recv("pair_op",None,talkTo,tag=_PAIR_OP_) print("(",rank,talkTo,") done in",time.time()-t1,"s")
[docs]def migrateBodies(ids,origin,destination): ''' Reassign bodies from origin to destination. The function has to be called by both origin (send) and destination (recv). Note: subD.completeSendBodies() will have to be called after a series of reassignement since subD.sendBodies() is non-blocking ''' ts = time.time() if rank==origin: if USE_CPP_REALLOC: O.subD.migrateBodiesSend(ids, destination) else: thisSubD = O.subD.subdomains[rank-1] for id in ids: if not O.bodies[id]: mprint("reassignBodies failed,",id," is not in subdomain ",rank) O.bodies[id].subdomain = destination createInteraction(thisSubD,id,virtualI=True) # link translated body to subdomain, since there is initially no interaction with local bodies #for k in O.subD.intersections[rank]: #if k==0: continue #if id in O.subD.intersections[k]: #O.subD.intersections[k].remove(id) # so we don't send the same body to multiple domains... O.subD.sendBodies(destination,ids) elif rank==destination: O.subD.receiveBodies(origin) te = time.time()
#mprint("time in migrateBodies--> ", te-ts, " rank = ", rank)
[docs]def projectedBounds(i,j): ''' Returns sorted list of projections of bounds on a given axis, with bounds taken in i->j and j->i intersections ''' useAABB=False #using center of subdomain AABB is a bit unstable since the movement of on single body can change it greatly if (useAABB): pt1 = 0.5*(O.bodies[O.subD.subdomains[i-1]].bound.min+O.bodies[O.subD.subdomains[i-1]].bound.max) pt2 = 0.5*(O.bodies[O.subD.subdomains[j-1]].bound.min+O.bodies[O.subD.subdomains[j-1]].bound.max) else: #use center of mass pt1 = O.subD._centers_of_mass[i] pt2 = O.subD._centers_of_mass[j] axis=pt2-pt1 axis.normalize() pos = [[O.subD.boundOnAxis(O.bodies[k].bound,axis,True),i,k] for k in O.subD.intersections[j]]+[[O.subD.boundOnAxis(O.bodies[k].bound,axis,False),j,k] for k in O.subD.mirrorIntersections[j]] pos.sort(key= lambda x: x[0]) return pos
[docs]def medianFilter(i,j,giveAway): ''' Returns bodies in "i" to be assigned to "j" based on median split between the center points of subdomain's AABBs If giveAway!=0, positive or negative, "i" will give/acquire this number to "j" with nothing in return (for load balancing purposes) ''' bodiesToSend=[] if USE_CPP_REALLOC: useAABB = False; otherSubDCM = O.subD._centers_of_mass[j] subDCM = O.subD._centers_of_mass[i] bodiesToSend= O.subD.medianFilterCPP(j, otherSubDCM, subDCM, giveAway, useAABB) else: pos = projectedBounds(i,j) firstJ, lastI = len(pos),0 for n in range(len(pos)): if pos[n][1] == i: lastI = n elif n<firstJ: firstJ = n finalSize = min( max(0, len(O.subD.intersections[j]) - giveAway) , len(pos)) if finalSize > lastI: finalSize=lastI+1 if finalSize < firstJ: finalSize=firstJ+1 bodiesToSend= [x[2] for x in pos[finalSize:] if x[1]==i] #bodiesToRecv= [x[2] for x in pos[:finalSize] if x[1]==j] #for debugging only return bodiesToSend
REALLOCATE_FILTER=medianFilter #that's currently default and only option
[docs]def reallocateBodiesToSubdomains(_filter=REALLOCATE_FILTER,blocking=True): ''' Re-assign bodies to subdomains based on '_filter' argument. Requirement: '_filter' is a function taking ranks of origin and destination and returning the list of bodies (by index) to be moved. That's where the decomposition strategy is defined. See example medianFilter (used by default). This function must be called in parallel, hence if ran interactively the command needs to be sent explicitely: mp.sendCommand("all","reallocateBodiesToSubdomains(medianFilter)",True) ''' O.subD._centers_of_mass=[Vector3(0,0,0) for k in range(numThreads)] O.subD._centers_of_mass[rank]=O.subD.centerOfMass() if blocking: # the filter will be applied sequentially to each other domain. It can include blocking communications in subdomain pairs _functor = lambda x : reallocateBodiesPairWiseBlocking(_filter,x) runOnSynchronouslPairs(O.subD.intersections[rank],_functor) else: # non-blocking method, migrated bodies are decided unilateraly by each subdomain # doesn't work with medianFilter if rank>0: for worker in O.subD.intersections[rank]: if worker==0: continue candidates = _filter(rank,worker) wprint("sending to ",worker,": ",len(candidates)) migrateBodies(candidates,rank,worker) #send migrateBodies(None,worker,rank) #recv O.subD.completeSendBodies() ts = time.time() if USE_CPP_REALLOC: O.subD.updateLocalIds(ERASE_REMOTE_MASTER) if not ERASE_REMOTE_MASTER: if rank == 0 : if (AUTO_COLOR) : colorDomains() else: O.subD.ids = [b.id for b in O.bodies if (b.subdomain==rank and not b.isSubdomain)] #update local ids if not ERASE_REMOTE_MASTER: # update remote ids in master if rank>0: req = comm.isend(O.subD.ids,dest=0,tag=_ASSIGNED_IDS_) req.wait() else: #master will update subdomains for correct display (besides, keeping 'ids' updated for remote subdomains may not be a strict requirement) for k in range(1,numThreads): ids=comm.recv(source=k,tag=_ASSIGNED_IDS_) O.bodies[O.subD.subdomains[k-1]].shape.ids=ids for i in ids: O.bodies[i].subdomain=k if (AUTO_COLOR): colorDomains() # update intersections and mirror #mprint("updating all intersections") updateAllIntersections() #triggers communication
[docs]def reallocateBodiesPairWiseBlocking(_filter,otherDomain): ''' Re-assign bodies from/to otherDomain based on '_filter' argument. Requirement: '_filter' is a function taking ranks of origin and destination and returning the list of bodies (by index) to be moved. That's where the decomposition strategy is defined. See example medianFilter (used by default). ''' #if rank==0: return ts = time.time() if True: #clean intersections, remove bodies already moved to other domain if USE_CPP_REALLOC: O.subD.cleanIntersections(otherDomain) else: ints = [ii for ii in O.subD.intersections[otherDomain] if O.bodies[ii].subdomain==rank] #make sure we don't send ids of already moved bodies O.subD.intersections=O.subD.intersections[:otherDomain]+[ints]+O.subD.intersections[otherDomain+1:] te = time.time() numSubscribedHere = len(O.subD.ids) #total weigth of this domain, for load balancing after some erase req = comm.isend([O.subD.intersections[otherDomain],O.subD._centers_of_mass[rank],numSubscribedHere],dest=otherDomain,tag=_MIRROR_INTERSECTIONS_) newMirror = comm.recv(source=otherDomain,tag=_MIRROR_INTERSECTIONS_) req.wait() ts = time.time() if USE_CPP_REALLOC: O.subD.updateNewMirrorIntrs(otherDomain, newMirror[0]) else: O.subD.mirrorIntersections=O.subD.mirrorIntersections[:otherDomain]+[newMirror[0]]+O.subD.mirrorIntersections[otherDomain+1:] te = time.time() O.subD._centers_of_mass[otherDomain]=newMirror[1] numSubscribedOther = newMirror[2] giveAway = int(0.5* (numSubscribedHere-numSubscribedOther)) candidates = _filter(rank,otherDomain,giveAway) migrateBodies(candidates,rank,otherDomain) #send migrateBodies(None,otherDomain,rank) #recv