diff --git a/Cassiopee/Apps/Apps/Fast/WindTunnelOutPres.py b/Cassiopee/Apps/Apps/Fast/WindTunnelOutPres.py new file mode 100644 index 000000000..248ad7a5d --- /dev/null +++ b/Cassiopee/Apps/Apps/Fast/WindTunnelOutPres.py @@ -0,0 +1,224 @@ +import Converter.Internal as Internal +import Converter.Mpi as Cmpi +import Converter.PyTree as C +import Post.Probe as Probe +import Post.PyTree as Post +import Geom.IBM as D_IBM +import sys, os +import numpy, scipy.optimize + +R_GAZ = 287.05 # J/kg/K +GAMMA = 1.4 + +def tauFunction__(mach, gamma=GAMMA): + """Return :math:`1 + \\frac{\\gamma - 1}{2} M^2` + + Parameters + ---------- + mach : array_like + Value of Mach number :math:`M` + gamma : float + Specific heat ratio :math:`\\gamma` + """ + return 1. + 0.5*(gamma - 1.)*mach**2 + + +def sigmaFunction__(mach, gamma=GAMMA): + """Return :math:`\\Sigma(M)` + + Parameters + ---------- + mach : array_like + Value of Mach number :math:`M` + gamma : float + Specific heat ratio :math:`\\gamma` + """ + return ( 2./(gamma + 1.) + (gamma - 1.)/(gamma + 1.) * mach**2 )**(0.5*(gamma + 1.)/(gamma - 1.)) / mach + + +def sigmaFunctionDerivate__(mach, gamma=GAMMA): + """Return :math:`\\Sigma'(M)`, first derivative of :math:`\\Sigma(M)` + + Parameters + ---------- + mach : array_like + Value of Mach number :math:`M` + gamma : float + Specific heat ratio :math:`\\gamma` + """ + return ( 2./(gamma + 1.) + (gamma - 1.)/(gamma + 1.) * mach**2 )**(0.5*(gamma + 1.)/(gamma - 1.)) \ + * (-1/mach**2 + 1./(2./(gamma + 1.) + (gamma - 1.)/(gamma + 1.) * mach**2)) + + +def _scalarSigmaFunctionInv__(s, gamma=GAMMA, range='subsonic'): + eps = numpy.finfo(float).eps + if range == 'subsonic': + sol = scipy.optimize.root_scalar(lambda M: sigmaFunction__(M, gamma) - s, + x0=0.5, bracket=(2*eps, 1. - 2*eps), method='brentq') + elif range == 'supersonic': + sol = scipy.optimize.root_scalar(lambda M: sigmaFunction__(M, gamma) - s, + x0=1.5, bracket=(1. + 2*eps, 1e3), method='brentq') # it is unlikely that user require Mach number above 1000. + else: + raise RuntimeError("Unexpected value for `range`: {:s}".format(str(range))) + return sol.root + + +def sigmaFunctionInv__(s, gamma=1.4, range='subsonic'): + # This method vectorizes _scalarSigmaFunctionInv__ + """Return the inverse of the function :math:`\\Sigma(M)` + + Parameters + ---------- + s : array_like + Value of :math:`\\Sigma` + gamma : float + Specific heat ratio :math:`\\gamma` + range : ['subsonic', 'supersonic'] + Range over which the inverse is to be looked for + """ + with numpy.nditer([s, None], + op_flags = [['readonly'], ['writeonly', 'allocate', 'no_broadcast']], + op_dtypes=['float64', 'float64']) as it: + for x, y in it: + y[...] = _scalarSigmaFunctionInv__(s, gamma=gamma) + return it.operands[1] + + + +## ∆p_(2n) = G∆M-(1n) + D[∆M_(1n) − ∆M_(n−ncontrol) ] +## G = -0.7 * ∂p2/∂M1 || G: proportional coefficeint +## 0.7 :: damping coefficient for stability +## ∂p2/∂M1 ≈ −p_(2,is) γ M_(2,is) A_2 Σ′(M_1)/(τ(M_(2,is)) A_1 Σ′(M_(2,is))) +## D = G(n_charac/n_control) +## n_charac = approx. characteristic response time of the CFD of the wind tunnel - a prior run to determine this +## τ :: tau_function +## Σ′:: sigma_function_inv + +def getInfo(tcase,familyName): + familyNameExists=False + for z in Internal.getZones(tcase): + FamNode = Internal.getNodeFromType(z, 'FamilyName_t') + if FamNode is not None: + FamName = Internal.getValue(FamNode) + if FamName == familyName: + familyNameExists = True + solDef = Internal.getNodeFromName(z, '.Solver#define') + controlProbeName = Internal.getNodeFromName(solDef, 'probeName') + controlProbeName = Internal.getValue(controlProbeName) + A1 = Internal.getNodeFromName(solDef, 'AtestSection')[1][0] + A2 = Internal.getNodeFromName(solDef, 'AOutPress')[1][0] + m1 = Internal.getNodeFromName(solDef, 'machTarget')[1][0] + pi1 = Internal.getNodeFromName(solDef, 'pStatTarget')[1][0] + ti1 = Internal.getNodeFromName(solDef, 'tStatTarget')[1][0] + lbda = Internal.getNodeFromName(solDef, 'lmbd')[1][0] + cxSupport = Internal.getNodeFromName(solDef, 'cxSupport')[1][0] + sSupport = Internal.getNodeFromName(solDef, 'sSupport')[1][0] + itExtrctPrb = Internal.getNodeFromName(solDef, 'itExtrctPrb')[1][0] + break + if not familyNameExists: + print('Family name %s is not an outlet pressure family name'%(familyName)) + print('Exiting') + exit() + + + + ## Calculated Values + _tau = tauFunction__(m1) # Eqn. (10) || τ(M)= 1 + (γ-1)/2 M² + m2is = sigmaFunctionInv__(A2/A1 * sigmaFunction__(m1)) # Eqn. (11) || M_2,is=Σ⁻¹(A_2 /A_1 Σ(M_1)) + p2is = pi1 * tauFunction__(m2is)**(-GAMMA/(GAMMA - 1.)) # Eqn. (12) || p_2,is=p_i1 τ(M_2,is)^(−γ∕(γ−1)) + + # coefficient de perte de charge entre l'entrée et la sortie du domaine, + # ici uniquement du au support, et calculé à partir d'une estimation de la traînée du support + if lbda <0: lbda = cxSupport * sSupport/A1 + + p0 = pi1 / _tau **(GAMMA/(GAMMA - 1.)) # Pa - static pressure for m1 + q0 = 0.5*p0*GAMMA*m1**2 # Pa - dynamic pressure for m1 + p2 = p2is - q0*lbda # Eqn. (13) || Pa - outlet static pressure || λ = (p_2,is-p_2)/q_1 + + values4gain =[p2, + m1, + p2is*GAMMA*m2is, + tauFunction__(m2is), + A2/A1, + sigmaFunctionDerivate__(m1), + sigmaFunctionDerivate__(m2is), + 0, + 0] + + return values4gain,controlProbeName,itExtrctPrb + + +def _setUpOutletPressure(values4gain,itValues4gain): + # G : proportional gain Eqn. (15) || G=-0.7* (∂p_2/∂M_1) + # The 0.7 constant is for stability of the algorithm and is independent of the WT + values4gain[7]=-0.7 * (-values4gain[2] / values4gain[3] * values4gain[4] * values4gain[5] / values4gain[6]) + + # D: derivative coefficient || D=G(n_charac/n_control) + values4gain[8]=values4gain[7] * float(itValues4gain[2]) / itValues4gain[1] + return None + + +def getPointsFromTree(tree): + dct_points = {} + for z in Internal.getZones(tree): + name = z[0] + xnode = Internal.getNodeFromName(z,'CoordinateX') + ynode = Internal.getNodeFromName(z,'CoordinateY') + znode = Internal.getNodeFromName(z,'CoordinateZ') + point = [Internal.getValue(xnode),Internal.getValue(ynode),Internal.getValue(znode)] + dct_points[name] = point + return dct_points + + +def setupMachProbe(t,buffer_size,isRestart,DIRECTORY_PROBES): + Post._computeVariables(t, ['centers:Mach']) + + dct_probe_point = {} + dct_points_for_probes = getPointsFromTree(C.convertFile2PyTree(os.path.join(DIRECTORY_PROBES, "probes.cgns"))) + + for name, point in dct_points_for_probes.items(): + print(name,flush=True) + probe = Probe.Probe(os.path.join(DIRECTORY_PROBES, "probe_{:s}.cgns".format(name)), t, X=point, fields=['centers:Mach'], bufferSize=buffer_size, append=isRestart) + dct_probe_point[name] = probe + + C._rmVars(t, ['centers:Mach']) + return dct_points_for_probes,dct_probe_point + + +def recordDataMach(t,dct_probe_point,it): + Post._computeVariables(t, ['centers:Mach']) + for name, probe in dct_probe_point.items(): probe.extract(t, time=it) + C._rmVars(t, ['centers:Mach']) + return dct_probe_point + + +def _controlOutletPressureMachProbe(tc,dctProbes,controlProbeName,DIRECTORY_PROBES,itValues4gain,values4gain,itExtractProbe,it,familyName): + for name, probe in dctProbes.items(): + probe.flush() + Cmpi.barrier() + + if Cmpi.rank == 0: + print(" iteration {:06d}: adjusting back pressure...".format(it), end='') + probe_tmp = C.convertFile2PyTree(os.path.join(DIRECTORY_PROBES, "probe_{:s}.cgns".format(controlProbeName))) + time = numpy.concatenate([node[1] for node in Internal.getNodesFromName(probe_tmp,'time')]) + mach = numpy.concatenate([node[1] for node in Internal.getNodesFromName(probe_tmp,'Mach')]) + select = (time >= 0.) # select values that have been recorded: array is prefilled with -1 + time = time[select] + mach = mach[select] + + index_current = -1 + index_previous = -2 #-1 - (itValues4gain[1] // itExtractProbe) + current_it = time[index_current] + current_mach = mach[index_current] + + previous_it = time[index_previous] + previous_mach = mach[index_previous] + #print(current_it, previous_it, current_mach, previous_mach,flush=True) + oldPressure = values4gain[0] + values4gain[0] += values4gain[7] * (current_mach - values4gain[1]) + values4gain[8] * (current_mach - previous_mach) + print("Control of Output Pressure:: target Ma = {:.5f}, current Ma = {:.5f}|| old Pout = {:.1f} Pa, new Pout = {:.1f} Pa, ".format(values4gain[1], current_mach, oldPressure, values4gain[0]), flush=True) + Cmpi.bcast(values4gain[0], root=0) + D_IBM._initOutflow(tc, familyName, values4gain[0]) + Cmpi.barrier() + return None + diff --git a/Cassiopee/Connector/Connector/getIBMPtsWithFront.cpp b/Cassiopee/Connector/Connector/getIBMPtsWithFront.cpp index 5438ded36..2d7384ca4 100644 --- a/Cassiopee/Connector/Connector/getIBMPtsWithFront.cpp +++ b/Cassiopee/Connector/Connector/getIBMPtsWithFront.cpp @@ -631,6 +631,14 @@ PyObject* K_CONNECTOR::getIBMPtsWithFront(PyObject* self, PyObject* args) // BILAN }// ind in zone // printf(" ZONE %d Nb de type 1 : %d, type2 = %d, type3 = %d, type4 = %d\n", noz, nType1, nType2,nType3,nType4); + if (nType3 > 0 || nType4 > 0) + { + printf("Warning getIBMPtsWithFront: #Zone " SF_D_ " has " SF_D_ " pts of type 3 or 4 ! (Bilan: t1/t2/t3/t4 = %.2f%%/%.2f%%/%.2f%%/%.2f%%)\n", noz, nType2+nType3+nType4, + nType1/E_Float(nType1+nType2+nType3+nType4)*100., + nType2/E_Float(nType1+nType2+nType3+nType4)*100., + nType3/E_Float(nType1+nType2+nType3+nType4)*100., + nType4/E_Float(nType1+nType2+nType3+nType4)*100. ); + } PyObject* PyListIndicesByIBCTypeForZone = PyList_New(0); for (E_Int noibctype = 0; noibctype < nbodies; noibctype++) diff --git a/Cassiopee/Geom/Geom/IBM.py b/Cassiopee/Geom/Geom/IBM.py index deb436967..7f7608544 100644 --- a/Cassiopee/Geom/Geom/IBM.py +++ b/Cassiopee/Geom/Geom/IBM.py @@ -212,21 +212,21 @@ def _setFluidInside(t): #============================================================================== def setOutPressControlParam(t, probeName='pointOutPress', AtestSection=1, AOutPress=1, machTarget=0.1, pStatTarget=1e05, tStatTarget=298.15,lmbd=0.1, - cxSupport = 0.6, sSupport=0.1): + cxSupport = 0.6, sSupport=0.1, itExtrctPrb=10): """Set the user input parameters for the outpress control algorithm. - Usage: setOutPressControlParam(t, probeName='X', AtestSection=Y, AOutPress=Z, machTarget=XX,pStatTarget=YY,tStatTarget=ZZ,lmbd=XXX,cxSupport=YYY,sSupport=ZZZ)""" + Usage: setOutPressControlParam(t, probeName='X', AtestSection=Y, AOutPress=Z, machTarget=XX,pStatTarget=YY,tStatTarget=ZZ,lmbd=XXX,cxSupport=YYY,sSupport=ZZZ,itExtrctPrb=XXXX)""" tp = Internal.copyRef(t) _setOutPressControlParam(tp, probeName=probeName, AtestSection=AtestSection, AOutPress=AOutPress, machTarget=machTarget, pStatTarget=pStatTarget, tStatTarget=tStatTarget,lmbd=lmbd, - cxSupport = cxSupport, sSupport=sSupport) + cxSupport = cxSupport, sSupport=sSupport, itExtrctPrb=itExtrctPrb) return tp def _setOutPressControlParam(t, probeName='pointOutPress', AtestSection=1, AOutPress=1, machTarget=0.1, pStatTarget=1e05, tStatTarget=298.15,lmbd=0.1, - cxSupport = 0.6, sSupport=0.1): + cxSupport = 0.6, sSupport=0.1, itExtrctPrb=10): """Set the user input parameters for the outpress control algorithm. - Usage: _setOutPressControlParam(t, probeName='X', AtestSection=Y, AOutPress=Z, machTarget=XX,pStatTarget=YY,tStatTarget=ZZ,lmbd=XXX,cxSupport=YYY,sSupport=ZZZ)""" + Usage: _setOutPressControlParam(t, probeName='X', AtestSection=Y, AOutPress=Z, machTarget=XX,pStatTarget=YY,tStatTarget=ZZ,lmbd=XXX,cxSupport=YYY,sSupport=ZZZ,itExtrctPrb=XXXX)""" zones = Internal.getZones(t) for z in zones: Internal._createUniqueChild(z, '.Solver#define', 'UserDefinedData_t') @@ -240,6 +240,7 @@ def _setOutPressControlParam(t, probeName='pointOutPress', AtestSection=1, AOutP Internal._createUniqueChild(n, 'lmbd' , 'DataArray_t', lmbd) Internal._createUniqueChild(n, 'cxSupport' , 'DataArray_t', cxSupport) Internal._createUniqueChild(n, 'sSupport' , 'DataArray_t', sSupport) + Internal._createUniqueChild(n, 'itExtrctPrb' , 'DataArray_t', itExtrctPrb) return None #==============================================================================