Skip to content

Commit

Permalink
Merge pull request #44 from antjost/main
Browse files Browse the repository at this point in the history 
[IBM] Add rectilinear mesh capabilities to Connector/IBM.py IBM prep approach. Copy functionality from Apps/IBM.py.
  • Loading branch information
@benoit128
benoit128 authored Jun 17, 2024
2 parents 0bdf16f + c86309a commit 77578df
Show file tree
Hide file tree
Showing 2 changed files with 155 additions and 25 deletions.
89 changes: 69 additions & 20 deletions Cassiopee/Connector/Connector/IBM.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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}')

Expand Down
91 changes: 86 additions & 5 deletions Cassiopee/Generator/Generator/IBM.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -466,13 +468,50 @@ 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)
graph = Cmpi.computeGraph(tbb, type='bbox', intersectionsDict=interDict, reduction=False)
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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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

0 comments on commit 77578df

Please sign in to comment.