From 46f35c607cc976ca22df33849bb5991314f8a18c Mon Sep 17 00:00:00 2001 From: antjost Date: Wed, 26 Jun 2024 13:16:38 +0200 Subject: [PATCH 1/2] [IBM] Move extrude 2D to 3D IBM from Apps/IBM.py to IBM.py in each relevant module. --- Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py | 131 +++++++ Cassiopee/Apps/test/IBMrotationCompute_t1.py | 211 ++++++++++ Cassiopee/Apps/test/IBMrotation_t1.py | 4 + Cassiopee/Connector/Connector/IBM.py | 367 +++++++++++++++++- .../Connector/test/IBMExtrude2Dto3D_t1.py | 51 +++ Cassiopee/Generator/Generator/IBM.py | 188 +++++++++ 6 files changed, 942 insertions(+), 10 deletions(-) create mode 100644 Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py create mode 100644 Cassiopee/Apps/test/IBMrotationCompute_t1.py create mode 100644 Cassiopee/Connector/test/IBMExtrude2Dto3D_t1.py diff --git a/Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py b/Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py new file mode 100644 index 000000000..33f14cb78 --- /dev/null +++ b/Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py @@ -0,0 +1,131 @@ +# - Fast.IBM - +import Apps.Fast.IBM as App +import Connector.PyTree as X +import Connector.IBM as X_IBM +import Converter.Internal as Internal +import Converter.PyTree as C +import Fast.PyTree as Fast +import FastS.PyTree as FastS +import Generator.PyTree as G +import Geom.IBM as D_IBM +import Geom.PyTree as D +import Initiator.PyTree as Initiator +import KCore.test as test +import math +import numpy as np +test.TOLERANCE = 1.e-6 + +LOCAL = test.getLocal() + +#Fabrication case 2d pour maillage octree +lines=[] +lines.append( D.line((0.05,-0.15, 0), (0.1 ,-0.08,0), N=10) ) +lines.append( D.line((0.10,-0.03, 0), (0.1 ,-0.08,0), N=10) ) +lines.append( D.line((0.1 ,-0.03, 0), (0.1 , 0. ,0), N=8) ) + +lines.append( D.line((0.15,-0.08, 0), (0.15,-0.03,0), N=10) ) +lines.append( D.line((0.15,-0.08, 0), (0.2 ,-0.15,0), N=10) ) +lines.append( D.line((0.05,-0.15, 0), (0.2 ,-0.15,0), N=10) ) + +lines.append( D.line((0.15,-0.03, 0), (0.15, 0 ,0), N=8 ) ) +lines.append( D.line((0.10, 0.01, 0), (0.15, 0.01,0), N=10) ) +lines.append( D.line((0.15, 0.01, 0), (0.15, 0 ,0), N=3) ) +lines.append( D.line((0.10, 0.01, 0), (0.10, 0 ,0), N=3) ) +case = C.newPyTree(['BOBY', lines]) + +D_IBM._setSnear(case, 0.002) +D_IBM._setIBCType(case, "Musker") +D_IBM._setDfar(case, 0.) +D_IBM._setFluidInside(case) + +zones = Internal.getZones(case) +for z in zones: + if '8' in z[0] or '9' in z[0] or '6' in z[0]: + D_IBM._setIBCType(z, "overlap") + D_IBM._setSnear(z, 0.008) + if 'line.1' == z[0] or '10' in z[0]: + D_IBM._setSnear(z, 0.008) + if '2' in z[0] or '0' in z[0]: + D_IBM._setSnear(z, 0.004) + +zones[6][0]='line10' +zones[8][0]='line9' +U0= 120.0 +P0= 101325.0 +L0= 1. +R0= 1.18 +T0= 279.15 + +equation = 'NSLaminar' +C._addState(case, 'GoverningEquations', equation ) +C._addState(case, 'EquationDimension', 2) +C._addState(case, UInf=U0, RoInf=R0, PInf=P0, LInf=L0, alphaZ=0., adim='dim3') +#C.convertPyTree2File(case,'verifcase.cgns') + +# Prepare cas 2d IBM +t_2d, tc_2d = App.prepare1(case, None, None, frontType=42 ,cleanCellN=False) + +test.testT(tc_2d, 1) +test.testT(t_2d , 2) + +#extrusion 3D IBM +extrusion ='cart'; span = 0.01 + +t_3d, tb_3d = App.extrudeCartesian(t_2d, case, extrusion=extrusion, NPas=5, span=span) + +test.testT(tb_3d, 3) +test.testT(t_3d , 4) + +#calcul interp IBM 3D; celln, distnce paroi 3D,.... deja obtenue lors de l'extrusion +interpDataType = 1 # on suppose maillage non cartesion pour interp +order = 2 +t_3d, tc_3d = App.prepare1(tb_3d,None, None, t_in=t_3d, extrusion=extrusion, interpDataType=interpDataType, order=order) + +Internal._rmNodesFromType(tc_2d,'Rind_t') +Internal._rmNodesFromName(tc_2d,Internal.__GridCoordinates__) +test.testT(tc_3d, 5) +test.testT(t_3d , 6) + + +#Maillage plaque plane Curviligne sans ghost +a = G.cart((0,0,0), (0.005,0.005,0.01), (200,100,5)) +C._addBC2Zone(a, 'farfield', 'BCFarfield', 'jmax') +C._addBC2Zone(a, 'inflow', 'BCInflow', 'imin') +C._addBC2Zone(a, 'outflow','BCOutflow', 'imax') +C._addBC2Zone(a, 'wallIn','BCWallViscous', [1,22,1,1,1,5]) +C._addBC2Zone(a, 'overlap','BCOverlap', [22,30,1,1,1,5]) +C._addBC2Zone(a, 'wallOut','BCWallViscous', [30,200,1,1,1,5]) +t_curvi = C.newPyTree(['Base', a]) + +zones = Internal.getZones(t_curvi) +#les zones curviligne possedanr raccord chimere avec zone Cart IBC doit avoir la racine "joinIBC" dans leur nom +zones[0][0]='curvi_joinIBC' +t_curvi = X.connectMatchPeriodic(t_curvi, translation=[0.,0.,0.04]) +#stretch maillage plaque direction normal paroi +for z in zones: + coordy = Internal.getNodeFromName(z,'CoordinateY')[1] + sh = np.shape(coordy) + for j in range(1,sh[1]): + coordy[:,j,:]=coordy[:,j-1,:]+0.002*1.02**j + +C._addState(t_curvi, 'GoverningEquations', equation ) +C._addState(t_curvi, 'EquationDimension', 3) +C._addState(t_curvi, UInf=U0, RoInf=R0, PInf=P0, LInf=L0, alphaZ=0., adim='dim3') +C._initVars(t_curvi,"{centers:Density}=1.18") +C._initVars(t_curvi,"{centers:VelocityX}=0.01") +C._initVars(t_curvi,"{centers:VelocityY}=0") +C._initVars(t_curvi,"{centers:VelocityZ}=0") +C._initVars(t_curvi,"{centers:Temperature}=279.15") + +C._initVars(t_3d,"{centers:Density}=1.18") +C._initVars(t_3d,"{centers:VelocityX}=0.01") +C._initVars(t_3d,"{centers:VelocityY}=0") +C._initVars(t_3d,"{centers:VelocityZ}=0") +C._initVars(t_3d,"{centers:Temperature}=279.15") + +#Calcul interpolation entre maillage curviligne et cartesien IBC et ajout Ghost au maillage curvi +#t_final, tc_final= App.setInterpData_Hybride(t_3d, tc_3d, t_curvi) +t_final, tc_final= X_IBM.setInterpDataHybrid(t_3d, tc_3d, t_curvi) +test.testT(t_final , 7) +test.testT(tc_final, 8) + diff --git a/Cassiopee/Apps/test/IBMrotationCompute_t1.py b/Cassiopee/Apps/test/IBMrotationCompute_t1.py new file mode 100644 index 000000000..2ad88adad --- /dev/null +++ b/Cassiopee/Apps/test/IBMrotationCompute_t1.py @@ -0,0 +1,211 @@ +# case - cylinder - IBM rotation - compute +import Connector.IBM as X_IBM +import Connector.Mpi as Xmpi +import Connector.PyTree as X +import Converter.Internal as Internal +import Converter.Mpi as Cmpi +import Converter.PyTree as C +import Fast.PyTree as Fast +import FastS.PyTree as FastS +import RigidMotion.PyTree as R +import Transform.PyTree as T +import KCore.test as test +import math +import numpy + +LOCAL = test.getLocal() + +## ======================= +## ======================= +TInf = 298.15 +UInf = 0.0001 +rad = 0.0005 +diam = rad*2 +RoInf= 1.18 +PInf = 101325 + +# necessaire pour les IBM +Model = 'NSLaminar' +dimPb = 2 +modulo_verif = 5 + +rpm = 100 +OMG = rpm*2*math.pi/60 +LInf = diam*0.5 +TInf = 298.15 +aref = math.sqrt(1.4*287.058*298.15) +Vtip = OMG*LInf +Mtip = OMG*LInf/aref +Re = RoInf*Vtip*diam/1.85e-05 +print(OMG,LInf,Vtip,Re) + +time_step= 1/(OMG*180/math.pi)/600 +## ======================= +## ======================= + +NIT = 15 + +numb={"temporal_scheme": "implicit_local", + "omp_mode":0, + "modulo_verif":modulo_verif} + +numz={"time_step": time_step, + "time_step_nature": "global", + "epsi_newton": 0.1} + +numz1 = numz.copy() +numz1['scheme'] = 'ausmpred' + +numz2 = numz.copy() +numz2['scheme'] = 'ausmpred' + +# load t, tc +fileTIn ='/t.cgns' +fileTcIn='/tc.cgns' +t,tc,ts,graph = Fast.load(LOCAL+fileTIn,LOCAL+fileTcIn) +baseNameIBM ='CARTESIAN_NEARBODY' +baseNameBKGD ='CARTESIAN_OFFBODY' + + +Fast._setNum2Base(t, numb) +for n in [baseNameIBM]: + b = Internal.getNodeFromName1(t, n) + Fast._setNum2Zones(b, numz1) +for n in [baseNameBKGD]: + b = Internal.getNodeFromName1(t, n) + Fast._setNum2Zones(b, numz2) + + +# load bodies +tb = C.convertFile2PyTree(LOCAL+'/bodiesBlank.cgns') + +# Warmup +it0 = 0; time0 = 0. +first = Internal.getNodeFromName1(t, 'Iteration') +if first is not None: it0 = Internal.getValue(first) +first = Internal.getNodeFromName1(t, 'Time') +if first is not None: time0 = Internal.getValue(first) +time_step = Internal.getNodeFromName(t, 'time_step') +time_step = Internal.getValue(time_step) + +# Tag les zones pour le motion +zones = Internal.getZones(t) +for z in zones: + timeMotion = Internal.getNodeFromName(z, 'TimeMotion') + if timeMotion is not None: + define = Internal.getNodeFromName1(z, '.Solver#define') + Internal._createUniqueChild(define, 'motion', 'DataArray_t', value='deformation') + +# Set interpolated points sur les bases IBM +C._initVars(t, "{centers:cellN#Motion}=1.") +##IBM CART +X._applyBCOverlaps(Internal.getNodeFromName1(t,baseNameIBM), depth=2, val=2, cellNName='cellN#Motion') +C._initVars(t, "{centers:cellN#MotionInit}={centers:cellN#Motion}") + +# Data for chimera motion +nbases = len(Internal.getBases(t)) +nbodies = len(Internal.getBases(tb)) +BM = numpy.zeros((nbases, nbodies),dtype=numpy.int32) +BM[1,0] = 1. + +# Creation BB tree commun +tBB = Cmpi.createBBoxTree(t) + +# Attachement des donneurs en dur (id on all procs) +intersectionDict={} +for z1 in Internal.getZones(Internal.getNodeFromName1(tBB,baseNameIBM)): + for z2 in Internal.getZones(Internal.getNodeFromName1(tBB,baseNameBKGD)): + Fast._addPair(intersectionDict, z1[0], z2[0]) + +procDict = None; +graphX = {} +procDict = Cmpi.getProcDict(tBB) +graphX = Cmpi.computeGraph(tBB, type='bbox2', t2=None, procDict=procDict, reduction=False,intersectionsDict=intersectionDict) + +R._evalPositionIBC(tc,time0) + +(t, tc, metrics) = FastS.warmup(t, tc, graph) + +# Doivent etre apres le warmup +dictOfADT={} +(dictOfNobOfRcvZones,dictOfNozOfRcvZones) = Fast.getDictOfNobNozOfRcvZones(t, intersectionDict) +(dictOfNobOfRcvZonesC,dictOfNozOfRcvZonesC) = Fast.getDictOfNobNozOfRcvZones(tc, intersectionDict) +(dictOfNobOfDnrZones,dictOfNozOfDnrZones) = Fast.getDictOfNobNozOfDnrZones(tc, intersectionDict, dictOfADT,cartFilter='CARTESIAN',isIbmAle=True) + +# modifie la vitesse de reference (si M=0) +cont = Internal.getNodeFromType(t, 'ReferenceState_t') +Minf = Internal.getNodeFromName(cont, 'Mach') +zones = Internal.getZones(t) +for z in zones: + n = Internal.getNodeFromName2(z, 'Parameter_real')[1] + n[5] = max(30, UInf) + +timeiter = time0 +varType = 0 + +# Get models +eqs = Internal.getNodeFromType(Internal.getNodeFromName1(t,baseNameIBM),'GoverningEquations_t') +Model1 = Internal.getValue(eqs) +eqs = Internal.getNodeFromType(Internal.getNodeFromName1(t, baseNameBKGD),'GoverningEquations_t') +Model2 = Internal.getValue(eqs) + +if Model1 == 'NSTurbulent' and Model2 == 'NSTurbulent': varType = 1 + +RefState = Internal.getNodeFromType(t,'ReferenceState_t') +VInf = Internal.getNodeFromName1(RefState,'VelocityX')[1][0] +PInf = Internal.getNodeFromName1(RefState,'Pressure')[1][0] +rhoInf = Internal.getNodeFromName1(RefState,'Density')[1][0] + +for it in range(it0,NIT+it0): + if Cmpi.rank == 0: print("it=%d, time=%f, angle=%f, N_rot=%d"%(it, timeiter, (timeiter*OMG*180./numpy.pi) % 360 , (timeiter*OMG*180./numpy.pi)//360 ),flush=True) + timeiter += time_step + + # bouge tout + R._evalPosition(tb, timeiter) + R._evalPosition(t, timeiter) + R._evalPosition(tc, timeiter) + R._evalGridSpeed(t, timeiter) + + R._evalPositionIBC(tc,timeiter) + + FastS.copy_velocity_ale(t, metrics, it=it) #get motion velocities @ face centers + + # Masquage + C._initVars(t,"{centers:cellN#Motion}={centers:cellN#MotionInit}") + bodies=[] + for base in Internal.getBases(tb): + bodies.append(Internal.getZones(base)) + + XRAYDIM1 = 2000;RAYDIM2 = XRAYDIM1 + if dimPb == 2: XRAYDIM2 = 2 + X._blankCells(t, bodies, BM, cellNName='cellN#Motion',XRaydim1=XRAYDIM1, XRaydim2=XRAYDIM2, blankingType='cell_intersect',dim=dimPb) + X._setHoleInterpolatedPoints(Internal.getNodeFromName1(t,baseNameBKGD), depth=2, loc='centers',addGC=False, cellNName='cellN#Motion', dir=2) + C._initVars(t, "{centers:cellN#Motion}=({centers:cellN#Static}<2)*{centers:cellN#Motion}+({centers:cellN#Static}>1)*minimum(1,{centers:cellN#Motion})") + C._initVars(t,"{centers:cellN}=minimum({centers:cellN#Motion}*{centers:cellN#Static},2.)") + C._cpVars(t, 'centers:cellN',tc, 'cellN') + ucData = (graphX, intersectionDict, dictOfADT, + dictOfNobOfRcvZones, dictOfNozOfRcvZones, + dictOfNobOfDnrZones, dictOfNozOfDnrZones, + dictOfNobOfRcvZonesC, dictOfNozOfRcvZonesC, + timeiter, procDict, True, varType, 1, 2, 1) + + FastS._compute(t, metrics, it, tc, graph, layer="Python", ucData=ucData) + + if it%modulo_verif==0: + FastS.display_temporal_criteria(t, metrics, it, format='double') + + Cmpi.barrier() + +for adt in dictOfADT.values(): + if adt is not None: C.freeHook(adt) + +Internal.createUniqueChild(t, 'Iteration', 'DataArray_t', value=NIT+it0) +Internal.createUniqueChild(t, 'Time', 'DataArray_t', value=timeiter) + +test.testT(t,1) +test.testT(tc,2) + + +#Cmpi.convertPyTree2File(t , LOCAL+'/t_restart.cgns') +#Cmpi.convertPyTree2File(tc, LOCAL+'/tc_restart.cgns') + diff --git a/Cassiopee/Apps/test/IBMrotation_t1.py b/Cassiopee/Apps/test/IBMrotation_t1.py index 78bd9437a..0ae91fc61 100644 --- a/Cassiopee/Apps/test/IBMrotation_t1.py +++ b/Cassiopee/Apps/test/IBMrotation_t1.py @@ -65,6 +65,7 @@ for z in Internal.getZones(tb): dfars.append(dfar_nb) t_ibm = G_IBM.generateIBMMesh(tb, vmin=vmin, dfars=dfars, dimPb=2, expand=2, ext=3) + R._setPrescribedMotion3(t_ibm ,'rot', axis_pnt=(0.,0.,0.), axis_vct=(0,0,1),omega=-OMG) t_ibm[2][1][0]='CARTESIAN_NEARBODY' @@ -72,6 +73,7 @@ tc_out=None t_ibm, tc_ibm = X_IBM.prepareIBMData(tb, t_out, tc_out, t_ibm, frontType=frontType) + C._rmBCOfType(t_ibm,'BCFarfield') C._fillEmptyBCWith(t_ibm,'dummy','BCExtrapolate', dim=dimPb) @@ -100,6 +102,7 @@ for z in Internal.getZones(tb_off): dfars.append(dfar_ext) t_off = G_IBM.generateIBMMesh(tb_off, vmin=vmin, dfars=dfars, dimPb=2, expand=2, ext=3) + t_off[2][1][0]='CARTESIAN_OFFBODY' X._applyBCOverlaps(t_off,depth=2,loc='centers') tc_off = C.node2Center(t_off) @@ -154,6 +157,7 @@ vars = ['Density', 'MomentumX', 'MomentumY', 'MomentumZ', 'EnergyStagnationDensity'] C._initVars(t_off, "centers:TurbulentDistance", 1e3) + t = C.mergeTrees(t_ibm, t_off) tc = C.mergeTrees(tc_ibm, tc_off) diff --git a/Cassiopee/Connector/Connector/IBM.py b/Cassiopee/Connector/Connector/IBM.py index 24bab60b1..8557cf534 100644 --- a/Cassiopee/Connector/Connector/IBM.py +++ b/Cassiopee/Connector/Connector/IBM.py @@ -29,14 +29,14 @@ varsn = ['gradxTurbulentDistance','gradyTurbulentDistance','gradzTurbulentDistance'] varsnDouble = ['gradxTurbulentDistanceDouble','gradyTurbulentDistanceDouble','gradzTurbulentDistanceDouble'] -TOLDIST = 1.e-14 -SHIFTF = 1.e-10 -EPSCART = 1.e-6 -TOLCELLN = 0.01 +TOLDIST = 1.e-14 +SHIFTF = 1.e-10 +EPSCART = 1.e-6 +TOLCELLN = 0.01 LBM_IBC_NUM = 113 #from param_solver.h -TypesOfIBC = XOD.TypesOfIBC +TypesOfIBC = XOD.TypesOfIBC # computes the friction velocity def _computeFrictionVelocity(a): @@ -97,11 +97,11 @@ def _addOneOverLocally(FileName,oneOver): #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ## MACRO FUNCTIONS #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -def printTimeAndMemory__(message, time=-1): +def printTimeAndMemory__(message, time=-1, functionName='prepareIBMData'): Cmpi.barrier() - test.printMem("Info: prepareIBMData: %s [%s]"%(message, 'end' if time > 0 else 'start')) + test.printMem("Info: %s: %s [%s]"%(functionName,message, 'end' if time > 0 else 'start')) Cmpi.barrier() - if time > 0 and Cmpi.rank == 0: print("Info: prepareIBMData: %s running time = %4.4fs"%(message, time)) + if time > 0 and Cmpi.rank == 0: print("Info: %s: %s running time = %4.4fs"%(functionName,message, time)) Cmpi.barrier() return None @@ -193,6 +193,12 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N if isinstance(t_case, str): tb = C.convertFile2PyTree(t_case) else: tb = Internal.copyTree(t_case) + ##[AJ] keep for now...will delete in the near future + ##if isinstance(tc_out, str): + ## if '/' in tc_out: fileoutpre = tc_out.split('/') + ## else: fileoutpre = ['.','template.cgns'] + ##else: fileoutpre = ['.','template.cgns'] + refstate = Internal.getNodeFromName(tb, 'ReferenceState') flowEqn = Internal.getNodeFromName(tb, 'FlowEquationSet') @@ -339,6 +345,349 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N if tc2 is not None: return t, tc, tc2 else: return t, tc +def prepareIBMDataExtrude(t_case, t_out, tc_out, t, to=None, + depth=2, frontType=1, octreeMode=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., + tbox=None, extrusion='cart'): + + import Generator.IBM as G_IBM + import time as python_time + + if isinstance(t_case, str): tb = C.convertFile2PyTree(t_case) + else: tb = Internal.copyTree(t_case) + + tbF1 = None + if tbox: + tbF1 = Internal.getNodesFromNameAndType(tbox, '*KeepF1*', 'CGNSBase_t') + tbox = None + + ##if isinstance(tc_out, str): fileoutpre = tc_out.split('/') + ##else: fileoutpre = ['.','template.cgns'] + + refstate = Internal.getNodeFromName(tb, 'ReferenceState') + flowEqn = Internal.getNodeFromName(tb, 'FlowEquationSet') + + Reynolds = Internal.getNodeFromName(tb, 'Reynolds') + if Reynolds is not None: + Reynolds = Internal.getValue(Reynolds) + if Reynolds < 1.e5: frontType = 1 + else: Reynolds = 1.e6 + + expand = 3 if frontType != 42 else 4 + + dimPb = Internal.getNodeFromName(tb, 'EquationDimension') + if dimPb is None: raise ValueError('prepareIBMDataPara: EquationDimension is missing in input geometry tree.') + dimPb = Internal.getValue(dimPb) + if dimPb == 2: C._initVars(tb, 'CoordinateZ', 0.) + + model = Internal.getNodeFromName(tb, 'GoverningEquations') + if model is None: raise ValueError('prepareIBMDataPara: GoverningEquations is missing in input geometry tree.') + model = Internal.getValue(model) + + ibctypes = Internal.getNodesFromName(tb, 'ibctype') + if ibctypes is None: raise ValueError('prepareIBMDataPara: ibc type is missing in input geometry tree.') + ibctypes = list(set(Internal.getValue(ibc) for ibc in ibctypes)) + + if model == 'Euler': + if any(ibc in ['Musker', 'MuskerMob', 'Mafzal', 'Log', 'TBLE', 'TBLE_FULL'] for ibc in ibctypes): + raise ValueError("prepareIBMDataPara: governing equations (Euler) not consistent with ibc types %s"%(ibctypes)) + + #=================== + # STEP 0 : GET FILAMENT BODIES Modified + #=================== + [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]) + + # Modification needed for Extrude + # Keep F1 regions - for F1 & F42 synergy + if tbF1: + tbbBTmp = G.BB(tbF1) + interDict_scale = X.getIntersectingDomains(tbbBTmp, t) + for kk in interDict_scale: + for kkk in interDict_scale[kk]: + z=Internal.getNodeFromName(t, kkk) + Internal._createUniqueChild(z, '.Solver#defineTMP', 'UserDefinedData_t') + Internal._createUniqueChild(Internal.getNodeFromName1(z, '.Solver#defineTMP'), 'SaveF1', 'DataArray_t', value=1) + node=Internal.getNodeFromName(t, kkk) + + ## New to extrude Prep + if extrusion=='cyl': cartesian = False + + _redispatch__(t=t) + + C._initVars(t, 'centers:cellN', 1.) + X._applyBCOverlaps(t, depth=depth, loc='centers', val=2, cellNName='cellN') + C._initVars(t,'{centers:cellNChim}={centers:cellN}') + C._initVars(t, 'centers:cellN', 1.) + + #=================== + # STEP 1 : BLANKING IBM Modified + #=================== + if verbose: pt0 = python_time.time(); printTimeAndMemory__('blank by IBC bodies', time=-1, functionName='setInterpDataIBMExtrude') + _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) + + ##Modify needed for extrude & to replicate previous behavior + ##Ghost kmin et kmax donneuse potentiel + listvars_local =['cellNChim','cellNIBC'] + for z in Internal.getZones(t): + sol = Internal.getNodeFromName(z,'FlowSolution#Centers') + for var in listvars_local: + cellN = Internal.getNodeFromName(sol,var)[1] + sh = numpy.shape(cellN) + for k in [0,1, sh[2]-2, sh[2]-1]: + for j in range(sh[1]): + for i in range(sh[0]): + if cellN[i,j,k] != 0: cellN[i,j,k] =1 + C._initVars(t,'{centers:cellN}=maximum(0.,{centers:cellNChim})') # vaut -3, 0, 1, 2 initialement + + Cmpi.barrier() + _redispatch__(t=t) + if verbose: printTimeAndMemory__('blank by IBC bodies', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') + + ##Same as prepareIBMDataPara + #=================== + # STEP 2 : INTERP DATA CHIM + #=================== + ## cellN mush be correct here + ## _setInterpData uses cellN + if verbose: pt0 = python_time.time(); printTimeAndMemory__('compute interpolation data (Abutting & Chimera)', time=-1, functionName='setInterpDataIBMExtrude') + tc = C.node2Center(t) + + if Internal.getNodeFromType(t, "GridConnectivity1to1_t") is not None: + Xmpi._setInterpData(t, tc, nature=1, loc='centers', storage='inverse', sameName=1, dim=dimPb, itype='abutting', order=2, cartesian=cartesian) + Xmpi._setInterpData(t, tc, nature=1, loc='centers', storage='inverse', sameName=1, sameBase=1, dim=dimPb, itype='chimera', order=2, cartesian=cartesian) + if verbose: printTimeAndMemory__('compute interpolation data (Abutting & Chimera)', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') + + #=================== + # STEP 3 : BUILD FRONT + #=================== + if verbose: pt0 = python_time.time(); printTimeAndMemory__('build IBM front', time=-1, functionName='setInterpDataIBMExtrude') + t, tc, front, front2, frontWMM = buildFrontIBM(t, tc, dimPb=dimPb, frontType=frontType, + cartesian=cartesian, twoFronts=twoFronts, check=check, + isWireModel=isWireModel) + if verbose: printTimeAndMemory__('build IBM front', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') + + #=================== + # STEP 4 : INTERP DATA IBM + #=================== + if verbose: pt0 = python_time.time(); printTimeAndMemory__('compute interpolation data (IBM)', time=-1, functionName='setInterpDataIBMExtrude') + _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) + if verbose: printTimeAndMemory__('compute interpolation data (IBM)', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') + + + #=================== + # STEP 5 : INIT IBM Modified + #=================== + tc2 = Internal.copyTree(tc) if twoFronts or isWireModel else None + + if isWireModel: + Internal._rmNodesByName(tc2, 'IBCD*') + Internal._rmNodesByName(tc, '2_IBCD*') + D_IBM._transformTc2(tc2) + tc2 = Internal.rmNodesByName(tc2, 'ID*') + NewIBCD= 141 + _changeNameIBCD__(tc2,NewIBCD) + tc = Internal.merge([tc,tc2]) + tc2 = None + + ## New to extrude Prep + if extrusion == 'cyl': + T._cyl2Cart(t, (0,0,0),(1,0,0)) + T._cyl2Cart(tc,(0,0,0),(1,0,0)) + + # modif info mesh in zonesubregion_t + for z in Internal.getZones(tc): + for zsr in Internal.getNodesFromType(z, "ZoneSubRegion_t"): + zsrname = Internal.getName(zsr) + zsrname = zsrname.split('_') + if zsrname[0]=='IBCD': + for var in ['C','W','I']: + r = Internal.getNodeFromName(zsr,'CoordinateY_P'+var)[1] + theta = Internal.getNodeFromName(zsr,'CoordinateZ_P'+var)[1] + for l in range(numpy.size(r)): + yy = r[l]*numpy.cos( theta[l] ) + zz = r[l]*numpy.sin( theta[l] ) + r[l]= yy; theta[l] = zz + + if distribute and Cmpi.size > 1: _redispatch__(t=t, tc=tc, tc2=tc2, twoFronts=twoFronts) + + ## New to extrude Prep + vars = ['centers:TurbulentDistanceAllBC','centers:TurbulentDistanceWallBC', 'centers:cellNIBC_hole'] + C._rmVars(t, vars) + + if isinstance(tc_out, str): + tcp = Compressor.compressCartesian(tc) + Cmpi.convertPyTree2File(tcp, tc_out, ignoreProcNodes=True) + + if twoFronts: + tcp2 = Compressor.compressCartesian(tc2) + tc2_out = tc_out.replace('tc', 'tc2') if 'tc' in tc_out else 'tc2.cgns' + Cmpi.convertPyTree2File(tcp2, tc2_out, ignoreProcNodes=True) + + if isinstance(t_out, str): + tp = Compressor.compressCartesian(t) + Cmpi.convertPyTree2File(tp, t_out, ignoreProcNodes=True) + + _computeMeshInfo(t) + + if Cmpi.size > 1: Cmpi.barrier() + if verbose: printTimeAndMemory__('initialize and clean', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') + + if tc2 is not None: return t, tc, tc2 + else: return t, tc + +def setInterpDataHybrid(t_3d, tc_3d, t_curvi, extrusion=None, interpDataType=1): + import FastS.ToolboxChimera as TBX + #================================================================================ + # IBM prepare hybride cart + curvi + # IN: t_3d arbre cartesien (from prepare1) + # IN: tc_3d arbre cartesien connectivite (from prepare1) + # IN: t_curvi arbre curviligne (with bc and connectivity but without ghost) + #================================================================================ + overlap = '14' # Overlap IBC + InterpOrder = 2 + root_racHybride = 'joinIBC' + + #add 2 ghost cells to the curvilinear mesh + Internal._addGhostCells(t_curvi, t_curvi, 2, adaptBCs=1, fillCorner=0) + C._initVars(t_curvi,"centers:cellN",1.) + t_curvi = TBX.modifyBCOverlapsForGhostMesh(t_curvi,2) + + dimPb = Internal.getValue(Internal.getNodeFromName(t_3d, 'EquationDimension')) + + for z in Internal.getZones(t_curvi): + if root_racHybride in z[0]: + C._fillEmptyBCWith(z,'inactive','BCExtrapolated',dim=dimPb) + C._rmBCOfType(z,'BCOverlap') + C._fillEmptyBCWith(z,'trou','BCOverlap',dim=dimPb) + X._applyBCOverlaps(z,loc='centers',depth=2, val=0) + X._applyBCOverlaps(z,loc='centers',depth=4, val=2) + + if extrusion == 'cyl': + T._cart2Cyl(t_3d , (0,0,0),(1,0,0)) + T._cart2Cyl(tc_3d, (0,0,0),(1,0,0)) + T._cart2Cyl(t_curvi , (0,0,0),(1,0,0)) + + tBB2 = G.BB(t_curvi) + tBB = G.BB(t_3d) + + tc2 = C.node2Center(t_curvi) + + X._setInterpData(t_curvi, tc2, nature=1, loc='centers', storage='inverse', sameName=1, dim=dimPb, itype='abutting') + + C._initVars(t_3d,"{centers:cellN#Init}={centers:cellN}") + + for var in ['wall','racChimer']: + C._initVars(t_3d,"centers:cellN",1.) + if var == 'wall': itype ="3" + else: itype = overlap + for zc in Internal.getZones(tc_3d): + for zsr in Internal.getNodesFromType(zc, "ZoneSubRegion_t"): + zsrname = Internal.getName(zsr) + zsrname = zsrname.split('_') + if zsrname[0]=='IBCD' and zsrname[1] == itype: + zrname = Internal.getValue(zsr) + ptlistD= Internal.getNodeFromName(zsr,'PointListDonor')[1] + zloc = Internal.getNodeFromName(t_3d,zrname) + sol = Internal.getNodeFromName(zloc,'FlowSolution#Centers') + cellN = Internal.getNodeFromName(sol,'cellN')[1] + sh = numpy.shape(cellN) + ni= sh[0]; ninj = sh[0]*sh[1] + for l0 in range(numpy.size(ptlistD)): + l = ptlistD[ l0] + k = l//ninj + j = (l-k*ninj)//ni + i = l-k*ninj - j*ni + cellN[ i,j,k ]=2 + + if var == 'wall': C._initVars(t_3d,"{centers:cellN#wall}={centers:cellN}") + else: C._initVars(t_3d,"{centers:cellN#racChim}={centers:cellN}") + + Internal._rmNodesFromName(tc_3d,'IBCD_'+overlap+'_*') + + for z in Internal.getZones(t_3d): + sol = Internal.getNodeFromName(z,'FlowSolution#Centers') + cellNRac = Internal.getNodeFromName(sol,'cellN#racChim')[1] + cellN = Internal.getNodeFromName(sol,'cellN')[1] + C._initVars(z,"{centers:cellN}={centers:cellN#racChim}") + sh = numpy.shape(cellN) + + if sh[2]> 1: + for k in [0,1, sh[2]-2, sh[2]-1]: + for j in range(sh[1]): + for i in range(sh[0]): + if cellN[i,j,k] != 0: cellN[i,j,k] =1 + + interDict = X.getIntersectingDomains(tBB, t2=tBB2, method='AABB', taabb=tBB, taabb2=tBB2) + print(" Interp Cartesian from curvilinear") + zonelist=[] + for z in Internal.getZones(t_3d): + if C.getMaxValue(z,'centers:cellN')==2: + dnrZones = [] + for zdname in interDict[z[0]]: + zd = Internal.getNodeFromName(tc2,zdname) + dnrZones.append(zd) + X._setInterpData(z,dnrZones, nature=0,penalty=1,loc='centers',storage='inverse',sameName=1,\ + interpDataType=interpDataType, itype='chimera', order=InterpOrder) + z = X.getOversetInfo(z, dnrZones, loc='center',type='extrapolated') + zonelist.append(z) + + interDict_curvi = X.getIntersectingDomains(tBB2, t2=tBB, method='AABB', taabb=tBB2, taabb2=tBB) + print(" Interp curvi from Cartesian") + for z in Internal.getZones(t_curvi): + if C.getMaxValue(z,'centers:cellN')==2: + dnrZones = [] + for zdname in interDict_curvi[z[0]]: + zd = Internal.getNodeFromName(tc_3d,zdname) + dnrZones.append(zd) + + X._setInterpData(z,dnrZones,nature=0,penalty=1,loc='centers',storage='inverse',sameName=1,\ + interpDataType=interpDataType, itype='chimera', order=InterpOrder) + + # to check orphans + #z = X.getOversetInfo(z, dnrZones, loc='center',type='extrapolated') + #z = X.getOversetInfo(z, dnrZones, loc='center',type='orphan') + #zonelist.append(z) + #C.convertPyTree2File(z,"rcv2.cgns") + + C._initVars(t_3d, '{centers:cellN}={centers:cellN#Init}') + + if extrusion == 'cyl': + T._cyl2Cart(t_3d, (0,0,0),(1,0,0)) + T._cyl2Cart(t_curvi, (0,0,0),(1,0,0)) + T._cyl2Cart(tc_3d, (0,0,0),(1,0,0)) + T._cyl2Cart(tc2, (0,0,0),(1,0,0)) + + for z in Internal.getZones(t_curvi): + for bc in Internal.getNodesFromType(z,'BC_t'): + if 'inactive' in bc[0] and Internal.getValue(bc) == 'BCExtrapolated': + Internal._rmNodesByName(z, bc[0]) + + t = C.mergeTrees(t_3d ,t_curvi ) + tc = C.mergeTrees(tc_3d,tc2) + + return t, tc + def _redispatch__(t=None, tc=None, tc2=None): import Distributor2.Mpi as D2mpi @@ -537,7 +886,6 @@ def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1. # IN: wallAdaptF42 (cloud of IB target points with yplus information): # use a previous computation to adapt the positioning of IB target points around the geometry according to a target yplus (F42) # IN: heightMaxF42 (float): maximum modeling height for the location of IB target points around the geometry (F42) -# IN: listF1save: list of zones that will have an F1 front # IN: filamentBases: list of bases that are IBC filaments # IN: isFilamentOnly: boolean on whether there is only a IBC filament # IN: tbFilament: PyTree of the IBC filaments @@ -831,7 +1179,6 @@ def _blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, Internal.__FlowSolutionCenters__) C._initVars(t,'{centers:cellN}=maximum(0.,{centers:cellNChim})') # vaut -3, 0, 1, 2 initialement - return None #========================================================================= diff --git a/Cassiopee/Connector/test/IBMExtrude2Dto3D_t1.py b/Cassiopee/Connector/test/IBMExtrude2Dto3D_t1.py new file mode 100644 index 000000000..c1ba84205 --- /dev/null +++ b/Cassiopee/Connector/test/IBMExtrude2Dto3D_t1.py @@ -0,0 +1,51 @@ +# extrude 2D mesh to 3D with Cartesian approach +import Converter.PyTree as C +import Converter.Mpi as Cmpi +import Transform.PyTree as T +import Converter.Internal as Internal +import Generator.IBM as G_IBM +import Geom.IBM as Geom +import Connector.IBM as X_IBM +import KCore.test as test +import sys + +LOCAL = test.getLocal() +bodySurfaceFile = LOCAL + '/naca0012.cgns' + +Lcharac = 0.03362355 +Lz = 0.2*Lcharac +dfar = 10.*Lcharac +vmin = 16; +snears = 1; + +###Generate 2D Mesh +t2D,tc2D=X_IBM.prepareIBMDataPara(bodySurfaceFile , None , None , + snears=snears , dfar=dfar , vmin=vmin ) +test.testT(t2D ,1) +test.testT(tc2D,2) +#C.convertPyTree2File(t2D,LOCAL+'/t2D_checking.cgns') + +####Extrusion for 3D Mesh +bodySurface = C.convertFile2PyTree(bodySurfaceFile) + +extrusion = 'cart' +span = Lz +NPas = 10+1 #number of nodes +t3D, tb3D = G_IBM.extrudeCartesianZDir(t2D, bodySurface, extrusion=extrusion, NPas=NPas, span=span, dz=span/(NPas-1), isAutoPeriodic=True) + +for t in [t3D,tb3D]: + zmax = C.getMaxValue(t, 'CoordinateZ'); + zmin = C.getMinValue(t, 'CoordinateZ'); + zavg = (zmax+zmin)/2 + T._translate(t, (0,0,0-zavg)) +test.testT(t3D ,3) +test.testT(tb3D ,4) + +#C.convertPyTree2File(t3D,LOCAL+'/t3D_checking.cgns') +#C.convertPyTree2File(tb3D,LOCAL+'/tb3D_checking.cgns') + +#####Interpolation 3D +t3D, tc3D = X_IBM.prepareIBMDataExtrude(tb3D, None, None, t3D, extrusion=extrusion) +test.testT(t3D ,5) +test.testT(tc3D ,6) + diff --git a/Cassiopee/Generator/Generator/IBM.py b/Cassiopee/Generator/Generator/IBM.py index 6ce314525..5bf314dee 100644 --- a/Cassiopee/Generator/Generator/IBM.py +++ b/Cassiopee/Generator/Generator/IBM.py @@ -1118,3 +1118,191 @@ def _projectMeshSize(t, NPas=10, span=1, dictNz=None, isCartesianExtrude=False): 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 + +def extrudeCartesianZDir(t, tb, check=False, extrusion="cart", dz=0.01, NPas=10, span=1 , Ntranche=1, + dictNz=None, ific=2, isCartesianExtrude=False, isAutoPeriodic=False, nghost=0): + """Extrudes a 2D IBM grid and geoemtry. The extraction is done in the z-direction. + Usage: extrudeCartesianZDir(t, tb, check, extrusion, dz, NPas, span, Ntranche, dictNz, ific, isCartesianExtrude, isAutoPeriodic, nghost)""" + + ## isAutoPeriodic=True is prefered over isAutoPeriodic=False when extruding a Cartesian mesh. + ## nghost should be 2. It is currently 0 (default value) to pass the non-regression tests. This function was initially developed + ## an application and its default value of 0 is kept for reproduciability reasons for that original test case. nghost is only used when isAutoPeriodic=False. + ## isAutoPeridioc=True assumes 2 ghost cells. + + if extrusion == "cart": perio = span/float(Ntranche) + else: perio = span/180.*math.pi/float(Ntranche) + + #for z in Internal.getZones(t): + # cellN = Internal.getNodeFromName(z, "cellN")[1] + # sh = numpy.shape(cellN) + # # modif CellN pour filtrer les cellule solide de l'interp chimere + # # et autoriser les ghost comme donneuse + # for k in range(sh[2]): + # for j in range(sh[1]): + # for i in range(sh[0]): + # if cellN[i,j,k] != 0: cellN[i,j,k] =1 + + #Dim = 3D + c = 0 + for tree in [t,tb]: + for dim in Internal.getNodesFromName(tree,'EquationDimension'): dim[1]=3 + + dz_loc={}; Nk={} + if c ==0: + Nk_min = int(10e09) + for z in Internal.getZones(tree): + name_zone = z[0] + if not isCartesianExtrude: + if dictNz is not None: + Nk[z[0]] = int(dictNz[name_zone]) + dz_loc[z[0]] = span/float(Ntranche*Nk[z[0]]) + else: + Nk[z[0]] = NPas-1 + dz_loc[z[0]] = dz + else: + h = abs(C.getValue(z,'CoordinateX',0)-C.getValue(z,'CoordinateX',1)) + NPas_local = int(round(span/h)) + if NPas_local<4: + print("WARNING:: Zone %s has Nz=%d and is being clipped to Nz=4"%(z[0],NPas_local)) + NPas_local=4 + Nk[z[0]] = NPas_local + dz_loc[z[0]] = span/float(Ntranche*Nk[z[0]]) + if Nk[z[0]] < Nk_min: Nk_min =Nk[z[0]] + else: + for z in Internal.getZones(tree): + Nk[z[0]] = Nk_min + dz_loc[z[0]] = span/float(Ntranche*Nk_min) + + ##Adding ghost cells in K direction + for z in Internal.getZones(tree): + Nk[z[0]] += 2*ific-1+c # -1, for the mesh (already added in the mesh generation) || need to remove this as it is not the case for the geometry (tb) + + for z in Internal.getZones(tree): + yy_2d = Internal.getNodeFromName(z, "CoordinateY")[1] + zz_2d = Internal.getNodeFromName(z, "CoordinateZ")[1] + sh_2d = numpy.shape(yy_2d) + + T._addkplane(z ,N=Nk[z[0]]) + + zz = Internal.getNodeFromName(z, "CoordinateZ")[1] + yy = Internal.getNodeFromName(z, "CoordinateY")[1] + sh = numpy.shape(yy) + r = numpy.sqrt( zz**2+yy**2) + + perio_loc = perio + if c==1: period_loc= perio*1.5 #is periodic_loc a typo or a new variable that is not used anywhere else? + + if len(sh_2d) == 1: #NGON + if extrusion == "cart": + for k in range( sh[1]): + zz[:,k] = (k-ific)* dz_loc[z[0]] + else: + theta0 = numpy.arctan(zz_2d/yy_2d) + for k in range( sh[1]): + shift = perio_loc/float(sh[1]-5.)*(k-ific) + zz[:,k] = r[:,0]*numpy.sin(theta0[:] + shift) + yy[:,k] = r[:,0]*numpy.cos(theta0[:] + shift) + else: + if extrusion == "cart": + for k in range(sh[2]): + zz[:,:,k] = (k-ific)* dz_loc[z[0]] + else: + theta0 = numpy.arctan(zz_2d/yy_2d) + for k in range(sh[2]): + shift = perio_loc/float(sh[2]-5.)*(k-ific) + zz[:,:,k] = r[:,:,0]*numpy.sin(theta0[:,:,0] + shift) + yy[:,:,k] = r[:,:,0]*numpy.cos(theta0[:,:,0] + shift) + + if c==1: break + c += 1 + for z in Internal.getZones(tree): + zdim = Internal.getValue(z) + + # Modif rind cellule GH en kmin et kmax + rind = Internal.getNodeFromName1(z, "Rind")[1] + rind[4] = ific; rind[5] = ific + # Modif range BC + BCs = Internal.getNodesFromType(z, "BC_t") + for bc in BCs: + ptrg = Internal.getNodeFromName(bc, "PointRange")[1] + ptrg[2,0] = 3 + ptrg[2,1] = Nk[z[0]] + + ##Periodic boundary conditions in k direction + if not isAutoPeriodic: + # Creatioon connectivite perio dans t + for idir in ['_kmax','_kmin']: + if idir == '_kmax': + ktg = zdim[2,0] + ktgD = 1 + if extrusion =='cart': + angle = 0 + trans = -perio + else: + angle = -perio + trans = 0. + else: + ktgD = zdim[2,0] + ktg = 1 + if extrusion =='cart': + angle = 0 + trans = perio + else: + angle = perio + trans = 0. + + Conn = Internal.getNodeFromName(z, "ZoneGridConnectivity") + name = 'match_'+z[0]+idir + Internal.createUniqueChild(Conn, name, 'GridConnectivity1to1_t') + tmp1 = Internal.getNodeFromName(Conn, name) + Internal.setValue(tmp1, z[0]) + datap = numpy.empty( (3,2) , numpy.int32) + datap[0,0]=1+nghost;datap[1,0]=1+nghost;datap[2,0]=ktg + datap[0,1]=zdim[0,0]-nghost;datap[1,1]=zdim[1,0]-nghost;datap[2,1]= ktg + Internal.createUniqueChild(tmp1 ,'PointRange', 'IndexRange_t',datap) + datap = numpy.empty( (3,2) , numpy.int32) + datap[0,0]=1+nghost;datap[1,0]=1+nghost;datap[2,0]=ktgD + datap[0,1]=zdim[0,0]-nghost;datap[1,1]=zdim[1,0]-nghost;datap[2,1]= ktgD + Internal.createUniqueChild(tmp1 ,'PointRangeDonor', 'IndexRange_t',datap) + datap = numpy.empty( 3 , numpy.int32) + datap[0]=1;datap[1]=2;datap[2]=3 + Internal.createUniqueChild(tmp1 ,'Transform', '"int[IndexDimension]"',datap) + Internal.createUniqueChild(tmp1 ,'GridConnectivityProperty', 'GridConnectivityProperty_t') + + prop = Internal.getNodeFromName(tmp1, 'GridConnectivityProperty') + Internal.createUniqueChild(prop ,'Periodic', 'Periodic_t') + period = Internal.getNodeFromName(prop, 'Periodic') + datap = numpy.zeros( 3 , numpy.float64) + Internal.createUniqueChild(period ,'RotationCenter', 'DataArray_t',datap) + datap = numpy.zeros( 3 , numpy.float64) + datap[0]= angle + Internal.createUniqueChild(period ,'RotationAngle' ,'DataArray_t',datap) + datap = numpy.zeros( 3 , numpy.float64) + datap[2]= trans + Internal.createUniqueChild(period ,'Translation', 'DataArray_t',datap) + + rot = Internal.getNodeFromName(period, 'RotationAngle') + Units=['Kilogram','Meter','Second','Kelvin','Radian'] + Internal.createUniqueChild(rot ,'DimensionalUnits' ,'DimensionalUnits_t',Units) + + if isAutoPeriodic: + for node in Internal.getNodesFromName(t,'EquationDimension'): Internal.setValue(node,3) + for z in Internal.getZones(t): + C._addBC2Zone(z, z[0]+'periodKmin', 'BCautoperiod', 'kmin') + C._addBC2Zone(z, z[0]+'periodKmin', 'BCautoperiod', 'kmax') + BCs = Internal.getNodesFromType(t, "BC_t") + for bc in BCs: + if Internal.getValue(bc)=='BCautoperiod': + ptrg = Internal.getNodeFromName(bc, "PointRange")[1] + ptrg[0,0] = 3 + ptrg[0,1] = ptrg[0,1]-2 + + ptrg[1,0] = 3 + ptrg[1,1] = ptrg[1,1]-2 + + if extrusion == 'cyl': + T._cart2Cyl(t, (0,0,0),(1,0,0)) + T._cart2Cyl(tb, (0,0,0),(1,0,0)) + + return t, tb + From db6f8d11915fc2beff23032df8197f6236bc843e Mon Sep 17 00:00:00 2001 From: antjost Date: Thu, 29 Aug 2024 10:42:35 +0200 Subject: [PATCH 2/2] [IBM] Test case for Cartesian grid extrusion in Z direction. Splitting of previous test case into an extrude and an extrude and interpolate. The latter is need to make the interpolation test case self contained. --- Cassiopee/Apps/Apps/Fast/IBM.py | 2 +- Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py | 131 ----- Cassiopee/Apps/test/IBMrotationCompute_t1.py | 211 ------- Cassiopee/Apps/test/IBMrotation_t1.py | 4 - Cassiopee/Connector/Connector/IBM.py | 550 +++++++----------- ...BeCompletedIBMinterpolateExtrude2Dto3D.py} | 3 + Cassiopee/Generator/Generator/IBM.py | 11 +- .../Generator/test/IBMExtrude2Dto3D_t1.py | 67 +++ .../Generator/test/IBMExtrude2Dto3D_t2.py | 67 +++ 9 files changed, 353 insertions(+), 693 deletions(-) delete mode 100644 Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py delete mode 100644 Cassiopee/Apps/test/IBMrotationCompute_t1.py rename Cassiopee/Connector/test/{IBMExtrude2Dto3D_t1.py => ToBeCompletedIBMinterpolateExtrude2Dto3D.py} (95%) create mode 100644 Cassiopee/Generator/test/IBMExtrude2Dto3D_t1.py create mode 100644 Cassiopee/Generator/test/IBMExtrude2Dto3D_t2.py diff --git a/Cassiopee/Apps/Apps/Fast/IBM.py b/Cassiopee/Apps/Apps/Fast/IBM.py index eb8c99eb4..4a5303de1 100644 --- a/Cassiopee/Apps/Apps/Fast/IBM.py +++ b/Cassiopee/Apps/Apps/Fast/IBM.py @@ -614,7 +614,7 @@ def extrudeCartesian(t,tb, check=False, extrusion="cart", dz=0.01, NPas=10, span for node in Internal.getNodesFromName(t,'EquationDimension'): Internal.setValue(node,3) for z in Internal.getZones(t): C._addBC2Zone(z, z[0]+'periodKmin', 'BCautoperiod', 'kmin') - C._addBC2Zone(z, z[0]+'periodKmin', 'BCautoperiod', 'kmax') + C._addBC2Zone(z, z[0]+'periodKmax', 'BCautoperiod', 'kmax') BCs = Internal.getNodesFromType(t, "BC_t") for bc in BCs: if Internal.getValue(bc)=='BCautoperiod': diff --git a/Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py b/Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py deleted file mode 100644 index 33f14cb78..000000000 --- a/Cassiopee/Apps/test/FastIBMHybrideCurvi_t2.py +++ /dev/null @@ -1,131 +0,0 @@ -# - Fast.IBM - -import Apps.Fast.IBM as App -import Connector.PyTree as X -import Connector.IBM as X_IBM -import Converter.Internal as Internal -import Converter.PyTree as C -import Fast.PyTree as Fast -import FastS.PyTree as FastS -import Generator.PyTree as G -import Geom.IBM as D_IBM -import Geom.PyTree as D -import Initiator.PyTree as Initiator -import KCore.test as test -import math -import numpy as np -test.TOLERANCE = 1.e-6 - -LOCAL = test.getLocal() - -#Fabrication case 2d pour maillage octree -lines=[] -lines.append( D.line((0.05,-0.15, 0), (0.1 ,-0.08,0), N=10) ) -lines.append( D.line((0.10,-0.03, 0), (0.1 ,-0.08,0), N=10) ) -lines.append( D.line((0.1 ,-0.03, 0), (0.1 , 0. ,0), N=8) ) - -lines.append( D.line((0.15,-0.08, 0), (0.15,-0.03,0), N=10) ) -lines.append( D.line((0.15,-0.08, 0), (0.2 ,-0.15,0), N=10) ) -lines.append( D.line((0.05,-0.15, 0), (0.2 ,-0.15,0), N=10) ) - -lines.append( D.line((0.15,-0.03, 0), (0.15, 0 ,0), N=8 ) ) -lines.append( D.line((0.10, 0.01, 0), (0.15, 0.01,0), N=10) ) -lines.append( D.line((0.15, 0.01, 0), (0.15, 0 ,0), N=3) ) -lines.append( D.line((0.10, 0.01, 0), (0.10, 0 ,0), N=3) ) -case = C.newPyTree(['BOBY', lines]) - -D_IBM._setSnear(case, 0.002) -D_IBM._setIBCType(case, "Musker") -D_IBM._setDfar(case, 0.) -D_IBM._setFluidInside(case) - -zones = Internal.getZones(case) -for z in zones: - if '8' in z[0] or '9' in z[0] or '6' in z[0]: - D_IBM._setIBCType(z, "overlap") - D_IBM._setSnear(z, 0.008) - if 'line.1' == z[0] or '10' in z[0]: - D_IBM._setSnear(z, 0.008) - if '2' in z[0] or '0' in z[0]: - D_IBM._setSnear(z, 0.004) - -zones[6][0]='line10' -zones[8][0]='line9' -U0= 120.0 -P0= 101325.0 -L0= 1. -R0= 1.18 -T0= 279.15 - -equation = 'NSLaminar' -C._addState(case, 'GoverningEquations', equation ) -C._addState(case, 'EquationDimension', 2) -C._addState(case, UInf=U0, RoInf=R0, PInf=P0, LInf=L0, alphaZ=0., adim='dim3') -#C.convertPyTree2File(case,'verifcase.cgns') - -# Prepare cas 2d IBM -t_2d, tc_2d = App.prepare1(case, None, None, frontType=42 ,cleanCellN=False) - -test.testT(tc_2d, 1) -test.testT(t_2d , 2) - -#extrusion 3D IBM -extrusion ='cart'; span = 0.01 - -t_3d, tb_3d = App.extrudeCartesian(t_2d, case, extrusion=extrusion, NPas=5, span=span) - -test.testT(tb_3d, 3) -test.testT(t_3d , 4) - -#calcul interp IBM 3D; celln, distnce paroi 3D,.... deja obtenue lors de l'extrusion -interpDataType = 1 # on suppose maillage non cartesion pour interp -order = 2 -t_3d, tc_3d = App.prepare1(tb_3d,None, None, t_in=t_3d, extrusion=extrusion, interpDataType=interpDataType, order=order) - -Internal._rmNodesFromType(tc_2d,'Rind_t') -Internal._rmNodesFromName(tc_2d,Internal.__GridCoordinates__) -test.testT(tc_3d, 5) -test.testT(t_3d , 6) - - -#Maillage plaque plane Curviligne sans ghost -a = G.cart((0,0,0), (0.005,0.005,0.01), (200,100,5)) -C._addBC2Zone(a, 'farfield', 'BCFarfield', 'jmax') -C._addBC2Zone(a, 'inflow', 'BCInflow', 'imin') -C._addBC2Zone(a, 'outflow','BCOutflow', 'imax') -C._addBC2Zone(a, 'wallIn','BCWallViscous', [1,22,1,1,1,5]) -C._addBC2Zone(a, 'overlap','BCOverlap', [22,30,1,1,1,5]) -C._addBC2Zone(a, 'wallOut','BCWallViscous', [30,200,1,1,1,5]) -t_curvi = C.newPyTree(['Base', a]) - -zones = Internal.getZones(t_curvi) -#les zones curviligne possedanr raccord chimere avec zone Cart IBC doit avoir la racine "joinIBC" dans leur nom -zones[0][0]='curvi_joinIBC' -t_curvi = X.connectMatchPeriodic(t_curvi, translation=[0.,0.,0.04]) -#stretch maillage plaque direction normal paroi -for z in zones: - coordy = Internal.getNodeFromName(z,'CoordinateY')[1] - sh = np.shape(coordy) - for j in range(1,sh[1]): - coordy[:,j,:]=coordy[:,j-1,:]+0.002*1.02**j - -C._addState(t_curvi, 'GoverningEquations', equation ) -C._addState(t_curvi, 'EquationDimension', 3) -C._addState(t_curvi, UInf=U0, RoInf=R0, PInf=P0, LInf=L0, alphaZ=0., adim='dim3') -C._initVars(t_curvi,"{centers:Density}=1.18") -C._initVars(t_curvi,"{centers:VelocityX}=0.01") -C._initVars(t_curvi,"{centers:VelocityY}=0") -C._initVars(t_curvi,"{centers:VelocityZ}=0") -C._initVars(t_curvi,"{centers:Temperature}=279.15") - -C._initVars(t_3d,"{centers:Density}=1.18") -C._initVars(t_3d,"{centers:VelocityX}=0.01") -C._initVars(t_3d,"{centers:VelocityY}=0") -C._initVars(t_3d,"{centers:VelocityZ}=0") -C._initVars(t_3d,"{centers:Temperature}=279.15") - -#Calcul interpolation entre maillage curviligne et cartesien IBC et ajout Ghost au maillage curvi -#t_final, tc_final= App.setInterpData_Hybride(t_3d, tc_3d, t_curvi) -t_final, tc_final= X_IBM.setInterpDataHybrid(t_3d, tc_3d, t_curvi) -test.testT(t_final , 7) -test.testT(tc_final, 8) - diff --git a/Cassiopee/Apps/test/IBMrotationCompute_t1.py b/Cassiopee/Apps/test/IBMrotationCompute_t1.py deleted file mode 100644 index 2ad88adad..000000000 --- a/Cassiopee/Apps/test/IBMrotationCompute_t1.py +++ /dev/null @@ -1,211 +0,0 @@ -# case - cylinder - IBM rotation - compute -import Connector.IBM as X_IBM -import Connector.Mpi as Xmpi -import Connector.PyTree as X -import Converter.Internal as Internal -import Converter.Mpi as Cmpi -import Converter.PyTree as C -import Fast.PyTree as Fast -import FastS.PyTree as FastS -import RigidMotion.PyTree as R -import Transform.PyTree as T -import KCore.test as test -import math -import numpy - -LOCAL = test.getLocal() - -## ======================= -## ======================= -TInf = 298.15 -UInf = 0.0001 -rad = 0.0005 -diam = rad*2 -RoInf= 1.18 -PInf = 101325 - -# necessaire pour les IBM -Model = 'NSLaminar' -dimPb = 2 -modulo_verif = 5 - -rpm = 100 -OMG = rpm*2*math.pi/60 -LInf = diam*0.5 -TInf = 298.15 -aref = math.sqrt(1.4*287.058*298.15) -Vtip = OMG*LInf -Mtip = OMG*LInf/aref -Re = RoInf*Vtip*diam/1.85e-05 -print(OMG,LInf,Vtip,Re) - -time_step= 1/(OMG*180/math.pi)/600 -## ======================= -## ======================= - -NIT = 15 - -numb={"temporal_scheme": "implicit_local", - "omp_mode":0, - "modulo_verif":modulo_verif} - -numz={"time_step": time_step, - "time_step_nature": "global", - "epsi_newton": 0.1} - -numz1 = numz.copy() -numz1['scheme'] = 'ausmpred' - -numz2 = numz.copy() -numz2['scheme'] = 'ausmpred' - -# load t, tc -fileTIn ='/t.cgns' -fileTcIn='/tc.cgns' -t,tc,ts,graph = Fast.load(LOCAL+fileTIn,LOCAL+fileTcIn) -baseNameIBM ='CARTESIAN_NEARBODY' -baseNameBKGD ='CARTESIAN_OFFBODY' - - -Fast._setNum2Base(t, numb) -for n in [baseNameIBM]: - b = Internal.getNodeFromName1(t, n) - Fast._setNum2Zones(b, numz1) -for n in [baseNameBKGD]: - b = Internal.getNodeFromName1(t, n) - Fast._setNum2Zones(b, numz2) - - -# load bodies -tb = C.convertFile2PyTree(LOCAL+'/bodiesBlank.cgns') - -# Warmup -it0 = 0; time0 = 0. -first = Internal.getNodeFromName1(t, 'Iteration') -if first is not None: it0 = Internal.getValue(first) -first = Internal.getNodeFromName1(t, 'Time') -if first is not None: time0 = Internal.getValue(first) -time_step = Internal.getNodeFromName(t, 'time_step') -time_step = Internal.getValue(time_step) - -# Tag les zones pour le motion -zones = Internal.getZones(t) -for z in zones: - timeMotion = Internal.getNodeFromName(z, 'TimeMotion') - if timeMotion is not None: - define = Internal.getNodeFromName1(z, '.Solver#define') - Internal._createUniqueChild(define, 'motion', 'DataArray_t', value='deformation') - -# Set interpolated points sur les bases IBM -C._initVars(t, "{centers:cellN#Motion}=1.") -##IBM CART -X._applyBCOverlaps(Internal.getNodeFromName1(t,baseNameIBM), depth=2, val=2, cellNName='cellN#Motion') -C._initVars(t, "{centers:cellN#MotionInit}={centers:cellN#Motion}") - -# Data for chimera motion -nbases = len(Internal.getBases(t)) -nbodies = len(Internal.getBases(tb)) -BM = numpy.zeros((nbases, nbodies),dtype=numpy.int32) -BM[1,0] = 1. - -# Creation BB tree commun -tBB = Cmpi.createBBoxTree(t) - -# Attachement des donneurs en dur (id on all procs) -intersectionDict={} -for z1 in Internal.getZones(Internal.getNodeFromName1(tBB,baseNameIBM)): - for z2 in Internal.getZones(Internal.getNodeFromName1(tBB,baseNameBKGD)): - Fast._addPair(intersectionDict, z1[0], z2[0]) - -procDict = None; -graphX = {} -procDict = Cmpi.getProcDict(tBB) -graphX = Cmpi.computeGraph(tBB, type='bbox2', t2=None, procDict=procDict, reduction=False,intersectionsDict=intersectionDict) - -R._evalPositionIBC(tc,time0) - -(t, tc, metrics) = FastS.warmup(t, tc, graph) - -# Doivent etre apres le warmup -dictOfADT={} -(dictOfNobOfRcvZones,dictOfNozOfRcvZones) = Fast.getDictOfNobNozOfRcvZones(t, intersectionDict) -(dictOfNobOfRcvZonesC,dictOfNozOfRcvZonesC) = Fast.getDictOfNobNozOfRcvZones(tc, intersectionDict) -(dictOfNobOfDnrZones,dictOfNozOfDnrZones) = Fast.getDictOfNobNozOfDnrZones(tc, intersectionDict, dictOfADT,cartFilter='CARTESIAN',isIbmAle=True) - -# modifie la vitesse de reference (si M=0) -cont = Internal.getNodeFromType(t, 'ReferenceState_t') -Minf = Internal.getNodeFromName(cont, 'Mach') -zones = Internal.getZones(t) -for z in zones: - n = Internal.getNodeFromName2(z, 'Parameter_real')[1] - n[5] = max(30, UInf) - -timeiter = time0 -varType = 0 - -# Get models -eqs = Internal.getNodeFromType(Internal.getNodeFromName1(t,baseNameIBM),'GoverningEquations_t') -Model1 = Internal.getValue(eqs) -eqs = Internal.getNodeFromType(Internal.getNodeFromName1(t, baseNameBKGD),'GoverningEquations_t') -Model2 = Internal.getValue(eqs) - -if Model1 == 'NSTurbulent' and Model2 == 'NSTurbulent': varType = 1 - -RefState = Internal.getNodeFromType(t,'ReferenceState_t') -VInf = Internal.getNodeFromName1(RefState,'VelocityX')[1][0] -PInf = Internal.getNodeFromName1(RefState,'Pressure')[1][0] -rhoInf = Internal.getNodeFromName1(RefState,'Density')[1][0] - -for it in range(it0,NIT+it0): - if Cmpi.rank == 0: print("it=%d, time=%f, angle=%f, N_rot=%d"%(it, timeiter, (timeiter*OMG*180./numpy.pi) % 360 , (timeiter*OMG*180./numpy.pi)//360 ),flush=True) - timeiter += time_step - - # bouge tout - R._evalPosition(tb, timeiter) - R._evalPosition(t, timeiter) - R._evalPosition(tc, timeiter) - R._evalGridSpeed(t, timeiter) - - R._evalPositionIBC(tc,timeiter) - - FastS.copy_velocity_ale(t, metrics, it=it) #get motion velocities @ face centers - - # Masquage - C._initVars(t,"{centers:cellN#Motion}={centers:cellN#MotionInit}") - bodies=[] - for base in Internal.getBases(tb): - bodies.append(Internal.getZones(base)) - - XRAYDIM1 = 2000;RAYDIM2 = XRAYDIM1 - if dimPb == 2: XRAYDIM2 = 2 - X._blankCells(t, bodies, BM, cellNName='cellN#Motion',XRaydim1=XRAYDIM1, XRaydim2=XRAYDIM2, blankingType='cell_intersect',dim=dimPb) - X._setHoleInterpolatedPoints(Internal.getNodeFromName1(t,baseNameBKGD), depth=2, loc='centers',addGC=False, cellNName='cellN#Motion', dir=2) - C._initVars(t, "{centers:cellN#Motion}=({centers:cellN#Static}<2)*{centers:cellN#Motion}+({centers:cellN#Static}>1)*minimum(1,{centers:cellN#Motion})") - C._initVars(t,"{centers:cellN}=minimum({centers:cellN#Motion}*{centers:cellN#Static},2.)") - C._cpVars(t, 'centers:cellN',tc, 'cellN') - ucData = (graphX, intersectionDict, dictOfADT, - dictOfNobOfRcvZones, dictOfNozOfRcvZones, - dictOfNobOfDnrZones, dictOfNozOfDnrZones, - dictOfNobOfRcvZonesC, dictOfNozOfRcvZonesC, - timeiter, procDict, True, varType, 1, 2, 1) - - FastS._compute(t, metrics, it, tc, graph, layer="Python", ucData=ucData) - - if it%modulo_verif==0: - FastS.display_temporal_criteria(t, metrics, it, format='double') - - Cmpi.barrier() - -for adt in dictOfADT.values(): - if adt is not None: C.freeHook(adt) - -Internal.createUniqueChild(t, 'Iteration', 'DataArray_t', value=NIT+it0) -Internal.createUniqueChild(t, 'Time', 'DataArray_t', value=timeiter) - -test.testT(t,1) -test.testT(tc,2) - - -#Cmpi.convertPyTree2File(t , LOCAL+'/t_restart.cgns') -#Cmpi.convertPyTree2File(tc, LOCAL+'/tc_restart.cgns') - diff --git a/Cassiopee/Apps/test/IBMrotation_t1.py b/Cassiopee/Apps/test/IBMrotation_t1.py index 0ae91fc61..78bd9437a 100644 --- a/Cassiopee/Apps/test/IBMrotation_t1.py +++ b/Cassiopee/Apps/test/IBMrotation_t1.py @@ -65,7 +65,6 @@ for z in Internal.getZones(tb): dfars.append(dfar_nb) t_ibm = G_IBM.generateIBMMesh(tb, vmin=vmin, dfars=dfars, dimPb=2, expand=2, ext=3) - R._setPrescribedMotion3(t_ibm ,'rot', axis_pnt=(0.,0.,0.), axis_vct=(0,0,1),omega=-OMG) t_ibm[2][1][0]='CARTESIAN_NEARBODY' @@ -73,7 +72,6 @@ tc_out=None t_ibm, tc_ibm = X_IBM.prepareIBMData(tb, t_out, tc_out, t_ibm, frontType=frontType) - C._rmBCOfType(t_ibm,'BCFarfield') C._fillEmptyBCWith(t_ibm,'dummy','BCExtrapolate', dim=dimPb) @@ -102,7 +100,6 @@ for z in Internal.getZones(tb_off): dfars.append(dfar_ext) t_off = G_IBM.generateIBMMesh(tb_off, vmin=vmin, dfars=dfars, dimPb=2, expand=2, ext=3) - t_off[2][1][0]='CARTESIAN_OFFBODY' X._applyBCOverlaps(t_off,depth=2,loc='centers') tc_off = C.node2Center(t_off) @@ -157,7 +154,6 @@ vars = ['Density', 'MomentumX', 'MomentumY', 'MomentumZ', 'EnergyStagnationDensity'] C._initVars(t_off, "centers:TurbulentDistance", 1e3) - t = C.mergeTrees(t_ibm, t_off) tc = C.mergeTrees(tc_ibm, tc_off) diff --git a/Cassiopee/Connector/Connector/IBM.py b/Cassiopee/Connector/Connector/IBM.py index 8557cf534..35a85d819 100644 --- a/Cassiopee/Connector/Connector/IBM.py +++ b/Cassiopee/Connector/Connector/IBM.py @@ -345,348 +345,214 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N if tc2 is not None: return t, tc, tc2 else: return t, tc -def prepareIBMDataExtrude(t_case, t_out, tc_out, t, to=None, - depth=2, frontType=1, octreeMode=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., - tbox=None, extrusion='cart'): - - import Generator.IBM as G_IBM - import time as python_time - - if isinstance(t_case, str): tb = C.convertFile2PyTree(t_case) - else: tb = Internal.copyTree(t_case) - - tbF1 = None - if tbox: - tbF1 = Internal.getNodesFromNameAndType(tbox, '*KeepF1*', 'CGNSBase_t') - tbox = None - - ##if isinstance(tc_out, str): fileoutpre = tc_out.split('/') - ##else: fileoutpre = ['.','template.cgns'] - - refstate = Internal.getNodeFromName(tb, 'ReferenceState') - flowEqn = Internal.getNodeFromName(tb, 'FlowEquationSet') - - Reynolds = Internal.getNodeFromName(tb, 'Reynolds') - if Reynolds is not None: - Reynolds = Internal.getValue(Reynolds) - if Reynolds < 1.e5: frontType = 1 - else: Reynolds = 1.e6 - - expand = 3 if frontType != 42 else 4 - - dimPb = Internal.getNodeFromName(tb, 'EquationDimension') - if dimPb is None: raise ValueError('prepareIBMDataPara: EquationDimension is missing in input geometry tree.') - dimPb = Internal.getValue(dimPb) - if dimPb == 2: C._initVars(tb, 'CoordinateZ', 0.) - - model = Internal.getNodeFromName(tb, 'GoverningEquations') - if model is None: raise ValueError('prepareIBMDataPara: GoverningEquations is missing in input geometry tree.') - model = Internal.getValue(model) - - ibctypes = Internal.getNodesFromName(tb, 'ibctype') - if ibctypes is None: raise ValueError('prepareIBMDataPara: ibc type is missing in input geometry tree.') - ibctypes = list(set(Internal.getValue(ibc) for ibc in ibctypes)) - - if model == 'Euler': - if any(ibc in ['Musker', 'MuskerMob', 'Mafzal', 'Log', 'TBLE', 'TBLE_FULL'] for ibc in ibctypes): - raise ValueError("prepareIBMDataPara: governing equations (Euler) not consistent with ibc types %s"%(ibctypes)) - - #=================== - # STEP 0 : GET FILAMENT BODIES Modified - #=================== - [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]) - - # Modification needed for Extrude - # Keep F1 regions - for F1 & F42 synergy - if tbF1: - tbbBTmp = G.BB(tbF1) - interDict_scale = X.getIntersectingDomains(tbbBTmp, t) - for kk in interDict_scale: - for kkk in interDict_scale[kk]: - z=Internal.getNodeFromName(t, kkk) - Internal._createUniqueChild(z, '.Solver#defineTMP', 'UserDefinedData_t') - Internal._createUniqueChild(Internal.getNodeFromName1(z, '.Solver#defineTMP'), 'SaveF1', 'DataArray_t', value=1) - node=Internal.getNodeFromName(t, kkk) - - ## New to extrude Prep - if extrusion=='cyl': cartesian = False - - _redispatch__(t=t) - - C._initVars(t, 'centers:cellN', 1.) - X._applyBCOverlaps(t, depth=depth, loc='centers', val=2, cellNName='cellN') - C._initVars(t,'{centers:cellNChim}={centers:cellN}') - C._initVars(t, 'centers:cellN', 1.) - - #=================== - # STEP 1 : BLANKING IBM Modified - #=================== - if verbose: pt0 = python_time.time(); printTimeAndMemory__('blank by IBC bodies', time=-1, functionName='setInterpDataIBMExtrude') - _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) - - ##Modify needed for extrude & to replicate previous behavior - ##Ghost kmin et kmax donneuse potentiel - listvars_local =['cellNChim','cellNIBC'] - for z in Internal.getZones(t): - sol = Internal.getNodeFromName(z,'FlowSolution#Centers') - for var in listvars_local: - cellN = Internal.getNodeFromName(sol,var)[1] - sh = numpy.shape(cellN) - for k in [0,1, sh[2]-2, sh[2]-1]: - for j in range(sh[1]): - for i in range(sh[0]): - if cellN[i,j,k] != 0: cellN[i,j,k] =1 - C._initVars(t,'{centers:cellN}=maximum(0.,{centers:cellNChim})') # vaut -3, 0, 1, 2 initialement - - Cmpi.barrier() - _redispatch__(t=t) - if verbose: printTimeAndMemory__('blank by IBC bodies', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') - - ##Same as prepareIBMDataPara - #=================== - # STEP 2 : INTERP DATA CHIM - #=================== - ## cellN mush be correct here - ## _setInterpData uses cellN - if verbose: pt0 = python_time.time(); printTimeAndMemory__('compute interpolation data (Abutting & Chimera)', time=-1, functionName='setInterpDataIBMExtrude') - tc = C.node2Center(t) - - if Internal.getNodeFromType(t, "GridConnectivity1to1_t") is not None: - Xmpi._setInterpData(t, tc, nature=1, loc='centers', storage='inverse', sameName=1, dim=dimPb, itype='abutting', order=2, cartesian=cartesian) - Xmpi._setInterpData(t, tc, nature=1, loc='centers', storage='inverse', sameName=1, sameBase=1, dim=dimPb, itype='chimera', order=2, cartesian=cartesian) - if verbose: printTimeAndMemory__('compute interpolation data (Abutting & Chimera)', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') - - #=================== - # STEP 3 : BUILD FRONT - #=================== - if verbose: pt0 = python_time.time(); printTimeAndMemory__('build IBM front', time=-1, functionName='setInterpDataIBMExtrude') - t, tc, front, front2, frontWMM = buildFrontIBM(t, tc, dimPb=dimPb, frontType=frontType, - cartesian=cartesian, twoFronts=twoFronts, check=check, - isWireModel=isWireModel) - if verbose: printTimeAndMemory__('build IBM front', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') - - #=================== - # STEP 4 : INTERP DATA IBM - #=================== - if verbose: pt0 = python_time.time(); printTimeAndMemory__('compute interpolation data (IBM)', time=-1, functionName='setInterpDataIBMExtrude') - _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) - if verbose: printTimeAndMemory__('compute interpolation data (IBM)', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') - - - #=================== - # STEP 5 : INIT IBM Modified - #=================== - tc2 = Internal.copyTree(tc) if twoFronts or isWireModel else None - - if isWireModel: - Internal._rmNodesByName(tc2, 'IBCD*') - Internal._rmNodesByName(tc, '2_IBCD*') - D_IBM._transformTc2(tc2) - tc2 = Internal.rmNodesByName(tc2, 'ID*') - NewIBCD= 141 - _changeNameIBCD__(tc2,NewIBCD) - tc = Internal.merge([tc,tc2]) - tc2 = None - - ## New to extrude Prep - if extrusion == 'cyl': - T._cyl2Cart(t, (0,0,0),(1,0,0)) - T._cyl2Cart(tc,(0,0,0),(1,0,0)) - - # modif info mesh in zonesubregion_t - for z in Internal.getZones(tc): - for zsr in Internal.getNodesFromType(z, "ZoneSubRegion_t"): - zsrname = Internal.getName(zsr) - zsrname = zsrname.split('_') - if zsrname[0]=='IBCD': - for var in ['C','W','I']: - r = Internal.getNodeFromName(zsr,'CoordinateY_P'+var)[1] - theta = Internal.getNodeFromName(zsr,'CoordinateZ_P'+var)[1] - for l in range(numpy.size(r)): - yy = r[l]*numpy.cos( theta[l] ) - zz = r[l]*numpy.sin( theta[l] ) - r[l]= yy; theta[l] = zz - - if distribute and Cmpi.size > 1: _redispatch__(t=t, tc=tc, tc2=tc2, twoFronts=twoFronts) - - ## New to extrude Prep - vars = ['centers:TurbulentDistanceAllBC','centers:TurbulentDistanceWallBC', 'centers:cellNIBC_hole'] - C._rmVars(t, vars) - - if isinstance(tc_out, str): - tcp = Compressor.compressCartesian(tc) - Cmpi.convertPyTree2File(tcp, tc_out, ignoreProcNodes=True) - - if twoFronts: - tcp2 = Compressor.compressCartesian(tc2) - tc2_out = tc_out.replace('tc', 'tc2') if 'tc' in tc_out else 'tc2.cgns' - Cmpi.convertPyTree2File(tcp2, tc2_out, ignoreProcNodes=True) - - if isinstance(t_out, str): - tp = Compressor.compressCartesian(t) - Cmpi.convertPyTree2File(tp, t_out, ignoreProcNodes=True) - - _computeMeshInfo(t) - - if Cmpi.size > 1: Cmpi.barrier() - if verbose: printTimeAndMemory__('initialize and clean', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') - - if tc2 is not None: return t, tc, tc2 - else: return t, tc - -def setInterpDataHybrid(t_3d, tc_3d, t_curvi, extrusion=None, interpDataType=1): - import FastS.ToolboxChimera as TBX - #================================================================================ - # IBM prepare hybride cart + curvi - # IN: t_3d arbre cartesien (from prepare1) - # IN: tc_3d arbre cartesien connectivite (from prepare1) - # IN: t_curvi arbre curviligne (with bc and connectivity but without ghost) - #================================================================================ - overlap = '14' # Overlap IBC - InterpOrder = 2 - root_racHybride = 'joinIBC' - - #add 2 ghost cells to the curvilinear mesh - Internal._addGhostCells(t_curvi, t_curvi, 2, adaptBCs=1, fillCorner=0) - C._initVars(t_curvi,"centers:cellN",1.) - t_curvi = TBX.modifyBCOverlapsForGhostMesh(t_curvi,2) - - dimPb = Internal.getValue(Internal.getNodeFromName(t_3d, 'EquationDimension')) - - for z in Internal.getZones(t_curvi): - if root_racHybride in z[0]: - C._fillEmptyBCWith(z,'inactive','BCExtrapolated',dim=dimPb) - C._rmBCOfType(z,'BCOverlap') - C._fillEmptyBCWith(z,'trou','BCOverlap',dim=dimPb) - X._applyBCOverlaps(z,loc='centers',depth=2, val=0) - X._applyBCOverlaps(z,loc='centers',depth=4, val=2) - - if extrusion == 'cyl': - T._cart2Cyl(t_3d , (0,0,0),(1,0,0)) - T._cart2Cyl(tc_3d, (0,0,0),(1,0,0)) - T._cart2Cyl(t_curvi , (0,0,0),(1,0,0)) - - tBB2 = G.BB(t_curvi) - tBB = G.BB(t_3d) - - tc2 = C.node2Center(t_curvi) - - X._setInterpData(t_curvi, tc2, nature=1, loc='centers', storage='inverse', sameName=1, dim=dimPb, itype='abutting') - - C._initVars(t_3d,"{centers:cellN#Init}={centers:cellN}") - - for var in ['wall','racChimer']: - C._initVars(t_3d,"centers:cellN",1.) - if var == 'wall': itype ="3" - else: itype = overlap - for zc in Internal.getZones(tc_3d): - for zsr in Internal.getNodesFromType(zc, "ZoneSubRegion_t"): - zsrname = Internal.getName(zsr) - zsrname = zsrname.split('_') - if zsrname[0]=='IBCD' and zsrname[1] == itype: - zrname = Internal.getValue(zsr) - ptlistD= Internal.getNodeFromName(zsr,'PointListDonor')[1] - zloc = Internal.getNodeFromName(t_3d,zrname) - sol = Internal.getNodeFromName(zloc,'FlowSolution#Centers') - cellN = Internal.getNodeFromName(sol,'cellN')[1] - sh = numpy.shape(cellN) - ni= sh[0]; ninj = sh[0]*sh[1] - for l0 in range(numpy.size(ptlistD)): - l = ptlistD[ l0] - k = l//ninj - j = (l-k*ninj)//ni - i = l-k*ninj - j*ni - cellN[ i,j,k ]=2 - - if var == 'wall': C._initVars(t_3d,"{centers:cellN#wall}={centers:cellN}") - else: C._initVars(t_3d,"{centers:cellN#racChim}={centers:cellN}") - - Internal._rmNodesFromName(tc_3d,'IBCD_'+overlap+'_*') - - for z in Internal.getZones(t_3d): - sol = Internal.getNodeFromName(z,'FlowSolution#Centers') - cellNRac = Internal.getNodeFromName(sol,'cellN#racChim')[1] - cellN = Internal.getNodeFromName(sol,'cellN')[1] - C._initVars(z,"{centers:cellN}={centers:cellN#racChim}") - sh = numpy.shape(cellN) - - if sh[2]> 1: - for k in [0,1, sh[2]-2, sh[2]-1]: - for j in range(sh[1]): - for i in range(sh[0]): - if cellN[i,j,k] != 0: cellN[i,j,k] =1 - - interDict = X.getIntersectingDomains(tBB, t2=tBB2, method='AABB', taabb=tBB, taabb2=tBB2) - print(" Interp Cartesian from curvilinear") - zonelist=[] - for z in Internal.getZones(t_3d): - if C.getMaxValue(z,'centers:cellN')==2: - dnrZones = [] - for zdname in interDict[z[0]]: - zd = Internal.getNodeFromName(tc2,zdname) - dnrZones.append(zd) - X._setInterpData(z,dnrZones, nature=0,penalty=1,loc='centers',storage='inverse',sameName=1,\ - interpDataType=interpDataType, itype='chimera', order=InterpOrder) - z = X.getOversetInfo(z, dnrZones, loc='center',type='extrapolated') - zonelist.append(z) - - interDict_curvi = X.getIntersectingDomains(tBB2, t2=tBB, method='AABB', taabb=tBB2, taabb2=tBB) - print(" Interp curvi from Cartesian") - for z in Internal.getZones(t_curvi): - if C.getMaxValue(z,'centers:cellN')==2: - dnrZones = [] - for zdname in interDict_curvi[z[0]]: - zd = Internal.getNodeFromName(tc_3d,zdname) - dnrZones.append(zd) - - X._setInterpData(z,dnrZones,nature=0,penalty=1,loc='centers',storage='inverse',sameName=1,\ - interpDataType=interpDataType, itype='chimera', order=InterpOrder) - - # to check orphans - #z = X.getOversetInfo(z, dnrZones, loc='center',type='extrapolated') - #z = X.getOversetInfo(z, dnrZones, loc='center',type='orphan') - #zonelist.append(z) - #C.convertPyTree2File(z,"rcv2.cgns") - - C._initVars(t_3d, '{centers:cellN}={centers:cellN#Init}') - - if extrusion == 'cyl': - T._cyl2Cart(t_3d, (0,0,0),(1,0,0)) - T._cyl2Cart(t_curvi, (0,0,0),(1,0,0)) - T._cyl2Cart(tc_3d, (0,0,0),(1,0,0)) - T._cyl2Cart(tc2, (0,0,0),(1,0,0)) - - for z in Internal.getZones(t_curvi): - for bc in Internal.getNodesFromType(z,'BC_t'): - if 'inactive' in bc[0] and Internal.getValue(bc) == 'BCExtrapolated': - Internal._rmNodesByName(z, bc[0]) - - t = C.mergeTrees(t_3d ,t_curvi ) - tc = C.mergeTrees(tc_3d,tc2) - - return t, tc +#def prepareIBMDataExtrude(t_case, t_out, tc_out, t, to=None, +# depth=2, frontType=1, octreeMode=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., +# tbox=None, extrusion='cart'): +# STILL IN DEV - DO NOT USE +# import Generator.IBM as G_IBM +# import time as python_time +# +# if isinstance(t_case, str): tb = C.convertFile2PyTree(t_case) +# else: tb = Internal.copyTree(t_case) +# +# tbF1 = None +# if tbox: +# tbF1 = Internal.getNodesFromNameAndType(tbox, '*KeepF1*', 'CGNSBase_t') +# tbox = None +# +# ##if isinstance(tc_out, str): fileoutpre = tc_out.split('/') +# ##else: fileoutpre = ['.','template.cgns'] +# +# refstate = Internal.getNodeFromName(tb, 'ReferenceState') +# flowEqn = Internal.getNodeFromName(tb, 'FlowEquationSet') +# +# Reynolds = Internal.getNodeFromName(tb, 'Reynolds') +# if Reynolds is not None: +# Reynolds = Internal.getValue(Reynolds) +# if Reynolds < 1.e5: frontType = 1 +# else: Reynolds = 1.e6 +# +# expand = 3 if frontType != 42 else 4 +# +# dimPb = Internal.getNodeFromName(tb, 'EquationDimension') +# if dimPb is None: raise ValueError('prepareIBMDataPara: EquationDimension is missing in input geometry tree.') +# dimPb = Internal.getValue(dimPb) +# if dimPb == 2: C._initVars(tb, 'CoordinateZ', 0.) +# +# model = Internal.getNodeFromName(tb, 'GoverningEquations') +# if model is None: raise ValueError('prepareIBMDataPara: GoverningEquations is missing in input geometry tree.') +# model = Internal.getValue(model) +# +# ibctypes = Internal.getNodesFromName(tb, 'ibctype') +# if ibctypes is None: raise ValueError('prepareIBMDataPara: ibc type is missing in input geometry tree.') +# ibctypes = list(set(Internal.getValue(ibc) for ibc in ibctypes)) +# +# if model == 'Euler': +# if any(ibc in ['Musker', 'MuskerMob', 'Mafzal', 'Log', 'TBLE', 'TBLE_FULL'] for ibc in ibctypes): +# raise ValueError("prepareIBMDataPara: governing equations (Euler) not consistent with ibc types %s"%(ibctypes)) +# +# #=================== +# # STEP 0 : GET FILAMENT BODIES Modified +# #=================== +# [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]) +# +# # Modification needed for Extrude +# # Keep F1 regions - for F1 & F42 synergy +# if tbF1: +# tbbBTmp = G.BB(tbF1) +# interDict_scale = X.getIntersectingDomains(tbbBTmp, t) +# for kk in interDict_scale: +# for kkk in interDict_scale[kk]: +# z=Internal.getNodeFromName(t, kkk) +# Internal._createUniqueChild(z, '.Solver#defineTMP', 'UserDefinedData_t') +# Internal._createUniqueChild(Internal.getNodeFromName1(z, '.Solver#defineTMP'), 'SaveF1', 'DataArray_t', value=1) +# node=Internal.getNodeFromName(t, kkk) +# +# ## New to extrude Prep +# if extrusion=='cyl': cartesian = False +# +# C._initVars(t, 'centers:cellN', 1.) +# X._applyBCOverlaps(t, depth=depth, loc='centers', val=2, cellNName='cellN') +# C._initVars(t,'{centers:cellNChim}={centers:cellN}') +# C._initVars(t, 'centers:cellN', 1.) +# +# #=================== +# # STEP 1 : BLANKING IBM Modified +# #=================== +# if verbose: pt0 = python_time.time(); printTimeAndMemory__('blank by IBC bodies', time=-1, functionName='setInterpDataIBMExtrude') +# _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) +# +# ##Modify needed for extrude & to replicate previous behavior +# ##Ghost kmin et kmax donneuse potentiel +# listvars_local =['cellNChim','cellNIBC'] +# for z in Internal.getZones(t): +# sol = Internal.getNodeFromName(z,'FlowSolution#Centers') +# for var in listvars_local: +# cellN = Internal.getNodeFromName(sol,var)[1] +# sh = numpy.shape(cellN) +# for k in [0,1, sh[2]-2, sh[2]-1]: +# for j in range(sh[1]): +# for i in range(sh[0]): +# if cellN[i,j,k] != 0: cellN[i,j,k] =1 +# C._initVars(t,'{centers:cellN}=maximum(0.,{centers:cellNChim})') # vaut -3, 0, 1, 2 initialement +# +# Cmpi.barrier() +# _redispatch__(t=t) +# if verbose: printTimeAndMemory__('blank by IBC bodies', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') +# +# ##Same as prepareIBMDataPara +# #=================== +# # STEP 2 : INTERP DATA CHIM +# #=================== +# ## cellN mush be correct here +# ## _setInterpData uses cellN +# if verbose: pt0 = python_time.time(); printTimeAndMemory__('compute interpolation data (Abutting & Chimera)', time=-1, functionName='setInterpDataIBMExtrude') +# tc = C.node2Center(t) +# +# if Internal.getNodeFromType(t, "GridConnectivity1to1_t") is not None: +# Xmpi._setInterpData(t, tc, nature=1, loc='centers', storage='inverse', sameName=1, dim=dimPb, itype='abutting', order=2, cartesian=cartesian) +# Xmpi._setInterpData(t, tc, nature=1, loc='centers', storage='inverse', sameName=1, sameBase=1, dim=dimPb, itype='chimera', order=2, cartesian=cartesian) +# if verbose: printTimeAndMemory__('compute interpolation data (Abutting & Chimera)', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') +# +# #=================== +# # STEP 3 : BUILD FRONT +# #=================== +# if verbose: pt0 = python_time.time(); printTimeAndMemory__('build IBM front', time=-1, functionName='setInterpDataIBMExtrude') +# t, tc, front, front2, frontWMM = buildFrontIBM(t, tc, dimPb=dimPb, frontType=frontType, +# cartesian=cartesian, twoFronts=twoFronts, check=check, +# isWireModel=isWireModel) +# if verbose: printTimeAndMemory__('build IBM front', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') +# +# #=================== +# # STEP 4 : INTERP DATA IBM +# #=================== +# if verbose: pt0 = python_time.time(); printTimeAndMemory__('compute interpolation data (IBM)', time=-1, functionName='setInterpDataIBMExtrude') +# _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) +# if verbose: printTimeAndMemory__('compute interpolation data (IBM)', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') +# +# +# #=================== +# # STEP 5 : INIT IBM Modified +# #=================== +# tc2 = Internal.copyTree(tc) if twoFronts or isWireModel else None +# +# if isWireModel: +# Internal._rmNodesByName(tc2, 'IBCD*') +# Internal._rmNodesByName(tc, '2_IBCD*') +# D_IBM._transformTc2(tc2) +# tc2 = Internal.rmNodesByName(tc2, 'ID*') +# NewIBCD= 141 +# _changeNameIBCD__(tc2,NewIBCD) +# tc = Internal.merge([tc,tc2]) +# tc2 = None +# +# ## New to extrude Prep +# if extrusion == 'cyl': +# T._cyl2Cart(t, (0,0,0),(1,0,0)) +# T._cyl2Cart(tc,(0,0,0),(1,0,0)) +# +# # modif info mesh in zonesubregion_t +# for z in Internal.getZones(tc): +# for zsr in Internal.getNodesFromType(z, "ZoneSubRegion_t"): +# zsrname = Internal.getName(zsr) +# zsrname = zsrname.split('_') +# if zsrname[0]=='IBCD': +# for var in ['C','W','I']: +# r = Internal.getNodeFromName(zsr,'CoordinateY_P'+var)[1] +# theta = Internal.getNodeFromName(zsr,'CoordinateZ_P'+var)[1] +# for l in range(numpy.size(r)): +# yy = r[l]*numpy.cos( theta[l] ) +# zz = r[l]*numpy.sin( theta[l] ) +# r[l]= yy; theta[l] = zz +# +# if distribute and Cmpi.size > 1: _redispatch__(t=t, tc=tc, tc2=tc2, twoFronts=twoFronts) +# +# ## New to extrude Prep +# vars = ['centers:TurbulentDistanceAllBC','centers:TurbulentDistanceWallBC', 'centers:cellNIBC_hole'] +# C._rmVars(t, vars) +# +# if isinstance(tc_out, str): +# tcp = Compressor.compressCartesian(tc) +# Cmpi.convertPyTree2File(tcp, tc_out, ignoreProcNodes=True) +# +# if twoFronts: +# tcp2 = Compressor.compressCartesian(tc2) +# tc2_out = tc_out.replace('tc', 'tc2') if 'tc' in tc_out else 'tc2.cgns' +# Cmpi.convertPyTree2File(tcp2, tc2_out, ignoreProcNodes=True) +# +# if isinstance(t_out, str): +# tp = Compressor.compressCartesian(t) +# Cmpi.convertPyTree2File(tp, t_out, ignoreProcNodes=True) +# +# _computeMeshInfo(t) +# +# if Cmpi.size > 1: Cmpi.barrier() +# if verbose: printTimeAndMemory__('initialize and clean', time=python_time.time()-pt0, functionName='setInterpDataIBMExtrude') +# +# if tc2 is not None: return t, tc, tc2 +# else: return t, tc def _redispatch__(t=None, tc=None, tc2=None): import Distributor2.Mpi as D2mpi diff --git a/Cassiopee/Connector/test/IBMExtrude2Dto3D_t1.py b/Cassiopee/Connector/test/ToBeCompletedIBMinterpolateExtrude2Dto3D.py similarity index 95% rename from Cassiopee/Connector/test/IBMExtrude2Dto3D_t1.py rename to Cassiopee/Connector/test/ToBeCompletedIBMinterpolateExtrude2Dto3D.py index c1ba84205..674fb5d11 100644 --- a/Cassiopee/Connector/test/IBMExtrude2Dto3D_t1.py +++ b/Cassiopee/Connector/test/ToBeCompletedIBMinterpolateExtrude2Dto3D.py @@ -1,3 +1,6 @@ + +## DO NOT RUN THIS TEST CASE - IT IS IN PROGRESS AND WILL BE COMPLETED SHORTLY + # extrude 2D mesh to 3D with Cartesian approach import Converter.PyTree as C import Converter.Mpi as Cmpi diff --git a/Cassiopee/Generator/Generator/IBM.py b/Cassiopee/Generator/Generator/IBM.py index 5bf314dee..3e9bfbcd9 100644 --- a/Cassiopee/Generator/Generator/IBM.py +++ b/Cassiopee/Generator/Generator/IBM.py @@ -1125,13 +1125,16 @@ def extrudeCartesianZDir(t, tb, check=False, extrusion="cart", dz=0.01, NPas=10, Usage: extrudeCartesianZDir(t, tb, check, extrusion, dz, NPas, span, Ntranche, dictNz, ific, isCartesianExtrude, isAutoPeriodic, nghost)""" ## isAutoPeriodic=True is prefered over isAutoPeriodic=False when extruding a Cartesian mesh. - ## nghost should be 2. It is currently 0 (default value) to pass the non-regression tests. This function was initially developed + ## nghost should be 2. It is currently 0 (default value) to pass the non-regression tests. This function was initially developed for ## an application and its default value of 0 is kept for reproduciability reasons for that original test case. nghost is only used when isAutoPeriodic=False. ## isAutoPeridioc=True assumes 2 ghost cells. if extrusion == "cart": perio = span/float(Ntranche) else: perio = span/180.*math.pi/float(Ntranche) - + + ## Found in original function. Commented here as it is related to interpolations + ## and not geometric extrusion. This functionality is now in Connector/IBM.py - it is has been + ## deemed to be a more appropriate location. #for z in Internal.getZones(t): # cellN = Internal.getNodeFromName(z, "cellN")[1] # sh = numpy.shape(cellN) @@ -1289,7 +1292,7 @@ def extrudeCartesianZDir(t, tb, check=False, extrusion="cart", dz=0.01, NPas=10, for node in Internal.getNodesFromName(t,'EquationDimension'): Internal.setValue(node,3) for z in Internal.getZones(t): C._addBC2Zone(z, z[0]+'periodKmin', 'BCautoperiod', 'kmin') - C._addBC2Zone(z, z[0]+'periodKmin', 'BCautoperiod', 'kmax') + C._addBC2Zone(z, z[0]+'periodKmax', 'BCautoperiod', 'kmax') BCs = Internal.getNodesFromType(t, "BC_t") for bc in BCs: if Internal.getValue(bc)=='BCautoperiod': @@ -1303,6 +1306,6 @@ def extrudeCartesianZDir(t, tb, check=False, extrusion="cart", dz=0.01, NPas=10, if extrusion == 'cyl': T._cart2Cyl(t, (0,0,0),(1,0,0)) T._cart2Cyl(tb, (0,0,0),(1,0,0)) - + X_IBM._redispatch__(t=t) return t, tb diff --git a/Cassiopee/Generator/test/IBMExtrude2Dto3D_t1.py b/Cassiopee/Generator/test/IBMExtrude2Dto3D_t1.py new file mode 100644 index 000000000..891566cb0 --- /dev/null +++ b/Cassiopee/Generator/test/IBMExtrude2Dto3D_t1.py @@ -0,0 +1,67 @@ +# extrude 2D mesh to 3D with Cartesian approach +import Apps.Fast.IBM as App +import Converter.PyTree as C +import Transform.PyTree as T +import Converter.Internal as Internal +import Generator.IBM as G_IBM +import Connector.IBM as X_IBM +import KCore.test as test +import sys + +LOCAL = test.getLocal() +bodySurfaceFile = '../../Apps/test/naca1DNS.cgns' + +# Prepare +vmin = 42 +dfars = 5 +snears = 1 +t, tc = X_IBM.prepareIBMData(bodySurfaceFile, None , None , + snears=snears , dfars=dfars , vmin=vmin, + check=False , frontType=1) + +test.testT(t ,1) +test.testT(tc,2) +#C.convertPyTree2File(t,LOCAL+'/t2D_checking.cgns') + + +####Extrusion for 3D Mesh +bodySurface = C.convertFile2PyTree(bodySurfaceFile) + +t2 = Internal.copyTree(t) +bodySurface2 = Internal.copyTree(bodySurface) + +extrusion = 'cart' +span = 1 +NPas = 10+1 #number of nodes +t3D, tb3D = G_IBM.extrudeCartesianZDir(t, bodySurface, extrusion=extrusion, NPas=NPas, span=span, dz=span/(NPas-1), isAutoPeriodic=True) + +for t in [t3D,tb3D]: + zmax = C.getMaxValue(t, 'CoordinateZ'); + zmin = C.getMinValue(t, 'CoordinateZ'); + zavg = (zmax+zmin)/2 + T._translate(t, (0,0,0-zavg)) +C._rmVars(t3D, ['centers:cellN']) +test.testT(t3D ,3) +test.testT(tb3D ,4) + +#C.convertPyTree2File(t3D ,LOCAL+'/t3D_checking.cgns') +#C.convertPyTree2File(tb3D,LOCAL+'/tb3D_checking.cgns') +# + +t3Dorig, tb3Dorig = App.extrudeCartesian(t2, bodySurface2, extrusion=extrusion, NPas=NPas, span=span, dz=span/(NPas-1), isAutoPeriodic=True) +for t in [t3Dorig,tb3Dorig]: + zmax = C.getMaxValue(t, 'CoordinateZ'); + zmin = C.getMinValue(t, 'CoordinateZ'); + zavg = (zmax+zmin)/2 + T._translate(t, (0,0,0-zavg)) + +BCs = Internal.getNodesFromType(t3Dorig, "BC_t") +for bc in BCs: + if Internal.getValue(bc)=='BCautoperiod': + nameSplit=bc[0].split(".") + lenLocal=len(nameSplit) + bc[0]='.'.join(nameSplit[0:lenLocal-1]) +C._rmVars(t3Dorig, ['centers:cellN']) +test.testT(t3Dorig,3) + +#C.convertPyTree2File(t3Dorig ,LOCAL+'/t3Dorig_checking.cgns') diff --git a/Cassiopee/Generator/test/IBMExtrude2Dto3D_t2.py b/Cassiopee/Generator/test/IBMExtrude2Dto3D_t2.py new file mode 100644 index 000000000..583821144 --- /dev/null +++ b/Cassiopee/Generator/test/IBMExtrude2Dto3D_t2.py @@ -0,0 +1,67 @@ +# extrude 2D mesh to 3D with Cartesian approach +import Apps.Fast.IBM as App +import Converter.PyTree as C +import Transform.PyTree as T +import Converter.Internal as Internal +import Generator.IBM as G_IBM +import Connector.IBM as X_IBM +import KCore.test as test +import sys + +LOCAL = test.getLocal() +bodySurfaceFile = '../../Apps/test/naca1DNS.cgns' + +# Prepare +vmin = 42 +dfars = 5 +snears = 1 +t, tc = X_IBM.prepareIBMData(bodySurfaceFile, None , None , + snears=snears , dfars=dfars , vmin=vmin, + check=False , frontType=1) + +test.testT(t ,1) +test.testT(tc,2) +#C.convertPyTree2File(t,LOCAL+'/t2D_checking.cgns') + + +####Extrusion for 3D Mesh +bodySurface = C.convertFile2PyTree(bodySurfaceFile) + +t2 = Internal.copyTree(t) +bodySurface2 = Internal.copyTree(bodySurface) + +extrusion = 'cart' +span = 1 +NPas = 10+1 #number of nodes +t3D, tb3D = G_IBM.extrudeCartesianZDir(t, bodySurface, extrusion=extrusion, NPas=NPas, span=span, dz=span/(NPas-1)) + +for t in [t3D,tb3D]: + zmax = C.getMaxValue(t, 'CoordinateZ'); + zmin = C.getMinValue(t, 'CoordinateZ'); + zavg = (zmax+zmin)/2 + T._translate(t, (0,0,0-zavg)) +C._rmVars(t3D, ['centers:cellN']) +test.testT(t3D ,3) +test.testT(tb3D ,4) + +#C.convertPyTree2File(t3D ,LOCAL+'/t3D_checking.cgns') +#C.convertPyTree2File(tb3D,LOCAL+'/tb3D_checking.cgns') +# + +t3Dorig, tb3Dorig = App.extrudeCartesian(t2, bodySurface2, extrusion=extrusion, NPas=NPas, span=span, dz=span/(NPas-1)) +for t in [t3Dorig,tb3Dorig]: + zmax = C.getMaxValue(t, 'CoordinateZ'); + zmin = C.getMinValue(t, 'CoordinateZ'); + zavg = (zmax+zmin)/2 + T._translate(t, (0,0,0-zavg)) + +BCs = Internal.getNodesFromType(t3Dorig, "BC_t") +for bc in BCs: + if Internal.getValue(bc)=='BCautoperiod': + nameSplit=bc[0].split(".") + lenLocal=len(nameSplit) + bc[0]='.'.join(nameSplit[0:lenLocal-1]) +C._rmVars(t3Dorig, ['centers:cellN']) +test.testT(t3Dorig,3) + +#C.convertPyTree2File(t3Dorig ,LOCAL+'/t3Dorig_checking.cgns')