#!/usr/bin/env python # using bmiwrapper to communicate with 2 xbeach instances # using mpi4py to create two processes, each running a # bmi-aware xbeach # this program handles the "Nested-Island" case from mpi4py import MPI import sys import os from bmi.wrapper import BMIWrapper import numpy as np from interpol import interpoldln,interpolkd binroot = "/home/willem/projecten/parallel/xbeachbin" #xbeachso = binroot + "/willem/gnu/serial/no/both/debug/lib/libxbeach.so" xbeachso = binroot + "/willem/gnu/serial/no/both/opt/lib/libxbeach.so" exchange_interval = 1 # number of seconds between exchange maxtime = 1200 workdirs = ["/home/willem/tmp/xbeach/Nested_Island/inner", "/home/willem/tmp/xbeach/Nested_Island/outer"] using2 = True # use island2{grd,dep} let op: aanpassen in params.txt #using2 = False # use island{grd,dep} let op: aanpassen in params.txt in_inner = [True,False] excvars = ['zb','zs','uu','vv','uv','vu','ue','ve','ee','rr'] def parentcode(): comm = MPI.COMM_SELF.Spawn("wrap", args = 'wrap', maxprocs = 2) class xbchild: def __init__(self,xbeachso,comm,excvars,tinterval,tstop,inner): self.w = BMIWrapper(engine=xbeachso) self.comm = comm self.excvars = excvars self.tinterval = tinterval self.tstop = tstop self.inner = inner self.outer = not inner self.w.initialize("is-always-params.txt") x = self.w.get_var("x").copy() y = self.w.get_var("y").copy() print "x.shape",x.shape print "y.shape",y.shape print "xvalues",x print "yvalues",y self.nx, self.ny = x.shape self.nx -= 1 self.ny -= 1 print "nx,ny:",self.nx,self.ny self.xy = np.asarray((x.ravel(), y.ravel())).T print "xy.shape",self.xy.shape with open('xy','w') as f: for item in self.xy: print >>f,item[0],item[1] # coordinates of border: if self.inner: # border of inner # found by trial and error if using2: self.bxy = np.asarray(( # this is for island2.{grd,dep} np.concatenate([x[1:,0],x[1:,-1],x[1:,1],x[1:,-2]]), np.concatenate([y[1:,0],y[1:,-1],y[1:,1],y[1:,-2]]) )).T else: self.bxy = np.asarray(( # this is for island.{grd,dep} np.concatenate([x[:,0],x[:,-1],x[:,1],x[:,-2]]), np.concatenate([y[:,0],y[:,-1],y[:,1],y[:,-2]]) )).T if self.outer: # border of outer, is not used # found by guessing self.bxy = np.asarray(( np.concatenate([x[0,:], x[-1,:], x[1:-1,0],x[1:-1,-1], x[1,1:-1],x[-2,1:-1],x[2:-2,1],x[2:-2,-2]]), np.concatenate([y[0,:], y[-1,:], y[1:-1,0],y[1:-1,-1], y[1,1:-1],y[-2,1:-1],y[2:-2,1],y[2:-2,-2]]) )).T print "bxy.shape",self.bxy.shape with open("bxy","w") as f: for item in self.bxy: print >>f, item[0],item[1] if self.comm.rank ==1: self.othercomm = 0 else: self.othercomm = 1 # inner: # otherxy is all xy values of outer # otherbxy is all border xy values of outer (not used) # # outer: # otherxy is all xy values of inner # otherbxy is all border xy values of inner # if self.comm.rank == 0: self.otherxy = self.comm.recv(source=self.othercomm) self.otherbxy = self.comm.recv(source=self.othercomm) self.comm.send(self.xy,self.othercomm) self.comm.send(self.bxy,self.othercomm) else: self.comm.send(self.xy,self.othercomm) self.comm.send(self.bxy,self.othercomm) self.otherxy = self.comm.recv(source=self.othercomm) self.otherbxy = self.comm.recv(source=self.othercomm) # something like this to keep only the overlapping # coordinates self.oxy: if self.outer: distx = abs(self.xy[0,0]-self.xy[1,0]) disty = abs(self.xy[0,1]-self.xy[1,1]) maxdist = max(distx,disty) maxdist = np.sqrt(distx**2+disty**2) kd = interpolkd(self.otherxy, np.ones(len(self.otherxy)), fill_value = -10, max_dist = maxdist) z = kd(self.xy) print "z.shape:",z.shape # can also be done with interpoldln: #tri = interpoldln(self.otherxy, # np.ones(len(self.otherxy)), # fill_value = -10) #z = tri(self.xy) # but no: in the island example: also the middle area is filled in self.oxy_indices = np.array(np.where(z>0)).ravel() self.oxy= self.xy[self.oxy_indices] print "oxy.shape:",self.oxy.shape print "indices.shape:",self.oxy_indices.shape with open("oxy","w") as f: for item in self.oxy: print >>f,item[0],item[1] # create kd tree using xy values from inner (otherxy) self.okd = interpolkd(self.otherxy) if self.inner: self.okd = interpolkd(self.otherxy) # create kd tree for all xy values of other process #self.kd = interpolkd(self.otherxy,np.zeros(len(self.otherxy))) def run(self): curtime = 0 while True: extime = curtime + self.tinterval print 'extime before allreduce:',extime extime = self.comm.allreduce(extime,MPI.MAX) print 'extime after allreduce:',extime while curtime < extime: rc = self.w.update(exchange_interval) curtime = self.w.get_current_time() print 'curtime after update:',curtime print 'exchange time!' for a in self.excvars: if self.inner: self.exchange_inner(a) if self.outer: self.exchange_outer(a) if extime > self.tstop: print 'time to stop' self.w.finalize() break def exchange_inner(self,a): # receive innerborder as calculated by xbeach.outer # send all inner contents to outer print "inner handling",a # get all inner values from xbeach.inner var = self.w.get_var(a).copy() ndims = len(var.shape) print "var.shape:",var.shape # get all values from other process othervar = self.comm.recv(source=self.othercomm,tag=1) print "othervar.shape:",othervar.shape # send all inner contents to other process self.comm.send(var,self.othercomm,tag=2) # use 3-dim arrays, if (other)var is 2d, make it 3-d if ndims == 2: othervar.resize((1,)+othervar.shape) var.resize((1,)+var.shape) for i in range(len(othervar)): # len(a) is teh size of the forstdimension of a self.okd.setz(othervar[i].ravel()) z = self.okd(self.bxy) # test if border is indeed changed: if False: if a == "zb": z = z*0+100 print "z.shape after okd:",z.shape if using2: # using island2{grd,dep} p = 0 l = self.nx var[i,1:,0] = z[p:p+l] p += l var[i,1:,-1] = z[p:p+l] p += l var[i,1:,1] = z[p:p+l] p += l var[i,1:,-2] = z[p:p+l] var[i,0,0] = var[i,1,0] var[i,0,-1] = var[i,1,-1] var[i,0,-2] = var[i,1,-2] else: # using island{grd,dep} p = 0 l = self.nx+1 var[i,:,0] = z[p:p+l] p += l var[i,:,-1] = z[p:p+l] p += l var[i,:,1] = z[p:p+l] p += l var[i,:,-2] = z[p:p+l] if ndims == 2: # resize to original: var.resize(var.shape[1:]) self.w.set_var(a,var) def exchange_outer(self,a): # send innerborder # receive all inner contents print "outer handling",a var = self.w.get_var(a).copy() ndims = len(var.shape) print "var.shape:",var.shape self.comm.send(var,self.othercomm,tag=1) # would be more efficient to send only the border values, using otherbxy othervar = self.comm.recv(source=self.othercomm,tag=2) print "outer: othervar.shape:",othervar.shape print "outer: indices.shape:",self.oxy_indices.shape # apply received data and setvar it if ndims == 2: othervar.resize((1,)+othervar.shape) var.resize((1,)+var.shape) for i in range(len(othervar)): var1 = var[i].ravel() self.okd.setz(othervar[i].ravel()) z = self.okd(self.oxy) print "z.shape after okd:",z.shape # test to see the overlap if False: if a == "zb": z = z*0 + 100 var1[self.oxy_indices] = z #var1.resize(self.nx+1,self.ny+1) var[i]=var1.reshape(self.nx+1,self.ny+1) self.w.set_var(a,var) def childcode(): print 'childcode: size,rank:',comm.size,comm.rank os.chdir(workdirs[comm.rank]) print "Working in",os.getcwd() xb = xbchild(xbeachso,comm,excvars, exchange_interval,maxtime,in_inner[comm.rank]) xb.run() return comm = MPI.COMM_WORLD fcomm = np.array([comm.py2f()]) size = comm.size if comm.rank == 0: print 'MPI initialized: size = ',size if (MPI.Comm.Get_parent() == MPI.COMM_NULL): parentcode() else: childcode()