diff --git a/Cassiopee/Connector/Connector/IBM.py b/Cassiopee/Connector/Connector/IBM.py index 6dbcf1645..9147a279c 100644 --- a/Cassiopee/Connector/Connector/IBM.py +++ b/Cassiopee/Connector/Connector/IBM.py @@ -65,6 +65,33 @@ def getIBCDName(proposedName): (ibcname,__IBCNameServer__)=C.getUniqueName(proposedName, __IBCNameServer__) return ibcname + +#==================================================================================== +#Add .Solver#Define with dirx,diry, & dirz to the base of the tboneover. tboneover is the +#PyTree that defines the region in space wherein a one over n coarsening will be pursued +#during the automatic cartesian grid generator of FastIBC. +#IN: FileName: name of the file. tboneover is read in this function. +#IN: oneOver: list of list of dirx,diry,dirz for each base in tboneover. E.g. oneOver=[[1,1,2],[1,2,1],[2,1,1]] +# for a tboneover with 3 bases where the 1st base has dirx=1, diry=1, & dirz=2 +# 2nd base has dirx=1, diry=2, & dirz=1 +# 3rd base has dirx=2, diry=1, & dirz=1 +#OUT: Nothing. Rewrite tboneover with the same FileName as that original used +##NOTE # 1: To be run SEQUENTIALLY ONLY. This is ok as we are dealing with a surface geometry which tend to be +## relatively small. +##NOTE # 2: Generation of tboneover is similar to that used for tbox. +def _addOneOverLocally(FileName,oneOver): + count = 0 + t_local = C.convertFile2PyTree(FileName) + for b in Internal.getBases(t_local): + Internal._createUniqueChild(b, '.Solver#define', 'UserDefinedData_t') + n = Internal.getNodeFromName1(b, '.Solver#define') + Internal._createUniqueChild(n, 'dirx', 'DataArray_t', value=oneOver[count][0]) + Internal._createUniqueChild(n, 'diry', 'DataArray_t', value=oneOver[count][1]) + Internal._createUniqueChild(n, 'dirz', 'DataArray_t', value=oneOver[count][2]) + count+=1 + C.convertPyTree2File(t_local,FileName) + return None + #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ## MACRO FUNCTIONS #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -103,7 +130,11 @@ def computeMeshInfo__(z, dim): def _computeMeshInfo(t): np_total, nc_total, nf_total, nzones = 0,0,0,0 nc0, nc1, nc2 = 0,0,0 - + + NP = Cmpi.size + rank = Cmpi.rank + NPTS = numpy.zeros(NP) + NCELLS = numpy.zeros(NP) for z in Internal.getZones(t): cellN = Internal.getNodeFromName(z, 'cellN')[1] c0 = numpy.count_nonzero(cellN == 0); nc0 += c0 @@ -115,10 +146,22 @@ def _computeMeshInfo(t): np_total += np nc_total += nc nf_total += nf - nzones += 1 - - Cmpi.barrier() - print("Info: mesh info for rank {}: number of points: {:.2f}M / number of cells: {:.2f}M".format(Cmpi.rank, np_total/1.e6, nc_total/1.e6)) + nzones += 1 + + zNoGhost = C.rmGhostCells(z, z, 2) + NPTS[rank] += C.getNPts(zNoGhost) + NCELLS[rank] += C.getNCells(zNoGhost) + + NPTS = Cmpi.allreduce(NPTS ,op=Cmpi.SUM) + NCELLS = Cmpi.allreduce(NCELLS,op=Cmpi.SUM) + print("Info: mesh info for rank {}: number of points: {:.2f}M / number of cells: {:.2f}M".format(rank, np_total/1.e6, nc_total/1.e6),flush=True) + if rank ==0: + NcellsTot = numpy.sum(NCELLS) + ncells_percent= [] + for i in range(NP): + ncells_percent.append(NCELLS[i]/NcellsTot*100) + print('Info: Rank {} has {:.3f}e06 points & {:.3f}e06 cells & {} % of cells - no ghost cells'.format(i,int(NPTS[i])/1e06,int(NCELLS[i])/1e06,round(ncells_percent[i],2)),flush=True) + print('Info: Range of % of cells: {} - {}'.format(round(min(ncells_percent),2),round(max(ncells_percent),2)),flush=True) Cmpi.barrier() np_total = Cmpi.allreduce(np_total, op=Cmpi.SUM) @@ -130,18 +173,18 @@ def _computeMeshInfo(t): natures = [ncx/float(nc_total)*100. for ncx in [nc2, nc1, nc0]] if Cmpi.rank == 0: - print("Info: global mesh info: Interpolated cells (cellN 2): {:.2f}%, Computed cells (cellN 1): {:.2f}%, Blanked cells (cellN 0): {:.2f}%".format(*natures)) - print("Info: global mesh info: total number of points: {:.2f}M / total number of cells: {:.2f}M".format(np_total/1.e6, nc_total/1.e6)) - + print("Info: global mesh info: Interpolated cells (cellN 2): {:.2f}%, Computed cells (cellN 1): {:.2f}%, Blanked cells (cellN 0): {:.2f}%".format(*natures),flush=True) + print("Info: global mesh info: total number of points: {:.2f}M / total number of cells: {:.2f}M".format(np_total/1.e6, nc_total/1.e6),flush=True) Cmpi.barrier() return None def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=None, - snears=0.01, snearsf=None, dfar=10., dfarDir=0, dfarList=[], vmin=21, depth=2, frontType=1, mode=0, - IBCType=1, verbose=True, - check=False, balancing=False, distribute=False, twoFronts=False, cartesian=False, - yplus=100., Lref=1., correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1.): + snears=0.01, snearsf=None, dfar=10., dfarDir=0, dfarList=[], vmin=21, depth=2, frontType=1, mode=0, + IBCType=1, verbose=True, + check=False, balancing=False, distribute=False, twoFronts=False, cartesian=False, + yplus=100., Lref=1., correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1., + tbOneOver=None): import Generator.IBM as G_IBM import time as python_time @@ -180,12 +223,14 @@ def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tin #=================== # STEP 1 : GENERATE MESH #=================== + listF1save = [] if t_in is None: if verbose: pt0 = python_time.time(); printTimeAndMemory__('generate Cartesian mesh', time=-1) test.printMem("Info: prepareIBMDataPara: generate Cartesian mesh [start]") t = G_IBM.generateIBMMeshPara(tb, vmin=vmin, snears=snears, dimPb=dimPb, dfar=dfar, dfarList=dfarList, tbox=tbox, - snearsf=snearsf, check=check, to=to, ext=depth+1, - expand=expand, dfarDir=dfarDir, check_snear=False, mode=mode) + snearsf=snearsf, check=check, to=to, ext=depth+1, + expand=expand, dfarDir=dfarDir, check_snear=False, mode=mode, + tbOneOver=tbOneOver, listF1save = listF1save) Internal._rmNodesFromName(tb,"SYM") if balancing and Cmpi.size > 1: _redispatch__(t=t) @@ -213,7 +258,7 @@ def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tin _blankingIBM(t, tb, dimPb=dimPb, frontType=frontType, IBCType=IBCType, depth=depth, Reynolds=Reynolds, yplus=yplus, Lref=Lref, twoFronts=twoFronts, heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42, - wallAdaptF42=wallAdaptF42, blankingF42=blankingF42) + wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, listF1save = listF1save) Cmpi.barrier() _redispatch__(t=t) if verbose: printTimeAndMemory__('blank by IBC bodies', time=python_time.time()-pt0) @@ -426,7 +471,8 @@ def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1. # OUT: (optional) centers:cellNIBC_2, centers:cellNFront_2 fields for second image points #========================================================================= def _blankingIBM__(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, yplus=100, Lref=1., - correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1.): + correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1., + listF1save=[]): snear_min = 10e10 for z in Internal.getZones(tb): @@ -470,6 +516,8 @@ def _blankingIBM__(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e X._setHoleInterpolatedPoints(t,depth=depth,dir=0,loc='centers',cellNName='cellNMin',addGC=False) for z in Internal.getZones(t): + if z[0] in listF1save: continue + h = abs(C.getValue(z,'CoordinateX',0)-C.getValue(z,'CoordinateX',1)) if yplus > 0.: height = G_IBM_Height.computeModelisationHeight(Re=Reynolds, yplus=yplus, L=Lref) @@ -561,17 +609,18 @@ def findIsoFront(cellNFront, Dist_1, Dist_2): return None def blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, yplus=100, Lref=1., twoFronts=False, - correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1.): + correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1., listF1save=[]): """Blank the computational tree by IBC bodies for IBM pre-processing.""" tp = Internal.copyRef(t) _blankingIBM(t, tb, dimPb=dimPb, frontType=frontType, IBCType=IBCType, depth=depth, Reynolds=Reynolds, yplus=yplus, Lref=Lref, twoFronts=twoFronts, heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42, - wallAdaptF42=wallAdaptF42, blankingF42=blankingF42) + wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, listF1save=listF1save) return tp def _blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, yplus=100, Lref=1., twoFronts=False, - correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1.): + correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1., + listF1save=[]): """Blank the computational tree by IBC bodies for IBM pre-processing.""" if dimPb == 2: T._addkplane(tb) @@ -584,7 +633,7 @@ def _blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, _blankingIBM__(t, tb, dimPb=dimPb, frontType=frontType, IBCType=IBCType, depth=depth, Reynolds=Reynolds, yplus=yplus, Lref=Lref, heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42, - wallAdaptF42=wallAdaptF42, blankingF42=blankingF42) + wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, listF1save=listF1save) C._initVars(t, '{centers:cellNIBC}={centers:cellN}') diff --git a/Cassiopee/Generator/Generator/IBM.py b/Cassiopee/Generator/Generator/IBM.py index 3d8c5dafb..42f71019b 100644 --- a/Cassiopee/Generator/Generator/IBM.py +++ b/Cassiopee/Generator/Generator/IBM.py @@ -363,7 +363,9 @@ def _addBCsForSymmetry(t, bbox=None, dimPb=3, dir_sym=0, X_SYM=0., depth=2): def generateIBMMeshPara(tb, vmin=15, snears=None, dimPb=3, dfar=10., dfarList=[], tbox=None, snearsf=None, check=True, to=None, ext=2, - expand=3, dfarDir=0, check_snear=False, mode=0): + expand=3, dfarDir=0, check_snear=False, mode=0, + tbOneOver=None, listF1save = []): + import KCore.test as test # list of dfars if dfarList == []: zones = Internal.getZones(tb) @@ -455,7 +457,7 @@ def generateIBMMeshPara(tb, vmin=15, snears=None, dimPb=3, dfar=10., dfarList=[] del o # fill vmin + merge in parallel - res = octree2StructLoc__(p, vmin=vmin, ext=-1, optimized=0, parento=parento, sizeMax=1000000) + res = octree2StructLoc__(p, vmin=vmin, ext=-1, optimized=0, parento=parento, sizeMax=1000000, tbOneOver=tbOneOver) del p if parento is not None: for po in parento: del po @@ -466,6 +468,14 @@ def generateIBMMeshPara(tb, vmin=15, snears=None, dimPb=3, dfar=10., dfarList=[] C._addState(t, 'EquationDimension', dimPb) + ##Keep F1 regions - for F1 & F42 synergy + if tbOneOver: + tbF1 = Internal.getNodesFromNameAndType(tbOneOver, '*KeepF1*', 'CGNSBase_t') + tbbBTmp = G.BB(tbF1) + interDict_scale = X.getIntersectingDomains(tbbBTmp, t) + for kk in interDict_scale: + for kkk in interDict_scale[kk]: listF1save.append(kkk) + # Add xzones for ext tbb = Cmpi.createBBoxTree(t) interDict = X.getIntersectingDomains(tbb) @@ -473,6 +483,35 @@ def generateIBMMeshPara(tb, vmin=15, snears=None, dimPb=3, dfar=10., dfarList=[] del tbb Cmpi._addXZones(t, graph, variables=[], cartesian=True) + #Turn Cartesian grid into a rectilinear grid + test.printMem(">>> cart grids --> rectilinear grids [start]") + listDone = [] + if tbOneOver: + tbb = G.BB(t) + if dimPb==2: + T._addkplane(tbb) + T._contract(tbb, (0,0,0), (1,0,0), (0,1,0), 0.01) + + ##RECTILINEAR REGION + ##Select regions that need to be coarsened + tbbB = G.BB(tbOneOver) + interDict_scale = X.getIntersectingDomains(tbbB, tbb) + ##Avoid a zone to be coarsened twice + for i in interDict_scale: + (b,btmp) = Internal.getParentOfNode(tbOneOver,Internal.getNodeByName(tbOneOver,i)) + checkOneOver = Internal.getNodeByName(b,".Solver#define") ##Needed for F1 & F42 approach + if checkOneOver: + b = Internal.getNodeByName(b,".Solver#define") + oneoverX = int(Internal.getNodeByName(b, 'dirx')[1]) + oneoverY = int(Internal.getNodeByName(b, 'diry')[1]) + oneoverZ = int(Internal.getNodeByName(b, 'dirz')[1]) + for z in interDict_scale[i]: + if z not in listDone: + zLocal = Internal.getNodeFromName(t,z) + T._oneovern(zLocal, (oneoverX,oneoverY,oneoverZ)); + listDone.append(z) + test.printMem(">>> cart grids --> rectilinear grids [end]") + zones = Internal.getZones(t) coords = C.getFields(Internal.__GridCoordinates__, zones, api=2) if symmetry==0: extBnd = 0 @@ -758,7 +797,7 @@ def addRefinementZones(o, tb, tbox, snearsf, vmin, dim): return Internal.getNodeFromType2(to, 'Zone_t') # only in generateIBMMeshPara and generateCartMesh__ -def octree2StructLoc__(o, parento=None, vmin=21, ext=0, optimized=0, sizeMax=4e6,tbOneOver=None): +def octree2StructLoc__(o, parento=None, vmin=21, ext=0, optimized=0, sizeMax=4e6, tbOneOver=None): sizeMax=int(sizeMax) dim = Internal.getZoneDim(o) if dim[3] == 'QUAD': dimPb = 2 @@ -915,7 +954,7 @@ def octree2StructLoc__(o, parento=None, vmin=21, ext=0, optimized=0, sizeMax=4e6 for i in range(NBases): ZONEStbOneOverTmp[i] = T.mergeCart(ZONEStbOneOverTmp[i][0]+ZONEStbOneOverTmp[i][1]+ZONEStbOneOverTmp[i][2]+ \ ZONEStbOneOverTmp[i][3]+ZONEStbOneOverTmp[i][4]+ZONEStbOneOverTmp[i][5]+ \ - ZONEStbOneOverTmp[i][6]+ZONEStbOneOverTmp[i][7],sizeMax=sizeMax) + ZONEStbOneOverTmp[i][6]+ZONEStbOneOverTmp[i][7], sizeMax=sizeMax) zones +=ZONEStbOneOverTmp[i] del ZONEStbOneOver del ZONEStbOneOverTmp @@ -931,7 +970,7 @@ def octree2StructLoc__(o, parento=None, vmin=21, ext=0, optimized=0, sizeMax=4e6 if ZONEStbOneOver is not None: for i in range(NBases): ZONEStbOneOverTmp[i] = T.mergeCart(ZONEStbOneOverTmp[i][0]+ZONEStbOneOverTmp[i][1]+ \ - ZONEStbOneOverTmp[i][2]+ZONEStbOneOverTmp[i][3],sizeMax=sizeMax) + ZONEStbOneOverTmp[i][2]+ZONEStbOneOverTmp[i][3], sizeMax=sizeMax) zones +=ZONEStbOneOverTmp[i] del ZONEStbOneOver del ZONEStbOneOverTmp @@ -1050,3 +1089,45 @@ def buildParentOctrees__(o, tb, snears=None, snearFactor=4., dfar=10., dfarList= return OCTREEPARENTS + +def _projectMeshSize(t, NPas=10, span=1, dictNz=None, isCartesianExtrude=False): + """Predicts the final size of the mesh when extruding 2D to 3D in the z-direction. + Usage: loads(t, NPas, span, dictNz, isCartesianExtrude)""" + NP = Cmpi.size + rank = Cmpi.rank + NPTS = numpy.zeros(NP) + NCELLS = numpy.zeros(NP) + NPTS_noghost = numpy.zeros(NP) + NCELLS_noghost = numpy.zeros(NP) + if isinstance(t, str): + h = Filter.Handle(t) + t = h.loadFromProc(loadVariables=False) + h._loadVariables(t, var=['CoordinateX']) + + + for z in Internal.getZones(t): + name_zone = z[0] + if not isCartesianExtrude: + if dictNz is not None: Nk = int(dictNz[name_zone]) + else: Nk = NPas-1 + else: + h = abs(C.getValue(z,'CoordinateX',0)-C.getValue(z,'CoordinateX',1)) + NPas_local = int(round(span/h)/4) + if NPas_local < 2: + print("WARNING:: MPI rank %d || Zone %s has Nz=%d and is being clipped to Nz=2"%(rank,z[0],NPas_local)) + NPas_local = 2 + Nk = NPas_local + Nk += 1 + NPTS[rank] += C.getNPts(z)/2*(Nk+4) + NCELLS[rank] += C.getNCells(z)*(Nk+3) + NPTS_noghost[rank] += C.getNPts(C.rmGhostCells(z, z, 2))*Nk + NCELLS_noghost[rank] += C.getNCells(C.rmGhostCells(z, z, 2))*(Nk-1) + NPTS = Cmpi.allreduce(NPTS ,op=Cmpi.SUM) + NCELLS = Cmpi.allreduce(NCELLS,op=Cmpi.SUM) + NPTS_noghost = Cmpi.allreduce(NPTS_noghost ,op=Cmpi.SUM) + NCELLS_noghost = Cmpi.allreduce(NCELLS_noghost,op=Cmpi.SUM) + if rank ==0: + print('Projected mesh size with ghost: {} million points & {} million cells'.format(numpy.sum(NPTS)/1.e6,numpy.sum(NCELLS)/1.e6)) + print('Projected mesh size without ghost: {} million points & {} million cells'.format(numpy.sum(NPTS_noghost)/1.e6,numpy.sum(NCELLS_noghost)/1.e6)) + return None +