Skip to content

Commit

Permalink
Merge pull request #150 from antjost/devTbUnified
Browse files Browse the repository at this point in the history 
[IBM] Update split tb and tbFilament. Done as part of the IBM clean up of 2024.
  • Loading branch information
@vincentcasseau
vincentcasseau authored Aug 30, 2024
2 parents 6f53b17 + f7812c4 commit e0b5ad9
Show file tree
Hide file tree
Showing 3 changed files with 97 additions and 62 deletions.
1 change: 0 additions & 1 deletion Cassiopee/Apps/test/FastIBMWireMeshModel_t1.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,6 @@
t,tc = X_IBM.prepareIBMData(tb , None , None , tbox=tboffset,
snears=snears , dfars=dfars , vmin=vmin,
check=False , frontType=1)

test.testT(t , 1)
test.testT(tc, 2)

Expand Down
111 changes: 61 additions & 50 deletions Cassiopee/Connector/Connector/IBM.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -224,27 +224,23 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N
#===================
# STEP 0 : GET FILAMENT BODIES
#===================
[filamentBases, isFilamentOnly, isOrthoProjectFirst, tb, tbFilament]=D_IBM.determineClosedSolidFilament__(tb)
isWireModel=False
for z in Internal.getZones(t_case):
soldef = Internal.getNodeFromName(z,'.Solver#define')
if soldef is not None:
ibctype = Internal.getNodeFromName(soldef,'ibctype')
if ibctype is not None:
if Internal.getValue(ibctype) == "wiremodel":
isWireModel=True
break
if isFilamentOnly: tbLocal=tbFilament
else: tbLocal=Internal.merge([tb,tbFilament])
tb, tbFilament = D_IBM.determineClosedSolidFilament__(tb)
isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament)

#===================
# STEP 1 : GENERATE MESH
#===================
if isFilamentOnly: tbLocal=tbFilament
elif tbFilament: tbLocal=Internal.merge([tb,tbFilament])
else: tbLocal=tb


if t_in is None:
if verbose: pt0 = python_time.time(); printTimeAndMemory__('generate Cartesian mesh', time=-1)
test.printMem("Info: prepareIBMData: generate Cartesian mesh [start]")
t = G_IBM.generateIBMMesh(tbLocal, vmin=vmin, snears=snears, dimPb=dimPb, dfars=dfars, tbox=tbox,
snearsf=snearsf, check=check, to=to, ext=depth+1,
expand=expand, dfarDir=dfarDir, octreeMode=octreeMode)
snearsf=snearsf, check=check, to=to, ext=depth+1,
expand=expand, dfarDir=dfarDir, octreeMode=octreeMode)
Internal._rmNodesFromName(tb,"SYM")

if balancing and Cmpi.size > 1: _redispatch__(t=t)
Expand All @@ -263,7 +259,7 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N
if verbose: pt0 = python_time.time(); printTimeAndMemory__('compute wall distance', time=-1)
_dist2wallIBM(t, tb, dimPb=dimPb, frontType=frontType, Reynolds=Reynolds, yplus=yplus, Lref=Lref,
correctionMultiCorpsF42=correctionMultiCorpsF42, heightMaxF42=heightMaxF42,
filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament)
tbFilament=tbFilament)
if verbose: printTimeAndMemory__('compute wall distance', time=python_time.time()-pt0)

#===================
Expand All @@ -274,8 +270,10 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N
Reynolds=Reynolds, yplus=yplus, Lref=Lref, twoFronts=twoFronts,
heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42,
wallAdaptF42=wallAdaptF42, blankingF42=blankingF42,
filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament,
isWireModel=isWireModel)
tbFilament=tbFilament)
##Keep for now to check
#filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament,
#isWireModel=isWireModel)
Cmpi.barrier()
_redispatch__(t=t)
if verbose: printTimeAndMemory__('blank by IBC bodies', time=python_time.time()-pt0)
Expand All @@ -295,9 +293,9 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N
# STEP 4 : BUILD FRONT
#===================
if verbose: pt0 = python_time.time(); printTimeAndMemory__('build IBM front', time=-1)
t, tc, front, front2, frontWMM = buildFrontIBM(t, tc, dimPb=dimPb, frontType=frontType,
t, tc, front, front2, frontWMM = buildFrontIBM(t, tc, tb=tb, dimPb=dimPb, frontType=frontType,
cartesian=cartesian, twoFronts=twoFronts, check=check,
isWireModel=isWireModel)
tbFilament=tbFilament)
if verbose: printTimeAndMemory__('build IBM front', time=python_time.time()-pt0)

#===================
Expand All @@ -307,8 +305,7 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N
_setInterpDataIBM(t, tc, tb, front, front2=front2, dimPb=dimPb, frontType=frontType, IBCType=IBCType, depth=depth,
Reynolds=Reynolds, yplus=yplus, Lref=Lref,
cartesian=cartesian, twoFronts=twoFronts, check=check,
isWireModel=isWireModel, tbFilament=tbFilament, isOrthoProjectFirst=isOrthoProjectFirst, isFilamentOnly=isFilamentOnly,
frontWMM=frontWMM)
tbFilament=tbFilament, frontWMM=frontWMM)
if verbose: printTimeAndMemory__('compute interpolation data (IBM)', time=python_time.time()-pt0)

#===================
Expand All @@ -317,7 +314,7 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N
if verbose: pt0 = python_time.time(); printTimeAndMemory__('initialize and clean', time=-1)

t, tc, tc2 = initializeIBM(t, tc, tb, tinit=tinit, tbCurvi=tbCurvi, dimPb=dimPb, twoFronts=twoFronts,
isWireModel=isWireModel, tbFilament=tbFilament, filamentBases=filamentBases, isFilamentOnly=isFilamentOnly)
tbFilament=tbFilament)

if distribute and Cmpi.size > 1: _redispatch__(t=t, tc=tc, tc2=tc2, twoFronts=twoFronts)

Expand Down Expand Up @@ -393,13 +390,12 @@ def _redispatch__(t=None, tc=None, tc2=None, twoFronts=False):
# OUT: (optional) centers:TurbulentDistance_body%i fields
#=========================================================================
def dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1.,
correctionMultiCorpsF42=False, heightMaxF42=-1., filamentBases=[], isFilamentOnly=False,
tbFilament=None):
correctionMultiCorpsF42=False, heightMaxF42=-1., tbFilament=None):
"""Compute the wall distance for IBM pre-processing."""
tp = Internal.copyRef(t)
_dist2wallIBM(t, tb, dimPb=dimPb, frontType=frontType, Reynolds=Reynolds, yplus=yplus, Lref=Lref,
correctionMultiCorpsF42=correctionMultiCorpsF42, heightMaxF42=heightMaxF42,
filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament)
tbFilament=tbFilament)
return tp

def _dist2wallIBMFilamentWMM__(t, tb2, tbsave, dimPb):
Expand Down Expand Up @@ -435,9 +431,11 @@ def _dist2wallIBMFilamentWMM__(t, tb2, tbsave, dimPb):
return None

def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1.,
correctionMultiCorpsF42=False, heightMaxF42=-1.,
filamentBases=[], isFilamentOnly=False, tbFilament=None):
correctionMultiCorpsF42=False, heightMaxF42=-1., tbFilament=None):
"""Compute the wall distance for IBM pre-processing."""

isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament)

if dimPb == 2:
# Set CoordinateZ to a fixed value
z0 = Internal.getNodeFromType2(t, "Zone_t")
Expand All @@ -452,7 +450,7 @@ def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1.

tbsave = tb2

if filamentBases and not isFilamentOnly:
if tbFilament and not isFilamentOnly:
C._initVars(t,'{centers:TurbulentDistanceSolid}={centers:TurbulentDistance}')
if dimPb ==2: tb2 = C.initVars(tbFilament, 'CoordinateZ', dz)
_dist2wallIBMFilamentWMM__(t, tb2, tbsave, dimPb)
Expand Down Expand Up @@ -548,7 +546,9 @@ def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1.
#=========================================================================
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.,
filamentBases=[], isFilamentOnly=False, tbFilament=None):
tbFilament=None):

isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament)

snear_min = 10e10
for z in Internal.getZones(tb):
Expand All @@ -568,7 +568,7 @@ def _blankingIBM__(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e

C._initVars(t,'{centers:cellN}={centers:cellNIBC}')

if filamentBases:
if tbFilament:
C._initVars(t,'{centers:cellNFil}={centers:cellN}')
C._initVars(t,'{centers:cellNFilWMM}={centers:cellN}')

Expand Down Expand Up @@ -687,7 +687,7 @@ def findIsoFront(cellNFront, Dist_1, Dist_2):
if wallAdaptF42 is None: X._maximizeBlankedCells(t, depth=2, addGC=False)
else: print("Info: blankingIBM: blankingF42 cannot operate with a local modeling height")

if filamentBases:
if tbFilament:
tbFilamentWMM = []
for z in Internal.getZones(tbFilament):
soldef = Internal.getNodeFromName(z,'.Solver#define')
Expand Down Expand Up @@ -754,28 +754,31 @@ def findIsoFront(cellNFront, Dist_1, Dist_2):
cellNFil[i,j]=2
cellN[i,j]=2
C._rmVars(t,['centers:TurbulentDistanceFilament','centers:TurbulentDistanceFilamentWMM'])
if filamentBases: C._rmVars(t,['centers:TurbulentDistanceSolid'])
if tbFilament: C._rmVars(t,['centers:TurbulentDistanceSolid'])
if isFilamentOnly: C._initVars(t,'{centers:cellNFilWMM}={centers:cellNFil}')

if not isFilamentOnly: _removeBlankedGrids(t, loc='centers')
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.,
filamentBases=[], isFilamentOnly=False, tbFilament=None, isWireModel=False):
tbFilament=None):
"""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,
filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament, isWireModel=isWireModel)
tbFilament=tbFilament)
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.,
filamentBases=[], isFilamentOnly=False, tbFilament=None, isWireModel=False):
tbFilament=None):
"""Blank the computational tree by IBC bodies for IBM pre-processing."""

isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament)

if dimPb == 2:
T._addkplane(tb)
T._contract(tb, (0,0,0), (1,0,0), (0,1,0), 0.01)
Expand All @@ -791,7 +794,7 @@ def _blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6,
Reynolds=Reynolds, yplus=yplus, Lref=Lref,
heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42,
wallAdaptF42=wallAdaptF42, blankingF42=blankingF42,
filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament)
tbFilament=tbFilament)

C._initVars(t, '{centers:cellNIBC}={centers:cellN}')

Expand Down Expand Up @@ -945,9 +948,12 @@ def _pushBackImageFront2__(t, tc, tbbc, cartesian=False):

return None

def buildFrontIBM(t, tc, dimPb=3, frontType=1, cartesian=False, twoFronts=False, check=False,
isWireModel=False):
def buildFrontIBM(t, tc, tb=None, dimPb=3, frontType=1, cartesian=False, twoFronts=False, check=False,
tbFilament=None):
"""Build the IBM front for IBM pre-processing."""

isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament)

tbbc = Cmpi.createBBoxTree(tc)

interpDataType = 0 if cartesian else 1
Expand Down Expand Up @@ -1023,20 +1029,22 @@ def buildFrontIBM(t, tc, dimPb=3, frontType=1, cartesian=False, twoFronts=False,
#=========================================================================
def setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6,
yplus=100, Lref=1., cartesian=False, twoFronts=False, check=False,
isWireModel=False, tbFilament=None, isOrthoProjectFirst=False, frontWMM=None):
tbFilament=None, frontWMM=None):
"""Compute the transfer coefficients and data for IBM pre-processing."""
tp = Internal.copyRef(t)
_setInterpDataIBM(t, tc, tb, front, front2=front2, dimPb=dimPb, frontType=frontType, IBCType=IBCType, depth=depth,
Reynolds=Reynolds, yplus=yplus, Lref=Lref,
cartesian=cartesian, twoFronts=twoFronts, check=check,
isWireModel=isWireModel, tbFilament=tbFilament, isOrthoProjectFirst=isOrthoProjectFirst,
frontWMM=frontWMM)
tbFilament=tbFilament, frontWMM=frontWMM)
return tp

def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6,
yplus=100, Lref=1., cartesian=False, twoFronts=False, check=False,
isWireModel=False, tbFilament=None, isOrthoProjectFirst=False, isFilamentOnly=False, frontWMM=None):
tbFilament=None, frontWMM=None):
"""Compute the transfer coefficients and data for IBM pre-processing."""

isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament)

tbbc = Cmpi.createBBoxTree(tc)

interpDataType = 0 if cartesian else 1
Expand All @@ -1057,7 +1065,7 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy
else: tb_local= Internal.merge([tb,tbFilament])
res = getAllIBMPoints(zonesRIBC, loc='centers',tb=tb_local, tfront=front, frontType=frontType,
cellNName='cellNIBC', depth=depth, IBCType=IBCType, Reynolds=Reynolds, yplus=yplus, Lref=Lref,
isOrthoFirst=isOrthoProjectFirst)
isOrthoFirst=isFilamentOnly)
if twoFronts:
res2 = getAllIBMPoints(zonesRIBC, loc='centers',tb=tb, tfront=front2, frontType=frontType,
cellNName='cellNIBC', depth=depth, IBCType=IBCType, Reynolds=Reynolds, yplus=yplus, Lref=Lref)
Expand All @@ -1070,11 +1078,11 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy
Internal._rmNode(tb_filament_localWMM,zlocal)
res2 = getAllIBMPoints(zonesRIBC, loc='centers',tb=tb_filament_localWMM, tfront=frontWMM, frontType=frontType,
cellNName='cellNFilWMM', depth=depth, IBCType=IBCType, Reynolds=Reynolds, yplus=yplus, Lref=Lref,
isWireModel=isWireModel, isOrthoFirst=isOrthoProjectFirst)
isWireModel=isWireModel, isOrthoFirst=isFilamentOnly)

restmp = getAllIBMPoints(zonesRIBC, loc='centers',tb=tb_filament_localWMM, tfront=frontWMM, frontType=frontType,
cellNName='cellNFilWMM', depth=depth, IBCType=IBCType, Reynolds=Reynolds,
yplus=yplus, Lref=Lref, isOrthoFirst=isOrthoProjectFirst)
yplus=yplus, Lref=Lref, isOrthoFirst=isFilamentOnly)

for j in range(3):
##delete in res
Expand Down Expand Up @@ -1325,8 +1333,10 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy
# OUT: updated t, tc
# OUT: (optional) new connectivity tree tc2
#=========================================================================
def _recomputeDistForViscousWall__(t, tb, tbCurvi=None, dimPb=3, tbFilament=None, filamentBases=[], isFilamentOnly=False):
def _recomputeDistForViscousWall__(t, tb, tbCurvi=None, dimPb=3, tbFilament=None):

isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament)

filteredBC=['outpress', 'inj', 'slip', 'overlap']
recompute=False
tb2 = Internal.copyRef(tb)
Expand All @@ -1351,7 +1361,7 @@ def _recomputeDistForViscousWall__(t, tb, tbCurvi=None, dimPb=3, tbFilament=None

tbsave = tb2

if filamentBases and not isFilamentOnly:
if tbFilament and not isFilamentOnly:
if dimPb == 2: tb2 = C.initVars(tbFilament, 'CoordinateZ', dz)
tbFilamentnoWMM = []
tbFilamentWMM = []
Expand Down Expand Up @@ -1440,10 +1450,11 @@ def _tInitialize__(t, tinit=None, model='NSTurbulent', isWireModel=False):
C._initVars(z,'{centers:'+v_local+'_WM}=0.')
return None

def initializeIBM(t, tc, tb, tinit=None, tbCurvi=None, dimPb=3, twoFronts=False, isWireModel=False,
tbFilament=None, filamentBases=[], isFilamentOnly=False):
def initializeIBM(t, tc, tb, tinit=None, tbCurvi=None, dimPb=3, twoFronts=False, tbFilament=None):
"""Initialize the computational and connectivity trees for IBM pre-processing."""

isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament)

model = Internal.getNodeFromName(tb, 'GoverningEquations')
if model is None: raise ValueError('initializeIBM: GoverningEquations is missing in input geometry tree.')
model = Internal.getValue(model)
Expand All @@ -1453,7 +1464,7 @@ def initializeIBM(t, tc, tb, tinit=None, tbCurvi=None, dimPb=3, twoFronts=False,
ibctypes = list(set(Internal.getValue(ibc) for ibc in ibctypes))

if model != 'Euler':
_recomputeDistForViscousWall__(t, tb, tbCurvi=tbCurvi, dimPb=dimPb, tbFilament=tbFilament, filamentBases=filamentBases, isFilamentOnly=isFilamentOnly )
_recomputeDistForViscousWall__(t, tb, tbCurvi=tbCurvi, dimPb=dimPb, tbFilament=tbFilament)

tc2 = Internal.copyTree(tc) if twoFronts or isWireModel else None

Expand Down
Loading

0 comments on commit e0b5ad9

Please sign in to comment.