SA rough

.gitignore

@@ -25,6 +25,12 @@ src/build/libadflow-f2pywrappers.f
 src/build/libadflowmodule.c
 src/build/importTest.py
 
+# unused fast_b routines generated by Tapenade
+src/adjoint/outputReverseFast/adjointextra_fast_b.f90
+src/adjoint/outputReverseFast/bcdata_fast_b.f90
+src/adjoint/outputReverseFast/oversetutilities_fast_b.f90
+src/adjoint/outputReverseFast/zipperintegrations_fast_b.f90
+
 # Regression test auxiliary files:
 input_files.tar.gz
 adflow_input_files.tar.gz

adflow/pyADflow.py

@@ -141,6 +141,9 @@ def __init__(self, comm=None, options=None, debug=False, dtype="d"):
         # Set all internal adflow default options before we set anything from python
         self.adflow.inputparamroutines.setdefaultvalues()
 
+        # surface-roughness overwrite-dict
+        self.surfaceRoughness = dict()
+
         defSetupTime = time.time()
 
         # Initialize the inherited AeroSolver
@@ -1162,6 +1165,50 @@ def setRotationRate(self, rotCenter, rotRate, cgnsBlocks=None):
         self.adflow.preprocessingapi.updaterotationrate(rotCenter, rotations, cgnsBlocks)
         self._updateVelInfo = True
 
+    def setSurfaceRoughness(self, ks, groupName=None):
+        """
+        Sets the equivalent sandgrain roughness on the specified surface family.
+        It may only be used when the rough SA-variant (useRoughSA = True) is
+        active.
+
+        This function takes precendence over all boundary conditions set in the
+        mesh.
+
+        Parameters:
+        ----------
+        ks : float
+            Equivalent sandgrain roughness
+        groupName : str
+            Family group to use. Default to all walls if not given (None)
+        """
+
+        if not self.options["useroughsa"]:
+            raise Error(
+                "It is not possible to set a surface roughness value without "
+                "using the rough SA-variant (useRoughSA = False)"
+            )
+
+        if groupName is None:
+            groupName = self.allWallsGroup
+
+        if groupName in self.surfaceRoughness:
+            raise Error(
+                f"The roughness value for surface group '{groupName}' has "
+                f"allready been set to {self.surfaceRoughness[groupName]}"
+            )
+
+        self.surfaceRoughness[groupName] = ks
+
+    def _setSurfaceRoughness(self):
+        """
+        Set the actual roughness values. It must be postponed for as long as
+        possible. Otherwise it might get overwritten.
+        """
+
+        for groupName, ks in self.surfaceRoughness.items():
+            famList = self._getFamilyList(groupName)
+            self.adflow.preprocessingapi.updatesurfaceroughness(ks, famList)
+
     def checkPartitioning(self, nprocs):
         """This function determines the potential load balancing for
         nprocs. The intent is this function can be run in serial
@@ -3706,6 +3753,13 @@ def _setAeroProblemData(self, aeroProblem, firstCall=False):
                 self.adflow.preprocessingapi.updatemetricsalllevels()
             self.adflow.preprocessingapi.updategridvelocitiesalllevels()
 
+            # overwrite BC surface roughness
+            if len(self.surfaceRoughness) > 0:
+                self._setSurfaceRoughness()
+
+            # Propagate roughness values through volume
+            self.adflow.walldistance.updatewallroughness()
+
     def _getBCDataFromAeroProblem(self, AP):
         variables = []
         dataArray = []
@@ -5675,6 +5729,7 @@ def _getDefaultOptions():
             "useQCR": [bool, False],
             "useRotationSA": [bool, False],
             "useft2SA": [bool, True],
+            "useRoughSA": [bool, False],
             "eddyVisInfRatio": [float, 0.009],
             "useWallFunctions": [bool, False],
             "useApproxWallDistance": [bool, True],
@@ -6070,6 +6125,7 @@ def _getOptionMap(self):
             "useqcr": ["physics", "useqcr"],
             "userotationsa": ["physics", "userotationsa"],
             "useft2sa": ["physics", "useft2sa"],
+            "useroughsa": ["physics", "useroughsa"],
             "eddyvisinfratio": ["physics", "eddyvisinfratio"],
             "usewallfunctions": ["physics", "wallfunctions"],
             "walldistcutoff": ["physics", "walldistcutoff"],

doc/options.yaml

@@ -237,6 +237,15 @@ useft2SA:
   desc: >
     Include the ft2 term in the Spalart-Allmaras turbulence model.
 
+useRoughSA:
+  desc: >
+    Activates the rough SA variant described `here
+    <https://turbmodels.larc.nasa.gov/spalart.html#sarough>`__. To use it,
+    blockettes must be disabled ``useBlockettes: False``. The roughness value on
+    the boundary may be prescribed like a regular boundary value using the name
+    ``SandGrainRoughness``. Additionally, the ADflow method ``setSurfaceRoughness``
+    may be used.
+
 eddyVisInfRatio:
   desc: >
     Free stream value of eddy viscosity.
@@ -1613,4 +1622,3 @@ meshMaxSkewness:
 useSkewnessCheck:
   desc: >
     When set to true, ADflow computes the `skewness` of each cell and throws an error if it is above ``meshMaxSkewness``. See also ``printBadlySkewedCells``.
-

src/NKSolver/blockette.F90

@@ -86,7 +86,7 @@ subroutine blocketteRes(useDissApprox, useViscApprox, useUpdateIntermed, useFlow
         use initializeFlow, only: referenceState
         use section, only: sections, nSections
         use iteration, only: rFil, currentLevel
-        use haloExchange, only: exchangeCoor, whalo2
+        use haloExchange, only: exchangeCoor, whalo2, exchanged2Wall
         use wallDistance, only: updateWallDistancesQuickly
         use utils, only: setPointers, EChk
         use turbUtils, only: computeEddyViscosity
@@ -181,12 +181,20 @@ subroutine blocketteRes(useDissApprox, useViscApprox, useUpdateIntermed, useFlow
                 do nn = 1, nDom
                     call setPointers(nn, currentLevel, sps)
                     call xhalo_block()
+
+                    if (equations == RANSEquations .and. useApproxWallDistance) then
+                        call updateWallDistancesQuickly(nn, 1, sps)
+                    end if
                 end do
             end do
 
             ! Now exchange the coordinates (fine level only)
             call exchangecoor(1)
 
+            if (equations == RANSEquations .and. useApproxWallDistance) then
+                call exchanged2Wall(1)
+            end if
+
             do sps = 1, nTimeIntervalsSpectral
                 ! Update overset connectivity if necessary
                 if (oversetPresent .and. oversetUpdateMode == updateFast) then

src/adjoint/Makefile_tapenade

@@ -132,17 +132,17 @@ flowUtils%computeEtotBlock(w, p) > \
 \
 turbUtils%computeEddyViscosity(w,rlv)>(w,rlv,rev) \
 \
-BCRoutines%bcSymm1stHalo(ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2, bcData%norm) > \
-                        (ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2, bcData%norm) \
+BCRoutines%bcSymm1stHalo(ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2, d2wall1, d2wall2, bcData%norm) > \
+                        (ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2, d2wall1, d2wall2, bcData%norm) \
 \
-BCRoutines%bcSymm2ndHalo(ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, bcData%norm) > \
-                        (ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, bcData%norm) \
+BCRoutines%bcSymm2ndHalo(ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, bd2wall0, d2wall3, cData%norm) > \
+                        (ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, bd2wall0, d2wall3, bcData%norm) \
 \
-BCRoutines%bcSymmPolar1stHalo(xx, ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2) > \
-                             (xx, ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2) \
+BCRoutines%bcSymmPolar1stHalo(xx, ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2, d2wall1, d2wall2) > \
+                             (xx, ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2, d2wall1, d2wall2) \
 \
-BCRoutines%bcSymmPolar2ndHalo(xx, ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3) > \
-                             (xx, ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3) \
+BCRoutines%bcSymmPolar2ndHalo(xx, ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, d2wall0, d2wall3) > \
+                             (xx, ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3, d2wall0, d2wall3) \
 \
 BCRoutines%bcNSWallAdiabatic(ww0, ww1, ww2, pp0, pp1, pp2, pp3, rlv0, rlv1, rlv2, rev0, rev1, rev2, bcData%uSlip) > \
                             (ww0, ww1, ww2, pp0, pp1, pp2, pp3, rlv0, rlv1, rlv2, rev0, rev1, rev2, bcData%uSlip) \
@@ -162,11 +162,11 @@ BCRoutines%bcSubsonicInflow(ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0
 BCRoutines%bcSubsonicOutflow(ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, bcData%norm, bcData%Ps) > \
                             (ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, bcData%norm, bcData%Ps) \
 \
-turbBCRoutines%applyAllTurbBCThisBlock(rev, w, bvtj1, bvtj2, bvtk1, bvtk2, bvti1, bvti2) > \
-                                      (rev, w) \
+turbBCRoutines%applyAllTurbBCThisBlock(rev, w, bvtj1, bvtj2, bvtk1, bvtk2, bvti1, bvti2, bmti1, bmti2, bmtj1, bmtj2, bmtk1, bmtk2, d2wall) > \
+                                      (rev, w, bvtj1, bvtj2, bvtk1, bvtk2, bvti1, bvti2, bmti1, bmti2, bmtj1, bmtj2, bmtk1, bmtk2, d2wall) \
 \
-turbBCRoutines%bcTurbTreatment(w, rlv, d2wall, winf) > \
-                              (w, rlv, d2wall, winf, bvtj1, bvtj2, bvtk1, bvtk2, bvti1, bvti2) \
+turbBCRoutines%bcTurbTreatment(w, rlv, d2wall, winf, bvtj1, bvtj2, bvtk1, bvtk2, bvti1, bvti2, bmti1, bmti2, bmtj1, bmtj2, bmtk1, bmtk2) > \
+                              (w, rlv, d2wall, winf, bvtj1, bvtj2, bvtk1, bvtk2, bvti1, bvti2, bmti1, bmti2, bmtj1, bmtj2, bmtk1, bmtk2) \
 \
 solverUtils%timeStep_block(w, pInfCorr, rhoInf, si, sj, sk, sFaceI, sFaceJ, sFaceK, p, radi, radj, radk, dtl, rlv, rev, vol) > \
               (w, pInfCorr, rhoInf, si, sj, sk, sFaceI, sFaceJ, sFaceK, p, radi, radj, radk, dtl, rlv, rev, vol) \

src/adjoint/adjointUtils.F90

@@ -772,8 +772,10 @@ subroutine allocDerivativeValues(level)
 
                 ! Allocate d2wall if not already done so
                 if (.not. associated(flowDoms(nn, 1, sps)%d2wall)) then
-                    allocate (flowDoms(nn, 1, sps)%d2wall(2:il, 2:jl, 2:kl))
+                    allocate (flowDoms(nn, 1, sps)%d2wall(0:ib, 0:jb, 0:kb))
                     call EChk(ierr, __FILE__, __LINE__)
+
+                    flowDoms(nn, level, sps)%d2Wall = 0.01
                 end if
 
                 ! Now allocate all valus that have a differentiable
@@ -829,7 +831,7 @@ subroutine allocDerivativeValues(level)
                     flowDomsd(nn, level, sps)%bvtj2(ie, ke, nt1:nt2), &
                     flowDomsd(nn, level, sps)%bvtk1(ie, je, nt1:nt2), &
                     flowDomsd(nn, level, sps)%bvtk2(ie, je, nt1:nt2), &
-                    flowDomsd(nn, level, sps)%d2Wall(2:il, 2:jl, 2:kl), &
+                    flowDomsd(nn, level, sps)%d2Wall(0:ib, 0:jb, 0:kb), &
                     stat=ierr)
                 call EChk(ierr, __FILE__, __LINE__)
 

src/adjoint/outputForward/BCData_d.f90

@@ -13,11 +13,28 @@ module bcdata_d
   subroutine setbcvarnamesisothermalwall()
     use cgnsnames
     use constants
+    use inputphysics, only : useroughsa
     implicit none
     nbcvar = nbcvarisothermalwall
     bcvarnames(1) = cgnstemp
+    if (useroughsa) then
+      nbcvar = nbcvar + 1
+      bcvarnames(2) = cgnssandgrainroughness
+    end if
   end subroutine setbcvarnamesisothermalwall
 
+  subroutine setbcvarnamesadiabaticwall()
+    use cgnsnames
+    use constants
+    use inputphysics, only : useroughsa
+    implicit none
+    nbcvar = nbcvaradiabaticwall
+    if (useroughsa) then
+      nbcvar = nbcvar + 1
+      bcvarnames(1) = cgnssandgrainroughness
+    end if
+  end subroutine setbcvarnamesadiabaticwall
+
   subroutine setbcvarnamessubsonicinflow()
     use constants
     use cgnsnames
@@ -329,10 +346,10 @@ subroutine unitvectorscylsystem(boco)
   end subroutine unitvectorscylsystem
 
 !  differentiation of bcdataisothermalwall in forward (tangent) mode (with options i4 dr8 r8):
-!   variations   of useful results: *(*bcdata.tns_wall)
+!   variations   of useful results: *(*bcdata.tns_wall) bcvararray
 !   with respect to varying inputs: tref bcvararray
 !   rw status of diff variables: tref:in *(*bcdata.tns_wall):out
-!                bcvararray:in
+!                bcvararray:in-out
 !   plus diff mem management of: bcdata:in *bcdata.tns_wall:in
 ! ---------------------------------------------------------------
 ! routines that set the actual bcdata values from the cgns data set
@@ -350,6 +367,7 @@ subroutine bcdataisothermalwall_d(boco, bcvararray, bcvararrayd, ibeg&
     use blockpointers, only : bcfaceid, bcdata, bcdatad, nbkglobal
     use utils_d, only : terminate, sitemperature
     use flowvarrefstate, only : tref, trefd
+    use inputphysics, only : useroughsa
     use diffsizes
 !  hint: isize1ofdrfbcdata should be the size of dimension 1 of array *bcdata
     implicit none
@@ -388,6 +406,21 @@ subroutine bcdataisothermalwall_d(boco, bcvararray, bcvararrayd, ibeg&
         bcdata(boco)%tns_wall(i, j) = temp0
       end do
     end do
+! set a value of 0 if it was not possible to determine the
+! sand grain roughness
+    if (useroughsa) then
+! set a value of 0 if it was not possible to determine the
+! sand grain roughness
+      if (.not.bcvarpresent(1)) then
+        bcvararrayd(:, :, 2) = 0.0_8
+        bcvararray(:, :, 2) = zero
+      end if
+      do j=jbeg,jend
+        do i=ibeg,iend
+          bcdata(boco)%ksns_wall(i, j) = bcvararray(i, j, 2)
+        end do
+      end do
+    end if
   end subroutine bcdataisothermalwall_d
 
 ! ---------------------------------------------------------------
@@ -406,6 +439,7 @@ subroutine bcdataisothermalwall(boco, bcvararray, ibeg, iend, jbeg, &
     use blockpointers, only : bcfaceid, bcdata, nbkglobal
     use utils_d, only : terminate, sitemperature
     use flowvarrefstate, only : tref
+    use inputphysics, only : useroughsa
     implicit none
 !
 !      subroutine arguments.
@@ -431,8 +465,54 @@ subroutine bcdataisothermalwall(boco, bcvararray, ibeg, iend, jbeg, &
 &         tref
       end do
     end do
+! set a value of 0 if it was not possible to determine the
+! sand grain roughness
+    if (useroughsa) then
+! set a value of 0 if it was not possible to determine the
+! sand grain roughness
+      if (.not.bcvarpresent(1)) bcvararray(:, :, 2) = zero
+      do j=jbeg,jend
+        do i=ibeg,iend
+          bcdata(boco)%ksns_wall(i, j) = bcvararray(i, j, 2)
+        end do
+      end do
+    end if
   end subroutine bcdataisothermalwall
 
+  subroutine bcdataadiabaticwall(boco, bcvararray, ibeg, iend, jbeg, &
+&   jend)
+!
+!       tries to extract the equivalent sand grain roughness. it sets
+!       a default value of 0.0
+!
+    use constants
+    use cgnsnames
+    use inputphysics, only : useroughsa
+    use blockpointers, only : bcfaceid, bcdata, nbkglobal
+    implicit none
+!
+!      subroutine arguments.
+!
+    integer(kind=inttype) :: boco
+    integer(kind=inttype) :: ibeg, iend, jbeg, jend
+    real(kind=realtype), dimension(ibeg:iend, jbeg:jend, nbcvarmax) :: &
+&   bcvararray
+!
+!      local variables.
+!
+    integer(kind=inttype) :: i, j
+    if (useroughsa) then
+! set a value of 0 if it was not possible to determine the
+! sand grain roughness
+      if (.not.bcvarpresent(1)) bcvararray(:, :, 1) = zero
+      do j=jbeg,jend
+        do i=ibeg,iend
+          bcdata(boco)%ksns_wall(i, j) = bcvararray(i, j, 1)
+        end do
+      end do
+    end if
+  end subroutine bcdataadiabaticwall
+
 !  differentiation of bcdatasubsonicinflow in forward (tangent) mode (with options i4 dr8 r8):
 !   variations   of useful results: *(*bcdata.ptinlet) *(*bcdata.ttinlet)
 !                *(*bcdata.htinlet) *(*bcdata.turbinlet)