diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 932a46ad4..bf1e91486 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -5766,8 +5766,15 @@ def _getDefaultOptions(): "verifySpatial": [bool, True], "verifyExtra": [bool, True], # Function parmeters + "computeSepSensorKs": [bool, False], + "sepSensorMaxRho": [float, 1000.0], "sepSensorOffset": [float, 0.0], + "sepSensorOffsetOne": [float, 0.0], + "sepSensorOffsetTwo": [float, 0.0], "sepSensorSharpness": [float, 10.0], + "sepSensorSharpnessOne": [float, 25.0], + "sepSensorSharpnessTwo": [float, 25.0], + "sepAngleDeviation": [float, 90.0], "cavSensorOffset": [float, 0.0], "cavSensorSharpness": [float, 10.0], "cavExponent": [int, 0], @@ -6181,8 +6188,15 @@ def _getOptionMap(self): "verifyextra": ["adjoint", "verifyextra"], "usematrixfreedrdw": ["adjoint", "usematrixfreedrdw"], # Parameters for functions + "computesepsensorks": ["cost", "computesepsensorks"], + "sepsensormaxrho": ["physics", "sepsenmaxrho"], "sepsensoroffset": ["cost", "sepsensoroffset"], "sepsensorsharpness": ["cost", "sepsensorsharpness"], + "sepsensoroffsetone": ["cost", "sepsensoroffsetone"], + "sepsensorsharpnessone": ["cost", "sepsensorsharpnessone"], + "sepsensoroffsettwo": ["cost", "sepsensoroffsettwo"], + "sepsensorsharpnesstwo": ["cost", "sepsensorsharpnesstwo"], + "sepangledeviation": ["cost", "sepangledeviation"], "cavsensoroffset": ["cost", "cavsensoroffset"], "cavsensorsharpness": ["cost", "cavsensorsharpness"], "cavexponent": ["cost", "cavexponent"], @@ -6326,6 +6340,8 @@ def _getObjectivesAndDVs(self): "clqdot": self.adflow.constants.costfuncclqdot, "cbend": self.adflow.constants.costfuncbendingcoef, "sepsensor": self.adflow.constants.costfuncsepsensor, + "sepsensorks": self.adflow.constants.costfuncsepsensorks, + "sepsensorarea": self.adflow.constants.costfuncsepsensorarea, "sepsensoravgx": self.adflow.constants.costfuncsepsensoravgx, "sepsensoravgy": self.adflow.constants.costfuncsepsensoravgy, "sepsensoravgz": self.adflow.constants.costfuncsepsensoravgz, diff --git a/doc/costFunctions.yaml b/doc/costFunctions.yaml index 5d833ea00..1ed4f3ba0 100644 --- a/doc/costFunctions.yaml +++ b/doc/costFunctions.yaml @@ -307,6 +307,22 @@ sepsensor: The separation values for the given surface is provided by this cost function. See :cite:t:`Kenway2017b` for more details. +sepsensorks: + desc: > + The separation sensor value based on the KS aggregation for the given surface is provided by this cost function. + We first compute the deviation of the local velocity from the projected freestream on the desired surface. + Then, we use a trigonometric function to compute the sensor by providing an allowable flow deviation angle from this projected vector. + As a result, the sensor provides values form ``1`` to ``-.inf``. + Any values that are greater than 0 are out of the alowable flow deviation. + Thus, to constraint the separation, we use KS aggregation to find the maximum value in the sensor and can constraint it to be less than or equal to ``0``. + +sepsensorarea: + desc: > + The area separated based on the KS aggregation approach.. + This is evaluated by the considering two separated heaviside smoothing with KS aggregation approach. + It is recommended to start the initial angle of attack in the separated region to compute this value because of the sensor's highly nonlinear behaviour. + Please choose a smaller major steplimit parameter in the optimizer or smaller inexact step in Schur complement solver. + sepsensoravgx: desc: > The separation sensor average in x direction. diff --git a/doc/options.yaml b/doc/options.yaml index d993ae409..34b6fd768 100644 --- a/doc/options.yaml +++ b/doc/options.yaml @@ -1527,16 +1527,58 @@ verifyExtra: This option is for debugging the adjoint only. It is used to verify dIda. +sepSensorMaxRho: + desc: > + The rho parameter used for the KS aggregation used for computing maximum separation sensor value in ``sepSensorKs`` cost function. + KS aggregation is used to compute the differentiable max separation sensor within the given family group. + A higher rho value will approach the actual max separation sensor value more accurately, at the cost of a more nonlinear function. + sepSensorOffset: desc: > - The offset value used for the separation sensor. + The offset value used for the separation sensor in ``sepSensor`` cost function. See "Buffet-Onset Constraint Formulation for Aerodynamic Shape Optimization" (Kenway2017b) for more details. sepSensorSharpness: desc: > - The sharpness parameter for the separation sensor. + The sharpness parameter for the separation sensor in ``sepSensor`` cost function. + See "Buffet-Onset Constraint Formulation for Aerodynamic Shape Optimization" (Kenway2017b) for more details. + +sepSensorOffsetOne: + desc: > + The offset value used for the separation sensor in ``sepSensor`` cost function. + See "Buffet-Onset Constraint Formulation for Aerodynamic Shape Optimization" (Kenway2017b) for more details. + +sepSensorSharpnessOne: + desc: > + The sharpness parameter for the separation sensor in ``sepSensor`` cost function. + See "Buffet-Onset Constraint Formulation for Aerodynamic Shape Optimization" (Kenway2017b) for more details. + +sepSensorOffsetTwo: + desc: > + The offset value used for the separation sensor in ``sepSensor`` cost function. + See "Buffet-Onset Constraint Formulation for Aerodynamic Shape Optimization" (Kenway2017b) for more details. + +sepSensorSharpnessTwo: + desc: > + The sharpness parameter for the separation sensor in ``sepSensor`` cost function. See "Buffet-Onset Constraint Formulation for Aerodynamic Shape Optimization" (Kenway2017b) for more details. +sepAngleDeviation: + desc: > + This option is for KS based separation sensor, ``sepSensorKs`` cost function. + This angle is used to compute the allowable flow deviation region with respect to the freestream projected direction on the desired surface. + For example, on a wing, the streamlines on the upper surface flows inboard direction. + And in some regions, it goes in the outboard directions. + Therefore, we provide an allowable flow deviation from the projected freestream on the surface to compute the sensor. + For instance, when the angle is ``90`` degrees, this approximately closes to the verge of separation. + This angle (in degrees) is for 3D cases only. + For airfoils, by default, this is ``0``. + +computeSepSensorKs: + desc: > + Whether or not to compute KS based separation related cost functions, which is ``sepSensorKs`` cost function. + If this option is not set to ``True``, the code will return zero for this cost function. + computeCavitation: desc: > Whether or not to compute cavitation related cost functions, which are `cavitation` and `cpMin`. diff --git a/src/adjoint/adjointUtils.F90 b/src/adjoint/adjointUtils.F90 index bc46bc30c..8d248d478 100644 --- a/src/adjoint/adjointUtils.F90 +++ b/src/adjoint/adjointUtils.F90 @@ -2087,6 +2087,14 @@ subroutine setDiffSizes ISIZE1OFDrfDrfbcdata_sepSensor = 0 ISIZE2OFDrfDrfbcdata_sepSensor = 0 + ! sepSensorKs + ISIZE1OFDrfDrfbcdata_sepSensorKs = 0 + ISIZE2OFDrfDrfbcdata_sepSensorKs = 0 + + ! sepSensorArea + ISIZE1OFDrfDrfbcdata_sepSensorArea = 0 + ISIZE2OFDrfDrfbcdata_sepSensorArea = 0 + ! Cavitation ISIZE1OFDrfDrfbcdata_Cavitation = 0 ISIZE2OFDrfDrfbcdata_Cavitation = 0 diff --git a/src/adjoint/outputForward/surfaceIntegrations_d.f90 b/src/adjoint/outputForward/surfaceIntegrations_d.f90 index b9169b906..2049be283 100644 --- a/src/adjoint/outputForward/surfaceIntegrations_d.f90 +++ b/src/adjoint/outputForward/surfaceIntegrations_d.f90 @@ -20,8 +20,9 @@ subroutine getcostfunctions_d(globalvals, globalvalsd, funcvalues, & use inputphysics, only : liftdirection, liftdirectiond, & & dragdirection, dragdirectiond, surfaceref, machcoef, machcoefd, & & lengthref, alpha, alphad, beta, betad, liftindex, cpmin_family, & -& cpmin_rho - use inputcostfunctions, only : computecavitation +& cpmin_rho, sepsenmaxfamily, sepsenmaxrho + use inputcostfunctions, only : computecavitation, computesepsensorks& +& , sepsensorsharpnessone, sepsensoroffsetone use inputtsstabderiv, only : tsstability use utils_d, only : computetsderivatives use flowutils_d, only : getdirvector @@ -47,14 +48,17 @@ subroutine getcostfunctions_d(globalvals, globalvalsd, funcvalues, & real(kind=realtype) :: mavgptotd, mavgttotd, mavgrhod, mavgpsd, & & mflowd, mavgmnd, mavgad, mavgvxd, mavgvyd, mavgvzd, garead, mavgvid& & , fxliftd, fyliftd, fzliftd - real(kind=realtype) :: vdotn, mag, u, v, w + real(kind=realtype) :: vdotn, mag, u, v, w, ks_comp + real(kind=realtype) :: ks_compd integer(kind=inttype) :: sps real(kind=realtype), dimension(8) :: dcdq, dcdqdot real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot real(kind=realtype), dimension(8) :: coef0 intrinsic log + intrinsic exp intrinsic sqrt real(kind=realtype) :: arg1 + real(kind=realtype) :: arg1d real(kind=realtype) :: temp real(kind=realtype) :: temp0 real(kind=realtype), dimension(3) :: temp1 @@ -296,6 +300,30 @@ subroutine getcostfunctions_d(globalvals, globalvalsd, funcvalues, & & ovrnts*cmomentd(3, sps) funcvalues(costfuncmomzcoef) = funcvalues(costfuncmomzcoef) + & & ovrnts*cmoment(3, sps) +! final part of the ks computation + if (computesepsensorks) then +! only calculate the log part if we are actually computing for separation for ks method. + ks_compd = ovrnts*globalvalsd(isepsensorks, sps)/(sepsenmaxrho*& +& globalvals(isepsensorks, sps)) + ks_comp = ovrnts*(sepsenmaxfamily(sps)+log(globalvals(& +& isepsensorks, sps))/sepsenmaxrho) + funcvaluesd(costfuncsepsensorks) = funcvaluesd(& +& costfuncsepsensorks) + ks_compd + funcvalues(costfuncsepsensorks) = funcvalues(costfuncsepsensorks& +& ) + ks_comp + arg1d = sepsensorsharpnessone*2*ks_compd + arg1 = 2*sepsensorsharpnessone*(ks_comp+sepsensoroffsetone) + temp0 = one + exp(arg1) + temp = globalvals(isepsensorareaks, sps)*ks_comp/temp0 + funcvaluesd(costfuncsepsensorarea) = funcvaluesd(& +& costfuncsepsensorarea) + ovrnts*one*(ks_comp*globalvalsd(& +& isepsensorareaks, sps)+globalvals(isepsensorareaks, sps)*& +& ks_compd-temp*exp(arg1)*arg1d)/temp0 + ovrnts*globalvalsd(& +& isepsensorarea, sps) + funcvalues(costfuncsepsensorarea) = funcvalues(& +& costfuncsepsensorarea) + ovrnts*one*temp + ovrnts*globalvals(& +& isepsensorarea, sps) + end if funcvaluesd(costfuncsepsensor) = funcvaluesd(costfuncsepsensor) + & & ovrnts*globalvalsd(isepsensor, sps) funcvalues(costfuncsepsensor) = funcvalues(costfuncsepsensor) + & @@ -701,8 +729,10 @@ subroutine getcostfunctions(globalvals, funcvalues) use flowvarrefstate, only : pref, rhoref, tref, lref, gammainf, & & pinf, uref, uinf use inputphysics, only : liftdirection, dragdirection, surfaceref,& -& machcoef, lengthref, alpha, beta, liftindex, cpmin_family, cpmin_rho - use inputcostfunctions, only : computecavitation +& machcoef, lengthref, alpha, beta, liftindex, cpmin_family, cpmin_rho& +& , sepsenmaxfamily, sepsenmaxrho + use inputcostfunctions, only : computecavitation, computesepsensorks& +& , sepsensorsharpnessone, sepsensoroffsetone use inputtsstabderiv, only : tsstability use utils_d, only : computetsderivatives use flowutils_d, only : getdirvector @@ -719,12 +749,13 @@ subroutine getcostfunctions(globalvals, funcvalues) real(kind=realtype) :: mavgptot, mavgttot, mavgrho, mavgps, mflow, & & mavgmn, mavga, mavgvx, mavgvy, mavgvz, garea, mavgvi, fxlift, fylift& & , fzlift - real(kind=realtype) :: vdotn, mag, u, v, w + real(kind=realtype) :: vdotn, mag, u, v, w, ks_comp integer(kind=inttype) :: sps real(kind=realtype), dimension(8) :: dcdq, dcdqdot real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot real(kind=realtype), dimension(8) :: coef0 intrinsic log + intrinsic exp intrinsic sqrt real(kind=realtype) :: arg1 ! factor used for time-averaged quantities. @@ -855,6 +886,19 @@ subroutine getcostfunctions(globalvals, funcvalues) & ovrnts*cmoment(2, sps) funcvalues(costfuncmomzcoef) = funcvalues(costfuncmomzcoef) + & & ovrnts*cmoment(3, sps) +! final part of the ks computation + if (computesepsensorks) then +! only calculate the log part if we are actually computing for separation for ks method. + ks_comp = ovrnts*(sepsenmaxfamily(sps)+log(globalvals(& +& isepsensorks, sps))/sepsenmaxrho) + funcvalues(costfuncsepsensorks) = funcvalues(costfuncsepsensorks& +& ) + ks_comp + arg1 = 2*sepsensorsharpnessone*(ks_comp+sepsensoroffsetone) + funcvalues(costfuncsepsensorarea) = funcvalues(& +& costfuncsepsensorarea) + ovrnts*globalvals(isepsensorareaks, & +& sps)*ks_comp*one/(one+exp(arg1)) + ovrnts*globalvals(& +& isepsensorarea, sps) + end if funcvalues(costfuncsepsensor) = funcvalues(costfuncsepsensor) + & & ovrnts*globalvals(isepsensor, sps) funcvalues(costfunccavitation) = funcvalues(costfunccavitation) + & @@ -1067,7 +1111,8 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) use inputcostfunctions use inputphysics, only : machcoef, machcoefd, pointref, pointrefd,& & veldirfreestream, veldirfreestreamd, equations, momentaxis, & -& cpmin_family, cpmin_rho, cavitationnumber +& cpmin_family, cpmin_rho, cavitationnumber, sepsenmaxfamily, & +& sepsenmaxrho use bcpointers_d implicit none ! input/output variables @@ -1081,13 +1126,18 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) real(kind=realtype), dimension(3) :: fpd, fvd, mpd, mvd real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz real(kind=realtype), dimension(3) :: cofsumfxd, cofsumfyd, cofsumfzd - real(kind=realtype) :: yplusmax, sepsensor, sepsensoravg(3), & -& cavitation, cpmin_ks_sum - real(kind=realtype) :: sepsensord, sepsensoravgd(3), cavitationd, & -& cpmin_ks_sumd + real(kind=realtype) :: yplusmax, sepsensorks, sepsensor, & +& sepsensoravg(3), sepsensorarea, cavitation, cpmin_ks_sum, & +& sepsensorareaks + real(kind=realtype) :: sepsensorksd, sepsensord, sepsensoravgd(3), & +& sepsensoraread, cavitationd, cpmin_ks_sumd, sepsensorareaksd integer(kind=inttype) :: i, j, ii, blk real(kind=realtype) :: pm1, fx, fy, fz, fn real(kind=realtype) :: pm1d, fxd, fyd, fzd + real(kind=realtype) :: vecttangential(3) + real(kind=realtype) :: vecttangentiald(3) + real(kind=realtype) :: vectdotproductfsnormal + real(kind=realtype) :: vectdotproductfsnormald real(kind=realtype) :: xc, xco, yc, yco, zc, zco, qf(3), r(3), n(3)& & , l real(kind=realtype) :: xcd, xcod, ycd, ycod, zcd, zcod, rd(3) @@ -1112,14 +1162,16 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) intrinsic mod intrinsic max intrinsic sqrt + intrinsic cos intrinsic exp real(kind=realtype) :: arg1 real(kind=realtype) :: arg1d real(kind=realtype) :: result1 real(kind=realtype) :: result1d real(kind=realtype) :: temp - real(kind=realtype) :: tempd real(kind=realtype) :: temp0 + real(kind=realtype) :: temp1 + real(kind=realtype) :: tempd select case (bcfaceid(mm)) case (imin, jmin, kmin) fact = -one @@ -1153,6 +1205,9 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) cofsumfz = zero yplusmax = zero sepsensor = zero + sepsensorks = zero + sepsensorarea = zero + sepsensorareaks = zero cavitation = zero cpmin_ks_sum = zero sepsensoravg = zero @@ -1167,8 +1222,12 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) cpmin_ks_sumd = 0.0_8 vd = 0.0_8 mpaxisd = 0.0_8 + vecttangentiald = 0.0_8 + sepsensorksd = 0.0_8 + sepsensorareaksd = 0.0_8 cperror2d = 0.0_8 fpd = 0.0_8 + sepsensoraread = 0.0_8 mpd = 0.0_8 cavitationd = 0.0_8 sepsensord = 0.0_8 @@ -1368,6 +1427,72 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) result1 = temp vd = (vd-v*result1d/(result1+1e-16))/(result1+1e-16) v = v/(result1+1e-16) + if (computesepsensorks) then +! freestream projection over the surface. + temp = bcdata(mm)%norm(i, j, 1) + temp0 = bcdata(mm)%norm(i, j, 2) + temp1 = bcdata(mm)%norm(i, j, 3) + vectdotproductfsnormald = temp*veldirfreestreamd(1) + temp0*& +& veldirfreestreamd(2) + temp1*veldirfreestreamd(3) + vectdotproductfsnormal = temp*veldirfreestream(1) + temp0*& +& veldirfreestream(2) + temp1*veldirfreestream(3) +! tangential vector on the surface, which is the freestream projected vector + temp1 = bcdata(mm)%norm(i, j, 1) + vecttangentiald(1) = veldirfreestreamd(1) - temp1*& +& vectdotproductfsnormald + vecttangential(1) = veldirfreestream(1) - temp1*& +& vectdotproductfsnormal + temp1 = bcdata(mm)%norm(i, j, 2) + vecttangentiald(2) = veldirfreestreamd(2) - temp1*& +& vectdotproductfsnormald + vecttangential(2) = veldirfreestream(2) - temp1*& +& vectdotproductfsnormal + temp1 = bcdata(mm)%norm(i, j, 3) + vecttangentiald(3) = veldirfreestreamd(3) - temp1*& +& vectdotproductfsnormald + vecttangential(3) = veldirfreestream(3) - temp1*& +& vectdotproductfsnormal + arg1d = 2*vecttangential(1)*vecttangentiald(1) + 2*& +& vecttangential(2)*vecttangentiald(2) + 2*vecttangential(3)*& +& vecttangentiald(3) + arg1 = vecttangential(1)**2 + vecttangential(2)**2 + & +& vecttangential(3)**2 + temp1 = sqrt(arg1) + if (arg1 .eq. 0.0_8) then + result1d = 0.0_8 + else + result1d = arg1d/(2.0*temp1) + end if + result1 = temp1 + vecttangentiald = (vecttangentiald-vecttangential*result1d/(& +& result1+1e-16))/(result1+1e-16) + vecttangential = vecttangential/(result1+1e-16) +! computing separation sensor +! velocity dot products + sensord = vecttangential(1)*vd(1) + v(1)*vecttangentiald(1) + & +& vecttangential(2)*vd(2) + v(2)*vecttangentiald(2) + & +& vecttangential(3)*vd(3) + v(3)*vecttangentiald(3) + sensor = v(1)*vecttangential(1) + v(2)*vecttangential(2) + v(3)*& +& vecttangential(3) +! sepsensor value + temp1 = cos(zero) - cos(degtorad*sepangledeviation) + 1e-16 + sensord = -(sensord/temp1) + sensor = (cos(degtorad*sepangledeviation)-sensor)/temp1 +! also do the ks-based spensenor max computation + call ksaggregationfunction_d(sensor, sensord, sepsenmaxfamily(& +& spectralsol), sepsenmaxrho, ks_exponent, & +& ks_exponentd) + sepsensorareaksd = sepsensorareaksd + blk*cellaread + sepsensorareaks = sepsensorareaks + blk*cellarea + arg1d = -(sepsensorsharpnesstwo*2*sensord) + arg1 = -(2*sepsensorsharpnesstwo*(sensor+sepsensoroffsettwo)) + temp1 = one + exp(arg1) + sepsensoraread = blk*one*(cellaread-cellarea*exp(arg1)*arg1d/& +& temp1)/temp1 + sepsensoraread + sepsensorarea = blk*one*(cellarea/temp1) + sepsensorarea + sepsensorksd = sepsensorksd + blk*ks_exponentd + sepsensorks = sepsensorks + ks_exponent*blk + end if ! dot product with free stream sensord = -(veldirfreestream(1)*vd(1)) - v(1)*veldirfreestreamd(1)& & - veldirfreestream(2)*vd(2) - v(2)*veldirfreestreamd(2) - & @@ -1377,9 +1502,9 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) !now run through a smooth heaviside function: arg1d = -(sepsensorsharpness*2*sensord) arg1 = -(2*sepsensorsharpness*(sensor-sepsensoroffset)) - temp = one/(one+exp(arg1)) - sensord = -(temp*exp(arg1)*arg1d/(one+exp(arg1))) - sensor = temp + temp1 = one/(one+exp(arg1)) + sensord = -(temp1*exp(arg1)*arg1d/(one+exp(arg1))) + sensor = temp1 ! and integrate over the area of this cell and save, blanking as we go. sensord = blk*(cellarea*sensord+sensor*cellaread) sensor = sensor*cellarea*blk @@ -1396,33 +1521,33 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) if (computecavitation) then plocald = pp2d(i, j) plocal = pp2(i, j) - temp = two/(gammainf*(machcoef*machcoef)) - tmpd = -(temp*2*machcoefd/machcoef) - tmp = temp + temp1 = two/(gammainf*(machcoef*machcoef)) + tmpd = -(temp1*2*machcoefd/machcoef) + tmp = temp1 cpd = (plocal-pinf)*tmpd + tmp*(plocald-pinfd) cp = tmp*(plocal-pinf) sensor1d = -cpd sensor1 = -cp - cavitationnumber arg1d = -(cavsensorsharpness*2*sensor1d) arg1 = 2*cavsensorsharpness*(-sensor1+cavsensoroffset) - temp = one + exp(arg1) - temp0 = sensor1**cavexponent/temp + temp1 = one + exp(arg1) + temp = sensor1**cavexponent/temp1 if (sensor1 .le. 0.0_8 .and. (cavexponent .eq. 0.0_8 .or. & & cavexponent .ne. int(cavexponent))) then tempd = 0.0_8 else tempd = cavexponent*sensor1**(cavexponent-1)*sensor1d end if - sensor1d = (tempd-temp0*exp(arg1)*arg1d)/temp - sensor1 = temp0 + sensor1d = (tempd-temp*exp(arg1)*arg1d)/temp1 + sensor1 = temp sensor1d = blk*(cellarea*sensor1d+sensor1*cellaread) sensor1 = sensor1*cellarea*blk cavitationd = cavitationd + sensor1d cavitation = cavitation + sensor1 ! also do the ks-based cpmin computation - temp0 = cpmin_rho*(cpmin_family(spectralsol)-cp) - ks_exponentd = -(exp(temp0)*cpmin_rho*cpd) - ks_exponent = exp(temp0) + temp1 = cpmin_rho*(cpmin_family(spectralsol)-cp) + ks_exponentd = -(exp(temp1)*cpmin_rho*cpd) + ks_exponent = exp(temp1) cpmin_ks_sumd = cpmin_ks_sumd + blk*ks_exponentd cpmin_ks_sum = cpmin_ks_sum + ks_exponent*blk end if @@ -1467,24 +1592,24 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) tauyz = viscsubface(mm)%tau(i, j, 6) ! compute the viscous force on the face. a minus sign ! is now present, due to the definition of this force. - temp0 = tauxx*ssi(i, j, 1) + tauxy*ssi(i, j, 2) + tauxz*ssi(i, j& + temp1 = tauxx*ssi(i, j, 1) + tauxy*ssi(i, j, 2) + tauxz*ssi(i, j& & , 3) fxd = -(fact*(pref*(ssi(i, j, 1)*tauxxd+tauxx*ssid(i, j, 1)+ssi(& & i, j, 2)*tauxyd+tauxy*ssid(i, j, 2)+ssi(i, j, 3)*tauxzd+tauxz*& -& ssid(i, j, 3))+temp0*prefd)) - fx = -(fact*(temp0*pref)) - temp0 = tauxy*ssi(i, j, 1) + tauyy*ssi(i, j, 2) + tauyz*ssi(i, j& +& ssid(i, j, 3))+temp1*prefd)) + fx = -(fact*(temp1*pref)) + temp1 = tauxy*ssi(i, j, 1) + tauyy*ssi(i, j, 2) + tauyz*ssi(i, j& & , 3) fyd = -(fact*(pref*(ssi(i, j, 1)*tauxyd+tauxy*ssid(i, j, 1)+ssi(& & i, j, 2)*tauyyd+tauyy*ssid(i, j, 2)+ssi(i, j, 3)*tauyzd+tauyz*& -& ssid(i, j, 3))+temp0*prefd)) - fy = -(fact*(temp0*pref)) - temp0 = tauxz*ssi(i, j, 1) + tauyz*ssi(i, j, 2) + tauzz*ssi(i, j& +& ssid(i, j, 3))+temp1*prefd)) + fy = -(fact*(temp1*pref)) + temp1 = tauxz*ssi(i, j, 1) + tauyz*ssi(i, j, 2) + tauzz*ssi(i, j& & , 3) fzd = -(fact*(pref*(ssi(i, j, 1)*tauxzd+tauxz*ssid(i, j, 1)+ssi(& & i, j, 2)*tauyzd+tauyz*ssid(i, j, 2)+ssi(i, j, 3)*tauzzd+tauzz*& -& ssid(i, j, 3))+temp0*prefd)) - fz = -(fact*(temp0*pref)) +& ssid(i, j, 3))+temp1*prefd)) + fz = -(fact*(temp1*pref)) ! compute the coordinates of the centroid of the face ! relative from the moment reference point. due to the ! usage of pointers for xx and offset of 1 is present, @@ -1648,6 +1773,17 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm) & +2) + cofsumfz localvaluesd(isepsensor) = localvaluesd(isepsensor) + sepsensord localvalues(isepsensor) = localvalues(isepsensor) + sepsensor + localvaluesd(isepsensorks) = localvaluesd(isepsensorks) + & +& sepsensorksd + localvalues(isepsensorks) = localvalues(isepsensorks) + sepsensorks + localvaluesd(isepsensorareaks) = localvaluesd(isepsensorareaks) + & +& sepsensorareaksd + localvalues(isepsensorareaks) = localvalues(isepsensorareaks) + & +& sepsensorareaks + localvaluesd(isepsensorarea) = localvaluesd(isepsensorarea) + & +& sepsensoraread + localvalues(isepsensorarea) = localvalues(isepsensorarea) + & +& sepsensorarea localvaluesd(icavitation) = localvaluesd(icavitation) + cavitationd localvalues(icavitation) = localvalues(icavitation) + cavitation localvaluesd(icpmin) = localvaluesd(icpmin) + cpmin_ks_sumd @@ -1681,7 +1817,8 @@ subroutine wallintegrationface(localvalues, mm) use flowvarrefstate use inputcostfunctions use inputphysics, only : machcoef, pointref, veldirfreestream, & -& equations, momentaxis, cpmin_family, cpmin_rho, cavitationnumber +& equations, momentaxis, cpmin_family, cpmin_rho, cavitationnumber, & +& sepsenmaxfamily, sepsenmaxrho use bcpointers_d implicit none ! input/output variables @@ -1691,10 +1828,13 @@ subroutine wallintegrationface(localvalues, mm) ! local variables. real(kind=realtype), dimension(3) :: fp, fv, mp, mv real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz - real(kind=realtype) :: yplusmax, sepsensor, sepsensoravg(3), & -& cavitation, cpmin_ks_sum + real(kind=realtype) :: yplusmax, sepsensorks, sepsensor, & +& sepsensoravg(3), sepsensorarea, cavitation, cpmin_ks_sum, & +& sepsensorareaks integer(kind=inttype) :: i, j, ii, blk real(kind=realtype) :: pm1, fx, fy, fz, fn + real(kind=realtype) :: vecttangential(3) + real(kind=realtype) :: vectdotproductfsnormal real(kind=realtype) :: xc, xco, yc, yco, zc, zco, qf(3), r(3), n(3)& & , l real(kind=realtype) :: fact, rho, mul, yplus, dwall @@ -1710,6 +1850,7 @@ subroutine wallintegrationface(localvalues, mm) intrinsic mod intrinsic max intrinsic sqrt + intrinsic cos intrinsic exp real(kind=realtype) :: arg1 real(kind=realtype) :: result1 @@ -1742,6 +1883,9 @@ subroutine wallintegrationface(localvalues, mm) cofsumfz = zero yplusmax = zero sepsensor = zero + sepsensorks = zero + sepsensorarea = zero + sepsensorareaks = zero cavitation = zero cpmin_ks_sum = zero sepsensoravg = zero @@ -1870,6 +2014,37 @@ subroutine wallintegrationface(localvalues, mm) arg1 = v(1)**2 + v(2)**2 + v(3)**2 result1 = sqrt(arg1) v = v/(result1+1e-16) + if (computesepsensorks) then +! freestream projection over the surface. + vectdotproductfsnormal = veldirfreestream(1)*bcdata(mm)%norm(i, & +& j, 1) + veldirfreestream(2)*bcdata(mm)%norm(i, j, 2) + & +& veldirfreestream(3)*bcdata(mm)%norm(i, j, 3) +! tangential vector on the surface, which is the freestream projected vector + vecttangential(1) = veldirfreestream(1) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 1) + vecttangential(2) = veldirfreestream(2) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 2) + vecttangential(3) = veldirfreestream(3) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 3) + arg1 = vecttangential(1)**2 + vecttangential(2)**2 + & +& vecttangential(3)**2 + result1 = sqrt(arg1) + vecttangential = vecttangential/(result1+1e-16) +! computing separation sensor +! velocity dot products + sensor = v(1)*vecttangential(1) + v(2)*vecttangential(2) + v(3)*& +& vecttangential(3) +! sepsensor value + sensor = (cos(degtorad*sepangledeviation)-sensor)/(-cos(degtorad& +& *sepangledeviation)+cos(zero)+1e-16) +! also do the ks-based spensenor max computation + call ksaggregationfunction(sensor, sepsenmaxfamily(spectralsol)& +& , sepsenmaxrho, ks_exponent) + sepsensorareaks = sepsensorareaks + blk*cellarea + arg1 = -(2*sepsensorsharpnesstwo*(sensor+sepsensoroffsettwo)) + sepsensorarea = cellarea*blk*one/(one+exp(arg1)) + sepsensorarea + sepsensorks = sepsensorks + ks_exponent*blk + end if ! dot product with free stream sensor = -(v(1)*veldirfreestream(1)+v(2)*veldirfreestream(2)+v(3)*& & veldirfreestream(3)) @@ -2041,6 +2216,11 @@ subroutine wallintegrationface(localvalues, mm) localvalues(icoforcez:icoforcez+2) = localvalues(icoforcez:icoforcez& & +2) + cofsumfz localvalues(isepsensor) = localvalues(isepsensor) + sepsensor + localvalues(isepsensorks) = localvalues(isepsensorks) + sepsensorks + localvalues(isepsensorareaks) = localvalues(isepsensorareaks) + & +& sepsensorareaks + localvalues(isepsensorarea) = localvalues(isepsensorarea) + & +& sepsensorarea localvalues(icavitation) = localvalues(icavitation) + cavitation localvalues(icpmin) = localvalues(icpmin) + cpmin_ks_sum localvalues(isepavg:isepavg+2) = localvalues(isepavg:isepavg+2) + & @@ -2050,6 +2230,30 @@ subroutine wallintegrationface(localvalues, mm) localvalues(icperror2) = localvalues(icperror2) + cperror2 end subroutine wallintegrationface +! differentiation of ksaggregationfunction in forward (tangent) mode (with options i4 dr8 r8): +! variations of useful results: ks_g +! with respect to varying inputs: g + subroutine ksaggregationfunction_d(g, gd, max_g, g_rho, ks_g, ks_gd) + use precision + implicit none + real(kind=realtype), intent(inout) :: ks_g + real(kind=realtype), intent(inout) :: ks_gd + real(kind=realtype), intent(in) :: g, max_g, g_rho + real(kind=realtype), intent(in) :: gd + intrinsic exp + ks_gd = exp(g_rho*(g-max_g))*g_rho*gd + ks_g = exp(g_rho*(g-max_g)) + end subroutine ksaggregationfunction_d + + subroutine ksaggregationfunction(g, max_g, g_rho, ks_g) + use precision + implicit none + real(kind=realtype), intent(inout) :: ks_g + real(kind=realtype), intent(in) :: g, max_g, g_rho + intrinsic exp + ks_g = exp(g_rho*(g-max_g)) + end subroutine ksaggregationfunction + ! differentiation of flowintegrationface in forward (tangent) mode (with options i4 dr8 r8): ! variations of useful results: localvalues ! with respect to varying inputs: timeref tref pref rgas rhoref diff --git a/src/adjoint/outputReverse/surfaceIntegrations_b.f90 b/src/adjoint/outputReverse/surfaceIntegrations_b.f90 index 937c195db..008907a65 100644 --- a/src/adjoint/outputReverse/surfaceIntegrations_b.f90 +++ b/src/adjoint/outputReverse/surfaceIntegrations_b.f90 @@ -21,8 +21,9 @@ subroutine getcostfunctions_b(globalvals, globalvalsd, funcvalues, & use inputphysics, only : liftdirection, liftdirectiond, & & dragdirection, dragdirectiond, surfaceref, machcoef, machcoefd, & & lengthref, alpha, alphad, beta, betad, liftindex, cpmin_family, & -& cpmin_rho - use inputcostfunctions, only : computecavitation +& cpmin_rho, sepsenmaxfamily, sepsenmaxrho + use inputcostfunctions, only : computecavitation, computesepsensorks& +& , sepsensorsharpnessone, sepsensoroffsetone use inputtsstabderiv, only : tsstability use utils_b, only : computetsderivatives use flowutils_b, only : getdirvector @@ -48,15 +49,18 @@ subroutine getcostfunctions_b(globalvals, globalvalsd, funcvalues, & real(kind=realtype) :: mavgptotd, mavgttotd, mavgrhod, mavgpsd, & & mflowd, mavgmnd, mavgad, mavgvxd, mavgvyd, mavgvzd, garead, mavgvid& & , fxliftd, fyliftd, fzliftd - real(kind=realtype) :: vdotn, mag, u, v, w + real(kind=realtype) :: vdotn, mag, u, v, w, ks_comp + real(kind=realtype) :: ks_compd integer(kind=inttype) :: sps real(kind=realtype), dimension(8) :: dcdq, dcdqdot real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot real(kind=realtype), dimension(8) :: coef0 intrinsic log + intrinsic exp intrinsic sqrt real(kind=realtype) :: temp real(kind=realtype) :: temp0 + real(kind=realtype) :: temp1 real(kind=realtype) :: tempd real(kind=realtype) :: tmp real(kind=realtype) :: tmpd @@ -97,6 +101,7 @@ subroutine getcostfunctions_b(globalvals, globalvalsd, funcvalues, & real(kind=realtype) :: tmpd16 real(kind=realtype) :: tmp17 real(kind=realtype) :: tmpd17 + real(kind=realtype) :: tempd1 integer :: branch ! factor used for time-averaged quantities. ovrnts = one/ntimeintervalsspectral @@ -230,6 +235,18 @@ subroutine getcostfunctions_b(globalvals, globalvalsd, funcvalues, & & ovrnts*cmoment(2, sps) funcvalues(costfuncmomzcoef) = funcvalues(costfuncmomzcoef) + & & ovrnts*cmoment(3, sps) +! final part of the ks computation + if (computesepsensorks) then +! only calculate the log part if we are actually computing for separation for ks method. + ks_comp = ovrnts*(sepsenmaxfamily(sps)+log(globalvals(& +& isepsensorks, sps))/sepsenmaxrho) + funcvalues(costfuncsepsensorks) = funcvalues(costfuncsepsensorks& +& ) + ks_comp + funcvalues(costfuncsepsensorarea) = funcvalues(& +& costfuncsepsensorarea) + ovrnts*globalvals(isepsensorareaks, & +& sps)*ks_comp*one/(one+exp(2*sepsensorsharpnessone*(ks_comp+& +& sepsensoroffsetone))) + ovrnts*globalvals(isepsensorarea, sps) + end if funcvalues(costfuncsepsensor) = funcvalues(costfuncsepsensor) + & & ovrnts*globalvals(isepsensor, sps) funcvalues(costfunccavitation) = funcvalues(costfunccavitation) + & @@ -779,6 +796,15 @@ subroutine getcostfunctions_b(globalvals, globalvalsd, funcvalues, & ! fy ! fz ! ------------ +! final part of the ks computation + if (computesepsensorks) then +! only calculate the log part if we are actually computing for separation for ks method. + ks_comp = ovrnts*(sepsenmaxfamily(sps)+log(globalvals(& +& isepsensorks, sps))/sepsenmaxrho) + call pushcontrol1b(0) + else + call pushcontrol1b(1) + end if ! final part of the ks computation if (computecavitation) then call pushcontrol1b(0) @@ -878,6 +904,23 @@ subroutine getcostfunctions_b(globalvals, globalvalsd, funcvalues, & & ovrnts*funcvaluesd(costfunccavitation) globalvalsd(isepsensor, sps) = globalvalsd(isepsensor, sps) + & & ovrnts*funcvaluesd(costfuncsepsensor) + call popcontrol1b(branch) + if (branch .eq. 0) then + temp1 = 2*sepsensorsharpnessone*(sepsensoroffsetone+ks_comp) + temp0 = one + exp(temp1) + tempd1 = ovrnts*one*funcvaluesd(costfuncsepsensorarea)/temp0 + globalvalsd(isepsensorarea, sps) = globalvalsd(isepsensorarea& +& , sps) + ovrnts*funcvaluesd(costfuncsepsensorarea) + globalvalsd(isepsensorareaks, sps) = globalvalsd(& +& isepsensorareaks, sps) + ks_comp*tempd1 + ks_compd = (globalvals(isepsensorareaks, sps)-& +& sepsensorsharpnessone*2*exp(temp1)*globalvals(& +& isepsensorareaks, sps)*ks_comp/temp0)*tempd1 + funcvaluesd(& +& costfuncsepsensorks) + globalvalsd(isepsensorks, sps) = globalvalsd(isepsensorks, sps& +& ) + ovrnts*ks_compd/(globalvals(isepsensorks, sps)*& +& sepsenmaxrho) + end if cmomentd(3, sps) = cmomentd(3, sps) + ovrnts*funcvaluesd(& & costfuncmomzcoef) cmomentd(2, sps) = cmomentd(2, sps) + ovrnts*funcvaluesd(& @@ -1023,8 +1066,10 @@ subroutine getcostfunctions(globalvals, funcvalues) use flowvarrefstate, only : pref, rhoref, tref, lref, gammainf, & & pinf, uref, uinf use inputphysics, only : liftdirection, dragdirection, surfaceref,& -& machcoef, lengthref, alpha, beta, liftindex, cpmin_family, cpmin_rho - use inputcostfunctions, only : computecavitation +& machcoef, lengthref, alpha, beta, liftindex, cpmin_family, cpmin_rho& +& , sepsenmaxfamily, sepsenmaxrho + use inputcostfunctions, only : computecavitation, computesepsensorks& +& , sepsensorsharpnessone, sepsensoroffsetone use inputtsstabderiv, only : tsstability use utils_b, only : computetsderivatives use flowutils_b, only : getdirvector @@ -1041,12 +1086,13 @@ subroutine getcostfunctions(globalvals, funcvalues) real(kind=realtype) :: mavgptot, mavgttot, mavgrho, mavgps, mflow, & & mavgmn, mavga, mavgvx, mavgvy, mavgvz, garea, mavgvi, fxlift, fylift& & , fzlift - real(kind=realtype) :: vdotn, mag, u, v, w + real(kind=realtype) :: vdotn, mag, u, v, w, ks_comp integer(kind=inttype) :: sps real(kind=realtype), dimension(8) :: dcdq, dcdqdot real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot real(kind=realtype), dimension(8) :: coef0 intrinsic log + intrinsic exp intrinsic sqrt ! factor used for time-averaged quantities. ovrnts = one/ntimeintervalsspectral @@ -1176,6 +1222,18 @@ subroutine getcostfunctions(globalvals, funcvalues) & ovrnts*cmoment(2, sps) funcvalues(costfuncmomzcoef) = funcvalues(costfuncmomzcoef) + & & ovrnts*cmoment(3, sps) +! final part of the ks computation + if (computesepsensorks) then +! only calculate the log part if we are actually computing for separation for ks method. + ks_comp = ovrnts*(sepsenmaxfamily(sps)+log(globalvals(& +& isepsensorks, sps))/sepsenmaxrho) + funcvalues(costfuncsepsensorks) = funcvalues(costfuncsepsensorks& +& ) + ks_comp + funcvalues(costfuncsepsensorarea) = funcvalues(& +& costfuncsepsensorarea) + ovrnts*globalvals(isepsensorareaks, & +& sps)*ks_comp*one/(one+exp(2*sepsensorsharpnessone*(ks_comp+& +& sepsensoroffsetone))) + ovrnts*globalvals(isepsensorarea, sps) + end if funcvalues(costfuncsepsensor) = funcvalues(costfuncsepsensor) + & & ovrnts*globalvals(isepsensor, sps) funcvalues(costfunccavitation) = funcvalues(costfunccavitation) + & @@ -1389,7 +1447,8 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) use inputcostfunctions use inputphysics, only : machcoef, machcoefd, pointref, pointrefd,& & veldirfreestream, veldirfreestreamd, equations, momentaxis, & -& cpmin_family, cpmin_rho, cavitationnumber +& cpmin_family, cpmin_rho, cavitationnumber, sepsenmaxfamily, & +& sepsenmaxrho use bcpointers_b implicit none ! input/output variables @@ -1403,13 +1462,18 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) real(kind=realtype), dimension(3) :: fpd, fvd, mpd, mvd real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz real(kind=realtype), dimension(3) :: cofsumfxd, cofsumfyd, cofsumfzd - real(kind=realtype) :: yplusmax, sepsensor, sepsensoravg(3), & -& cavitation, cpmin_ks_sum - real(kind=realtype) :: sepsensord, sepsensoravgd(3), cavitationd, & -& cpmin_ks_sumd + real(kind=realtype) :: yplusmax, sepsensorks, sepsensor, & +& sepsensoravg(3), sepsensorarea, cavitation, cpmin_ks_sum, & +& sepsensorareaks + real(kind=realtype) :: sepsensorksd, sepsensord, sepsensoravgd(3), & +& sepsensoraread, cavitationd, cpmin_ks_sumd, sepsensorareaksd integer(kind=inttype) :: i, j, ii, blk real(kind=realtype) :: pm1, fx, fy, fz, fn real(kind=realtype) :: pm1d, fxd, fyd, fzd + real(kind=realtype) :: vecttangential(3) + real(kind=realtype) :: vecttangentiald(3) + real(kind=realtype) :: vectdotproductfsnormal + real(kind=realtype) :: vectdotproductfsnormald real(kind=realtype) :: xc, xco, yc, yco, zc, zco, qf(3), r(3), n(3)& & , l real(kind=realtype) :: xcd, xcod, ycd, ycod, zcd, zcod, rd(3) @@ -1434,13 +1498,18 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) intrinsic mod intrinsic max intrinsic sqrt + intrinsic cos intrinsic exp real(kind=realtype) :: temp real(kind=realtype) :: tempd real(kind=realtype), dimension(3) :: tmp0 real(kind=realtype), dimension(3) :: tmpd0 real(kind=realtype) :: temp0 + real(kind=realtype), dimension(3) :: tmp1 + real(kind=realtype), dimension(3) :: tmpd1 real(kind=realtype) :: tempd0 + real(kind=realtype) :: temp1 + real(kind=realtype) :: tempd1 integer :: branch select case (bcfaceid(mm)) case (imin, jmin, kmin) @@ -1465,6 +1534,7 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) call pushreal8array(n, 3) call pushreal8array(r, 3) call pushreal8array(v, 3) + call pushreal8array(vecttangential, 3) !$fwd-of ii-loop ! ! integrate the inviscid contribution over the solid walls, @@ -1584,6 +1654,36 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) v(2) = ww2(i, j, ivy) v(3) = ww2(i, j, ivz) v = v/(sqrt(v(1)**2+v(2)**2+v(3)**2)+1e-16) + if (computesepsensorks) then +! freestream projection over the surface. + vectdotproductfsnormal = veldirfreestream(1)*bcdata(mm)%norm(i, & +& j, 1) + veldirfreestream(2)*bcdata(mm)%norm(i, j, 2) + & +& veldirfreestream(3)*bcdata(mm)%norm(i, j, 3) +! tangential vector on the surface, which is the freestream projected vector + vecttangential(1) = veldirfreestream(1) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 1) + vecttangential(2) = veldirfreestream(2) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 2) + vecttangential(3) = veldirfreestream(3) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 3) + vecttangential = vecttangential/(sqrt(vecttangential(1)**2+& +& vecttangential(2)**2+vecttangential(3)**2)+1e-16) +! computing separation sensor +! velocity dot products + sensor = v(1)*vecttangential(1) + v(2)*vecttangential(2) + v(3)*& +& vecttangential(3) +! sepsensor value + sensor = (cos(degtorad*sepangledeviation)-sensor)/(-cos(degtorad& +& *sepangledeviation)+cos(zero)+1e-16) +! also do the ks-based spensenor max computation + call ksaggregationfunction(sensor, sepsenmaxfamily(spectralsol)& +& , sepsenmaxrho, ks_exponent) + sepsensorareaks = sepsensorareaks + blk*cellarea + sepsensorarea = cellarea*blk*one/(one+exp(-(2*& +& sepsensorsharpnesstwo*(sensor+sepsensoroffsettwo)))) + & +& sepsensorarea + sepsensorks = sepsensorks + ks_exponent*blk + end if ! dot product with free stream sensor = -(v(1)*veldirfreestream(1)+v(2)*veldirfreestream(2)+v(3)*& & veldirfreestream(3)) @@ -1629,6 +1729,9 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) sepsensoravgd = localvaluesd(isepavg:isepavg+2) cpmin_ks_sumd = localvaluesd(icpmin) cavitationd = localvaluesd(icavitation) + sepsensoraread = localvaluesd(isepsensorarea) + sepsensorareaksd = localvaluesd(isepsensorareaks) + sepsensorksd = localvaluesd(isepsensorks) sepsensord = localvaluesd(isepsensor) cofsumfzd = 0.0_8 cofsumfzd = localvaluesd(icoforcez:icoforcez+2) @@ -1728,10 +1831,10 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) mxd = blk*mvd(1) myd = blk*mvd(2) mzd = blk*mvd(3) - tempd0 = blk*mvaxisd - m0xd = n(1)*tempd0 - m0yd = n(2)*tempd0 - m0zd = n(3)*tempd0 + tempd1 = blk*mvaxisd + m0xd = n(1)*tempd1 + m0yd = n(2)*tempd1 + m0zd = n(3)*tempd1 fzd = bcdatad(mm)%fv(i, j, 3) + r(2)*m0xd - r(1)*m0yd + zco*blk*& & cofsumfzd(3) + yco*blk*cofsumfzd(2) + xco*blk*cofsumfzd(1) + & & yc*mxd - xc*myd + blk*fvd(3) @@ -1747,91 +1850,91 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) rd(1) = rd(1) + fy*m0zd - fz*m0yd rd(2) = rd(2) + fz*m0xd - fx*m0zd rd(3) = rd(3) + fx*m0yd - fy*m0xd - tempd0 = fourth*rd(3) + tempd1 = fourth*rd(3) rd(3) = 0.0_8 - xxd(i, j, 3) = xxd(i, j, 3) + tempd0 - xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd0 - xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd0 - xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd0 - tempd0 = fourth*rd(2) + xxd(i, j, 3) = xxd(i, j, 3) + tempd1 + xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd1 + xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd1 + xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd1 + tempd1 = fourth*rd(2) rd(2) = 0.0_8 - xxd(i, j, 2) = xxd(i, j, 2) + tempd0 - xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd0 - xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd0 - xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd0 - tempd0 = fourth*rd(1) + xxd(i, j, 2) = xxd(i, j, 2) + tempd1 + xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd1 + xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd1 + xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd1 + tempd1 = fourth*rd(1) rd(1) = 0.0_8 - xxd(i, j, 1) = xxd(i, j, 1) + tempd0 - xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd0 - xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd0 - xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd0 + xxd(i, j, 1) = xxd(i, j, 1) + tempd1 + xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd1 + xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd1 + xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd1 zcod = fz*blk*cofsumfzd(3) + fy*blk*cofsumfyd(3) + fx*blk*& & cofsumfxd(3) ycod = fz*blk*cofsumfzd(2) + fy*blk*cofsumfyd(2) + fx*blk*& & cofsumfxd(2) xcod = fz*blk*cofsumfzd(1) + fy*blk*cofsumfyd(1) + fx*blk*& & cofsumfxd(1) - tempd0 = fourth*zcod - xxd(i, j, 3) = xxd(i, j, 3) + tempd0 - xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd0 - xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd0 - xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd0 - tempd0 = fourth*ycod - xxd(i, j, 2) = xxd(i, j, 2) + tempd0 - xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd0 - xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd0 - xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd0 - tempd0 = fourth*xcod - xxd(i, j, 1) = xxd(i, j, 1) + tempd0 - xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd0 - xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd0 - xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd0 + tempd1 = fourth*zcod + xxd(i, j, 3) = xxd(i, j, 3) + tempd1 + xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd1 + xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd1 + xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd1 + tempd1 = fourth*ycod + xxd(i, j, 2) = xxd(i, j, 2) + tempd1 + xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd1 + xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd1 + xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd1 + tempd1 = fourth*xcod + xxd(i, j, 1) = xxd(i, j, 1) + tempd1 + xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd1 + xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd1 + xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd1 xcd = fy*mzd - fz*myd ycd = fz*mxd - fx*mzd zcd = fx*myd - fy*mxd - tempd0 = fourth*zcd + tempd1 = fourth*zcd refpointd(3) = refpointd(3) - zcd - xxd(i, j, 3) = xxd(i, j, 3) + tempd0 - xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd0 - xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd0 - xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd0 - tempd0 = fourth*ycd + xxd(i, j, 3) = xxd(i, j, 3) + tempd1 + xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd1 + xxd(i, j+1, 3) = xxd(i, j+1, 3) + tempd1 + xxd(i+1, j+1, 3) = xxd(i+1, j+1, 3) + tempd1 + tempd1 = fourth*ycd refpointd(2) = refpointd(2) - ycd - xxd(i, j, 2) = xxd(i, j, 2) + tempd0 - xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd0 - xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd0 - xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd0 - tempd0 = fourth*xcd + xxd(i, j, 2) = xxd(i, j, 2) + tempd1 + xxd(i+1, j, 2) = xxd(i+1, j, 2) + tempd1 + xxd(i, j+1, 2) = xxd(i, j+1, 2) + tempd1 + xxd(i+1, j+1, 2) = xxd(i+1, j+1, 2) + tempd1 + tempd1 = fourth*xcd refpointd(1) = refpointd(1) - xcd - xxd(i, j, 1) = xxd(i, j, 1) + tempd0 - xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd0 - xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd0 - xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd0 - tempd0 = -(pref*fact*fzd) + xxd(i, j, 1) = xxd(i, j, 1) + tempd1 + xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd1 + xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd1 + xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd1 + tempd1 = -(pref*fact*fzd) prefd = prefd - (tauxz*ssi(i, j, 1)+tauyz*ssi(i, j, 2)+tauzz*ssi& & (i, j, 3))*fact*fzd - (tauxy*ssi(i, j, 1)+tauyy*ssi(i, j, 2)+& & tauyz*ssi(i, j, 3))*fact*fyd - (tauxx*ssi(i, j, 1)+tauxy*ssi(i& & , j, 2)+tauxz*ssi(i, j, 3))*fact*fxd - tauxzd = ssi(i, j, 1)*tempd0 - ssid(i, j, 1) = ssid(i, j, 1) + tauxz*tempd0 - tauyzd = ssi(i, j, 2)*tempd0 - ssid(i, j, 2) = ssid(i, j, 2) + tauyz*tempd0 - tauzzd = ssi(i, j, 3)*tempd0 - ssid(i, j, 3) = ssid(i, j, 3) + tauzz*tempd0 - tempd0 = -(pref*fact*fyd) - tauxyd = ssi(i, j, 1)*tempd0 - ssid(i, j, 1) = ssid(i, j, 1) + tauxy*tempd0 - tauyyd = ssi(i, j, 2)*tempd0 - ssid(i, j, 2) = ssid(i, j, 2) + tauyy*tempd0 - tauyzd = tauyzd + ssi(i, j, 3)*tempd0 - ssid(i, j, 3) = ssid(i, j, 3) + tauyz*tempd0 - tempd0 = -(pref*fact*fxd) - tauxxd = ssi(i, j, 1)*tempd0 - ssid(i, j, 1) = ssid(i, j, 1) + tauxx*tempd0 - tauxyd = tauxyd + ssi(i, j, 2)*tempd0 - ssid(i, j, 2) = ssid(i, j, 2) + tauxy*tempd0 - tauxzd = tauxzd + ssi(i, j, 3)*tempd0 - ssid(i, j, 3) = ssid(i, j, 3) + tauxz*tempd0 + tauxzd = ssi(i, j, 1)*tempd1 + ssid(i, j, 1) = ssid(i, j, 1) + tauxz*tempd1 + tauyzd = ssi(i, j, 2)*tempd1 + ssid(i, j, 2) = ssid(i, j, 2) + tauyz*tempd1 + tauzzd = ssi(i, j, 3)*tempd1 + ssid(i, j, 3) = ssid(i, j, 3) + tauzz*tempd1 + tempd1 = -(pref*fact*fyd) + tauxyd = ssi(i, j, 1)*tempd1 + ssid(i, j, 1) = ssid(i, j, 1) + tauxy*tempd1 + tauyyd = ssi(i, j, 2)*tempd1 + ssid(i, j, 2) = ssid(i, j, 2) + tauyy*tempd1 + tauyzd = tauyzd + ssi(i, j, 3)*tempd1 + ssid(i, j, 3) = ssid(i, j, 3) + tauyz*tempd1 + tempd1 = -(pref*fact*fxd) + tauxxd = ssi(i, j, 1)*tempd1 + ssid(i, j, 1) = ssid(i, j, 1) + tauxx*tempd1 + tauxyd = tauxyd + ssi(i, j, 2)*tempd1 + ssid(i, j, 2) = ssid(i, j, 2) + tauxy*tempd1 + tauxzd = tauxzd + ssi(i, j, 3)*tempd1 + ssid(i, j, 3) = ssid(i, j, 3) + tauxz*tempd1 viscsubfaced(mm)%tau(i, j, 6) = viscsubfaced(mm)%tau(i, j, 6) + & & tauyzd viscsubfaced(mm)%tau(i, j, 5) = viscsubfaced(mm)%tau(i, j, 5) + & @@ -1851,6 +1954,8 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) refpointd = 0.0_8 end if vd = 0.0_8 + vecttangentiald = 0.0_8 + call popreal8array(vecttangential, 3) call popreal8array(v, 3) call popreal8array(r, 3) call popreal8array(n, 3) @@ -1931,7 +2036,36 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) tmp0 = v/(sqrt(v(1)**2+v(2)**2+v(3)**2)+1e-16) call pushreal8array(v, 3) v = tmp0 + if (computesepsensorks) then +! freestream projection over the surface. + vectdotproductfsnormal = veldirfreestream(1)*bcdata(mm)%norm(i, & +& j, 1) + veldirfreestream(2)*bcdata(mm)%norm(i, j, 2) + & +& veldirfreestream(3)*bcdata(mm)%norm(i, j, 3) +! tangential vector on the surface, which is the freestream projected vector + vecttangential(1) = veldirfreestream(1) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 1) + vecttangential(2) = veldirfreestream(2) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 2) + vecttangential(3) = veldirfreestream(3) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 3) + tmp1 = vecttangential/(sqrt(vecttangential(1)**2+vecttangential(& +& 2)**2+vecttangential(3)**2)+1e-16) + call pushreal8array(vecttangential, 3) + vecttangential = tmp1 +! computing separation sensor +! velocity dot products + sensor = v(1)*vecttangential(1) + v(2)*vecttangential(2) + v(3)*& +& vecttangential(3) +! sepsensor value + sensor = (cos(degtorad*sepangledeviation)-sensor)/(-cos(degtorad& +& *sepangledeviation)+cos(zero)+1e-16) +! also do the ks-based spensenor max computation + call pushcontrol1b(0) + else + call pushcontrol1b(1) + end if ! dot product with free stream + call pushreal8(sensor) sensor = -(v(1)*veldirfreestream(1)+v(2)*veldirfreestream(2)+v(3)*& & veldirfreestream(3)) !now run through a smooth heaviside function: @@ -1956,16 +2090,16 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) cellaread = sensor1*blk*sensor1d sensor1d = cellarea*blk*sensor1d call popreal8(sensor1) - temp0 = 2*cavsensorsharpness*(cavsensoroffset-sensor1) - temp = one + exp(temp0) + temp1 = 2*cavsensorsharpness*(cavsensoroffset-sensor1) + temp0 = one + exp(temp1) if (sensor1 .le. 0.0_8 .and. (cavexponent .eq. 0.0_8 .or. & & cavexponent .ne. int(cavexponent))) then - sensor1d = cavsensorsharpness*2*exp(temp0)*sensor1**& -& cavexponent*sensor1d/temp**2 + sensor1d = cavsensorsharpness*2*exp(temp1)*sensor1**& +& cavexponent*sensor1d/temp0**2 else - sensor1d = (cavexponent*sensor1**(cavexponent-1)/temp+& -& cavsensorsharpness*2*exp(temp0)*sensor1**cavexponent/temp**2& -& )*sensor1d + sensor1d = (cavexponent*sensor1**(cavexponent-1)/temp0+& +& cavsensorsharpness*2*exp(temp1)*sensor1**cavexponent/temp0**& +& 2)*sensor1d end if ks_exponentd = blk*cpmin_ks_sumd cpd = -(cpmin_rho*exp(cpmin_rho*(cpmin_family(spectralsol)-cp))*& @@ -1973,38 +2107,86 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) tmpd = (plocal-pinf)*cpd plocald = tmp*cpd pinfd = pinfd - tmp*cpd - temp0 = gammainf*(machcoef*machcoef) - machcoefd = machcoefd - 2*machcoef*gammainf*two*tmpd/temp0**2 + temp1 = gammainf*(machcoef*machcoef) + machcoefd = machcoefd - 2*machcoef*gammainf*two*tmpd/temp1**2 tmp = two/(gammainf*pinf*machcoef*machcoef) pp2d(i, j) = pp2d(i, j) + plocald else cellaread = 0.0_8 end if - mxd = blk*mpd(1) - myd = blk*mpd(2) - mzd = blk*mpd(3) sensord = zco*sepsensoravgd(3) + yco*sepsensoravgd(2) + xco*& & sepsensoravgd(1) + sepsensord - zcod = sensor*sepsensoravgd(3) + fz*blk*cofsumfzd(3) + fy*blk*& -& cofsumfyd(3) + fx*blk*cofsumfxd(3) - ycod = sensor*sepsensoravgd(2) + fz*blk*cofsumfzd(2) + fy*blk*& -& cofsumfyd(2) + fx*blk*cofsumfxd(2) - xcod = sensor*sepsensoravgd(1) + fz*blk*cofsumfzd(1) + fy*blk*& -& cofsumfyd(1) + fx*blk*cofsumfxd(1) + zcod = sensor*sepsensoravgd(3) + ycod = sensor*sepsensoravgd(2) + xcod = sensor*sepsensoravgd(1) call popreal8(sensor) - cellaread = cellaread + sensor*blk*sensord + bcdatad(mm)%area(i, j& -& ) + cellaread = cellaread + sensor*blk*sensord sensord = cellarea*blk*sensord call popreal8(sensor) - temp0 = -(2*sepsensorsharpness*(sensor-sepsensoroffset)) - temp = one + exp(temp0) - sensord = sepsensorsharpness*2*exp(temp0)*one*sensord/temp**2 + temp1 = -(2*sepsensorsharpness*(sensor-sepsensoroffset)) + temp0 = one + exp(temp1) + sensord = sepsensorsharpness*2*exp(temp1)*one*sensord/temp0**2 + call popreal8(sensor) vd(1) = vd(1) - veldirfreestream(1)*sensord veldirfreestreamd(1) = veldirfreestreamd(1) - v(1)*sensord vd(2) = vd(2) - veldirfreestream(2)*sensord veldirfreestreamd(2) = veldirfreestreamd(2) - v(2)*sensord vd(3) = vd(3) - veldirfreestream(3)*sensord veldirfreestreamd(3) = veldirfreestreamd(3) - v(3)*sensord + call popcontrol1b(branch) + if (branch .eq. 0) then + ks_exponentd = blk*sepsensorksd + temp0 = -(2*sepsensorsharpnesstwo*(sepsensoroffsettwo+sensor)) + temp = one + exp(temp0) + tempd1 = blk*one*sepsensoraread/temp + cellaread = cellaread + tempd1 + blk*sepsensorareaksd + sensord = sepsensorsharpnesstwo*2*exp(temp0)*cellarea*tempd1/& +& temp + call ksaggregationfunction_b(sensor, sensord, sepsenmaxfamily(& +& spectralsol), sepsenmaxrho, ks_exponent, & +& ks_exponentd) + sensord = -(sensord/(cos(zero)-cos(degtorad*sepangledeviation)+& +& 1e-16)) + vd(1) = vd(1) + vecttangential(1)*sensord + vecttangentiald(1) = vecttangentiald(1) + v(1)*sensord + vd(2) = vd(2) + vecttangential(2)*sensord + vecttangentiald(2) = vecttangentiald(2) + v(2)*sensord + vd(3) = vd(3) + vecttangential(3)*sensord + vecttangentiald(3) = vecttangentiald(3) + v(3)*sensord + call popreal8array(vecttangential, 3) + tmpd1 = vecttangentiald + temp0 = vecttangential(1)*vecttangential(1) + vecttangential(2)*& +& vecttangential(2) + vecttangential(3)*vecttangential(3) + temp = sqrt(temp0) + vecttangentiald = tmpd1/(temp+1e-16) + if (temp0 .eq. 0.0_8) then + tempd0 = 0.0_8 + else + tempd0 = -(sum(vecttangential*tmpd1)/(2.0*temp*(temp+1e-16)**2& +& )) + end if + vecttangentiald(1) = vecttangentiald(1) + 2*vecttangential(1)*& +& tempd0 + vecttangentiald(2) = vecttangentiald(2) + 2*vecttangential(2)*& +& tempd0 + vecttangentiald(3) = vecttangentiald(3) + 2*vecttangential(3)*& +& tempd0 + vectdotproductfsnormald = -(bcdata(mm)%norm(i, j, 3)*& +& vecttangentiald(3)) - bcdata(mm)%norm(i, j, 2)*vecttangentiald& +& (2) - bcdata(mm)%norm(i, j, 1)*vecttangentiald(1) + veldirfreestreamd(3) = veldirfreestreamd(3) + vecttangentiald(3)& +& + bcdata(mm)%norm(i, j, 3)*vectdotproductfsnormald + vecttangentiald(3) = 0.0_8 + veldirfreestreamd(2) = veldirfreestreamd(2) + vecttangentiald(2)& +& + bcdata(mm)%norm(i, j, 2)*vectdotproductfsnormald + vecttangentiald(2) = 0.0_8 + veldirfreestreamd(1) = veldirfreestreamd(1) + vecttangentiald(1)& +& + bcdata(mm)%norm(i, j, 1)*vectdotproductfsnormald + vecttangentiald(1) = 0.0_8 + end if + mxd = blk*mpd(1) + myd = blk*mpd(2) + mzd = blk*mpd(3) call popreal8array(v, 3) tmpd0 = vd temp = v(1)*v(1) + v(2)*v(2) + v(3)*v(3) @@ -2024,6 +2206,7 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) vd(2) = 0.0_8 ww2d(i, j, ivx) = ww2d(i, j, ivx) + vd(1) vd(1) = 0.0_8 + cellaread = cellaread + bcdatad(mm)%area(i, j) bcdatad(mm)%area(i, j) = 0.0_8 if (ssi(i, j, 1)**2 + ssi(i, j, 2)**2 + ssi(i, j, 3)**2 .eq. 0.0_8& & ) then @@ -2072,6 +2255,12 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm) xxd(i+1, j, 1) = xxd(i+1, j, 1) + tempd xxd(i, j+1, 1) = xxd(i, j+1, 1) + tempd xxd(i+1, j+1, 1) = xxd(i+1, j+1, 1) + tempd + zcod = zcod + fz*blk*cofsumfzd(3) + fy*blk*cofsumfyd(3) + fx*blk*& +& cofsumfxd(3) + ycod = ycod + fz*blk*cofsumfzd(2) + fy*blk*cofsumfyd(2) + fx*blk*& +& cofsumfxd(2) + xcod = xcod + fz*blk*cofsumfzd(1) + fy*blk*cofsumfyd(1) + fx*blk*& +& cofsumfxd(1) tempd = fourth*zcod xxd(i, j, 3) = xxd(i, j, 3) + tempd xxd(i+1, j, 3) = xxd(i+1, j, 3) + tempd @@ -2153,7 +2342,8 @@ subroutine wallintegrationface(localvalues, mm) use flowvarrefstate use inputcostfunctions use inputphysics, only : machcoef, pointref, veldirfreestream, & -& equations, momentaxis, cpmin_family, cpmin_rho, cavitationnumber +& equations, momentaxis, cpmin_family, cpmin_rho, cavitationnumber, & +& sepsenmaxfamily, sepsenmaxrho use bcpointers_b implicit none ! input/output variables @@ -2163,10 +2353,13 @@ subroutine wallintegrationface(localvalues, mm) ! local variables. real(kind=realtype), dimension(3) :: fp, fv, mp, mv real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz - real(kind=realtype) :: yplusmax, sepsensor, sepsensoravg(3), & -& cavitation, cpmin_ks_sum + real(kind=realtype) :: yplusmax, sepsensorks, sepsensor, & +& sepsensoravg(3), sepsensorarea, cavitation, cpmin_ks_sum, & +& sepsensorareaks integer(kind=inttype) :: i, j, ii, blk real(kind=realtype) :: pm1, fx, fy, fz, fn + real(kind=realtype) :: vecttangential(3) + real(kind=realtype) :: vectdotproductfsnormal real(kind=realtype) :: xc, xco, yc, yco, zc, zco, qf(3), r(3), n(3)& & , l real(kind=realtype) :: fact, rho, mul, yplus, dwall @@ -2182,6 +2375,7 @@ subroutine wallintegrationface(localvalues, mm) intrinsic mod intrinsic max intrinsic sqrt + intrinsic cos intrinsic exp select case (bcfaceid(mm)) case (imin, jmin, kmin) @@ -2212,6 +2406,9 @@ subroutine wallintegrationface(localvalues, mm) cofsumfz = zero yplusmax = zero sepsensor = zero + sepsensorks = zero + sepsensorarea = zero + sepsensorareaks = zero cavitation = zero cpmin_ks_sum = zero sepsensoravg = zero @@ -2337,6 +2534,36 @@ subroutine wallintegrationface(localvalues, mm) v(2) = ww2(i, j, ivy) v(3) = ww2(i, j, ivz) v = v/(sqrt(v(1)**2+v(2)**2+v(3)**2)+1e-16) + if (computesepsensorks) then +! freestream projection over the surface. + vectdotproductfsnormal = veldirfreestream(1)*bcdata(mm)%norm(i, & +& j, 1) + veldirfreestream(2)*bcdata(mm)%norm(i, j, 2) + & +& veldirfreestream(3)*bcdata(mm)%norm(i, j, 3) +! tangential vector on the surface, which is the freestream projected vector + vecttangential(1) = veldirfreestream(1) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 1) + vecttangential(2) = veldirfreestream(2) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 2) + vecttangential(3) = veldirfreestream(3) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 3) + vecttangential = vecttangential/(sqrt(vecttangential(1)**2+& +& vecttangential(2)**2+vecttangential(3)**2)+1e-16) +! computing separation sensor +! velocity dot products + sensor = v(1)*vecttangential(1) + v(2)*vecttangential(2) + v(3)*& +& vecttangential(3) +! sepsensor value + sensor = (cos(degtorad*sepangledeviation)-sensor)/(-cos(degtorad& +& *sepangledeviation)+cos(zero)+1e-16) +! also do the ks-based spensenor max computation + call ksaggregationfunction(sensor, sepsenmaxfamily(spectralsol)& +& , sepsenmaxrho, ks_exponent) + sepsensorareaks = sepsensorareaks + blk*cellarea + sepsensorarea = cellarea*blk*one/(one+exp(-(2*& +& sepsensorsharpnesstwo*(sensor+sepsensoroffsettwo)))) + & +& sepsensorarea + sepsensorks = sepsensorks + ks_exponent*blk + end if ! dot product with free stream sensor = -(v(1)*veldirfreestream(1)+v(2)*veldirfreestream(2)+v(3)*& & veldirfreestream(3)) @@ -2507,6 +2734,11 @@ subroutine wallintegrationface(localvalues, mm) localvalues(icoforcez:icoforcez+2) = localvalues(icoforcez:icoforcez& & +2) + cofsumfz localvalues(isepsensor) = localvalues(isepsensor) + sepsensor + localvalues(isepsensorks) = localvalues(isepsensorks) + sepsensorks + localvalues(isepsensorareaks) = localvalues(isepsensorareaks) + & +& sepsensorareaks + localvalues(isepsensorarea) = localvalues(isepsensorarea) + & +& sepsensorarea localvalues(icavitation) = localvalues(icavitation) + cavitation localvalues(icpmin) = localvalues(icpmin) + cpmin_ks_sum localvalues(isepavg:isepavg+2) = localvalues(isepavg:isepavg+2) + & @@ -2516,6 +2748,29 @@ subroutine wallintegrationface(localvalues, mm) localvalues(icperror2) = localvalues(icperror2) + cperror2 end subroutine wallintegrationface +! differentiation of ksaggregationfunction in reverse (adjoint) mode (with options noisize i4 dr8 r8): +! gradient of useful results: ks_g g +! with respect to varying inputs: g + subroutine ksaggregationfunction_b(g, gd, max_g, g_rho, ks_g, ks_gd) + use precision + implicit none + real(kind=realtype), intent(inout) :: ks_g + real(kind=realtype), intent(inout) :: ks_gd + real(kind=realtype), intent(in) :: g, max_g, g_rho + real(kind=realtype) :: gd + intrinsic exp + gd = gd + g_rho*exp(g_rho*(g-max_g))*ks_gd + end subroutine ksaggregationfunction_b + + subroutine ksaggregationfunction(g, max_g, g_rho, ks_g) + use precision + implicit none + real(kind=realtype), intent(inout) :: ks_g + real(kind=realtype), intent(in) :: g, max_g, g_rho + intrinsic exp + ks_g = exp(g_rho*(g-max_g)) + end subroutine ksaggregationfunction + ! differentiation of flowintegrationface in reverse (adjoint) mode (with options noisize i4 dr8 r8): ! gradient of useful results: timeref tref pref rgas rhoref ! *xx *pp1 *pp2 *ssi *ww1 *ww2 pointref localvalues diff --git a/src/adjoint/outputReverseFast/surfaceIntegrations_fast_b.f90 b/src/adjoint/outputReverseFast/surfaceIntegrations_fast_b.f90 index b5274cb37..4b28bd072 100644 --- a/src/adjoint/outputReverseFast/surfaceIntegrations_fast_b.f90 +++ b/src/adjoint/outputReverseFast/surfaceIntegrations_fast_b.f90 @@ -11,8 +11,10 @@ subroutine getcostfunctions(globalvals, funcvalues) use flowvarrefstate, only : pref, rhoref, tref, lref, gammainf, pinf& & , uref, uinf use inputphysics, only : liftdirection, dragdirection, surfaceref, & -& machcoef, lengthref, alpha, beta, liftindex, cpmin_family, cpmin_rho - use inputcostfunctions, only : computecavitation +& machcoef, lengthref, alpha, beta, liftindex, cpmin_family, cpmin_rho& +& , sepsenmaxfamily, sepsenmaxrho + use inputcostfunctions, only : computecavitation, computesepsensorks& +& , sepsensorsharpnessone, sepsensoroffsetone use inputtsstabderiv, only : tsstability use utils_fast_b, only : computetsderivatives use flowutils_fast_b, only : getdirvector @@ -29,12 +31,13 @@ subroutine getcostfunctions(globalvals, funcvalues) real(kind=realtype) :: mavgptot, mavgttot, mavgrho, mavgps, mflow, & & mavgmn, mavga, mavgvx, mavgvy, mavgvz, garea, mavgvi, fxlift, fylift& & , fzlift - real(kind=realtype) :: vdotn, mag, u, v, w + real(kind=realtype) :: vdotn, mag, u, v, w, ks_comp integer(kind=inttype) :: sps real(kind=realtype), dimension(8) :: dcdq, dcdqdot real(kind=realtype), dimension(8) :: dcdalpha, dcdalphadot real(kind=realtype), dimension(8) :: coef0 intrinsic log + intrinsic exp intrinsic sqrt ! factor used for time-averaged quantities. ovrnts = one/ntimeintervalsspectral @@ -164,6 +167,18 @@ subroutine getcostfunctions(globalvals, funcvalues) & ovrnts*cmoment(2, sps) funcvalues(costfuncmomzcoef) = funcvalues(costfuncmomzcoef) + & & ovrnts*cmoment(3, sps) +! final part of the ks computation + if (computesepsensorks) then +! only calculate the log part if we are actually computing for separation for ks method. + ks_comp = ovrnts*(sepsenmaxfamily(sps)+log(globalvals(& +& isepsensorks, sps))/sepsenmaxrho) + funcvalues(costfuncsepsensorks) = funcvalues(costfuncsepsensorks& +& ) + ks_comp + funcvalues(costfuncsepsensorarea) = funcvalues(& +& costfuncsepsensorarea) + ovrnts*globalvals(isepsensorareaks, & +& sps)*ks_comp*one/(one+exp(2*sepsensorsharpnessone*(ks_comp+& +& sepsensoroffsetone))) + ovrnts*globalvals(isepsensorarea, sps) + end if funcvalues(costfuncsepsensor) = funcvalues(costfuncsepsensor) + & & ovrnts*globalvals(isepsensor, sps) funcvalues(costfunccavitation) = funcvalues(costfunccavitation) + & @@ -361,7 +376,8 @@ subroutine wallintegrationface(localvalues, mm) use flowvarrefstate use inputcostfunctions use inputphysics, only : machcoef, pointref, veldirfreestream, & -& equations, momentaxis, cpmin_family, cpmin_rho, cavitationnumber +& equations, momentaxis, cpmin_family, cpmin_rho, cavitationnumber, & +& sepsenmaxfamily, sepsenmaxrho use bcpointers implicit none ! input/output variables @@ -371,10 +387,13 @@ subroutine wallintegrationface(localvalues, mm) ! local variables. real(kind=realtype), dimension(3) :: fp, fv, mp, mv real(kind=realtype), dimension(3) :: cofsumfx, cofsumfy, cofsumfz - real(kind=realtype) :: yplusmax, sepsensor, sepsensoravg(3), & -& cavitation, cpmin_ks_sum + real(kind=realtype) :: yplusmax, sepsensorks, sepsensor, & +& sepsensoravg(3), sepsensorarea, cavitation, cpmin_ks_sum, & +& sepsensorareaks integer(kind=inttype) :: i, j, ii, blk real(kind=realtype) :: pm1, fx, fy, fz, fn + real(kind=realtype) :: vecttangential(3) + real(kind=realtype) :: vectdotproductfsnormal real(kind=realtype) :: xc, xco, yc, yco, zc, zco, qf(3), r(3), n(3)& & , l real(kind=realtype) :: fact, rho, mul, yplus, dwall @@ -390,6 +409,7 @@ subroutine wallintegrationface(localvalues, mm) intrinsic mod intrinsic max intrinsic sqrt + intrinsic cos intrinsic exp select case (bcfaceid(mm)) case (imin, jmin, kmin) @@ -420,6 +440,9 @@ subroutine wallintegrationface(localvalues, mm) cofsumfz = zero yplusmax = zero sepsensor = zero + sepsensorks = zero + sepsensorarea = zero + sepsensorareaks = zero cavitation = zero cpmin_ks_sum = zero sepsensoravg = zero @@ -545,6 +568,36 @@ subroutine wallintegrationface(localvalues, mm) v(2) = ww2(i, j, ivy) v(3) = ww2(i, j, ivz) v = v/(sqrt(v(1)**2+v(2)**2+v(3)**2)+1e-16) + if (computesepsensorks) then +! freestream projection over the surface. + vectdotproductfsnormal = veldirfreestream(1)*bcdata(mm)%norm(i, & +& j, 1) + veldirfreestream(2)*bcdata(mm)%norm(i, j, 2) + & +& veldirfreestream(3)*bcdata(mm)%norm(i, j, 3) +! tangential vector on the surface, which is the freestream projected vector + vecttangential(1) = veldirfreestream(1) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 1) + vecttangential(2) = veldirfreestream(2) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 2) + vecttangential(3) = veldirfreestream(3) - vectdotproductfsnormal& +& *bcdata(mm)%norm(i, j, 3) + vecttangential = vecttangential/(sqrt(vecttangential(1)**2+& +& vecttangential(2)**2+vecttangential(3)**2)+1e-16) +! computing separation sensor +! velocity dot products + sensor = v(1)*vecttangential(1) + v(2)*vecttangential(2) + v(3)*& +& vecttangential(3) +! sepsensor value + sensor = (cos(degtorad*sepangledeviation)-sensor)/(-cos(degtorad& +& *sepangledeviation)+cos(zero)+1e-16) +! also do the ks-based spensenor max computation + call ksaggregationfunction(sensor, sepsenmaxfamily(spectralsol)& +& , sepsenmaxrho, ks_exponent) + sepsensorareaks = sepsensorareaks + blk*cellarea + sepsensorarea = cellarea*blk*one/(one+exp(-(2*& +& sepsensorsharpnesstwo*(sensor+sepsensoroffsettwo)))) + & +& sepsensorarea + sepsensorks = sepsensorks + ks_exponent*blk + end if ! dot product with free stream sensor = -(v(1)*veldirfreestream(1)+v(2)*veldirfreestream(2)+v(3)*& & veldirfreestream(3)) @@ -715,6 +768,11 @@ subroutine wallintegrationface(localvalues, mm) localvalues(icoforcez:icoforcez+2) = localvalues(icoforcez:icoforcez& & +2) + cofsumfz localvalues(isepsensor) = localvalues(isepsensor) + sepsensor + localvalues(isepsensorks) = localvalues(isepsensorks) + sepsensorks + localvalues(isepsensorareaks) = localvalues(isepsensorareaks) + & +& sepsensorareaks + localvalues(isepsensorarea) = localvalues(isepsensorarea) + & +& sepsensorarea localvalues(icavitation) = localvalues(icavitation) + cavitation localvalues(icpmin) = localvalues(icpmin) + cpmin_ks_sum localvalues(isepavg:isepavg+2) = localvalues(isepavg:isepavg+2) + & @@ -724,6 +782,15 @@ subroutine wallintegrationface(localvalues, mm) localvalues(icperror2) = localvalues(icperror2) + cperror2 end subroutine wallintegrationface + subroutine ksaggregationfunction(g, max_g, g_rho, ks_g) + use precision + implicit none + real(kind=realtype), intent(inout) :: ks_g + real(kind=realtype), intent(in) :: g, max_g, g_rho + intrinsic exp + ks_g = exp(g_rho*(g-max_g)) + end subroutine ksaggregationfunction + subroutine flowintegrationface(isinflow, localvalues, mm) use constants use blockpointers, only : bctype, bcfaceid, bcdata, & diff --git a/src/f2py/adflow.pyf b/src/f2py/adflow.pyf index c666338fc..000596b84 100644 --- a/src/f2py/adflow.pyf +++ b/src/f2py/adflow.pyf @@ -185,6 +185,8 @@ python module libadflow integer(kind=inttype) parameter,optional :: costfunccoflifty=100 integer(kind=inttype) parameter,optional :: costfunccofliftz=101 integer(kind=inttype) parameter,optional :: costfuncmavgvi=102 + integer(kind=inttype) parameter,optional :: costfuncsepsensorks=103 + integer(kind=inttype) parameter,optional :: costfuncsepsensorarea=104 end module constants module utils @@ -1196,6 +1198,7 @@ python module libadflow integer(kind=inttype) :: liftindex real(kind=realtype) :: cavitationnumber real(kind=realtype) :: cpmin_rho + real(kind=realtype) :: sepsenmaxrho end module inputphysics module inputadjoint ! in :adflow:../modules/inputParam.f90 @@ -1377,8 +1380,14 @@ python module libadflow module inputcostfunctions ! in :test:costFunctions.F90 use constants + logical::computesepsensorks real(kind=realtype) sepsensoroffset + real(kind=realtype) sepsensoroffsetone + real(kind=realtype) sepsensoroffsettwo real(kind=realtype) sepsensorsharpness + real(kind=realtype) sepsensorsharpnessone + real(kind=realtype) sepsensorsharpnesstwo + real(kind=realtype) sepangledeviation real(kind=realtype) cavsensoroffset real(kind=realtype) cavsensorsharpness integer(kind=inttype) cavexponent @@ -1508,4 +1517,4 @@ python module libadflow end python module libadflow_cs #else end python module libadflow -#endif +#endif \ No newline at end of file diff --git a/src/inputParam/inputParamRoutines.F90 b/src/inputParam/inputParamRoutines.F90 index fc1eb8101..e254a7cf6 100644 --- a/src/inputParam/inputParamRoutines.F90 +++ b/src/inputParam/inputParamRoutines.F90 @@ -298,6 +298,9 @@ subroutine checkMonitor case (cgnsAxisMoment) sortNumber(i) = 116 + case (cgnsSepSensorArea) + sortNumber(i) = 117 + case (cgnsHdiffMax) sortNumber(i) = 201 @@ -1594,6 +1597,10 @@ subroutine monitorVariables(variables) nMon = nMon + 1; nMonSum = nMonSum + 1 tmpNames(nMon) = cgnsSepSensor + case ("sepsensorarea") + nMon = nMon + 1; nMonSum = nMonSum + 1 + tmpNames(nMon) = cgnsSepSensorArea + case ("cavitation") nMon = nMon + 1; nMonSum = nMonSum + 1 tmpNames(nMon) = cgnsCavitation @@ -2352,6 +2359,8 @@ subroutine surfaceVariables(variables) surfWriteBlank = .false. surfWriteSepSensor = .false. + surfWriteSepSensorKs = .false. + surfWriteSepSensorArea = .false. surfWriteCavitation = .false. surfWriteAxisMoment = .false. surfWriteGC = .false. @@ -2470,6 +2479,14 @@ subroutine surfaceVariables(variables) surfWriteSepSensor = .true. nVarSpecified = nVarSpecified + 1 + case ("sepsensorks") + surfWriteSepSensorKs = .true. + nVarSpecified = nVarSpecified + 1 + + case ("sepsensorarea") + surfWriteSepSensorArea = .true. + nVarSpecified = nVarSpecified + 1 + case ("cavitation") surfWriteCavitation = .true. nVarSpecified = nVarSpecified + 1 @@ -3697,6 +3714,17 @@ subroutine checkInputParam cpmin_family = zero end if + ! Allocate the memory for sepsenmaxfamily. We had to wait until + ! nTimeIntervalsSpectral was set. + if (.not. allocated(sepSenMaxFamily)) then + allocate (sepSenMaxFamily(nTimeIntervalsSpectral), stat=ierr) + if (ierr /= 0) & + call terminate("checkInputParam", & + "Memory allocation failure for & + &sepSenMaxFamily") + sepSenMaxFamily = zero + end if + end subroutine checkInputParam subroutine setDefaultValues ! @@ -4046,8 +4074,6 @@ subroutine setDefaultValues adjointPETScVarsAllocated = .False. adjointPETScPreProcVarsAllocated = .False. usematrixfreedrdw = .False. - sepSensorOffset = zero - sepSensorSharpness = 10_realType end subroutine setDefaultValues subroutine initializeIsoSurfaceVariables(values, nValues) diff --git a/src/modules/cgnsNames.f90 b/src/modules/cgnsNames.f90 index c3c67549e..4e4f0743c 100644 --- a/src/modules/cgnsNames.f90 +++ b/src/modules/cgnsNames.f90 @@ -272,6 +272,8 @@ module cgnsNames ! Names of the "lift" force, separation sensor and cavitation ! character(len=maxCGNSNameLen), parameter :: cgnsSepSensor = "SepSensor" + character(len=maxCGNSNameLen), parameter :: cgnsSepSensorKs = "SepSensorKs" + character(len=maxCGNSNameLen), parameter :: cgnsSepSensorArea = "SepSensorArea" character(len=maxCGNSNameLen), parameter :: cgnsCavitation = "Cavitation" character(len=maxCGNSNameLen), parameter :: cgnsAxisMoment = "AxisMoment" ! diff --git a/src/modules/constants.F90 b/src/modules/constants.F90 index 086240822..bfd6af353 100644 --- a/src/modules/constants.F90 +++ b/src/modules/constants.F90 @@ -85,6 +85,9 @@ module constants real(kind=realType), parameter :: threefourth = 0.75_realType real(kind=realType), parameter :: sqrtthree = 1.7320508075688772_realType + ! radian to degree conversion + real(kind=realType), parameter :: degtorad = pi / 180.0_realType + ! String constants CHARACTER(*), PARAMETER :: LOWER_CASE = 'abcdefghijklmnopqrstuvwxyz' CHARACTER(*), PARAMETER :: UPPER_CASE = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' @@ -341,7 +344,7 @@ module constants integer(kind=intType), parameter :: iTotal = 16 ! Cost functions. - integer(kind=intType), parameter :: nCostFunction = 102 + integer(kind=intType), parameter :: nCostFunction = 104 integer(kind=intType), parameter :: & costFuncLift = 1, & costFuncDrag = 2, & @@ -444,9 +447,11 @@ module constants costfuncCofLiftX = 99, & costfuncCofLiftY = 100, & costfuncCofLiftZ = 101, & - costfuncmavgvi = 102 + costfuncmavgvi = 102, & + costFuncSepSensorKs = 103, & + costFuncSepSensorArea = 104 - integer(kind=intType), parameter :: nLocalValues = 60 + integer(kind=intType), parameter :: nLocalValues = 63 integer(kind=intType), parameter :: & iFp = 1, & iFv = 4, & @@ -486,7 +491,10 @@ module constants iCoForceX = 51, & iCoForceY = 54, & iCoForceZ = 57, & - iMassVi = 60 + iMassVi = 60, & + iSepSensorKs = 61, & + iSepSensorArea = 62, & + iSepSensorAreaKs = 63 ! Constants for zipper comm diff --git a/src/modules/diffSizes.f90 b/src/modules/diffSizes.f90 index b765c08de..c3634e8b3 100644 --- a/src/modules/diffSizes.f90 +++ b/src/modules/diffSizes.f90 @@ -75,6 +75,8 @@ module diffSizes integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_norm, ISIZE2OFDrfDrfbcdata_norm, ISIZE3OFDrfDrfbcdata_norm integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_oarea, ISIZE2OFDrfDrfbcdata_oarea, ISIZE3OFDrfDrfbcdata_oarea integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_sepsensor, ISIZE2OFDrfDrfbcdata_sepSensor + integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_sepsensorks, ISIZE2OFDrfDrfbcdata_sepSensorKs + integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_sepsensorarea, ISIZE2OFDrfDrfbcdata_sepSensorArea integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_Cavitation, ISIZE2OFDrfDrfbcdata_Cavitation integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_axisMoment, ISIZE2OFDrfDrfbcdata_axisMoment integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_rface, ISIZE2OFDrfDrfbcdata_rface, ISIZE3OFDrfDrfbcdata_rface diff --git a/src/modules/extraOutput.f90 b/src/modules/extraOutput.f90 index 8a6e33e1e..bf75d99e1 100644 --- a/src/modules/extraOutput.f90 +++ b/src/modules/extraOutput.f90 @@ -19,7 +19,7 @@ module extraOutput logical :: surfWriteRMach logical :: surfWriteCf, surfWriteCh, surfWriteYPlus logical :: surfWriteCfx, surfWriteCfy, surfWriteCfz - logical :: surfWriteBlank, surfWriteSepSensor + logical :: surfWriteBlank, surfWriteSepSensor, surfWriteSepSensorKs, surfWriteSepSensorArea logical :: surfWriteCavitation, surfWriteGC, surfWriteAxisMoment ! ! The logical variables, which define the extra volume variables diff --git a/src/modules/inputParam.F90 b/src/modules/inputParam.F90 index 2d994711c..2dca9d9fd 100644 --- a/src/modules/inputParam.F90 +++ b/src/modules/inputParam.F90 @@ -299,8 +299,14 @@ end module inputIteration module inputCostFunctions use constants - real(kind=realtype) :: sepSensorOffset = zero - real(kind=realtype) :: sepSensorSharpness = 10.0_realType + logical :: computeSepSensorKs + real(kind=realtype) :: sepSensorOffset + real(kind=realtype) :: sepSensorOffsetOne + real(kind=realtype) :: sepSensorOffsetTwo + real(kind=realtype) :: sepAngleDeviation + real(kind=realtype) :: sepSensorSharpness + real(kind=realtype) :: sepSensorSharpnessOne + real(kind=realtype) :: sepSensorSharpnessTwo real(kind=realtype) :: cavSensorOffset real(kind=realtype) :: cavSensorSharpness integer(kind=inttype) :: cavExponent @@ -571,6 +577,9 @@ module inputPhysics ! cpmin_rho The rho parameter used with the KS-based cavitation sensor. ! cpmin_family The cpmin for a given surface family that does not use ! KS-aggregation, but rather an exact min computation. + ! sepSenMaxRho The rho parameter used with the KS-based separation sensor. + ! sepSenMaxFamily The maximum sepsensor value for a given surface family that does not use + ! KS-aggregation, but rather an exact max computation. integer(kind=intType) :: equations, equationMode, flowType integer(kind=intType) :: turbModel, cpModel, turbProd @@ -597,6 +606,8 @@ module inputPhysics real(kind=realType) :: cavitationnumber real(kind=realType) :: cpmin_rho real(kind=realType), dimension(:), allocatable :: cpmin_family + real(kind=realType) :: sepSenMaxRho + real(kind=realType), dimension(:), allocatable :: sepSenMaxFamily #ifndef USE_TAPENADE real(kind=realType) :: alphad, betad diff --git a/src/output/outputMod.F90 b/src/output/outputMod.F90 index 1a954888b..c6ff293a6 100644 --- a/src/output/outputMod.F90 +++ b/src/output/outputMod.F90 @@ -165,8 +165,10 @@ subroutine numberOfSurfSolVariables(nSolVar) if (surfWriteCfz) nSolVar = nSolVar + 1 if (surfWriteBlank) nSolVar = nSolVar + 1 if (surfWriteSepSensor) nSolVar = nSolVar + 1 - if (surfWriteCavitation) nsolVar = nsolVar + 1 - if (surfWriteGC) nsolVar = nsolVar + 1 + if (surfWriteSepSensorKs) nSolVar = nSolVar + 1 + if (surfWriteSepSensorArea) nSolVar = nSolVar + 1 + if (surfWriteCavitation) nSolVar = nSolVar + 1 + if (surfWriteGC) nSolVar = nSolVar + 1 end subroutine numberOfSurfSolVariables @@ -699,6 +701,16 @@ subroutine surfSolNames(solNames) solNames(nn) = cgnsSepSensor end if + if (surfWriteSepSensorKs) then + nn = nn + 1 + solNames(nn) = cgnsSepSensorKs + end if + + if (surfWriteSepSensorArea) then + nn = nn + 1 + solNames(nn) = cgnsSepSensorArea + end if + if (surfWriteCavitation) then nn = nn + 1 solNames(nn) = cgnsCavitation @@ -1424,6 +1436,8 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, & real(kind=realType) :: tauxx, tauyy, tauzz real(kind=realType) :: tauxy, tauxz, tauyz real(kind=realType) :: pm1, a, sensor, plocal, sensor1 + real(kind=realType) :: vectTangential(3) + real(kind=realType) :: vectDotProductFsNormal real(kind=realType), dimension(3) :: norm, V real(kind=realType), dimension(:, :, :), pointer :: ww1, ww2 @@ -1471,8 +1485,8 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, & select case (solName) - case (cgnsSkinFmag, cgnsStanton, cgnsYplus, & - cgnsSkinFx, cgnsSkinFy, cgnsSkinFz) + case (cgnsSkinFmag, cgnsStanton, cgnsYplus, cgnsSkinFx, & + cgnsSkinFy, cgnsSkinFz, cgnsSepSensor, cgnsSepSensorKs, cgnsSepSensorArea) ! Update the counter and set this entry of buffer to 0. @@ -2165,12 +2179,154 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, & end do end do - case (cgnsSepSensor) + case (cgnsSepSensorKs) do j = rangeFace(2, 1), rangeFace(2, 2) + if (present(jBeg) .and. present(jEnd) .and. (useRindLayer)) then + jor = j + jBegOr - 1 + if (jor == jBeg) then + jj = j + 1 + else if (jor == jEnd + 1) then + jj = j - 1 + else + jj = j + end if + else + jj = j + + end if + do i = rangeFace(1, 1), rangeFace(1, 2) + + if (present(iBeg) .and. present(iEnd) .and. (useRindLayer)) then + ior = i + iBegor - 1 + if (ior == iBeg) then + ii = i + 1 + else if (ior == iEnd + 1) then + ii = i - 1 + else + ii = i + end if + else + ii = i + end if + nn = nn + 1 + ! Get normalized surface velocity: + v(1) = ww2(ii, jj, ivx) + v(2) = ww2(ii, jj, ivy) + v(3) = ww2(ii, jj, ivz) + + ! Normalize + v = v / (sqrt(v(1)**2 + v(2)**2 + v(3)**2) + 1e-16) + mm = viscPointer(ii, jj) + + norm(1) = BCData(mm)%norm(ii, jj, 1) + norm(2) = BCData(mm)%norm(ii, jj, 2) + norm(3) = BCData(mm)%norm(ii, jj, 3) + + vectDotProductFsNormal = velDirFreeStream(1) * norm(1) + & + velDirFreeStream(2) * norm(2) + & + velDirFreeStream(3) * norm(3) + + vectTangential(1) = velDirFreeStream(1) - vectDotProductFsNormal * norm(1) + vectTangential(2) = velDirFreeStream(2) - vectDotProductFsNormal * norm(2) + vectTangential(3) = velDirFreeStream(3) - vectDotProductFsNormal * norm(3) + + vectTangential = vectTangential / (sqrt(vectTangential(1)**2 + vectTangential(2)**2 + & + vectTangential(3)**2) + 1e-16) + + ! computing separation sensor + ! velocity dot products + sensor = (v(1) * vectTangential(1) + v(2) * vectTangential(2) + & + v(3) * vectTangential(3)) + + ! sepsensor value + sensor = (cos(degtorad * sepAngleDeviation) - sensor) / & + (-cos(degtorad * sepAngleDeviation) + cos(zero) + 1e-16) + + buffer(nn) = sensor + end do + end do + + case (cgnsSepSensorArea) + + do j = rangeFace(2, 1), rangeFace(2, 2) + if (present(jBeg) .and. present(jEnd) .and. (useRindLayer)) then + jor = j + jBegOr - 1 + if (jor == jBeg) then + jj = j + 1 + else if (jor == jEnd + 1) then + jj = j - 1 + else + jj = j + end if + else + jj = j + + end if + + do i = rangeFace(1, 1), rangeFace(1, 2) + + if (present(iBeg) .and. present(iEnd) .and. (useRindLayer)) then + ior = i + iBegor - 1 + if (ior == iBeg) then + ii = i + 1 + else if (ior == iEnd + 1) then + ii = i - 1 + else + ii = i + end if + else + ii = i + end if + + nn = nn + 1 + + ! Get normalized surface velocity: + v(1) = ww2(ii, jj, ivx) + v(2) = ww2(ii, jj, ivy) + v(3) = ww2(ii, jj, ivz) + + ! Normalize + v = v / (sqrt(v(1)**2 + v(2)**2 + v(3)**2) + 1e-16) + mm = viscPointer(ii, jj) + + norm(1) = BCData(mm)%norm(ii, jj, 1) + norm(2) = BCData(mm)%norm(ii, jj, 2) + norm(3) = BCData(mm)%norm(ii, jj, 3) + + vectDotProductFsNormal = velDirFreeStream(1) * norm(1) + & + velDirFreeStream(2) * norm(2) + & + velDirFreeStream(3) * norm(3) + + vectTangential(1) = velDirFreeStream(1) - vectDotProductFsNormal * norm(1) + vectTangential(2) = velDirFreeStream(2) - vectDotProductFsNormal * norm(2) + vectTangential(3) = velDirFreeStream(3) - vectDotProductFsNormal * norm(3) + + vectTangential = vectTangential / (sqrt(vectTangential(1)**2 + vectTangential(2)**2 + & + vectTangential(3)**2) + 1e-16) + + ! computing separation sensor + ! velocity dot products + sensor = (v(1) * vectTangential(1) + v(2) * vectTangential(2) + & + v(3) * vectTangential(3)) + + ! sepsensor value + sensor = (cos(degtorad * sepAngleDeviation) - sensor) / & + (-cos(degtorad * sepAngleDeviation) + cos(zero) + 1e-16) + + sensor = one / (one + exp(-2 * sepSensorSharpnessTwo * (sensor + sepSensorOffsetTwo))) + + buffer(nn) = sensor + end do + end do + + case (cgnsSepSensor) + do j = rangeFace(2, 1), rangeFace(2, 2) + do i = rangeFace(1, 1), rangeFace(1, 2) + nn = nn + 1 ! Get normalized surface velocity: v(1) = ww2(i, j, ivx) v(2) = ww2(i, j, ivy) diff --git a/src/solver/solvers.F90 b/src/solver/solvers.F90 index 40d254bbb..cb553148a 100644 --- a/src/solver/solvers.F90 +++ b/src/solver/solvers.F90 @@ -1440,6 +1440,9 @@ subroutine convergenceInfo case (cgnsSepSensor) monLoc(mm) = monLoc(mm) + localValues(isepSensor) + case (cgnsSepSensorArea) + monLoc(mm) = monLoc(mm) + localValues(isepSensorArea) + case (cgnsCavitation) monLoc(mm) = monLoc(mm) + localValues(iCavitation) diff --git a/src/solver/surfaceIntegrations.F90 b/src/solver/surfaceIntegrations.F90 index f5bcd81f2..96cc144da 100644 --- a/src/solver/surfaceIntegrations.F90 +++ b/src/solver/surfaceIntegrations.F90 @@ -9,8 +9,9 @@ subroutine getCostFunctions(globalVals, funcValues) use flowVarRefState, only: pRef, rhoRef, tRef, LRef, gammaInf, pInf, uRef, uInf use inputPhysics, only: liftDirection, dragDirection, surfaceRef, & machCoef, lengthRef, alpha, beta, liftIndex, cpmin_family, & - cpmin_rho - use inputCostFunctions, only: computeCavitation + cpmin_rho, sepSenMaxFamily, sepSenMaxRho + use inputCostFunctions, only: computeCavitation, computeSepSensorKs, sepSensorSharpnessOne, & + sepSensorOffsetOne use inputTSStabDeriv, only: TSstability use utils, only: computeTSDerivatives use flowUtils, only: getDirVector @@ -29,7 +30,7 @@ subroutine getCostFunctions(globalVals, funcValues) real(kind=realType) :: mAvgPtot, mAvgTtot, mAvgRho, mAvgPs, mFlow, mAvgMn, mAvga, & mAvgVx, mAvgVy, mAvgVz, gArea, mAvgVi, fxLift, fyLift, fzLift - real(kind=realType) :: vdotn, mag, u, v, w + real(kind=realType) :: vdotn, mag, u, v, w, ks_comp integer(kind=intType) :: sps real(kind=realType), dimension(8) :: dcdq, dcdqdot real(kind=realType), dimension(8) :: dcdalpha, dcdalphadot @@ -149,7 +150,24 @@ subroutine getCostFunctions(globalVals, funcValues) funcValues(costFuncMomYCoef) = funcValues(costFuncMomYCoef) + ovrNTS * cMoment(2, sps) funcValues(costFuncMomZCoef) = funcValues(costFuncMomZCoef) + ovrNTS * cMoment(3, sps) + ! final part of the KS computation + if (computeSepSensorKs) then + ! only calculate the log part if we are actually computing for separation for KS method. + ks_comp = ovrNTS * (sepSenMaxFamily(sps) + & + log(globalVals(iSepSensorKs, sps)) / sepSenMaxRho) + + funcValues(costFuncSepSensorKs) = funcValues(costFuncSepSensorKs) + ks_comp + + funcValues(costFuncSepSensorArea) = funcValues(costFuncSepSensorArea) + & + ovrNTS * globalVals(iSepSensorAreaKs, sps) * ks_comp * & + one / (one + exp(2 * sepSensorSharpnessOne & + * (ks_comp + sepSensorOffsetOne))) + & + ovrNTS * globalVals(iSepSensorArea, sps) + + end if + funcValues(costFuncSepSensor) = funcValues(costFuncSepSensor) + ovrNTS * globalVals(iSepSensor, sps) + funcValues(costFuncCavitation) = funcValues(costFuncCavitation) + ovrNTS * globalVals(iCavitation, sps) ! final part of the KS computation if (computeCavitation) then @@ -402,7 +420,8 @@ subroutine wallIntegrationFace(localValues, mm) use flowVarRefState use inputCostFunctions use inputPhysics, only: MachCoef, pointRef, velDirFreeStream, & - equations, momentAxis, cpmin_family, cpmin_rho, cavitationnumber + equations, momentAxis, cpmin_family, cpmin_rho, & + cavitationnumber, sepSenMaxFamily, sepSenMaxRho use BCPointers implicit none @@ -413,10 +432,13 @@ subroutine wallIntegrationFace(localValues, mm) ! Local variables. real(kind=realType), dimension(3) :: Fp, Fv, Mp, Mv real(kind=realType), dimension(3) :: COFSumFx, COFSumFy, COFSumFz - real(kind=realType) :: yplusMax, sepSensor, sepSensorAvg(3), Cavitation, cpmin_ks_sum + real(kind=realType) :: yplusMax, sepSensorKs, sepSensor, sepSensorAvg(3), & + sepSensorArea, Cavitation, cpmin_ks_sum, sepSensorAreaKs integer(kind=intType) :: i, j, ii, blk real(kind=realType) :: pm1, fx, fy, fz, fn + real(kind=realType) :: vectTangential(3) + real(kind=realType) :: vectDotProductFsNormal real(kind=realType) :: xc, xco, yc, yco, zc, zco, qf(3), r(3), n(3), L real(kind=realType) :: fact, rho, mul, yplus, dwall real(kind=realType) :: V(3), sensor, sensor1, Cp, tmp, plocal, ks_exponent @@ -458,6 +480,9 @@ subroutine wallIntegrationFace(localValues, mm) COFSumFx = zero; COFSumFy = zero; COFSumFz = zero yplusMax = zero sepSensor = zero + sepSensorKs = zero + sepSensorArea = zero + sepSensorAreaKs = zero Cavitation = zero cpmin_ks_sum = zero sepSensorAvg = zero @@ -601,6 +626,39 @@ subroutine wallIntegrationFace(localValues, mm) v(3) = ww2(i, j, ivz) v = v / (sqrt(v(1)**2 + v(2)**2 + v(3)**2) + 1e-16) + if (computeSepSensorKs) then + ! Freestream projection over the surface. + vectDotProductFsNormal = velDirFreeStream(1) * BCData(mm)%norm(i, j, 1) + & + velDirFreeStream(2) * BCData(mm)%norm(i, j, 2) + & + velDirFreeStream(3) * BCData(mm)%norm(i, j, 3) + ! Tangential Vector on the surface, which is the freestream projected vector + vectTangential(1) = velDirFreeStream(1) - vectDotProductFsNormal * BCData(mm)%norm(i, j, 1) + vectTangential(2) = velDirFreeStream(2) - vectDotProductFsNormal * BCData(mm)%norm(i, j, 2) + vectTangential(3) = velDirFreeStream(3) - vectDotProductFsNormal * BCData(mm)%norm(i, j, 3) + + vectTangential = vectTangential / (sqrt(vectTangential(1)**2 + vectTangential(2)**2 + & + vectTangential(3)**2) + 1e-16) + + ! computing separation sensor + ! velocity dot products + sensor = (v(1) * vectTangential(1) + v(2) * vectTangential(2) + & + v(3) * vectTangential(3)) + + ! sepsensor value + sensor = (cos(degtorad * sepangledeviation) - sensor) / & + (-cos(degtorad * sepangledeviation) + cos(zero) + 1e-16) + + ! also do the ks-based spensenor max computation + call KSaggregationFunction(sensor, sepSenMaxFamily(spectralSol), sepSenMaxRho, ks_exponent) + + sepSensorAreaKs = sepSensorAreaKs + blk * cellArea + sepSensorArea = cellArea * blk * one / (one + exp(-2 * sepSensorSharpnessTwo & + * (sensor + sepSensorOffsetTwo))) + sepSensorArea + + sepSensorKs = sepSensorKs + ks_exponent * blk + + end if + ! Dot product with free stream sensor = -(v(1) * velDirFreeStream(1) + v(2) * velDirFreeStream(2) + & v(3) * velDirFreeStream(3)) @@ -808,6 +866,9 @@ subroutine wallIntegrationFace(localValues, mm) localValues(iCoForceY:iCoForceY + 2) = localValues(iCoForceY:iCoForceY + 2) + COFSumFy localValues(iCoForceZ:iCoForceZ + 2) = localValues(iCoForceZ:iCoForceZ + 2) + COFSumFz localValues(iSepSensor) = localValues(iSepSensor) + sepSensor + localValues(iSepSensorKs) = localValues(iSepSensorKs) + sepSensorKs + localValues(iSepSensorAreaKs) = localValues(iSepSensorAreaKs) + sepSensorAreaKs + localValues(iSepSensorArea) = localValues(iSepSensorArea) + sepSensorArea localValues(iCavitation) = localValues(iCavitation) + cavitation localValues(iCpMin) = localValues(iCpMin) + cpmin_ks_sum localValues(iSepAvg:iSepAvg + 2) = localValues(iSepAvg:iSepAvg + 2) + sepSensorAvg @@ -819,6 +880,17 @@ subroutine wallIntegrationFace(localValues, mm) #endif end subroutine wallIntegrationFace + subroutine KSaggregationFunction(g, max_g, g_rho, ks_g) + + use precision + + real(kind=realType), intent(inout) :: ks_g + real(kind=realType), intent(in) :: g, max_g, g_rho + + ks_g = exp(g_rho * (g - max_g)) + + end subroutine KSaggregationFunction + subroutine flowIntegrationFace(isInflow, localValues, mm) use constants @@ -1164,7 +1236,7 @@ subroutine getSolution(famLists, funcValues, globalValues) use zipperIntegrations, only: integrateZippers use userSurfaceIntegrations, only: integrateUserSurfaces use actuatorRegion, only: integrateActuatorRegions - use inputCostFunctions, only: computeCavitation + use inputCostFunctions, only: computeCavitation, computeSepSensorKs implicit none ! Input/Output Variables @@ -1191,6 +1263,11 @@ subroutine getSolution(famLists, funcValues, globalValues) call computeCpMinFamily(famList) end if + ! compute the current sepsensor max value for the separation computation with KS aggregation + if (computeSepSensorKs) then + call computeSepSenMaxFamily(famList) + end if + do sps = 1, nTimeIntervalsSpectral ! Integrate the normal block surfaces. do nn = 1, nDom @@ -1359,6 +1436,111 @@ subroutine computeCpMinFamily(famList) end subroutine computeCpMinFamily + subroutine computeSepSenMaxFamily(famList) + + use constants + use inputTimeSpectral, only: nTimeIntervalsSpectral + use communication, only: ADflow_comm_world, myID + use blockPointers, only: nDom + use inputPhysics, only: sepSenMaxFamily, velDirFreeStream + use blockPointers + use flowVarRefState + use inputCostFunctions, only: sepangledeviation + use BCPointers + use utils, only: setPointers, setBCPointers, isWallType, EChk + use sorting, only: famInList + + implicit none + + integer(kind=intType), dimension(:), intent(in) :: famList + integer(kind=intType) :: mm, nn, sps + integer(kind=intType) :: i, j, ii, blk, ierr + real(kind=realType) :: sepsensor_local + real(kind=realType) :: vectTangential(3), v(3) + real(kind=realType) :: vectDotProductFsNormal, sensor + + ! this routine loops over the surface cells in the given family + ! and computes the true maximum sepsensor value. + ! this is then used in the surface integration routine to compute + ! the maximum sepsensor value using KS aggregation. + ! the goal is to get a differentiable maximum sepsensor output. + + ! loop over the TS instances and compute sepSenMaxFamily for each TS instance + do sps = 1, nTimeIntervalsSpectral + ! set the local sepsensor to a smaller value so that we get the actual max + sepsensor_local = -10000.0_realType + do nn = 1, nDom + call setPointers(nn, 1, sps) + + ! Loop over all possible boundary conditions + bocos: do mm = 1, nBocos + ! Determine if this boundary condition is to be incldued in the + ! currently active group + famInclude: if (famInList(BCData(mm)%famID, famList)) then + + ! Set a bunch of pointers depending on the face id to make + ! a generic treatment possible. + call setBCPointers(mm, .True.) + + ! no post gathered integrations currently + isWall: if (isWallType(BCType(mm))) then + + !$AD II-LOOP + do ii = 0, (BCData(mm)%jnEnd - bcData(mm)%jnBeg) * (bcData(mm)%inEnd - bcData(mm)%inBeg) - 1 + i = mod(ii, (bcData(mm)%inEnd - bcData(mm)%inBeg)) + bcData(mm)%inBeg + 1 + j = ii / (bcData(mm)%inEnd - bcData(mm)%inBeg) + bcData(mm)%jnBeg + 1 + + ! only take this if its a compute cell + if (BCData(mm)%iblank(i, j) .eq. one) then + ! Get normalized surface velocity: + v(1) = ww2(i, j, ivx) + v(2) = ww2(i, j, ivy) + v(3) = ww2(i, j, ivz) + v = v / (sqrt(v(1)**2 + v(2)**2 + v(3)**2) + 1e-16) + + vectDotProductFsNormal = velDirFreeStream(1) * BCData(mm)%norm(i, j, 1) + & + velDirFreeStream(2) * BCData(mm)%norm(i, j, 2) + & + velDirFreeStream(3) * BCData(mm)%norm(i, j, 3) + ! Tangential Vector on the surface, which is the freestream projected vector + vectTangential(1) = velDirFreeStream(1) - & + vectDotProductFsNormal * BCData(mm)%norm(i, j, 1) + vectTangential(2) = velDirFreeStream(2) - & + vectDotProductFsNormal * BCData(mm)%norm(i, j, 2) + vectTangential(3) = velDirFreeStream(3) - & + vectDotProductFsNormal * BCData(mm)%norm(i, j, 3) + + vectTangential = vectTangential / & + (sqrt(vectTangential(1)**2 + vectTangential(2)**2 + & + vectTangential(3)**2) + 1e-16) + + ! computing separation sensor + ! velocity dot products + sensor = (v(1) * vectTangential(1) + v(2) * vectTangential(2) + & + v(3) * vectTangential(3)) + + ! sepsensor value + sensor = (cos(degtorad * sepangledeviation) - sensor) / & + (-cos(degtorad * sepangledeviation) + cos(zero) + 1e-16) + + ! compare it against the current value on this proc + sepsensor_local = max(sepsensor_local, sensor) + end if + + ! end if + end do + end if isWall + + end if famInclude + end do bocos + end do + ! finally communicate across all processors for this time spectral instance + call mpi_allreduce(sepsensor_local, sepSenMaxFamily(sps), 1, MPI_DOUBLE, & + MPI_MAX, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + end do + + end subroutine computeSepSenMaxFamily + #ifndef USE_COMPLEX subroutine integrateSurfaces_d(localValues, localValuesd, famList) !------------------------------------------------------------------------ @@ -1468,7 +1650,7 @@ subroutine getSolution_d(famLists, funcValues, funcValuesd) use zipperIntegrations, only: integrateZippers_d use userSurfaceIntegrations, only: integrateUserSurfaces_d use actuatorRegion, only: integrateActuatorRegions_d - use inputCostFunctions, only: computeCavitation + use inputCostFunctions, only: computeCavitation, computeSepSensorKs implicit none ! Input/Output Variables @@ -1492,6 +1674,11 @@ subroutine getSolution_d(famLists, funcValues, funcValuesd) call computeCpMinFamily(famList) end if + ! compute the current sepsensor max value for the separation computation with KS aggregation + if (computeSepSensorKs) then + call computeSepSenMaxFamily(famList) + end if + localVal = zero localVald = zero do sps = 1, nTimeIntervalsSpectral @@ -1550,7 +1737,7 @@ subroutine getSolution_b(famLists, funcValues, funcValuesd) use zipperIntegrations, only: integrateZippers_b use userSurfaceIntegrations, only: integrateUserSurfaces_b use actuatorRegion, only: integrateActuatorRegions_b - use inputCostFunctions, only: computeCavitation + use inputCostFunctions, only: computeCavitation, computeSepSensorKs implicit none ! Input/Output Variables @@ -1577,6 +1764,11 @@ subroutine getSolution_b(famLists, funcValues, funcValuesd) call computeCpMinFamily(famList) end if + ! compute the current sepsensor max value for the separation computation with KS aggregation + if (computeSepSensorKs) then + call computeSepSenMaxFamily(famList) + end if + localVal = zero localVald = zero diff --git a/src/utils/utils.F90 b/src/utils/utils.F90 index b50cc8c17..b36fc0cbf 100644 --- a/src/utils/utils.F90 +++ b/src/utils/utils.F90 @@ -4737,7 +4737,7 @@ subroutine releaseMemoryPart2 ! use block use inputTimeSpectral - use inputPhysics, only: cpmin_family + use inputPhysics, only: cpmin_family, sepSenMaxFamily use ADjointPETSc use cgnsGrid implicit none @@ -4769,6 +4769,10 @@ subroutine releaseMemoryPart2 if (allocated(cpmin_family)) & deallocate (cpmin_family) + ! deallocate the sepSenMaxFamily array allocated in inputParamRoutines + if (allocated(sepSenMaxFamily)) & + deallocate (sepSenMaxFamily) + ! Destroy variables allocated in preprocessingAdjoint if (adjointPETScPreProcVarsAllocated) then call vecDestroy(w_like1, PETScIerr) @@ -6168,6 +6172,9 @@ subroutine convergenceHeader case (cgnsSepSensor) write (*, "(a)", advance="no") " SepSensor |" + case (cgnsSepSensorArea) + write (*, "(a)", advance="no") " SepSensorArea |" + case (cgnsCavitation) write (*, "(a)", advance="no") " Cavitation |" @@ -6270,6 +6277,9 @@ subroutine convergenceHeader case (cgnsSepSensor) write (*, "(a)", advance="no") " SepSensor |" + case (cgnsSepSensorArea) + write (*, "(a)", advance="no") " SepSensorArea |" + case (cgnsCavitation) write (*, "(a)", advance="no") " Cavitation |" diff --git a/tests/reg_tests/refs/separation_tests_sepsensor.json b/tests/reg_tests/refs/separation_tests_sepsensor.json new file mode 100644 index 000000000..2fec4c3cb --- /dev/null +++ b/tests/reg_tests/refs/separation_tests_sepsensor.json @@ -0,0 +1,275 @@ +{ + "Dot product test for w -> sepsensor_wingup": -3.413578625244386e-05, + "Dot product test for xV -> sepsensor_wingup": -0.1708149227129096, + "metadata": { + "ADPC": false, + "ANKADPC": false, + "ANKAMGLevels": 2, + "ANKAMGNSmooth": 1, + "ANKASMOverlap": 1, + "ANKCFL0": 5.0, + "ANKCFLCutback": 0.5, + "ANKCFLExponent": 0.5, + "ANKCFLFactor": 10.0, + "ANKCFLLimit": 100000.0, + "ANKCFLMin": 1.0, + "ANKCFLReset": true, + "ANKCharTimeStepType": "None", + "ANKConstCFLStep": 0.4, + "ANKCoupledSwitchTol": 0.0001, + "ANKGlobalPreconditioner": "additive Schwarz", + "ANKInnerPreconIts": 1, + "ANKJacobianLag": 10, + "ANKLinResMax": 0.1, + "ANKLinearSolveBuffer": 0.01, + "ANKLinearSolveTol": 0.05, + "ANKMaxIter": 40, + "ANKNSubiterTurb": 1, + "ANKOuterPreconIts": 1, + "ANKPCILUFill": 2, + "ANKPCUpdateCutoff": 1e-16, + "ANKPCUpdateTol": 0.5, + "ANKPCUpdateTolAfterCutoff": 0.0001, + "ANKPhysicalLSTol": 0.2, + "ANKPhysicalLSTolTurb": 0.99, + "ANKSecondOrdSwitchTol": 0.001, + "ANKStepFactor": 1.0, + "ANKStepMin": 0.01, + "ANKSubspaceSize": -1, + "ANKSwitchTol": 1000.0, + "ANKTurbCFLScale": 1.0, + "ANKTurbKSPDebug": false, + "ANKUnsteadyLSTol": 1.0, + "ANKUseApproxSA": false, + "ANKUseFullVisc": true, + "ANKUseMatrixFree": true, + "ANKUseTurbDADI": true, + "ASMOverlap": 1, + "CFL": 1.7, + "CFLCoarse": 1.0, + "CFLLimit": 1.5, + "GMRESOrthogonalizationType": "modified Gram-Schmidt", + "ILUFill": 2, + "L2Convergence": 1e-15, + "L2ConvergenceCoarse": 0.01, + "L2ConvergenceRel": 1e-16, + "MGCycle": "sg", + "MGStartLevel": -1, + "NKADPC": false, + "NKAMGLevels": 2, + "NKAMGNSmooth": 1, + "NKASMOverlap": 1, + "NKFixedStep": 0.25, + "NKGlobalPreconditioner": "additive Schwarz", + "NKInnerPreconIts": 1, + "NKJacobianLag": 20, + "NKLS": "cubic", + "NKLinearSolveTol": 0.3, + "NKOuterPreconIts": 1, + "NKPCILUFill": 2, + "NKSubspaceSize": 60, + "NKSwitchTol": 1e-05, + "NKUseEW": true, + "NKViscPC": false, + "RKReset": false, + "TSStability": false, + "acousticScaleFactor": 1.0, + "adjointAMGLevels": 2, + "adjointAMGNSmooth": 1, + "adjointDivTol": 100000.0, + "adjointL2Convergence": 1e-15, + "adjointL2ConvergenceAbs": 1e-16, + "adjointL2ConvergenceRel": 1e-16, + "adjointMaxIter": 500, + "adjointMaxL2DeviationFactor": 1.0, + "adjointMonitorStep": 10, + "adjointSolver": "GMRES", + "adjointSubspaceSize": 100, + "alphaFollowing": true, + "alphaMode": false, + "altitudeMode": false, + "applyAdjointPCSubspaceSize": 20, + "applyPCSubspaceSize": 10, + "approxPC": true, + "backgroundVolScale": 1.0, + "betaMode": false, + "blockSplitting": true, + "cavExponent": 0, + "cavSensorOffset": 0.0, + "cavSensorSharpness": 10.0, + "cavitationNumber": 1.4, + "closedSurfaceFamilies": null, + "coarseDiscretization": "central plus scalar dissipation", + "computeCavitation": false, + "computeSepSensorKs": false, + "coupledSolution": false, + "cpMinRho": 100.0, + "cutCallback": null, + "debugZipper": false, + "deltaT": 0.01, + "designSurfaceFamily": null, + "discretization": "central plus scalar dissipation", + "dissContMagnitude": 1.0, + "dissContMidpoint": 3.0, + "dissContSharpness": 3.0, + "dissipationLumpingParameter": 6.0, + "dissipationScalingExponent": 0.67, + "eddyVisInfRatio": 0.009, + "equationMode": "steady", + "equationType": "RANS", + "eulerWallTreatment": "linear pressure extrapolation", + "explicitSurfaceCallback": null, + "firstRun": true, + "flowType": "external", + "forcesAsTractions": true, + "frozenTurbulence": false, + "globalPreconditioner": "additive Schwarz", + "gridFile": "input_files/naca0012_L3_SEP.cgns", + "gridPrecision": "double", + "gridPrecisionSurface": "single", + "infChangeCorrection": true, + "infChangeCorrectionTol": 1e-12, + "infChangeCorrectionType": "offset", + "innerPreconIts": 1, + "isoVariables": [], + "isosurface": {}, + "liftIndex": 2, + "limiter": "van Albada", + "loadBalanceIter": 10, + "loadImbalance": 0.1, + "localPreconditioner": "ILU", + "lowSpeedPreconditioner": false, + "machMode": false, + "matrixOrdering": "RCM", + "maxL2DeviationFactor": 1.0, + "meshMaxSkewness": 1.0, + "meshSurfaceFamily": null, + "monitorVariables": [ + "cl", + "cd", + "sepsensor" + ], + "nCycles": 5000, + "nCyclesCoarse": 500, + "nFloodIter": -1, + "nRKReset": 5, + "nRefine": 10, + "nSaveSurface": 1, + "nSaveVolume": 1, + "nSubiter": 1, + "nSubiterTurb": 35, + "nTimeStepsCoarse": 48, + "nTimeStepsFine": 400, + "nearWallDist": 0.1, + "numberSolutions": true, + "outerPreconIts": 3, + "outputDirectory": "tests/output_files", + "outputSurfaceFamily": "allSurfaces", + "overlapFactor": 0.9, + "oversetDebugPrint": false, + "oversetLoadBalance": true, + "oversetPriority": {}, + "oversetProjTol": 1e-12, + "oversetUpdateMode": "frozen", + "pMode": false, + "partitionLikeNProc": -1, + "partitionOnly": false, + "preconditionerSide": "right", + "printAllOptions": true, + "printBadlySkewedCells": false, + "printIntro": true, + "printIterations": true, + "printNegativeVolumes": false, + "printTiming": true, + "printWarnings": true, + "qMode": false, + "rMode": false, + "resAveraging": "alternate", + "restartAdjoint": true, + "restartFile": null, + "restrictionRelaxation": 0.8, + "selfZipCutoff": 120.0, + "sepAngleDeviation": 90.0, + "sepSensorMaxRho": 1000.0, + "sepSensorOffset": 0.0, + "sepSensorOffsetOne": 0.0, + "sepSensorOffsetTwo": 0.0, + "sepSensorSharpness": 10.0, + "sepSensorSharpnessOne": 25.0, + "sepSensorSharpnessTwo": 25.0, + "setMonitor": true, + "skipAfterFailedAdjoint": true, + "smoothParameter": 1.5, + "smoother": "Runge-Kutta", + "solutionPrecision": "single", + "solutionPrecisionSurface": "single", + "storeConvHist": true, + "storeRindLayer": true, + "surfaceVariables": [ + "cp", + "vx", + "vy", + "vz", + "mach" + ], + "timeAccuracy": 2, + "timeIntegrationScheme": "BDF", + "timeIntervals": 1, + "timeLimit": -1.0, + "turbResScale": 10000.0, + "turbulenceModel": "SA", + "turbulenceOrder": "first order", + "turbulenceProduction": "strain", + "useALE": true, + "useANKSolver": true, + "useApproxWallDistance": true, + "useBlockettes": true, + "useDiagTSPC": true, + "useDissContinuation": false, + "useExternalDynamicMesh": false, + "useGridMotion": false, + "useLinResMonitor": false, + "useMatrixFreedrdw": true, + "useNKSolver": false, + "useOversetWallScaling": false, + "useQCR": false, + "useRotationSA": false, + "useSkewnessCheck": false, + "useTSInterpolatedGridVelocity": false, + "useWallFunctions": false, + "useZipperMesh": true, + "useft2SA": true, + "verifyExtra": true, + "verifySpatial": true, + "verifyState": true, + "vis2": 0.25, + "vis2Coarse": 0.5, + "vis4": 0.0156, + "viscPC": false, + "viscWallTreatment": "constant pressure extrapolation", + "viscousSurfaceVelocities": true, + "volumeVariables": [ + "temp", + "mach", + "resrho", + "cp" + ], + "wallDistCutoff": 1e+20, + "windAxis": false, + "writeSolutionDigits": 3, + "writeSolutionEachIter": false, + "writeSurfaceSolution": false, + "writeTecplotSurfaceSolution": false, + "writeVolumeSolution": false, + "zipperSurfaceFamily": null + }, + "separation outputs": { + "naca0012_rans_2D_sepsensor_sepsensor_wingup": 0.9436658909967628 + }, + "separation totals": { + "naca0012_rans_2D_sepsensor_sepsensor_wingup": { + "alpha": -2.0453157612660585e-09 + } + }, + "total sepsensor_wingup derivative wrt random volume perturbation": -0.17081860114602293 +} \ No newline at end of file diff --git a/tests/reg_tests/refs/separation_tests_sepsensorarea.json b/tests/reg_tests/refs/separation_tests_sepsensorarea.json new file mode 100644 index 000000000..f3a88b25c --- /dev/null +++ b/tests/reg_tests/refs/separation_tests_sepsensorarea.json @@ -0,0 +1,289 @@ +{ + "Dot product test for w -> sepsensor_wingup": -3.413578625244386e-05, + "Dot product test for w -> sepsensorarea_wingup": 3.156991882291173e-22, + "Dot product test for w -> sepsensorks_wingup": -0.012285733905643411, + "Dot product test for xV -> sepsensor_wingup": -0.1708149227129096, + "Dot product test for xV -> sepsensorarea_wingup": -0.17081501916514247, + "Dot product test for xV -> sepsensorks_wingup": 0.0, + "metadata": { + "ADPC": false, + "ANKADPC": false, + "ANKAMGLevels": 2, + "ANKAMGNSmooth": 1, + "ANKASMOverlap": 1, + "ANKCFL0": 5.0, + "ANKCFLCutback": 0.5, + "ANKCFLExponent": 0.5, + "ANKCFLFactor": 10.0, + "ANKCFLLimit": 100000.0, + "ANKCFLMin": 1.0, + "ANKCFLReset": true, + "ANKCharTimeStepType": "None", + "ANKConstCFLStep": 0.4, + "ANKCoupledSwitchTol": 0.0001, + "ANKGlobalPreconditioner": "additive Schwarz", + "ANKInnerPreconIts": 1, + "ANKJacobianLag": 10, + "ANKLinResMax": 0.1, + "ANKLinearSolveBuffer": 0.01, + "ANKLinearSolveTol": 0.05, + "ANKMaxIter": 40, + "ANKNSubiterTurb": 1, + "ANKOuterPreconIts": 1, + "ANKPCILUFill": 2, + "ANKPCUpdateCutoff": 1e-16, + "ANKPCUpdateTol": 0.5, + "ANKPCUpdateTolAfterCutoff": 0.0001, + "ANKPhysicalLSTol": 0.2, + "ANKPhysicalLSTolTurb": 0.99, + "ANKSecondOrdSwitchTol": 0.001, + "ANKStepFactor": 1.0, + "ANKStepMin": 0.01, + "ANKSubspaceSize": -1, + "ANKSwitchTol": 1000.0, + "ANKTurbCFLScale": 1.0, + "ANKTurbKSPDebug": false, + "ANKUnsteadyLSTol": 1.0, + "ANKUseApproxSA": false, + "ANKUseFullVisc": true, + "ANKUseMatrixFree": true, + "ANKUseTurbDADI": true, + "ASMOverlap": 1, + "CFL": 1.7, + "CFLCoarse": 1.0, + "CFLLimit": 1.5, + "GMRESOrthogonalizationType": "modified Gram-Schmidt", + "ILUFill": 2, + "L2Convergence": 1e-15, + "L2ConvergenceCoarse": 0.01, + "L2ConvergenceRel": 1e-16, + "MGCycle": "sg", + "MGStartLevel": -1, + "NKADPC": false, + "NKAMGLevels": 2, + "NKAMGNSmooth": 1, + "NKASMOverlap": 1, + "NKFixedStep": 0.25, + "NKGlobalPreconditioner": "additive Schwarz", + "NKInnerPreconIts": 1, + "NKJacobianLag": 20, + "NKLS": "cubic", + "NKLinearSolveTol": 0.3, + "NKOuterPreconIts": 1, + "NKPCILUFill": 2, + "NKSubspaceSize": 60, + "NKSwitchTol": 1e-05, + "NKUseEW": true, + "NKViscPC": false, + "RKReset": false, + "TSStability": false, + "acousticScaleFactor": 1.0, + "adjointAMGLevels": 2, + "adjointAMGNSmooth": 1, + "adjointDivTol": 100000.0, + "adjointL2Convergence": 1e-15, + "adjointL2ConvergenceAbs": 1e-16, + "adjointL2ConvergenceRel": 1e-16, + "adjointMaxIter": 500, + "adjointMaxL2DeviationFactor": 1.0, + "adjointMonitorStep": 10, + "adjointSolver": "GMRES", + "adjointSubspaceSize": 100, + "alphaFollowing": true, + "alphaMode": false, + "altitudeMode": false, + "applyAdjointPCSubspaceSize": 20, + "applyPCSubspaceSize": 10, + "approxPC": true, + "backgroundVolScale": 1.0, + "betaMode": false, + "blockSplitting": true, + "cavExponent": 0, + "cavSensorOffset": 0.0, + "cavSensorSharpness": 10.0, + "cavitationNumber": 1.4, + "closedSurfaceFamilies": null, + "coarseDiscretization": "central plus scalar dissipation", + "computeCavitation": false, + "computeSepSensorKs": true, + "coupledSolution": false, + "cpMinRho": 100.0, + "cutCallback": null, + "debugZipper": false, + "deltaT": 0.01, + "designSurfaceFamily": null, + "discretization": "central plus scalar dissipation", + "dissContMagnitude": 1.0, + "dissContMidpoint": 3.0, + "dissContSharpness": 3.0, + "dissipationLumpingParameter": 6.0, + "dissipationScalingExponent": 0.67, + "eddyVisInfRatio": 0.009, + "equationMode": "steady", + "equationType": "RANS", + "eulerWallTreatment": "linear pressure extrapolation", + "explicitSurfaceCallback": null, + "firstRun": true, + "flowType": "external", + "forcesAsTractions": true, + "frozenTurbulence": false, + "globalPreconditioner": "additive Schwarz", + "gridFile": "input_files/naca0012_L3_SEP.cgns", + "gridPrecision": "double", + "gridPrecisionSurface": "single", + "infChangeCorrection": true, + "infChangeCorrectionTol": 1e-12, + "infChangeCorrectionType": "offset", + "innerPreconIts": 1, + "isoVariables": [], + "isosurface": {}, + "liftIndex": 2, + "limiter": "van Albada", + "loadBalanceIter": 10, + "loadImbalance": 0.1, + "localPreconditioner": "ILU", + "lowSpeedPreconditioner": false, + "machMode": false, + "matrixOrdering": "RCM", + "maxL2DeviationFactor": 1.0, + "meshMaxSkewness": 1.0, + "meshSurfaceFamily": null, + "monitorVariables": [ + "cl", + "cd", + "sepsensor" + ], + "nCycles": 5000, + "nCyclesCoarse": 500, + "nFloodIter": -1, + "nRKReset": 5, + "nRefine": 10, + "nSaveSurface": 1, + "nSaveVolume": 1, + "nSubiter": 1, + "nSubiterTurb": 35, + "nTimeStepsCoarse": 48, + "nTimeStepsFine": 400, + "nearWallDist": 0.1, + "numberSolutions": true, + "outerPreconIts": 3, + "outputDirectory": "tests/output_files", + "outputSurfaceFamily": "allSurfaces", + "overlapFactor": 0.9, + "oversetDebugPrint": false, + "oversetLoadBalance": true, + "oversetPriority": {}, + "oversetProjTol": 1e-12, + "oversetUpdateMode": "frozen", + "pMode": false, + "partitionLikeNProc": -1, + "partitionOnly": false, + "preconditionerSide": "right", + "printAllOptions": true, + "printBadlySkewedCells": false, + "printIntro": true, + "printIterations": true, + "printNegativeVolumes": false, + "printTiming": true, + "printWarnings": true, + "qMode": false, + "rMode": false, + "resAveraging": "alternate", + "restartAdjoint": true, + "restartFile": null, + "restrictionRelaxation": 0.8, + "selfZipCutoff": 120.0, + "sepAngleDeviation": 90.0, + "sepSensorMaxRho": 1000.0, + "sepSensorOffset": 0.0, + "sepSensorOffsetOne": 0.0, + "sepSensorOffsetTwo": 0.0, + "sepSensorSharpness": 10.0, + "sepSensorSharpnessOne": 25.0, + "sepSensorSharpnessTwo": 25.0, + "setMonitor": true, + "skipAfterFailedAdjoint": true, + "smoothParameter": 1.5, + "smoother": "Runge-Kutta", + "solutionPrecision": "single", + "solutionPrecisionSurface": "single", + "storeConvHist": true, + "storeRindLayer": true, + "surfaceVariables": [ + "cp", + "vx", + "vy", + "vz", + "mach" + ], + "timeAccuracy": 2, + "timeIntegrationScheme": "BDF", + "timeIntervals": 1, + "timeLimit": -1.0, + "turbResScale": 10000.0, + "turbulenceModel": "SA", + "turbulenceOrder": "first order", + "turbulenceProduction": "strain", + "useALE": true, + "useANKSolver": true, + "useApproxWallDistance": true, + "useBlockettes": true, + "useDiagTSPC": true, + "useDissContinuation": false, + "useExternalDynamicMesh": false, + "useGridMotion": false, + "useLinResMonitor": false, + "useMatrixFreedrdw": true, + "useNKSolver": false, + "useOversetWallScaling": false, + "useQCR": false, + "useRotationSA": false, + "useSkewnessCheck": false, + "useTSInterpolatedGridVelocity": false, + "useWallFunctions": false, + "useZipperMesh": true, + "useft2SA": true, + "verifyExtra": true, + "verifySpatial": true, + "verifyState": true, + "vis2": 0.25, + "vis2Coarse": 0.5, + "vis4": 0.0156, + "viscPC": false, + "viscWallTreatment": "constant pressure extrapolation", + "viscousSurfaceVelocities": true, + "volumeVariables": [ + "temp", + "mach", + "resrho", + "cp" + ], + "wallDistCutoff": 1e+20, + "windAxis": false, + "writeSolutionDigits": 3, + "writeSolutionEachIter": false, + "writeSurfaceSolution": false, + "writeTecplotSurfaceSolution": false, + "writeVolumeSolution": false, + "zipperSurfaceFamily": null + }, + "separation outputs": { + "naca0012_rans_2D_sepsensorarea_sepsensor_wingup": 0.9436658909967628, + "naca0012_rans_2D_sepsensorarea_sepsensorarea_wingup": 0.943665896753898, + "naca0012_rans_2D_sepsensorarea_sepsensorks_wingup": 1.0039701459789896 + }, + "separation totals": { + "naca0012_rans_2D_sepsensorarea_sepsensor_wingup": { + "alpha": -2.0453157612660585e-09 + }, + "naca0012_rans_2D_sepsensorarea_sepsensorarea_wingup": { + "alpha": -4.80058469332402e-40 + }, + "naca0012_rans_2D_sepsensorarea_sepsensorks_wingup": { + "alpha": 3.7782394013417314e-06 + } + }, + "total sepsensor_wingup derivative wrt random volume perturbation": -0.17081860114602293, + "total sepsensorarea_wingup derivative wrt random volume perturbation": -0.17081501916514247, + "total sepsensorks_wingup derivative wrt random volume perturbation": 0.1637025596134265 +} \ No newline at end of file diff --git a/tests/reg_tests/refs/separation_tests_sepsensorks.json b/tests/reg_tests/refs/separation_tests_sepsensorks.json new file mode 100644 index 000000000..06705788d --- /dev/null +++ b/tests/reg_tests/refs/separation_tests_sepsensorks.json @@ -0,0 +1,282 @@ +{ + "Dot product test for w -> sepsensor_wingup": -3.413578625244386e-05, + "Dot product test for w -> sepsensorks_wingup": -0.012285733905643411, + "Dot product test for xV -> sepsensor_wingup": -0.1708149227129096, + "Dot product test for xV -> sepsensorks_wingup": 0.0, + "metadata": { + "ADPC": false, + "ANKADPC": false, + "ANKAMGLevels": 2, + "ANKAMGNSmooth": 1, + "ANKASMOverlap": 1, + "ANKCFL0": 5.0, + "ANKCFLCutback": 0.5, + "ANKCFLExponent": 0.5, + "ANKCFLFactor": 10.0, + "ANKCFLLimit": 100000.0, + "ANKCFLMin": 1.0, + "ANKCFLReset": true, + "ANKCharTimeStepType": "None", + "ANKConstCFLStep": 0.4, + "ANKCoupledSwitchTol": 0.0001, + "ANKGlobalPreconditioner": "additive Schwarz", + "ANKInnerPreconIts": 1, + "ANKJacobianLag": 10, + "ANKLinResMax": 0.1, + "ANKLinearSolveBuffer": 0.01, + "ANKLinearSolveTol": 0.05, + "ANKMaxIter": 40, + "ANKNSubiterTurb": 1, + "ANKOuterPreconIts": 1, + "ANKPCILUFill": 2, + "ANKPCUpdateCutoff": 1e-16, + "ANKPCUpdateTol": 0.5, + "ANKPCUpdateTolAfterCutoff": 0.0001, + "ANKPhysicalLSTol": 0.2, + "ANKPhysicalLSTolTurb": 0.99, + "ANKSecondOrdSwitchTol": 0.001, + "ANKStepFactor": 1.0, + "ANKStepMin": 0.01, + "ANKSubspaceSize": -1, + "ANKSwitchTol": 1000.0, + "ANKTurbCFLScale": 1.0, + "ANKTurbKSPDebug": false, + "ANKUnsteadyLSTol": 1.0, + "ANKUseApproxSA": false, + "ANKUseFullVisc": true, + "ANKUseMatrixFree": true, + "ANKUseTurbDADI": true, + "ASMOverlap": 1, + "CFL": 1.7, + "CFLCoarse": 1.0, + "CFLLimit": 1.5, + "GMRESOrthogonalizationType": "modified Gram-Schmidt", + "ILUFill": 2, + "L2Convergence": 1e-15, + "L2ConvergenceCoarse": 0.01, + "L2ConvergenceRel": 1e-16, + "MGCycle": "sg", + "MGStartLevel": -1, + "NKADPC": false, + "NKAMGLevels": 2, + "NKAMGNSmooth": 1, + "NKASMOverlap": 1, + "NKFixedStep": 0.25, + "NKGlobalPreconditioner": "additive Schwarz", + "NKInnerPreconIts": 1, + "NKJacobianLag": 20, + "NKLS": "cubic", + "NKLinearSolveTol": 0.3, + "NKOuterPreconIts": 1, + "NKPCILUFill": 2, + "NKSubspaceSize": 60, + "NKSwitchTol": 1e-05, + "NKUseEW": true, + "NKViscPC": false, + "RKReset": false, + "TSStability": false, + "acousticScaleFactor": 1.0, + "adjointAMGLevels": 2, + "adjointAMGNSmooth": 1, + "adjointDivTol": 100000.0, + "adjointL2Convergence": 1e-15, + "adjointL2ConvergenceAbs": 1e-16, + "adjointL2ConvergenceRel": 1e-16, + "adjointMaxIter": 500, + "adjointMaxL2DeviationFactor": 1.0, + "adjointMonitorStep": 10, + "adjointSolver": "GMRES", + "adjointSubspaceSize": 100, + "alphaFollowing": true, + "alphaMode": false, + "altitudeMode": false, + "applyAdjointPCSubspaceSize": 20, + "applyPCSubspaceSize": 10, + "approxPC": true, + "backgroundVolScale": 1.0, + "betaMode": false, + "blockSplitting": true, + "cavExponent": 0, + "cavSensorOffset": 0.0, + "cavSensorSharpness": 10.0, + "cavitationNumber": 1.4, + "closedSurfaceFamilies": null, + "coarseDiscretization": "central plus scalar dissipation", + "computeCavitation": false, + "computeSepSensorKs": true, + "coupledSolution": false, + "cpMinRho": 100.0, + "cutCallback": null, + "debugZipper": false, + "deltaT": 0.01, + "designSurfaceFamily": null, + "discretization": "central plus scalar dissipation", + "dissContMagnitude": 1.0, + "dissContMidpoint": 3.0, + "dissContSharpness": 3.0, + "dissipationLumpingParameter": 6.0, + "dissipationScalingExponent": 0.67, + "eddyVisInfRatio": 0.009, + "equationMode": "steady", + "equationType": "RANS", + "eulerWallTreatment": "linear pressure extrapolation", + "explicitSurfaceCallback": null, + "firstRun": true, + "flowType": "external", + "forcesAsTractions": true, + "frozenTurbulence": false, + "globalPreconditioner": "additive Schwarz", + "gridFile": "input_files/naca0012_L3_SEP.cgns", + "gridPrecision": "double", + "gridPrecisionSurface": "single", + "infChangeCorrection": true, + "infChangeCorrectionTol": 1e-12, + "infChangeCorrectionType": "offset", + "innerPreconIts": 1, + "isoVariables": [], + "isosurface": {}, + "liftIndex": 2, + "limiter": "van Albada", + "loadBalanceIter": 10, + "loadImbalance": 0.1, + "localPreconditioner": "ILU", + "lowSpeedPreconditioner": false, + "machMode": false, + "matrixOrdering": "RCM", + "maxL2DeviationFactor": 1.0, + "meshMaxSkewness": 1.0, + "meshSurfaceFamily": null, + "monitorVariables": [ + "cl", + "cd", + "sepsensor" + ], + "nCycles": 5000, + "nCyclesCoarse": 500, + "nFloodIter": -1, + "nRKReset": 5, + "nRefine": 10, + "nSaveSurface": 1, + "nSaveVolume": 1, + "nSubiter": 1, + "nSubiterTurb": 35, + "nTimeStepsCoarse": 48, + "nTimeStepsFine": 400, + "nearWallDist": 0.1, + "numberSolutions": true, + "outerPreconIts": 3, + "outputDirectory": "tests/output_files", + "outputSurfaceFamily": "allSurfaces", + "overlapFactor": 0.9, + "oversetDebugPrint": false, + "oversetLoadBalance": true, + "oversetPriority": {}, + "oversetProjTol": 1e-12, + "oversetUpdateMode": "frozen", + "pMode": false, + "partitionLikeNProc": -1, + "partitionOnly": false, + "preconditionerSide": "right", + "printAllOptions": true, + "printBadlySkewedCells": false, + "printIntro": true, + "printIterations": true, + "printNegativeVolumes": false, + "printTiming": true, + "printWarnings": true, + "qMode": false, + "rMode": false, + "resAveraging": "alternate", + "restartAdjoint": true, + "restartFile": null, + "restrictionRelaxation": 0.8, + "selfZipCutoff": 120.0, + "sepAngleDeviation": 90.0, + "sepSensorMaxRho": 1000.0, + "sepSensorOffset": 0.0, + "sepSensorOffsetOne": 0.0, + "sepSensorOffsetTwo": 0.0, + "sepSensorSharpness": 10.0, + "sepSensorSharpnessOne": 25.0, + "sepSensorSharpnessTwo": 25.0, + "setMonitor": true, + "skipAfterFailedAdjoint": true, + "smoothParameter": 1.5, + "smoother": "Runge-Kutta", + "solutionPrecision": "single", + "solutionPrecisionSurface": "single", + "storeConvHist": true, + "storeRindLayer": true, + "surfaceVariables": [ + "cp", + "vx", + "vy", + "vz", + "mach" + ], + "timeAccuracy": 2, + "timeIntegrationScheme": "BDF", + "timeIntervals": 1, + "timeLimit": -1.0, + "turbResScale": 10000.0, + "turbulenceModel": "SA", + "turbulenceOrder": "first order", + "turbulenceProduction": "strain", + "useALE": true, + "useANKSolver": true, + "useApproxWallDistance": true, + "useBlockettes": true, + "useDiagTSPC": true, + "useDissContinuation": false, + "useExternalDynamicMesh": false, + "useGridMotion": false, + "useLinResMonitor": false, + "useMatrixFreedrdw": true, + "useNKSolver": false, + "useOversetWallScaling": false, + "useQCR": false, + "useRotationSA": false, + "useSkewnessCheck": false, + "useTSInterpolatedGridVelocity": false, + "useWallFunctions": false, + "useZipperMesh": true, + "useft2SA": true, + "verifyExtra": true, + "verifySpatial": true, + "verifyState": true, + "vis2": 0.25, + "vis2Coarse": 0.5, + "vis4": 0.0156, + "viscPC": false, + "viscWallTreatment": "constant pressure extrapolation", + "viscousSurfaceVelocities": true, + "volumeVariables": [ + "temp", + "mach", + "resrho", + "cp" + ], + "wallDistCutoff": 1e+20, + "windAxis": false, + "writeSolutionDigits": 3, + "writeSolutionEachIter": false, + "writeSurfaceSolution": false, + "writeTecplotSurfaceSolution": false, + "writeVolumeSolution": false, + "zipperSurfaceFamily": null + }, + "separation outputs": { + "naca0012_rans_2D_sepsensorks_sepsensor_wingup": 0.9436658909967628, + "naca0012_rans_2D_sepsensorks_sepsensorks_wingup": 1.0039701459789896 + }, + "separation totals": { + "naca0012_rans_2D_sepsensorks_sepsensor_wingup": { + "alpha": -2.0453157612660585e-09 + }, + "naca0012_rans_2D_sepsensorks_sepsensorks_wingup": { + "alpha": 3.7782394013417314e-06 + } + }, + "total sepsensor_wingup derivative wrt random volume perturbation": -0.17081860114602293, + "total sepsensorks_wingup derivative wrt random volume perturbation": 0.1637025596134265 +} \ No newline at end of file diff --git a/tests/reg_tests/reg_aeroproblems.py b/tests/reg_tests/reg_aeroproblems.py index 2804df415..2d557243c 100644 --- a/tests/reg_tests/reg_aeroproblems.py +++ b/tests/reg_tests/reg_aeroproblems.py @@ -234,3 +234,13 @@ yRef=0.0, zRef=0.0, ) + +ap_naca0012_separation = AeroProblem( + name="0012separation", + alpha=20.0, + mach=0.35, + altitude=3048, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["sepsensor_wingup", "sepsensorks_wingup", "sepsensorarea_wingup"], +) diff --git a/tests/reg_tests/test_separation.py b/tests/reg_tests/test_separation.py new file mode 100644 index 000000000..740afa8d7 --- /dev/null +++ b/tests/reg_tests/test_separation.py @@ -0,0 +1,347 @@ +# built-ins +import unittest +import numpy as np +import os +import copy + +# MACH classes +from adflow import ADFLOW +from adflow import ADFLOW_C + +# import the testing utilities +from parameterized import parameterized_class +from reg_default_options import adflowDefOpts +from reg_aeroproblems import ap_naca0012_separation +import reg_test_classes + +baseDir = os.path.dirname(os.path.abspath(__file__)) +gridFile = os.path.join(baseDir, "../../input_files/naca0012_L3_SEP.cgns") +test_params = [ + { + "name": "naca0012_rans_2D_sepsensor", + "ref_file": "separation_tests_sepsensor.json", + "aero_prob": ap_naca0012_separation, + "eval_funcs": ["sepsensor_wingup"], + "N_PROCS": 2, + "options": { + "monitorvariables": ["cl", "cd", "sepsensor"], + }, + }, + { + "name": "naca0012_rans_2D_sepsensorks", + "ref_file": "separation_tests_sepsensorks.json", + "aero_prob": ap_naca0012_separation, + "eval_funcs": ["sepsensor_wingup", "sepsensorks_wingup"], + "N_PROCS": 2, + "options": { + "computeSepSensorKs": True, + # "sepangledeviation": 0.0, + "monitorvariables": ["cl", "cd", "sepsensor"], + }, + }, + { + "name": "naca0012_rans_2D_sepsensorarea", + "ref_file": "separation_tests_sepsensorarea.json", + "aero_prob": ap_naca0012_separation, + "eval_funcs": ["sepsensor_wingup", "sepsensorks_wingup", "sepsensorarea_wingup"], + "N_PROCS": 2, + "options": { + "computeSepSensorKs": True, + # "sepangledeviation": 0.0, + "monitorvariables": ["cl", "cd", "sepsensor"], + }, + }, +] + + +@parameterized_class(test_params) +class SeparationBasicTests(reg_test_classes.RegTest): + """ + Tests for the separation metrics. + """ + + N_PROCS = 2 + # ref_file = None + options = None + + def setUp(self): + if not hasattr(self, "name"): + # return immediately when the setup method is being called on the based class and NOT the + # classes created using parametrized + # this will happen when training, but will hopefully be fixed down the line + return + super().setUp() + + options = copy.copy(adflowDefOpts) + options.update( + { + "gridfile": gridFile, + "outputdirectory": os.path.join(baseDir, "../output_files"), + "writevolumesolution": False, + "writesurfacesolution": False, + "writetecplotsurfacesolution": False, + "mgcycle": "sg", + "ncycles": 5000, + "useanksolver": True, + # "usenksolver": True, + "nsubiterturb": 35, + "anksecondordswitchtol": 1e-3, + "ankcoupledswitchtol": 1e-4, + # "nkswitchtol": 1e-6, + "volumevariables": ["temp", "mach", "resrho", "cp"], + "equationType": "RANS", + "l2convergence": 1e-15, + "adjointl2convergence": 1e-15, + } + ) + options.update(self.options) + + # Setup aeroproblem + self.ap = copy.deepcopy(self.aero_prob) + # change the name + self.ap.name = self.name + + # add aoa DV + self.ap.addDV("alpha", name="alpha") + + CFDSolver = ADFLOW(options=options) + + self.CFDSolver = CFDSolver + self.CFDSolver.addFamilyGroup("wingup", ["upper_skin"]) + + self.CFDSolver.addFunction("sepsensor", "wingup") + self.CFDSolver.addFunction("sepsensorks", "wingup") + self.CFDSolver.addFunction("sepsensorarea", "wingup") + + def test_separation_metrics_and_derivatives(self): + evalFuncs = self.eval_funcs + + self.CFDSolver(self.ap) + + # check if solution failed + self.assert_solution_failure() + + ################ + # TEST ALL FUNCS + ################ + funcs = {} + self.CFDSolver.evalFunctions(self.ap, funcs, evalFuncs=evalFuncs) + self.handler.root_print("separation outputs") + self.handler.root_add_dict("separation outputs", funcs, rtol=5e-10, atol=5e-10) + + ############# + # TEST TOTALS + ############# + funcsSens = {} + self.CFDSolver.evalFunctionsSens(self.ap, funcsSens, evalFuncs=evalFuncs) + self.handler.root_print("separation totals") + self.handler.root_add_dict("separation totals", funcsSens, rtol=5e-10, atol=5e-10) + + for funcName in evalFuncs: + ################## + # GEOMETRIC TOTALS + ################## + + # random state and volume coordinate perturbations for tests below + wDot = self.CFDSolver.getStatePerturbation(314) + xVDot = self.CFDSolver.getSpatialPerturbation(314) + + # Now that we have the adjoint solution for both functionals + # also get a total derivative wrt a random spatial perturbation DV + psi = -self.CFDSolver.getAdjoint(funcName) + funcsBar = self.CFDSolver._getFuncsBar(funcName) + + # this is the reverse seed up to the volume coordinates + xVBar = self.CFDSolver.computeJacobianVectorProductBwd(resBar=psi, funcsBar=funcsBar, xVDeriv=True) + + # dot product these two vectors to get a total derivative + dotLocal = np.dot(xVDot, xVBar) + # this is not the best test out there; the final answer does get affected quite a bit + # by the processor count or architecture. We just have it here to have some test on the + # separation functionals' derivatives w.r.t. spatial changes in a lazy way. + self.handler.par_add_sum( + f"total {funcName} derivative wrt random volume perturbation", + dotLocal, + rtol=1e-3, + ) + + ################## + # DOT PRODUCT TEST + ################## + fDot_w = self.CFDSolver.computeJacobianVectorProductFwd(wDot=wDot, funcDeriv=True) + fDot_xv = self.CFDSolver.computeJacobianVectorProductFwd(xVDot=xVDot, funcDeriv=True) + funcsBar = {funcName: 1.0} + wBar, xVBar = self.CFDSolver.computeJacobianVectorProductBwd(funcsBar=funcsBar, wDeriv=True, xVDeriv=True) + + # do the dot product. we test both state partials and volume coordinate partials + # state + dotLocal1 = np.dot(wDot, wBar) + dotLocal2 = fDot_w[funcName] / self.CFDSolver.comm.size + + self.handler.par_add_sum(f"Dot product test for w -> {funcName}", dotLocal1, rtol=1e-8, atol=1e-8) + self.handler.par_add_sum( + f"Dot product test for w -> {funcName}", + dotLocal2, + rtol=1e-8, + atol=1e-8, + compare=True, + ) + + # volume coords + dotLocal1 = np.dot(xVDot, xVBar) + dotLocal2 = fDot_xv[funcName] / self.CFDSolver.comm.size + self.handler.par_add_sum(f"Dot product test for xV -> {funcName}", dotLocal1, rtol=1e-8, atol=1e-8) + self.handler.par_add_sum( + f"Dot product test for xV -> {funcName}", + dotLocal2, + rtol=5e-10, + compare=True, + ) + + +@parameterized_class(test_params) +class SeparationCmplxTests(reg_test_classes.CmplxRegTest): + """ + Complex step tests for the separation functions. + """ + + N_PROCS = 2 + ref_file = None + h = 1e-40 + + def setUp(self): + if not hasattr(self, "name"): + # return immediately when the setup method is being called on the based class and NOT the + # classes created using parametrized + # this will happen when training, but will hopefully be fixed down the line + return + super().setUp() + + options = copy.copy(adflowDefOpts) + options.update( + { + "gridfile": gridFile, + "outputdirectory": os.path.join(baseDir, "../output_files"), + "writevolumesolution": False, + "writesurfacesolution": False, + "writetecplotsurfacesolution": False, + "mgcycle": "sg", + "ncycles": 5000, + "useanksolver": True, + # "usenksolver": True, + "nsubiterturb": 35, + "anksecondordswitchtol": 1e-3, + "ankcoupledswitchtol": 1e-4, + # "nkswitchtol": 1e-6, + "volumevariables": ["temp", "mach", "resrho", "cp"], + "equationType": "RANS", + "l2convergence": 1e-15, + "adjointl2convergence": 1e-15, + } + ) + options.update(self.options) + + # Setup aeroproblem + self.ap = copy.deepcopy(ap_naca0012_separation) + # change the name + self.ap.name = self.name + + # add aoa DV + self.ap.addDV("alpha", name="alpha") + + CFDSolver = ADFLOW_C(options=options, debug=True) + + self.CFDSolver = CFDSolver + + self.CFDSolver.addFamilyGroup("wingup", ["upper_skin"]) + + self.CFDSolver.addFunction("sepsensor", "wingup") + self.CFDSolver.addFunction("sepsensorks", "wingup") + self.CFDSolver.addFunction("sepsensorarea", "wingup") + + def cmplx_test_separation_adjoints(self): + if not hasattr(self, "name"): + # return immediately when the setup method is being called on the based class and NOT the + # classes created using parametrized + # this will happen when training, but will hopefully be fixed down the line + return + aDV = {"alpha": self.ap.alpha} + funcs = {} + funcsSensCS = {} + + # run the baseline evaluation + self.CFDSolver(self.ap) + # check if solution failed + self.assert_solution_failure() + self.CFDSolver.evalFunctions(self.ap, funcs) + + evalFuncs = self.eval_funcs + funcs_plus = {} + for dv in ["alpha", "vol_perturbation"]: + if dv == "alpha": + # save the old alpha + dvsave = aDV["alpha"] + # perturb + aDV["alpha"] += self.h * 1j + self.ap.setDesignVars(aDV) + elif dv == "vol_perturbation": + xVDot = self.CFDSolver.getSpatialPerturbation(314) + # gridSave = np.zeros_like(xVDot) + # get the original grid + gridSave = self.CFDSolver.adflow.warping.getgrid(self.CFDSolver.getSpatialSize()) + # perturb using the random seed + # this is very intrusive but it works and its testing these functions sensitivities + # wrt geometric DVs w/o actually having a mesh or dvgeo object + self.CFDSolver.adflow.warping.setgrid(gridSave + self.h * 1j * xVDot) + self.CFDSolver.adflow.preprocessingapi.updatecoordinatesalllevels() + self.CFDSolver.adflow.walldistance.updatewalldistancealllevels() + self.CFDSolver.adflow.preprocessingapi.updatemetricsalllevels() + self.CFDSolver.adflow.preprocessingapi.updategridvelocitiesalllevels() + + # call solver again + # we can also not re-set the flow and call the solver a few times until complex parts converge + self.CFDSolver.resetFlow(self.ap) + self.CFDSolver(self.ap) # the complex residuals dont converge well for the coordinate perturbation + + # check if solution failed + self.assert_solution_failure() + + # save the new funcs in a dict + funcs_plus[dv] = {} + + # eval functions + self.CFDSolver.evalFunctions(self.ap, funcs_plus[dv]) + + # compute the sens + funcsSensCS[dv] = {} + for f in evalFuncs: + funcsSensCS[dv][self.ap[f]] = np.imag(funcs_plus[dv][self.ap[f]]) / self.h + + # reset the DV + if dv == "alpha": + aDV["alpha"] = dvsave + self.ap.setDesignVars(aDV) + elif dv == "vol_perturbation": + # set the original grid back + self.CFDSolver.adflow.warping.setgrid(gridSave) + + ################## + # TEST DERIVATIVES + ################## + + # here we compare the CS derivatives to the adjoint values computed by the real test + # we treat the CS value as the truth, so if this test passes, + # we assume the adjoint sensitivities are also true + + for funcName in evalFuncs: + fullName = f"{self.name}_{funcName}" + + refVal = self.handler.db["separation totals"][fullName]["alpha"] + np.testing.assert_allclose(funcsSensCS["alpha"][fullName], refVal, atol=1e-6, rtol=1e-6) + + refVal = self.handler.db[f"total {funcName} derivative wrt random volume perturbation"] + np.testing.assert_allclose(funcsSensCS["vol_perturbation"][fullName], refVal, rtol=9e-3, atol=5e-3) + + +if __name__ == "__main__": + unittest.main()