From f7812c45ae201356e574eb697678b89045537c76 Mon Sep 17 00:00:00 2001 From: antjost Date: Thu, 29 Aug 2024 13:41:31 +0200 Subject: [PATCH] [IBM] Move the split of tbOneOver & tbF1 from tbox into the macro IBM functions and not in the overarching prepareIBMData function. Done as part of the IBM clean up of 2024. --- .../Apps/test/FastIBMWireMeshModel_t1.py | 1 - Cassiopee/Connector/Connector/IBM.py | 111 ++++++++++-------- Cassiopee/Geom/Geom/IBM.py | 47 ++++++-- 3 files changed, 97 insertions(+), 62 deletions(-) diff --git a/Cassiopee/Apps/test/FastIBMWireMeshModel_t1.py b/Cassiopee/Apps/test/FastIBMWireMeshModel_t1.py index 008d773ff..76b6f079c 100644 --- a/Cassiopee/Apps/test/FastIBMWireMeshModel_t1.py +++ b/Cassiopee/Apps/test/FastIBMWireMeshModel_t1.py @@ -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) diff --git a/Cassiopee/Connector/Connector/IBM.py b/Cassiopee/Connector/Connector/IBM.py index 03ef64667..186a18e32 100644 --- a/Cassiopee/Connector/Connector/IBM.py +++ b/Cassiopee/Connector/Connector/IBM.py @@ -225,27 +225,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) @@ -264,7 +260,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) #=================== @@ -275,8 +271,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() if not skipRedispatchNonRegression: _redispatch__(t=t) if verbose: printTimeAndMemory__('blank by IBC bodies', time=python_time.time()-pt0) @@ -296,9 +294,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) #=================== @@ -308,8 +306,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) #=================== @@ -318,7 +315,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) @@ -394,13 +391,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): @@ -436,9 +432,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") @@ -453,7 +451,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) @@ -549,7 +547,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): @@ -569,7 +569,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}') @@ -688,7 +688,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') @@ -755,7 +755,7 @@ 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') @@ -763,20 +763,23 @@ def findIsoFront(cellNFront, Dist_1, Dist_2): 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) @@ -792,7 +795,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}') @@ -946,9 +949,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 @@ -1024,20 +1030,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 @@ -1058,7 +1066,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) @@ -1071,11 +1079,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 @@ -1326,8 +1334,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) @@ -1352,7 +1362,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 = [] @@ -1441,10 +1451,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) @@ -1454,7 +1465,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 diff --git a/Cassiopee/Geom/Geom/IBM.py b/Cassiopee/Geom/Geom/IBM.py index c9b948bc0..95be68b36 100644 --- a/Cassiopee/Geom/Geom/IBM.py +++ b/Cassiopee/Geom/Geom/IBM.py @@ -542,21 +542,46 @@ def determineClosedSolidFilament__(tb): filamentBases = [] isFilamentOnly= False - len_tb = len(Internal.getBases(tb)) for b in Internal.getBases(tb): if "IBCFil" in b[0]:filamentBases.append(b[0]) - if len(filamentBases) == len_tb:isFilamentOnly=True + if len(filamentBases) == len(Internal.getBases(tb)):isFilamentOnly=True isOrthoProjectFirst = isFilamentOnly - ## if tb has both a closed solid and filaments + ## assume isFilamentOnly=True tbFilament = Internal.copyTree(tb) + + ## if tb has 1) no filaments or 2) filaments & closed bodies if not isFilamentOnly: - tbFilament = [] - for b in filamentBases: - node_local = Internal.getNodeFromNameAndType(tb, b, 'CGNSBase_t') - tbFilament.append(node_local) - Internal._rmNode(tb,node_local) - isOrthoProjectFirst = True - tbFilament = C.newPyTree(tbFilament); - return [filamentBases, isFilamentOnly, isOrthoProjectFirst, tb, tbFilament] + if len(filamentBases)==0: + tbFilament=None + else: + ##tb : only closed bodies + ##tbFilament: filament bodies + tbFilament = [] + for b in filamentBases: + node_local = Internal.getNodeFromNameAndType(tb, b, 'CGNSBase_t') + tbFilament.append(node_local) + Internal._rmNode(tb,node_local) + isOrthoProjectFirst = True + tbFilament = C.newPyTree(tbFilament); + + return tb, tbFilament + + +def localWMMFlags__(tb,tbFilament): + isFilamentOnly=False + isWireModel =False + + if tbFilament: + if len(Internal.getBases(tbFilament))==len(Internal.getBases(tb)): + isFilamentOnly=True + for z in Internal.getZones(tbFilament): + 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 + return isFilamentOnly,isWireModel