Slice mods

adflow/airfoilClass.py

@@ -642,6 +642,31 @@ def plotAirfoil(self):
         fig.update_scenes(aspectmode="data")
         fig.show()
 
+    def addFunction(self, funcName, callBack):
+
+    def getNodeFiltering(self, begFrac, endFrac, surface="upper"):
+
+        return mask_array
+
+    def evalFuncs(self, )
+
+
+
+
+if __name__ == "__main__":
+    mask_array = airfoil_Instance.getNodeFiltering()
+
+    def callback(nodeData, func, mask_array=mask_array):
+        cp = nodeData["cp"]
+        x = nodeData["x"]
+
+        cp_local = cp[mask_array]
+        x_local = x[mask_array]
+
+    airfoil_instance.addFunc("fname", callback, active_data=["x", "cp"])
+
+
+
 
 class Wing:
     def __init__(self, list_data, list_conn, list_normal, list_points):

adflow/pyADflow.py

@@ -177,6 +177,9 @@ def __init__(self, comm=None, options=None, debug=False, dtype="d"):
         self.updateTime = 0.0
         self.nSlice = 0
         self.nLiftDist = 0
+        self.nActiveSlice = 0
+        # dictionary to hold the active slice name mapping to the fortran level slice data
+        self.activeSliceDict = OrderedDict()
 
         # Set default values
         self.adflow.inputio.autoparameterupdate = False
@@ -802,6 +805,85 @@ def addCylindricalSlices(
 
         self.nSlice += nSlice
 
+    def addActiveSlice(self, normal, point, sliceType="relative", groupName=None, sliceDir=None, name=None):
+        """
+        Add active slices that can be used during optimization
+
+        Parameters
+        ----------
+        normal : 3-d array
+            The normals of the slice directions. If an array of size 3 is passed,
+            we use the same normal for all slices. if an array of multiple normals
+            are passed, we use the individual normals for each point. in this case,
+            the numbers of points and normals must match.
+        point : 3-d array
+            Point coordinates that define a slicing plane along with the normals
+        sliceDir : 3d array
+            If this is provided, only the intersections that is in this direction
+            starting from the normal is kept. This is useful if you are slicing a
+            closed geometry and only want to keep one side of it. It is similar to
+            the functionality provided in :meth:`addCylindricalSlices <.addCylindricalSlices>`.
+        sliceType : str {'relative', 'absolute'}
+            Relative slices are 'sliced' at the beginning and then parametrically
+            move as the geometry deforms. As a result, the slice through the
+            geometry may not remain planar. An absolute slice is re-sliced for
+            every output so is always exactly planar and always at the initial
+            position the user indicated.
+        groupName : str
+             The family to use for the slices. Default is None corresponding to all
+             wall groups.
+        """
+
+        # call the preprocessing routine that avoids some code duplication across 3 types of slices
+        sliceType, famList = self._preprocessSliceInput(sliceType, groupName)
+
+        # check if a custom name is provided, if not, we name the slices sequentially
+        if name is None:
+            sliceName = f"active_slice_{self.nActiveSlice:03d}"
+        else:
+            # make sure this slice name is not already taken
+            if name in self.activeSliceDict:
+                raise Error(f"Active slice name {name} is already taken. Use another name for this slice.")
+            else:
+                sliceName = name
+
+        if sliceDir is None:
+            # if we dont have a direction vector to pick a projection direction, we can just set this
+            # array to an arbitrary direction. Because the useDir flag is false, the code won't actually
+            # use this array
+            dummySliceDir = [1.0, 0.0, 0.0]
+            useDir = False
+        else:
+            useDir = True
+
+        if useDir:
+            direction = sliceDir
+        else:
+            direction = dummySliceDir
+
+        sliceTitle = (
+            f"Slice_{self.nSlice + 1:04d} {sliceName} {groupName} {sliceType.capitalize()} Arbitrary "
+            f"Point = ({point[0]:.8f}, {point[1]:.8f}, {point[2]:.8f}) "
+            f"Normal = ({normal[0]:.8f}, {normal[1]:.8f}, {normal[2]:.8f})"
+        )
+        if sliceType == "relative":
+            sliceIndex = self.adflow.tecplotio.addparaslice(sliceTitle, point, normal, direction, useDir, famList)
+        else:
+            sliceIndex = self.adflow.tecplotio.addabsslice(sliceTitle, point, normal, direction, useDir, famList)
+
+        # save the fortran level data. specifically, we need the slice index and if this slice is parametric or absolute. Also, save the point and normal for convenience
+        self.activeSliceDict[sliceName] = {
+            "index": sliceIndex,
+            "type": sliceType,
+            "point": point,
+            "normal": normal,
+            "groupName": groupName,
+            "sliceDir": sliceDir,
+        }
+
+        self.nActiveSlice += 1
+        self.nSlice += 1
+
     def _preprocessSliceInput(self, sliceType, groupName):
         """
         Preprocessing routine that holds some of the duplicated code required for 3 types of slice methods.
@@ -821,6 +903,75 @@ def _preprocessSliceInput(self, sliceType, groupName):
 
         return sliceType, famList
 
+    def getActiveSliceData(self):
+        # returns a dictionary where the keys are the active slice names.
+        # values of the top level dictionary are dictionaries that hold the data
+        # for the slice itself. "conn" is the 0-based connectivity array
+        # each nodel data type is provided with its own key, i.e. keying CoordinateX returns
+        # the array of x coordinates of the nodes
+
+        # get the list of solution names
+        # TODO this can be done once and saved
+        # TODO the f2py wrapper here can be greatly improved...
+        nSolVar = self.adflow.outputmod.numberofsurfsolvariables()
+        surfSolNames = self.adflow.outputmod.surfsolnameswrap(nSolVar)
+        # convert to numpy array for easier processing
+        surfSolNames = numpy.array(surfSolNames).T
+
+        # save the values. We initialize the list with the coordinates, and the 6 tractions
+        surfSolNameList = [
+            "CoordinateX",
+            "CoordinateY",
+            "CoordinateZ",
+            "PressureTractionX",
+            "PressureTractionY",
+            "PressureTractionZ",
+            "ViscousTractionX",
+            "ViscousTractionY",
+            "ViscousTractionZ",
+        ]
+        for ii in range(nSolVar):
+            surfSolNameList.append(surfSolNames[ii, :].tostring().decode().strip())
+
+        # compute nodal values
+        # famList = self._getFamilyList(self.allWallsGroup)
+        self.adflow.tecplotio.computetempnodaldata()
+
+        # loop over each slice and call the fortran routine
+        sliceDataDict = OrderedDict()
+
+        for sliceName, sliceInfo in self.activeSliceDict.items():
+            sliceIndex = sliceInfo["index"]
+            sliceType = sliceInfo["type"]
+
+            # TODO also account for multiple time spectral instances
+            if sliceType == "relative":
+                nNode, nElem, nData = self.adflow.tecplotio.getparaslicedatasize(sliceIndex)
+
+                # allocate the numpy array to hold the data
+                sliceNodalData = numpy.zeros((nNode, nData), self.dtype)
+                sliceConn = numpy.zeros((nElem, 2), numpy.intc)
+
+                # the fortran layer broadcasts the global data. each proc gets the same data back in python
+                self.adflow.tecplotio.getparaslicedata(sliceIndex, sliceNodalData.T, sliceConn.T)
+
+            # convert conn to zero based
+            sliceConn -= 1
+
+            # initialize the data in the dict with conn array
+            sliceDataDict[sliceName] = {
+                "conn": sliceConn,
+            }
+
+            # loop over the nodal data and save one by one using the keys
+            for ii, dataName in enumerate(surfSolNameList):
+                sliceDataDict[sliceName][dataName] = sliceNodalData[:, ii]
+
+        # deallocate nodal values
+        self.adflow.tecplotio.deallocatetempnodaldata()
+
+        return sliceDataDict
+
     def addIntegrationSurface(self, fileName, familyName, isInflow=True, coordXfer=None):
         """Add a specific integration surface for performing massflow-like
         computations.
@@ -2780,7 +2931,7 @@ def writeSolution(self, outputDir=None, baseName=None, number=None, writeSlices=
         # Flag to write the tecplot surface solution or not
         writeSurf = self.getOption("writeTecplotSurfaceSolution")
 
-        # # Call fully compbined fortran routine.
+        # Call fully compbined fortran routine.
         self.adflow.tecplotio.writetecplot(sliceName, writeSlices, liftName, writeLift, surfName, writeSurf, famList)
 
     def writeMeshFile(self, fileName):

src/ADT/adtLocalSearch.F90

@@ -2327,11 +2327,14 @@ subroutine hessD2Hexa(xP, x1, x2, x3, x4, x5, x6, x7, x8, chi, hess, iErr)
         hess(3, 3) = 2.0_realType * ((dxdzeta * dxdzeta) + (dydzeta * dydzeta) + (dzdzeta * dzdzeta))
 
         hess(1, 2) = 2.0_realType * ((dxdksi * dxdeta) + (dydksi * dydeta) + (dzdksi * dzdeta) &
-                           - ((xP(1) - x0) * d2xdksideta) - ((xP(2) - y0) * d2ydksideta) - ((xP(3) - z0) * d2zdksideta))
+                                     - ((xP(1) - x0) * d2xdksideta) - &
+                                     ((xP(2) - y0) * d2ydksideta) - ((xP(3) - z0) * d2zdksideta))
         hess(1, 3) = 2.0_realType * ((dxdksi * dxdzeta) + (dydksi * dydzeta) + (dzdksi * dzdzeta) &
-                        - ((xP(1) - x0) * d2xdksidzeta) - ((xP(2) - y0) * d2ydksidzeta) - ((xP(3) - z0) * d2zdksidzeta))
+                                     - ((xP(1) - x0) * d2xdksidzeta) - &
+                                     ((xP(2) - y0) * d2ydksidzeta) - ((xP(3) - z0) * d2zdksidzeta))
         hess(2, 3) = 2.0_realType * ((dxdeta * dxdzeta) + (dydeta * dydzeta) + (dzdeta * dzdzeta) &
-                        - ((xP(1) - x0) * d2xdetadzeta) - ((xP(2) - y0) * d2ydetadzeta) - ((xP(3) - z0) * d2zdetadzeta))
+                                     - ((xP(1) - x0) * d2xdetadzeta) - &
+                                     ((xP(2) - y0) * d2ydetadzeta) - ((xP(3) - z0) * d2zdetadzeta))
 
         hess(2, 1) = hess(1, 2)
         hess(3, 1) = hess(1, 3)

src/NKSolver/NKSolvers.F90

@@ -3136,7 +3136,8 @@ subroutine physicalityCheckANK(lambdaP)
 #ifndef USE_COMPLEX
                                 ratio = (wvec_pointer(ii) / (dvec_pointer(ii) + eps)) * ANK_physLSTolTurb
 #else
-                               ratio = (real(wvec_pointer(ii)) / real(dvec_pointer(ii) + eps)) * real(ANK_physLSTolTurb)
+                                ratio = (real(wvec_pointer(ii)) &
+                                         / real(dvec_pointer(ii) + eps)) * real(ANK_physLSTolTurb)
 #endif
                                 ! if the ratio is less than min step, the update is either
                                 ! in the positive direction, therefore we do not clip it,

src/NKSolver/blockette.F90

@@ -3158,7 +3158,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw2 = w(i + 1, j, k, ivx) * w(i + 1, j, k, irho) - w(i, j, k, ivx) * w(i, j, k, irho)
                     fs = dis2 * ddw2 &
-       - dis4 * (w(i + 2, j, k, ivx) * w(i + 2, j, k, irho) - w(i - 1, j, k, ivx) * w(i - 1, j, k, irho) - three * ddw2)
+                         - dis4 * (w(i + 2, j, k, ivx) * w(i + 2, j, k, irho) - &
+                                   w(i - 1, j, k, ivx) * w(i - 1, j, k, irho) - three * ddw2)
 
                     fw(i + 1, j, k, imx) = fw(i + 1, j, k, imx) + fs
                     fw(i, j, k, imx) = fw(i, j, k, imx) - fs
@@ -3167,7 +3168,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw3 = w(i + 1, j, k, ivy) * w(i + 1, j, k, irho) - w(i, j, k, ivy) * w(i, j, k, irho)
                     fs = dis2 * ddw3 &
-       - dis4 * (w(i + 2, j, k, ivy) * w(i + 2, j, k, irho) - w(i - 1, j, k, ivy) * w(i - 1, j, k, irho) - three * ddw3)
+                         - dis4 * (w(i + 2, j, k, ivy) * w(i + 2, j, k, irho) - &
+                                   w(i - 1, j, k, ivy) * w(i - 1, j, k, irho) - three * ddw3)
 
                     fw(i + 1, j, k, imy) = fw(i + 1, j, k, imy) + fs
                     fw(i, j, k, imy) = fw(i, j, k, imy) - fs
@@ -3176,7 +3178,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw4 = w(i + 1, j, k, ivz) * w(i + 1, j, k, irho) - w(i, j, k, ivz) * w(i, j, k, irho)
                     fs = dis2 * ddw4 &
-       - dis4 * (w(i + 2, j, k, ivz) * w(i + 2, j, k, irho) - w(i - 1, j, k, ivz) * w(i - 1, j, k, irho) - three * ddw4)
+                         - dis4 * (w(i + 2, j, k, ivz) * w(i + 2, j, k, irho) - &
+                                   w(i - 1, j, k, ivz) * w(i - 1, j, k, irho) - three * ddw4)
 
                     fw(i + 1, j, k, imz) = fw(i + 1, j, k, imz) + fs
                     fw(i, j, k, imz) = fw(i, j, k, imz) - fs
@@ -3185,7 +3188,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw5 = (w(i + 1, j, k, irhoE) + P(i + 1, j, K)) - (w(i, j, k, irhoE) + P(i, j, k))
                     fs = dis2 * ddw5 &
-           - dis4 * ((w(i + 2, j, k, irhoE) + P(i + 2, j, k)) - (w(i - 1, j, k, irhoE) + P(i - 1, j, k)) - three * ddw5)
+                         - dis4 * ((w(i + 2, j, k, irhoE) + P(i + 2, j, k)) - &
+                                   (w(i - 1, j, k, irhoE) + P(i - 1, j, k)) - three * ddw5)
 
                     fw(i + 1, j, k, irhoE) = fw(i + 1, j, k, irhoE) + fs
                     fw(i, j, k, irhoE) = fw(i, j, k, irhoE) - fs
@@ -3223,7 +3227,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw2 = w(i, j + 1, k, ivx) * w(i, j + 1, k, irho) - w(i, j, k, ivx) * w(i, j, k, irho)
                     fs = dis2 * ddw2 &
-       - dis4 * (w(i, j + 2, k, ivx) * w(i, j + 2, k, irho) - w(i, j - 1, k, ivx) * w(i, j - 1, k, irho) - three * ddw2)
+                         - dis4 * (w(i, j + 2, k, ivx) * w(i, j + 2, k, irho) - &
+                                   w(i, j - 1, k, ivx) * w(i, j - 1, k, irho) - three * ddw2)
 
                     fw(i, j + 1, k, imx) = fw(i, j + 1, k, imx) + fs
                     fw(i, j, k, imx) = fw(i, j, k, imx) - fs
@@ -3232,7 +3237,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw3 = w(i, j + 1, k, ivy) * w(i, j + 1, k, irho) - w(i, j, k, ivy) * w(i, j, k, irho)
                     fs = dis2 * ddw3 &
-       - dis4 * (w(i, j + 2, k, ivy) * w(i, j + 2, k, irho) - w(i, j - 1, k, ivy) * w(i, j - 1, k, irho) - three * ddw3)
+                         - dis4 * (w(i, j + 2, k, ivy) * w(i, j + 2, k, irho) - &
+                                   w(i, j - 1, k, ivy) * w(i, j - 1, k, irho) - three * ddw3)
 
                     fw(i, j + 1, k, imy) = fw(i, j + 1, k, imy) + fs
                     fw(i, j, k, imy) = fw(i, j, k, imy) - fs
@@ -3241,7 +3247,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw4 = w(i, j + 1, k, ivz) * w(i, j + 1, k, irho) - w(i, j, k, ivz) * w(i, j, k, irho)
                     fs = dis2 * ddw4 &
-       - dis4 * (w(i, j + 2, k, ivz) * w(i, j + 2, k, irho) - w(i, j - 1, k, ivz) * w(i, j - 1, k, irho) - three * ddw4)
+                         - dis4 * (w(i, j + 2, k, ivz) * w(i, j + 2, k, irho) - &
+                                   w(i, j - 1, k, ivz) * w(i, j - 1, k, irho) - three * ddw4)
 
                     fw(i, j + 1, k, imz) = fw(i, j + 1, k, imz) + fs
                     fw(i, j, k, imz) = fw(i, j, k, imz) - fs
@@ -3250,7 +3257,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw5 = (w(i, j + 1, k, irhoE) + P(i, j + 1, k)) - (w(i, j, k, irhoE) + P(i, j, k))
                     fs = dis2 * ddw5 &
-           - dis4 * ((w(i, j + 2, k, irhoE) + P(i, j + 2, k)) - (w(i, j - 1, k, irhoE) + P(i, j - 1, k)) - three * ddw5)
+                         - dis4 * ((w(i, j + 2, k, irhoE) + P(i, j + 2, k)) - &
+                                   (w(i, j - 1, k, irhoE) + P(i, j - 1, k)) - three * ddw5)
 
                     fw(i, j + 1, k, irhoE) = fw(i, j + 1, k, irhoE) + fs
                     fw(i, j, k, irhoE) = fw(i, j, k, irhoE) - fs
@@ -3288,7 +3296,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw2 = w(i, j, k + 1, ivx) * w(i, j, k + 1, irho) - w(i, j, k, ivx) * w(i, j, k, irho)
                     fs = dis2 * ddw2 &
-       - dis4 * (w(i, j, k + 2, ivx) * w(i, j, k + 2, irho) - w(i, j, k - 1, ivx) * w(i, j, k - 1, irho) - three * ddw2)
+                         - dis4 * (w(i, j, k + 2, ivx) * w(i, j, k + 2, irho) - &
+                                   w(i, j, k - 1, ivx) * w(i, j, k - 1, irho) - three * ddw2)
 
                     fw(i, j, k + 1, imx) = fw(i, j, k + 1, imx) + fs
                     fw(i, j, k, imx) = fw(i, j, k, imx) - fs
@@ -3297,7 +3306,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw3 = w(i, j, k + 1, ivy) * w(i, j, k + 1, irho) - w(i, j, k, ivy) * w(i, j, k, irho)
                     fs = dis2 * ddw3 &
-       - dis4 * (w(i, j, k + 2, ivy) * w(i, j, k + 2, irho) - w(i, j, k - 1, ivy) * w(i, j, k - 1, irho) - three * ddw3)
+                         - dis4 * (w(i, j, k + 2, ivy) * w(i, j, k + 2, irho) - &
+                                   w(i, j, k - 1, ivy) * w(i, j, k - 1, irho) - three * ddw3)
 
                     fw(i, j, k + 1, imy) = fw(i, j, k + 1, imy) + fs
                     fw(i, j, k, imy) = fw(i, j, k, imy) - fs
@@ -3306,7 +3316,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw4 = w(i, j, k + 1, ivz) * w(i, j, k + 1, irho) - w(i, j, k, ivz) * w(i, j, k, irho)
                     fs = dis2 * ddw4 &
-       - dis4 * (w(i, j, k + 2, ivz) * w(i, j, k + 2, irho) - w(i, j, k - 1, ivz) * w(i, j, k - 1, irho) - three * ddw4)
+                         - dis4 * (w(i, j, k + 2, ivz) * w(i, j, k + 2, irho) - &
+                                   w(i, j, k - 1, ivz) * w(i, j, k - 1, irho) - three * ddw4)
 
                     fw(i, j, k + 1, imz) = fw(i, j, k + 1, imz) + fs
                     fw(i, j, k, imz) = fw(i, j, k, imz) - fs
@@ -3315,7 +3326,8 @@ subroutine inviscidDissFluxScalar
 
                     ddw5 = (w(i, j, k + 1, irhoE) + P(i, j, k + 1)) - (w(i, j, k, irhoE) + P(i, j, k))
                     fs = dis2 * ddw5 &
-           - dis4 * ((w(i, j, k + 2, irhoE) + P(i, j, k + 2)) - (w(i, j, k - 1, irhoE) + P(i, j, k - 1)) - three * ddw5)
+                         - dis4 * ((w(i, j, k + 2, irhoE) + P(i, j, k + 2)) - &
+                                   (w(i, j, k - 1, irhoE) + P(i, j, k - 1)) - three * ddw5)
 
                     fw(i, j, k + 1, irhoE) = fw(i, j, k + 1, irhoE) + fs
                     fw(i, j, k, irhoE) = fw(i, j, k, irhoE) - fs

src/adjoint/adjointAPI.F90

@@ -9,7 +9,8 @@ module adjointAPI
 contains
 #ifndef USE_COMPLEX
     subroutine computeMatrixFreeProductFwd(xvdot, extradot, wdot, bcDataValuesdot, useSpatial, &
-                    useState, famLists, bcDataNames, bcDataValues, bcDataFamLists, bcVarsEmpty, dwdot, funcsDot, fDot, &
+                                           useState, famLists, bcDataNames, bcDataValues, &
+                                           bcDataFamLists, bcVarsEmpty, dwdot, funcsDot, fDot, &
                                            costSize, fSize, nTime)
 
         ! This is the main matrix-free forward mode computation
@@ -901,7 +902,8 @@ subroutine setupPETScKsp
             ! linear system also serves as the preconditioning matrix. This
             ! is only valid if useMatrixFree is flase.
             if (useMatrixfreedRdw) then
-            call terminate("setupPETScKSP", "useMatrixFreedRdW option cannot be true when the approxPC option is False")
+                call terminate("setupPETScKSP", &
+                               "useMatrixFreedRdW option cannot be true when the approxPC option is False")
             end if
             call KSPSetOperators(adjointKSP, dRdWt, dRdWT, ierr)
             call EChk(ierr, __FILE__, __LINE__)