Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improved cavitation constraint #176

Merged
merged 10 commits into from
Jun 17, 2022
6 changes: 6 additions & 0 deletions adflow/pyADflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -5093,6 +5093,9 @@ def _getDefaultOptions():
# Function parmeters
"sepSensorOffset": [float, 0.0],
"sepSensorSharpness": [float, 10.0],
"cavSensorOffset": [float, 0.0],
"cavSensorSharpness": [float, 10.0],
"cavExponent": [int, 0],
"computeCavitation": [bool, False],
}

Expand Down Expand Up @@ -5468,6 +5471,9 @@ def _getOptionMap(self):
# Parameters for functions
"sepsensoroffset": ["cost", "sepsensoroffset"],
"sepsensorsharpness": ["cost", "sepsensorsharpness"],
"cavsensoroffset": ["cost", "cavsensoroffset"],
"cavsensorsharpness": ["cost", "cavsensorsharpness"],
"cavexponent": ["cost", "cavexponent"],
"computecavitation": ["cost", "computecavitation"],
"writesolutioneachiter": ["monitor", "writesoleachiter"],
"writesurfacesolution": ["monitor", "writesurface"],
Expand Down
24 changes: 20 additions & 4 deletions src/adjoint/outputForward/surfaceintegrations_d.f90
Original file line number Diff line number Diff line change
Expand Up @@ -875,6 +875,8 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm)
real(kind=realtype) :: arg1d
real(kind=realtype) :: result1
real(kind=realtype) :: result1d
real(kind=realtype) :: pwr1
real(kind=realtype) :: pwr1d
select case (bcfaceid(mm))
case (imin, jmin, kmin)
fact = -one
Expand Down Expand Up @@ -1116,9 +1118,20 @@ subroutine wallintegrationface_d(localvalues, localvaluesd, mm)
cp = tmp*(plocal-pinf)
sensor1d = -cpd
sensor1 = -cp - cavitationnumber
sensor1d = -((-(one*2*10*sensor1d*exp(-(2*10*sensor1))))/(one+&
& exp(-(2*10*sensor1)))**2)
sensor1 = one/(one+exp(-(2*10*sensor1)))
if (sensor1 .gt. 0.0_8 .or. (sensor1 .lt. 0.0_8 .and. &
& cavexponent .eq. int(cavexponent))) then
pwr1d = cavexponent*sensor1**(cavexponent-1)*sensor1d
else if (sensor1 .eq. 0.0_8 .and. cavexponent .eq. 1.0) then
pwr1d = sensor1d
else
pwr1d = 0.0_8
end if
pwr1 = sensor1**cavexponent
arg1d = -(2*cavsensorsharpness*sensor1d)
arg1 = -(2*cavsensorsharpness*(sensor1-cavsensoroffset))
sensor1d = (pwr1d*(one+exp(arg1))-pwr1*arg1d*exp(arg1))/(one+exp&
& (arg1))**2
sensor1 = pwr1/(one+exp(arg1))
sensor1d = blk*(sensor1d*cellarea+sensor1*cellaread)
sensor1 = sensor1*cellarea*blk
cavitationd = cavitationd + sensor1d
Expand Down Expand Up @@ -1357,6 +1370,7 @@ subroutine wallintegrationface(localvalues, mm)
intrinsic exp
real(kind=realtype) :: arg1
real(kind=realtype) :: result1
real(kind=realtype) :: pwr1
select case (bcfaceid(mm))
case (imin, jmin, kmin)
fact = -one
Expand Down Expand Up @@ -1507,7 +1521,9 @@ subroutine wallintegrationface(localvalues, mm)
tmp = two/(gammainf*machcoef*machcoef)
cp = tmp*(plocal-pinf)
sensor1 = -cp - cavitationnumber
sensor1 = one/(one+exp(-(2*10*sensor1)))
pwr1 = sensor1**cavexponent
arg1 = -(2*cavsensorsharpness*(sensor1-cavsensoroffset))
sensor1 = pwr1/(one+exp(arg1))
sensor1 = sensor1*cellarea*blk
cavitation = cavitation + sensor1
end if
Expand Down
21 changes: 16 additions & 5 deletions src/adjoint/outputReverse/surfaceintegrations_b.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1271,7 +1271,8 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm)
tmp = two/(gammainf*machcoef*machcoef)
cp = tmp*(plocal-pinf)
sensor1 = -cp - cavitationnumber
sensor1 = one/(one+exp(-(2*10*sensor1)))
sensor1 = sensor1**cavexponent/(one+exp(-(2*cavsensorsharpness*(&
& sensor1-cavsensoroffset))))
sensor1 = sensor1*cellarea*blk
cavitation = cavitation + sensor1
end if
Expand Down Expand Up @@ -1575,14 +1576,23 @@ subroutine wallintegrationface_b(localvalues, localvaluesd, mm)
cp = tmp*(plocal-pinf)
sensor1 = -cp - cavitationnumber
call pushreal8(sensor1)
sensor1 = one/(one+exp(-(2*10*sensor1)))
sensor1 = sensor1**cavexponent/(one+exp(-(2*cavsensorsharpness*(&
& sensor1-cavsensoroffset))))
sensor1d = cavitationd
cellaread = blk*sensor1*sensor1d
sensor1d = blk*cellarea*sensor1d
call popreal8(sensor1)
temp6 = -(10*2*sensor1)
temp6 = -(2*cavsensorsharpness*(sensor1-cavsensoroffset))
temp5 = one + exp(temp6)
sensor1d = exp(temp6)*one*10*2*sensor1d/temp5**2
if (sensor1 .le. 0.0_8 .and. (cavexponent .eq. 0.0_8 .or. &
& cavexponent .ne. int(cavexponent))) then
sensor1d = exp(temp6)*sensor1**cavexponent*cavsensorsharpness*&
& 2*sensor1d/temp5**2
else
sensor1d = (exp(temp6)*sensor1**cavexponent*cavsensorsharpness&
& *2/temp5**2+cavexponent*sensor1**(cavexponent-1)/temp5)*&
& sensor1d
end if
cpd = -sensor1d
tmpd = (plocal-pinf)*cpd
plocald = tmp*cpd
Expand Down Expand Up @@ -1940,7 +1950,8 @@ subroutine wallintegrationface(localvalues, mm)
tmp = two/(gammainf*machcoef*machcoef)
cp = tmp*(plocal-pinf)
sensor1 = -cp - cavitationnumber
sensor1 = one/(one+exp(-(2*10*sensor1)))
sensor1 = sensor1**cavexponent/(one+exp(-(2*cavsensorsharpness*(&
& sensor1-cavsensoroffset))))
sensor1 = sensor1*cellarea*blk
cavitation = cavitation + sensor1
end if
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -458,7 +458,8 @@ subroutine wallintegrationface(localvalues, mm)
tmp = two/(gammainf*machcoef*machcoef)
cp = tmp*(plocal-pinf)
sensor1 = -cp - cavitationnumber
sensor1 = one/(one+exp(-(2*10*sensor1)))
sensor1 = sensor1**cavexponent/(one+exp(-(2*cavsensorsharpness*(&
& sensor1-cavsensoroffset))))
sensor1 = sensor1*cellarea*blk
cavitation = cavitation + sensor1
end if
Expand Down
3 changes: 3 additions & 0 deletions src/f2py/adflow.pyf
Original file line number Diff line number Diff line change
Expand Up @@ -1385,6 +1385,9 @@ python module libadflow
use constants
real(kind=realtype) sepsensoroffset
real(kind=realtype) sepsensorsharpness
real(kind=realtype) cavsensoroffset
real(kind=realtype) cavsensorsharpness
integer(kind=inttype) cavexponent
logical::computecavitation
end module inputcostfunctions

Expand Down
3 changes: 3 additions & 0 deletions src/inputParam/inputParamRoutines.F90
Original file line number Diff line number Diff line change
Expand Up @@ -4042,6 +4042,9 @@ subroutine setDefaultValues
usematrixfreedrdw = .False.
sepSensorOffset = zero
sepSensorSharpness = 10_realType
cavSensorOffset = zero
cavSensorSharpness = 10_realType
cavExponent = 0
yqliaohk marked this conversation as resolved.
Show resolved Hide resolved
end subroutine setDefaultValues

subroutine initializeIsoSurfaceVariables(values, nValues)
Expand Down
7 changes: 5 additions & 2 deletions src/modules/inputParam.F90
Original file line number Diff line number Diff line change
Expand Up @@ -290,8 +290,11 @@ end module inputIteration

module inputCostFunctions
use constants
real(kind=realtype) :: sepSensorOffset= zero
real(kind=realtype) ::sepSensorSharpness=10.0_realType
real(kind=realtype) :: sepSensorOffset = zero
real(kind=realtype) :: sepSensorSharpness = 10.0_realType
real(kind=realtype) :: cavSensorOffset = zero
real(kind=realtype) :: cavSensorSharpness = 10.0_realType
integer(kind=inttype) :: cavExponent = 0
yqliaohk marked this conversation as resolved.
Show resolved Hide resolved
logical :: computeCavitation

end module inputCostFunctions
Expand Down
28 changes: 14 additions & 14 deletions src/output/outputMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1405,8 +1405,8 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
real(kind=realType), dimension(*), intent(out) :: buffer
character(len=*), intent(in) :: solName
logical, intent(in) :: viscousSubface, useRindLayer
! if useRindLayer is true, then iBeg, iEnd, jBeg, jEnd are use to determine

! if useRindLayer is true, then iBeg, iEnd, jBeg, jEnd are use to determine
! when the indices are in the rind layer.
integer(kind=intType), optional, intent(in) :: iBeg, iEnd, jBeg, jEnd
!
Expand Down Expand Up @@ -1999,7 +1999,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
else if (jor == jEnd +1 ) then
jj = j - 1
else
jj = j
jj = j
endif
else
jj = j
Expand All @@ -2014,12 +2014,12 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
else if (ior == iEnd + 1) then
ii = i - 1
else
ii = i
ii = i
endif
else
ii = i
endif

! Determine the viscous subface on which this
! face is located.

Expand Down Expand Up @@ -2105,29 +2105,29 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
! if statements are used to copy the value of the interior
! cell since the value isn't defined in the rind cell

if (present(jBeg) .and. present(jEnd) .and. (useRindLayer)) then
if (present(jBeg) .and. present(jEnd) .and. (useRindLayer)) then
jor = j + jBegOr - 1
if (jor == jBeg) then
jj = j + 1
if (jor == jBeg) then
jj = j + 1
else if (jor == jEnd + 1) then
jj = j - 1
else
jj = j
jj = j
endif
else
jj = j

end if

do i=rangeFace(1,1), rangeFace(1,2)
if (present(iBeg) .and. present( iEnd) .and. (useRindLayer)) then
if (present(iBeg) .and. present( iEnd) .and. (useRindLayer)) then
ior = i + iBegor - 1
if (ior == iBeg) then
ii = i + 1
if (ior == iBeg) then
ii = i + 1
else if (ior == iEnd + 1) then
ii = i - 1
else
ii = i
ii = i
endif
else
ii = i
Expand Down Expand Up @@ -2206,7 +2206,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
plocal = half*(pp1(i,j) + pp2(i,j))

sensor1 = (-(fact)*(plocal-pInf))- cavitationnumber
sensor1 = one/(one + exp(-2*10*sensor1))
sensor1 = (sensor1**cavExponent)/(one + exp(-2*cavSensorSharpness*(sensor1 - cavSensorOffset)))
yqliaohk marked this conversation as resolved.
Show resolved Hide resolved
buffer(nn) = sensor1
!print*, sensor
enddo
Expand Down
2 changes: 1 addition & 1 deletion src/solver/surfaceIntegrations.F90
Original file line number Diff line number Diff line change
Expand Up @@ -505,7 +505,7 @@ subroutine wallIntegrationFace(localValues, mm)
tmp = two/(gammaInf*MachCoef*MachCoef)
Cp = tmp*(plocal-pinf)
Sensor1 = -Cp - cavitationnumber
Sensor1 = one/(one+exp(-2*10*Sensor1))
Sensor1 = (Sensor1**cavExponent)/(one+exp(-2*cavSensorSharpness*(Sensor1-cavSensorOffset)))
yqliaohk marked this conversation as resolved.
Show resolved Hide resolved
Sensor1 = Sensor1 * cellArea * blk
Cavitation = Cavitation + Sensor1
end if
Expand Down