diff --git a/adflow/airfoilClass.py b/adflow/airfoilClass.py index 0b23bc9d3..de19bd83a 100644 --- a/adflow/airfoilClass.py +++ b/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): diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 6abe7481a..d13c5908d 100644 --- a/adflow/pyADflow.py +++ b/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): diff --git a/src/ADT/adtLocalSearch.F90 b/src/ADT/adtLocalSearch.F90 index d9a633400..7e2fc8646 100644 --- a/src/ADT/adtLocalSearch.F90 +++ b/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) diff --git a/src/NKSolver/NKSolvers.F90 b/src/NKSolver/NKSolvers.F90 index 74941a1a2..f2a2a2f77 100644 --- a/src/NKSolver/NKSolvers.F90 +++ b/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, diff --git a/src/NKSolver/blockette.F90 b/src/NKSolver/blockette.F90 index 319e294ef..7841e4942 100644 --- a/src/NKSolver/blockette.F90 +++ b/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 diff --git a/src/adjoint/adjointAPI.F90 b/src/adjoint/adjointAPI.F90 index 307b901d2..8f9ad758b 100644 --- a/src/adjoint/adjointAPI.F90 +++ b/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__) diff --git a/src/adjoint/adjointDebug.F90 b/src/adjoint/adjointDebug.F90 index 90aca5b51..24119b4d8 100644 --- a/src/adjoint/adjointDebug.F90 +++ b/src/adjoint/adjointDebug.F90 @@ -504,7 +504,8 @@ subroutine printADSeeds(nn, level, sps) maxval(commPatternOverset(level, sps)%sendList(i)%interpd), & norm2(commPatternOverset(level, sps)%sendList(i)%interpd) end do sends - write (*, *) 'internalOverset(level, sps)%donorInterpd ', minval(internalOverset(level, sps)%donorInterpd), & + write (*, *) 'internalOverset(level, sps)%donorInterpd ', & + minval(internalOverset(level, sps)%donorInterpd), & maxval(internalOverset(level, sps)%donorInterpd), & norm2(internalOverset(level, sps)%donorInterpd) end if @@ -546,7 +547,8 @@ subroutine printADSeeds(nn, level, sps) if (associated(cgnsDoms(iDom)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays)) then do iDirichlet = 1, size(cgnsDoms(iDom)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays) write (*, *) iDom, iBoco, iData, iDirichlet, 'dataArr(:) ' & - , cgnsDomsd(iDom)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays(iDirichlet)%dataArr(:) + , cgnsDomsd(iDom)%bocoInfo(iBoco)%dataSet(iData)% & + dirichletArrays(iDirichlet)%dataArr(:) end do end if end do diff --git a/src/adjoint/adjointUtils.F90 b/src/adjoint/adjointUtils.F90 index 6bed0e709..68b26f8f7 100644 --- a/src/adjoint/adjointUtils.F90 +++ b/src/adjoint/adjointUtils.F90 @@ -331,7 +331,8 @@ subroutine setupStateResidualMatrix(matrix, useAD, usePC, useTranspose, & do k = 0, kb do j = 0, jb do i = 0, ib - flowDoms(nn, level, sps2)%w(i, j, k, ll) = flowDoms(nn, 1, sps2)%wtmp(i, j, k, ll) + flowDoms(nn, level, sps2)%w(i, j, k, ll) = & + flowDoms(nn, 1, sps2)%wtmp(i, j, k, ll) end do end do end do @@ -452,7 +453,8 @@ subroutine setupStateResidualMatrix(matrix, useAD, usePC, useTranspose, & if (buildCoarseMats) then do lvl = 1, agmgLevels - 1 - coarseCols(m, lvl + 1) = coarseOversetIndices(nn, lvl)%arr(m, i, j, k) + coarseCols(m, lvl + 1) = & + coarseOversetIndices(nn, lvl)%arr(m, i, j, k) end do end if end do @@ -488,28 +490,35 @@ subroutine setupStateResidualMatrix(matrix, useAD, usePC, useTranspose, & if (buildCoarseMats) then do lvl = 1, agmgLevels - 1 - coarseRows(lvl + 1) = coarseIndices(nn, lvl)%arr(i + ii, j + jj, k + kk) + coarseRows(lvl + 1) = & + coarseIndices(nn, lvl)%arr(i + ii, j + jj, k + kk) end do end if - rowBlank: if (flowDoms(nn, level, sps)%iBlank(i + ii, j + jj, k + kk) == 1) then + rowBlank: if (flowDoms(nn, level, sps)% & + iBlank(i + ii, j + jj, k + kk) == 1) then centerCell: if (ii == 0 .and. jj == 0 .and. kk == 0) then useDiagPC: if (usePC .and. useDiagTSPC) then ! If we're doing the PC and we want ! to use TS diagonal form, only set ! values for on-time insintance - blk = flowDoms(nn, 1, sps)%dw_deriv(i + ii, j + jj, k + kk, & - lStart:lEnd, lStart:lEnd) + blk = flowDoms(nn, 1, sps)%dw_deriv(i + ii, & + j + jj, k + kk, & + lStart:lEnd, & + lStart:lEnd) call setBlock(blk) else ! Otherwise loop over spectral ! instances and set all. do sps2 = 1, nTimeIntervalsSpectral irow = flowDoms(nn, level, sps2)% & - globalCell(i + ii, j + jj, k + kk) - blk = flowDoms(nn, 1, sps2)%dw_deriv(i + ii, j + jj, k + kk, & - lStart:lEnd, lStart:lEnd) + globalCell(i + ii, & + j + jj, k + kk) + blk = flowDoms(nn, 1, sps2)% & + dw_deriv(i + ii, & + j + jj, k + kk, & + lStart:lEnd, lStart:lEnd) call setBlock(blk) end do end if useDiagPC @@ -892,9 +901,12 @@ subroutine allocDerivativeValues(level) allocate (cgnsDomsd(nn)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays(nDirichlet)) do iDirichlet = 1, nDirichlet - nArray = size(cgnsDoms(nn)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays(iDirichlet)%dataArr) - allocate (cgnsDomsd(nn)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays(iDirichlet)%dataArr(nArray)) - cgnsDomsd(nn)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays(iDirichlet)%dataArr(nArray) = zero + nArray = size(cgnsDoms(nn)%bocoInfo(iBoco)%dataSet(iData) & + %dirichletArrays(iDirichlet)%dataArr) + allocate (cgnsDomsd(nn)%bocoInfo(iBoco)% & + dataSet(iData)%dirichletArrays(iDirichlet)%dataArr(nArray)) + cgnsDomsd(nn)%bocoInfo(iBoco)%dataSet(iData)% & + dirichletArrays(iDirichlet)%dataArr(nArray) = zero end do end if end do @@ -1049,8 +1061,10 @@ subroutine zeroADSeeds(nn, level, sps) if (associated(cgnsDoms(iDom)%bocoInfo(iBoco)%dataSet)) then do iData = 1, size(cgnsDoms(iDom)%bocoInfo(iBoco)%dataSet) if (associated(cgnsDoms(iDom)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays)) then - do iDirichlet = 1, size(cgnsDoms(iDom)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays) - cgnsDomsd(iDom)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays(iDirichlet)%dataArr(:) = zero + do iDirichlet = 1, & + size(cgnsDoms(iDom)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays) + cgnsDomsd(iDom)%bocoInfo(iBoco)%dataSet(iData)%dirichletArrays(iDirichlet)%dataArr(:) & + = zero end do end if end do diff --git a/src/adjoint/fortranPC.F90 b/src/adjoint/fortranPC.F90 index e82512e94..351ebe228 100644 --- a/src/adjoint/fortranPC.F90 +++ b/src/adjoint/fortranPC.F90 @@ -268,7 +268,8 @@ subroutine setupPCMatrix(useAD, useTranspose, frozenTurb, level) do k = 0, kb do j = 0, jb do i = 0, ib - flowDoms(nn, level, sps2)%w(i, j, k, ll) = flowDoms(nn, 1, sps2)%wtmp(i, j, k, ll) + flowDoms(nn, level, sps2)%w(i, j, k, ll) = & + flowDoms(nn, 1, sps2)%wtmp(i, j, k, ll) end do end do end do @@ -370,31 +371,38 @@ subroutine setupPCMatrix(useAD, useTranspose, frozenTurb, level) ! Diagonal block is easy. if (onBlock(i, j, k)) then - PCMat(:, 1:nw, i, j, k) = flowDoms(nn, 1, sps)%dw_deriv(i, j, k, 1:nstate, 1:nstate) + PCMat(:, 1:nw, i, j, k) = & + flowDoms(nn, 1, sps)%dw_deriv(i, j, k, 1:nstate, 1:nstate) end if if (onBlock(i - 1, j, k)) then - PCMat(:, 2 * nw + 1:3 * nw, i - 1, j, k) = flowDoms(nn, 1, sps)%dw_deriv(i - 1, j, k, 1:nstate, 1:nstate) + PCMat(:, 2 * nw + 1:3 * nw, i - 1, j, k) = & + flowDoms(nn, 1, sps)%dw_deriv(i - 1, j, k, 1:nstate, 1:nstate) end if if (onBlock(i + 1, j, k)) then - PCMat(:, 1 * nw + 1:2 * nw, i + 1, j, k) = flowDoms(nn, 1, sps)%dw_deriv(i + 1, j, k, 1:nstate, 1:nstate) + PCMat(:, 1 * nw + 1:2 * nw, i + 1, j, k) = & + flowDoms(nn, 1, sps)%dw_deriv(i + 1, j, k, 1:nstate, 1:nstate) end if if (onBlock(i, j - 1, k)) then - PCMat(:, 4 * nw + 1:5 * nw, i, j - 1, k) = flowDoms(nn, 1, sps)%dw_deriv(i, j - 1, k, 1:nstate, 1:nstate) + PCMat(:, 4 * nw + 1:5 * nw, i, j - 1, k) = & + flowDoms(nn, 1, sps)%dw_deriv(i, j - 1, k, 1:nstate, 1:nstate) end if if (onBlock(i, j + 1, k)) then - PCMat(:, 3 * nw + 1:4 * nw, i, j + 1, k) = flowDoms(nn, 1, sps)%dw_deriv(i, j + 1, k, 1:nstate, 1:nstate) + PCMat(:, 3 * nw + 1:4 * nw, i, j + 1, k) = & + flowDoms(nn, 1, sps)%dw_deriv(i, j + 1, k, 1:nstate, 1:nstate) end if if (onBlock(i, j, k - 1)) then - PCMat(:, 6 * nw + 1:7 * nw, i, j, k - 1) = flowDoms(nn, 1, sps)%dw_deriv(i, j, k - 1, 1:nstate, 1:nstate) + PCMat(:, 6 * nw + 1:7 * nw, i, j, k - 1) = & + flowDoms(nn, 1, sps)%dw_deriv(i, j, k - 1, 1:nstate, 1:nstate) end if if (onBlock(i, j, k + 1)) then - PCMat(:, 5 * nw + 1:6 * nw, i, j, k + 1) = flowDoms(nn, 1, sps)%dw_deriv(i, j, k + 1, 1:nstate, 1:nstate) + PCMat(:, 5 * nw + 1:6 * nw, i, j, k + 1) = & + flowDoms(nn, 1, sps)%dw_deriv(i, j, k + 1, 1:nstate, 1:nstate) end if end if end do iLoop diff --git a/src/f2py/adflow.pyf b/src/f2py/adflow.pyf index f425ce9e8..fcd30981a 100644 --- a/src/f2py/adflow.pyf +++ b/src/f2py/adflow.pyf @@ -423,7 +423,7 @@ python module libadflow end module surfaceutils module tecplotio - subroutine addparaslice(slicename,pt,normal,dir_vec,use_dir,famlist,n) ! in :test:liftDistribution.F90 + subroutine addparaslice(slicename,pt,normal,dir_vec,use_dir,famlist,n,sliceindex) ! in :test:liftDistribution.F90 character*(*) intent(in) :: slicename real(kind=realtype) dimension(3),intent(in) :: pt real(kind=realtype) dimension(3),intent(in) :: normal @@ -431,9 +431,10 @@ python module libadflow logical, intent(in) :: use_dir integer(kind=inttype) dimension(n),intent(in) :: famlist integer(kind=inttype), optional,intent(in),check(len(famlist)>=n),depend(famlist) :: n=len(famlist) + integer(kind=inttype), intent(out) :: sliceindex end subroutine addparaslice - subroutine addabsslice(slicename,pt,normal,dir_vec,use_dir,famlist,n) ! in :test:liftDistribution.F90 + subroutine addabsslice(slicename,pt,normal,dir_vec,use_dir,famlist,n,sliceindex) ! in :test:liftDistribution.F90 character*(*) intent(in) :: slicename real(kind=realtype) dimension(3),intent(in) :: pt real(kind=realtype) dimension(3),intent(in) :: normal @@ -441,6 +442,7 @@ python module libadflow logical, intent(in) :: use_dir integer(kind=inttype) dimension(n),intent(in) :: famlist integer(kind=inttype), optional,intent(in),check(len(famlist)>=n),depend(famlist) :: n=len(famlist) + integer(kind=inttype), intent(out) :: sliceindex end subroutine addabsslice subroutine addliftdistribution(nsegments,normal,normal_ind, distname,famlist,n) ! in :test:liftDistribution.F90 @@ -465,8 +467,39 @@ python module libadflow subroutine initializeliftdistributiondata end subroutine initializeliftdistributiondata + + subroutine computetempnodaldata + end subroutine computetempnodaldata + + subroutine getparaslicedatasize(sliceindex,nnode,nelem,ndata) + integer(kind=inttype) intent(in) :: sliceindex + integer(kind=inttype) intent(out) :: nnode, nelem, ndata + end subroutine getparaslicedatasize + + subroutine getparaslicedata(sliceindex,slicenodaldata,nnode,ndata,sliceconn,nelem) + integer(kind=inttype) intent(in) :: sliceindex + real(kind=realtype), dimension(ndata, nnode),intent(inout) :: slicenodaldata + integer(kind=inttype), optional,intent(in),check(shape(slicenodaldata,1)==nnode),depend(slicenodaldata) :: nnode=shape(slicenodaldata,1) + integer(kind=inttype), optional,intent(in),check(shape(slicenodaldata,0)==ndata),depend(slicenodaldata) :: ndata=shape(slicenodaldata,0) + integer(kind=inttype), dimension(2, nelem),intent(inout) :: sliceconn + integer(kind=inttype), optional,intent(in),check(shape(sliceconn,1)==nelem),depend(sliceconn) :: nelem=shape(sliceconn,1) + end subroutine getparaslicedata + + subroutine deallocatetempnodaldata + end subroutine deallocatetempnodaldata end subroutine tecplotio + module outputmod + subroutine numberofsurfsolvariables(nsolvar) + integer(kind=inttype), intent(out) :: nsolvar + end subroutine numberofsurfsolvariables + + subroutine surfsolnameswrap(nsolvar, solnamesstring) + ! TODO fix the max stringlength parameter. hard coded to 256 for now. + character*maxstringlen, dimension(256, nsolvar), intent(out) :: solnamesstring + end subroutine surfsolnameswrap + end module outputmod + module surfaceintegrations subroutine getsolutionwrap(famlists,funcvalues,ncost,ngroups,nfammax) ! in :test:surfaceIntegrations.F90:surfaceintegrations integer(kind=inttype) dimension(ngroups,nfammax) :: famlists @@ -1498,4 +1531,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/output/outputMod.F90 b/src/output/outputMod.F90 index 6666ecc10..e6b72632b 100644 --- a/src/output/outputMod.F90 +++ b/src/output/outputMod.F90 @@ -716,6 +716,25 @@ subroutine surfSolNames(solNames) end subroutine surfSolNames + subroutine surfSolNamesWrap(nSolVar, surfSolNamesString) + use precision + use cgnsNames + use extraOutput + implicit none + ! + ! Subroutine argument. + ! + integer(kind=inttype), intent(in) ::nSolVar + character(len=maxStringLen), dimension(nSolVar), intent(out) :: surfSolNamesString + ! + ! Local variables. + ! + integer(kind=intType) :: nn, solNameLength + + call surfSolNames(surfSolNamesString) + + end subroutine surfSolNamesWrap + subroutine storeSolInBuffer(buffer, copyInBuffer, solName, & iBeg, iEnd, jBeg, jEnd, kBeg, kEnd) ! @@ -2309,7 +2328,8 @@ subroutine describeScheme(string) if (spaceDiscr == dissScalar) then if (dirScaling) then - write (string, "(2(A, 1X), ES12.5, A)") trim(string), "Directional scaling of dissipation with exponent", adis, "." + write (string, "(2(A, 1X), ES12.5, A)") trim(string), & + "Directional scaling of dissipation with exponent", adis, "." else write (string, stringSpace) trim(string), "No directional scaling of dissipation." end if @@ -2330,7 +2350,8 @@ subroutine describeScheme(string) write (string, stringSpace) trim(string), "Quadratic extrapolation of normal pressure gradIent", & "for inviscid wall boundary conditions." case (normalMomentum) - write (string, stringSpace) trim(string), "Normal momentum equation used to determine pressure gradient", & + write (string, stringSpace) trim(string), & + "Normal momentum equation used to determine pressure gradient", & "for inviscid wall boundary conditions." end select end if diff --git a/src/output/tecplotIO.F90 b/src/output/tecplotIO.F90 index 22a120972..b6dd91ac4 100644 --- a/src/output/tecplotIO.F90 +++ b/src/output/tecplotIO.F90 @@ -73,8 +73,11 @@ module tecplotIO character(len=maxCGNSNameLen), dimension(:), allocatable :: liftDistName integer(kind=intType), parameter :: nLiftDistVar = 26 + ! temp array to hold nodal values while we pass back slice data to python + real(kind=realType), dimension(:, :, :), allocatable :: tempNodalValues + contains - subroutine addParaSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n) + subroutine addParaSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n, sliceIndex) ! ! This subroutine is intended to be called from python. ! This routine will add a parametric slice to the list of user @@ -92,6 +95,9 @@ subroutine addParaSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n) logical, intent(in) :: use_dir integer(kind=intType), intent(in) :: n, famList(n) + ! Output parameters + integer(kind=intType), intent(out) :: sliceIndex + ! Working integer(kind=intType) :: sps, sizeNode, sizeCell integer(kind=intType), dimension(:), pointer :: wallList @@ -103,9 +109,12 @@ subroutine addParaSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n) allocate (paraSlices(nSliceMax, nTimeIntervalsSpectral)) end if + ! increment the slice counter. this also becomes the slice index + nParaSlices = nParaSlices + 1 + sliceIndex = nParaSlices + ! We have to add a slice for each spectral instance. do sps = 1, nTimeIntervalsSpectral - nParaSlices = nParaSlices + 1 if (nParaSlices > nSliceMax) then print *, 'Error: Exceeded the maximum number of slices. Increase nSliceMax' @@ -122,7 +131,7 @@ subroutine addParaSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n) call getSurfaceFamily(elemFam, sizeCell, wallList, size(wallList), .True.) ! Create actual slice - call createSlice(pts, conn, elemFam, paraSlices(nParaSlices, sps), pt, normal, dir_vec, & + call createSlice(pts, conn, elemFam, paraSlices(sliceIndex, sps), pt, normal, dir_vec, & use_dir, sliceName, famList) ! Clean up memory. @@ -131,7 +140,7 @@ subroutine addParaSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n) end subroutine addParaSlice - subroutine addAbsSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n) + subroutine addAbsSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n, sliceIndex) ! ! This subroutine is intended to be called from python. ! This routine will add an absolute slice to the list of user @@ -149,6 +158,9 @@ subroutine addAbsSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n) logical, intent(in) :: use_dir integer(kind=intType), intent(in) :: n, famList(n) + ! Output parameters + integer(kind=intType), intent(out) :: sliceIndex + ! Working integer(kind=intType) :: sps, sizeNode, sizeCell integer(kind=intType), dimension(:), pointer :: wallList @@ -160,8 +172,11 @@ subroutine addAbsSlice(sliceName, pt, normal, dir_vec, use_dir, famList, n) allocate (absSlices(nSliceMax, nTimeIntervalsSpectral)) end if + ! increment the slice counter. this also becomes the slice index + nAbsSlices = nAbsSlices + 1 + sliceIndex = nAbsSlices + do sps = 1, nTimeIntervalsSpectral - nAbsSlices = nAbsSlices + 1 if (nAbsSlices > nSliceMax) then print *, 'Error: Exceeded the maximum number of slices. Increase nSliceMax' @@ -217,6 +232,133 @@ subroutine addLiftDistribution(nSegments, normal, normal_ind, distName, famList, end subroutine addLiftDistribution + subroutine computeTempNodalData() + + use constants + use inputTimeSpectral, only: nTimeIntervalsSpectral + use surfaceFamilies, only: BCFamGroups, BCFamExchange + use surfaceUtils, only: getSurfaceSize + use oversetData, only: zipperMeshes + use outputMod, only: numberOfSurfSolVariables + implicit none + + ! Working + integer(kind=intType) :: sps, nSolVar, sizeNOde, sizeCell + integer(kind=intType), dimension(:), pointer :: wallLIST + + ! Determine the number of surface variables we have + call numberOfSurfSolVariables(nSolVar) + wallList => BCFamGroups(iBCGroupWalls)%famList + call getSurfaceSize(sizeNode, sizeCell, wallList, size(wallList), .True.) + + ! Allocate and compute the wall-based surface data for hte slices + ! and lift distributions. + allocate (tempNodalValues(max(sizeNode, 1), nSolVar + 6 + 3, nTimeIntervalsSpectral)) + + ! TODO modify this so we keep each sps instance separately + sps = 1 + call computeSurfaceOutputNodalData(BCFamExchange(iBCGroupWalls, sps), & + zipperMeshes(iBCGroupWalls), .True., tempNodalValues(:, :, sps)) + + + end subroutine computeTempNodalData + + subroutine deallocateTempNodalData() + implicit none + ! just deallocate the data we already computed + deallocate(tempNodalValues) + end subroutine deallocateTempNodalData + + subroutine getParaSliceDataSize(sliceIndex, nNode, nElem, nData) + use constants + use communication, only: adflow_comm_world + use outputMod, only: numberOfSurfSolVariables + use utils, only: EChk + implicit none + + integer(kind=intType), intent(in) :: sliceIndex + integer(kind=intType), intent(out) :: nNode, nElem, nData + + integer(kind=intType) :: sps, ierr, nElemLocal, nNodeLocal, nSolVar + + ! return slice node count and data count back to python so we can allocate the numpy array there + + sps = 1 + + ! sum up the number of nodes from all procs + nNodeLocal = paraSlices(sliceIndex, sps)%nNodes + call mpi_allreduce(nNodeLocal, nNode, 1, adflow_integer, MPI_SUM, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! add up the number of elements + nElemLocal = size(paraSlices(sliceIndex, sps)%conn, 2) + call mpi_allreduce(nElemLocal, nElem, 1, adflow_integer, MPI_SUM, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Determine the number of surface variables we have + call numberOfSurfSolVariables(nSolVar) + ! add 3 + 6 because 3 coordinates, 6 tractions (3 pressure and 3 viscous) + ! TODO we may want to change this; only pass a subset back + nData = nSolVar + 3 + 6 + + end subroutine + + subroutine getParaSliceData(sliceIndex, sliceNodalData, nNode, nData, sliceConn, nElem) + + use constants + ! use inputTimeSpectral, only: nTimeIntervalsSpectral + ! use surfaceFamilies, only: BCFamGroups, BCFamExchange + ! use surfaceUtils, only: getSurfaceSize + ! use oversetData, only: zipperMeshes + use communication, only: myid, adflow_comm_world + use outputMod, only: numberOfSurfSolVariables + use utils, only: EChk + implicit none + + ! Input/Output Params + integer(kind=intType), intent(in) :: sliceIndex + integer(kind=inttype), intent(in):: nNode, nData, nElem + real(kind=realtype) ,intent(inout), dimension(ndata, nnode) :: sliceNodalData + integer(kind=inttype) ,intent(inout), dimension(2, nElem) :: sliceConn + + ! Working + integer(kind=intType) :: sps, nSolVar, ii, jj, ierr + type(slice) :: globalSlice + + sps = 1 + call numberOfSurfSolVariables(nSolVar) + + ! assume the nodal data is already computed. + ! run the regular integration call + call integrateSlice(paraSlices(sliceIndex, sps), globalSlice, & + tempNodalValues(:, :, sps), nSolVar, .True.) + + ! copy the data we want + if (myid == 0) then + do ii = 1, nNode + do jj = 1, nData + sliceNodalData(jj, ii) = globalSlice%vars(jj, ii) + end do + end do + end if + + ! broadcast the nodal data + call mpi_bcast(sliceNodalData, nNode*nData, adflow_real, 0, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! do the same for the conn + if (myid == 0) then + do ii = 1, nElem + sliceConn(:, ii) = globalSlice%conn(:, ii) + end do + end if + + ! broadcast the nodal data + call mpi_bcast(sliceConn, nElem*2, adflow_integer, 0, adflow_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + end subroutine getParaSliceData + subroutine writeTecplot(sliceFile, writeSlices, liftFile, writeLift, & surfFile, writeSurf, famList, nFamList) ! @@ -1001,7 +1143,8 @@ subroutine writeTecplotSurfaceFile(fileName, famList) ! Gather values to the root proc. do i = 1, iSize call mpi_gatherv(nodalValues(:, i), sizeNode, & - adflow_real, vars(:, i), nodeSizes, nodeDisps, adflow_real, 0, adflow_comm_world, ierr) + adflow_real, vars(:, i), nodeSizes, nodeDisps, & + adflow_real, 0, adflow_comm_world, ierr) call EChk(ierr, __FILE__, __LINE__) end do deallocate (nodalValues) @@ -2207,19 +2350,19 @@ subroutine integrateSlice(lSlc, gSlc, nodalValues, nFields, doConnectivity) end if if (normal_ind == 1) then - M(1, 1) = one; M(1, 2) = zero; M(1, 3) = zero; - M(2, 1) = zero; M(2, 2) = one; M(2, 3) = zero; - M(3, 1) = zero; M(3, 2) = zero; M(3, 3) = one; + M(1, 1) = one; M(1, 2) = zero; M(1, 3) = zero; + M(2, 1) = zero; M(2, 2) = one; M(2, 3) = zero; + M(3, 1) = zero; M(3, 2) = zero; M(3, 3) = one; else if (normal_ind == 2) then ! Y-rotation matrix - M(1, 1) = cos(-theta); M(1, 2) = zero; M(1, 3) = sin(-theta); - M(2, 1) = zero; M(2, 2) = one; M(2, 3) = zero; - M(3, 1) = -sin(-theta); M(3, 2) = zero; M(3, 3) = cos(-theta); + M(1, 1) = cos(-theta); M(1, 2) = zero; M(1, 3) = sin(-theta); + M(2, 1) = zero; M(2, 2) = one; M(2, 3) = zero; + M(3, 1) = -sin(-theta); M(3, 2) = zero; M(3, 3) = cos(-theta); else ! Z rotation Matrix - M(1, 1) = cos(theta); M(1, 2) = -sin(theta); M(1, 3) = zero; - M(2, 1) = sin(theta); M(2, 2) = cos(theta); M(2, 3) = zero; - M(3, 1) = zero; M(3, 2) = zero; M(3, 3) = one; + M(1, 1) = cos(theta); M(1, 2) = -sin(theta); M(1, 3) = zero; + M(2, 1) = sin(theta); M(2, 2) = cos(theta); M(2, 3) = zero; + M(3, 1) = zero; M(3, 2) = zero; M(3, 3) = one; end if allocate (tempCoords(3, size(gSlc%vars, 2))) diff --git a/src/output/writeCGNSSurface.F90 b/src/output/writeCGNSSurface.F90 index ebbe0c622..8a6bcf3b6 100644 --- a/src/output/writeCGNSSurface.F90 +++ b/src/output/writeCGNSSurface.F90 @@ -1724,11 +1724,13 @@ subroutine writeIsoSurface(isoName, sps, nIsoSurfVar, isoSurfSolNames) call cg_coord_partial_write_f(cgnsInd, cgnsBase, cgnsZone, realDouble, & 'CoordinateY', cumNodes + 1, cumNodes + nPtsProc(iProc + 1), & - buffer(nPtsProc(iProc + 1) + 1:2 * nPtsProc(iProc + 1)), coordID, ierr) + buffer(nPtsProc(iProc + 1) + & + 1:2 * nPtsProc(iProc + 1)), coordID, ierr) call cg_coord_partial_write_f(cgnsInd, cgnsBase, cgnsZone, realDouble, & 'CoordinateZ', cumNodes + 1, cumNodes + nPtsProc(iProc + 1), & - buffer(2 * nPtsProc(iProc + 1) + 1:3 * nPtsProc(iProc + 1)), coordID, ierr) + buffer(2 * nPtsProc(iProc + 1) + & + 1:3 * nPtsProc(iProc + 1)), coordID, ierr) if (ierr /= CG_OK) & call terminate("writeIsoSurface", & @@ -1892,7 +1894,8 @@ subroutine writeIsoSurface(isoName, sps, nIsoSurfVar, isoSurfSolNames) ! received from iProc if (nPtsProc(iProc + 1) > 0) Then call cg_field_partial_write_f(cgnsInd, cgnsBase, cgnsZone, solID, realDouble, & - isoSurfSolNames(iVar), cumNodes + 1, cumNodes + nPtsProc(iProc + 1), & + isoSurfSolNames(iVar), cumNodes + 1, & + cumNodes + nPtsProc(iProc + 1), & buffer, fieldID, ierr) if (ierr /= CG_OK) & diff --git a/src/overset/cartMesh.F90 b/src/overset/cartMesh.F90 index ec41eac10..405191154 100644 --- a/src/overset/cartMesh.F90 +++ b/src/overset/cartMesh.F90 @@ -173,7 +173,8 @@ subroutine createCartMesh(level, sps) iDim = maxloc(abs(bcData(mm)%symNorm), 1) ! Now determine the average value in "iDim" dimension - coorAvg = sum(xx(iBeg + 1:iEnd + 1, jBeg + 1:jEnd + 1, iDim)) / ((iEnd - iBeg + 1) * (jEnd - jBeg + 1)) + coorAvg = sum(xx(iBeg + 1:iEnd + 1, jBeg + 1:jEnd + 1, iDim)) & + / ((iEnd - iBeg + 1) * (jEnd - jBeg + 1)) ! Check if it is sufficently close to the bounding box: err1 = abs(coorAvg - xMin(iDim)) / scaleSize @@ -919,7 +920,8 @@ subroutine writeCartMesh(blankVec, cellDims, xMin, h) call cg_base_write_f(cg, "Base#1", 3, 3, base, ierr) call cg_zone_write_f(cg, base, "cartblock", int((/cellDims(1) + 1, cellDims(2) + 1, cellDims(3) + 1, & - cellDims(1), cellDims(2), cellDims(3), 0, 0, 0/), cgsize_t), Structured, zoneID, ierr) + cellDims(1), cellDims(2), cellDims(3), 0, 0, 0/), & + cgsize_t), Structured, zoneID, ierr) allocate (xtmp(cellDims(1) + 1, cellDims(2) + 1, cellDims(3) + 1, 3)) diff --git a/src/overset/computeHolesInsideBody.F90 b/src/overset/computeHolesInsideBody.F90 index 463c662d7..327b41226 100644 --- a/src/overset/computeHolesInsideBody.F90 +++ b/src/overset/computeHolesInsideBody.F90 @@ -533,7 +533,8 @@ subroutine computeHolesInsideBody(level, sps) coor(4) = large call minDistancetreeSearchSinglePoint(walls(c)%ADT, coor, & - intInfo2, uvw2, dummy, 0, BB, frontLeaves, frontLeavesNew) + intInfo2, uvw2, dummy, 0, BB, & + frontLeaves, frontLeavesNew) cellID2 = intInfo2(3) if (uvw2(4) < nearWallDist**2) then @@ -588,7 +589,8 @@ subroutine computeHolesInsideBody(level, sps) ! dist cutoff. coor(4) = large call minDistancetreeSearchSinglePoint(fullWall%ADT, coor, & - intInfo, uvw, dummy, 0, BB, frontLeaves, frontLeavesNew) + intInfo, uvw, dummy, 0, BB, & + frontLeaves, frontLeavesNew) cellID = intInfo(3) ! Determine if it is inside: diff --git a/src/overset/floodInteriorCells.F90 b/src/overset/floodInteriorCells.F90 index faa896e80..458c9211c 100644 --- a/src/overset/floodInteriorCells.F90 +++ b/src/overset/floodInteriorCells.F90 @@ -157,7 +157,8 @@ subroutine floodInteriorCells(level, sps) k = stack(3, stackPointer) stackPointer = stackPointer - 1 - if (isCompute(status(i, j, k)) .and. .not. isReceiver(status(i, j, k)) .and. iblank(i, j, k) /= -4) then + if (isCompute(status(i, j, k)) .and. .not. isReceiver(status(i, j, k)) & + .and. iblank(i, j, k) /= -4) then ! Flag the cell (using changed) as being changed changed(i + 1, j + 1, k + 1) = 1 diff --git a/src/overset/makeBoundaryStrings.F90 b/src/overset/makeBoundaryStrings.F90 index a978fc86a..031974a02 100644 --- a/src/overset/makeBoundaryStrings.F90 +++ b/src/overset/makeBoundaryStrings.F90 @@ -437,11 +437,13 @@ subroutine makeGapBoundaryStrings(zipperFamList, level, sps, master) iSize = iEnd - iStart + 1 ! ----------- Node sized arrays ------------- - call MPI_Recv(globalStrings(c)%nodeData(:, iStart:iEnd), iSize * 10, adflow_real, iProc, iProc, & + call MPI_Recv(globalStrings(c)%nodeData(:, iStart:iEnd), iSize * 10, & + adflow_real, iProc, iProc, & adflow_comm_world, mpiStatus, ierr) call ECHK(ierr, __FILE__, __LINE__) - call MPI_Recv(globalStrings(c)%intNodeData(:, iStart:iEnd), iSize * 3, adflow_integer, iProc, iProc, & + call MPI_Recv(globalStrings(c)%intNodeData(:, iStart:iEnd), iSize * 3, & + adflow_integer, iProc, iProc, & adflow_comm_world, mpiStatus, ierr) call ECHK(ierr, __FILE__, __LINE__) diff --git a/src/overset/oversetAPI.F90 b/src/overset/oversetAPI.F90 index 824bd3c4c..698918c2d 100644 --- a/src/overset/oversetAPI.F90 +++ b/src/overset/oversetAPI.F90 @@ -34,10 +34,14 @@ subroutine oversetComm(level, firstTime, coarseLevel, closedFamList) exchangeFringes, sendOFringe, sendOBlock, setupFringeGlobalInd, & getFringeReturnSizes use oversetUtilities, only: isCompute, checkOverset, irregularCellCorrection, & - fringeReduction, transposeOverlap, setIBlankArray, deallocateOFringes, deallocateoBlocks, & - deallocateOSurfs, deallocateCSRMatrix, setIsCompute, getWorkArray, flagForcedRecv, & - qsortFringeType, isReceiver, setIsReceiver, addToFringeList, printOverlapMatrix, & - tic, toc, unwindIndex, windIndex, isFloodSeed, isFlooded, wallsOnBlock + fringeReduction, transposeOverlap, setIBlankArray, & + deallocateOFringes, deallocateoBlocks, & + deallocateOSurfs, deallocateCSRMatrix, & + setIsCompute, getWorkArray, flagForcedRecv, & + qsortFringeType, isReceiver, setIsReceiver, & + addToFringeList, printOverlapMatrix, & + tic, toc, unwindIndex, windIndex, isFloodSeed, & + isFlooded, wallsOnBlock use oversetPackingRoutines, only: packOFringe, packOBlock, unpackOFringe, unpackOBlock, & getOFringeBufferSizes, getOBlockBufferSizes, getOSurfBufferSizes implicit none @@ -847,11 +851,14 @@ subroutine oversetComm(level, firstTime, coarseLevel, closedFamList) aspect = one if (useOversetWallScaling) then if (CGNSDoms(nbkGlobal)%viscousDir(1)) & - aspect(1) = (half * (mynorm2(si(i - 1, j, k, :)) + mynorm2(si(i, j, k, :)))) / vol(i, j, k) + aspect(1) = (half * (mynorm2(si(i - 1, j, k, :)) + & + mynorm2(si(i, j, k, :)))) / vol(i, j, k) if (CGNSDoms(nbkGlobal)%viscousDir(2)) & - aspect(2) = (half * (mynorm2(sj(i, j - 1, k, :)) + mynorm2(sj(i, j, k, :)))) / vol(i, j, k) + aspect(2) = (half * (mynorm2(sj(i, j - 1, k, :)) + & + mynorm2(sj(i, j, k, :)))) / vol(i, j, k) if (CGNSDoms(nbkGlobal)%viscousDir(3)) & - aspect(3) = (half * (mynorm2(sk(i, j, k - 1, :)) + mynorm2(sk(i, j, k, :)))) / vol(i, j, k) + aspect(3) = (half * (mynorm2(sk(i, j, k - 1, :)) + & + mynorm2(sk(i, j, k, :)))) / vol(i, j, k) end if fact = min(aspect(1) * aspect(2) * aspect(3), 100.0_realType) @@ -2255,7 +2262,8 @@ subroutine flagCellsInSurface(pts, npts, conn, nconn, flag, ncell, blockids, nbl coor(4) = dStar intInfo(3) = 0 call minDistancetreeSearchSinglePoint(ADT, coor, intInfo, & - uvw, dummy, 0, BB, frontLeaves, frontLeavesNew) + uvw, dummy, 0, BB, & + frontLeaves, frontLeavesNew) cellID = intInfo(3) if (cellID > 0) then ! Now check if this was successful or not: diff --git a/src/overset/oversetCommUtilites.F90 b/src/overset/oversetCommUtilites.F90 index 35a0dd221..a94e78a29 100644 --- a/src/overset/oversetCommUtilites.F90 +++ b/src/overset/oversetCommUtilites.F90 @@ -150,7 +150,8 @@ subroutine getOSurfCommPattern(oMat, oMatT, sendList, nSend, & procsForThisRow(1:nnRow) = oMat%assignedProc(oMat%rowPtr(iDom):oMat%rowPtr(iDom + 1) - 1) nnRowT = oMatT%rowPtr(iDom + 1) - oMatT%rowPtr(iDom) - procsForThisRow(nnRow + 1:nnRow + nnRowT) = oMatT%assignedProc(oMatT%rowPtr(iDom):oMatT%rowPtr(iDom + 1) - 1) + procsForThisRow(nnRow + 1:nnRow + nnRowT) = & + oMatT%assignedProc(oMatT%rowPtr(iDom):oMatT%rowPtr(iDom + 1) - 1) call unique(procsForThisRow, nnRow + nnRowT, nUniqueProc, inverse) @@ -1948,7 +1949,8 @@ subroutine updateOversetConnectivity(level, sps) ! Compute new fraction frac0 = (/half, half, half/) call newtonUpdate(xCen, & - flowDoms(d2, level, sps)%x(i2 - 1:i2 + 1, j2 - 1:j2 + 1, k2 - 1:k2 + 1, :), frac0, frac) + flowDoms(d2, level, sps)%x(i2 - 1:i2 + 1, j2 - 1:j2 + 1, k2 - 1:k2 + 1, :), & + frac0, frac) ! Check if the fractions are between zero and one if (MAXVAL(frac) > one + fracTol .or. MINVAL(frac) < zero - fracTol) then @@ -2279,7 +2281,8 @@ subroutine updateOversetConnectivity_b(level, sps) ! Do newton update frac0 = (/half, half, half/) call newtonUpdate(xCen, & - flowDoms(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), frac0, frac) + flowDoms(d1, level, sps)%x(i1 - 1:i1 + 1, j1 - 1:j1 + 1, k1 - 1:k1 + 1, :), & + frac0, frac) ! Set the new weights call fracToWeights(frac, weight) diff --git a/src/overset/oversetInitialization.F90 b/src/overset/oversetInitialization.F90 index c83c39dc1..9fc99f346 100644 --- a/src/overset/oversetInitialization.F90 +++ b/src/overset/oversetInitialization.F90 @@ -147,11 +147,14 @@ subroutine initializeOBlock(oBlock, nn, level, sps) aspect = one if (useOversetWallScaling) then if (CGNSDoms(nbkGlobal)%viscousDir(1)) & - aspect(1) = (half * (mynorm2(si(i - 1, j, k, :)) + mynorm2(si(i, j, k, :)))) / vol(i, j, k) + aspect(1) = (half * (mynorm2(si(i - 1, j, k, :)) + & + mynorm2(si(i, j, k, :)))) / vol(i, j, k) if (CGNSDoms(nbkGlobal)%viscousDir(2)) & - aspect(2) = (half * (mynorm2(sj(i, j - 1, k, :)) + mynorm2(sj(i, j, k, :)))) / vol(i, j, k) + aspect(2) = (half * (mynorm2(sj(i, j - 1, k, :)) + & + mynorm2(sj(i, j, k, :)))) / vol(i, j, k) if (CGNSDoms(nbkGlobal)%viscousDir(3)) & - aspect(3) = (half * (mynorm2(sk(i, j, k - 1, :)) + mynorm2(sk(i, j, k, :)))) / vol(i, j, k) + aspect(3) = (half * (mynorm2(sk(i, j, k - 1, :)) + & + mynorm2(sk(i, j, k, :)))) / vol(i, j, k) end if fact = min(aspect(1) * aspect(2) * aspect(3), 100.0_realType) diff --git a/src/partitioning/gridChecking.F90 b/src/partitioning/gridChecking.F90 index c1b19bbe1..be3ef7bd3 100644 --- a/src/partitioning/gridChecking.F90 +++ b/src/partitioning/gridChecking.F90 @@ -165,7 +165,8 @@ subroutine checkFaces case (kMax) print strings, "Zone ", trim(cgnsDoms(nn)%zoneName), ": kMax block face not fully described" case default - print strings, "Zone ", trim(cgnsDoms(nn)%zoneName), ": Multiple block faces not fully described" + print strings, "Zone ", trim(cgnsDoms(nn)%zoneName), & + ": Multiple block faces not fully described" end select else @@ -184,7 +185,8 @@ subroutine checkFaces case (kMax) print strings, "Zone ", trim(cgnsDoms(nn)%zoneName), ": kMax block face multiple described" case default - print strings, "Zone ", trim(cgnsDoms(nn)%zoneName), ": Multiple block faces multiple described" + print strings, "Zone ", trim(cgnsDoms(nn)%zoneName), & + ": Multiple block faces multiple described" end select end if @@ -932,12 +934,14 @@ subroutine check1to1Subfaces if (nTimeIntervalsSpectral > 1) then if (badGlobal(4, nn) == 1) then - print spectralDevFormat, "# Spectral grid ", trim(intString), ", zone ", trim(cgnsDoms(i)%zoneName), & + print spectralDevFormat, "# Spectral grid ", trim(intString), & + ", zone ", trim(cgnsDoms(i)%zoneName), & ", periodic subface ", trim(cgnsDoms(i)%conn1to1(j)%connectName), & " does not match donor ", trim(cgnsDoms(i)%conn1to1(j)%donorName), & ". Maximum deviation: ", badDistGlobal(nn) else - print spectralDevFormat, "# Spectral grid ", trim(intString), ", zone ", trim(cgnsDoms(i)%zoneName), & + print spectralDevFormat, "# Spectral grid ", trim(intString), & + ", zone ", trim(cgnsDoms(i)%zoneName), & ", subface ", trim(cgnsDoms(i)%conn1to1(j)%connectName), & " does not match donor ", trim(cgnsDoms(i)%conn1to1(j)%donorName), & ". Maximum deviation: ", badDistGlobal(nn) diff --git a/src/partitioning/loadBalance.F90 b/src/partitioning/loadBalance.F90 index 72e39cb48..6b1fe989c 100644 --- a/src/partitioning/loadBalance.F90 +++ b/src/partitioning/loadBalance.F90 @@ -1720,10 +1720,12 @@ subroutine BCFacesSubblock(cgnsID, ii, jj) ! accordingly. if (flowType == internalFlow) then - write (errorMessage, strings) "Zone ", trim(zoneName), ", boundary subface ", trim(subName), & + write (errorMessage, strings) "Zone ", trim(zoneName), & + ", boundary subface ", trim(subName), & ": Not a valid boundary condition for internal flow" else - write (errorMessage, strings) "Zone ", trim(zoneName), ", boundary subface ", trim(subName), & + write (errorMessage, strings) "Zone ", trim(zoneName), & + ", boundary subface ", trim(subName), & ": Not a valid boundary condition for external flow" end if @@ -1868,7 +1870,8 @@ subroutine externalFacesSubblock(cgnsID, ii, jj, nSubPerCGNS, & if (myID == 0) then zoneName = cgnsDoms(cgnsID)%zoneName subName = cgnsDoms(cgnsID)%bocoInfo(j)%bocoName - write (errorMessage, strings) "Zone ", trim(zoneName), ", 1 to 1 block connectivity ", trim(subName), & + write (errorMessage, strings) "Zone ", trim(zoneName), & + ", 1 to 1 block connectivity ", trim(subName), & ": No constant index found for subface" call terminate("externalFacesSubblock", errorMessage) end if diff --git a/src/partitioning/readCGNSGrid.F90 b/src/partitioning/readCGNSGrid.F90 index 78436ad65..66dcf908a 100644 --- a/src/partitioning/readCGNSGrid.F90 +++ b/src/partitioning/readCGNSGrid.F90 @@ -1200,7 +1200,8 @@ subroutine read1to1Conn(cgnsInd, cgnsBase, nZone) cgnsDoms(nZone)%conn1to1(i)%translation = zero call cg_1to1_periodic_read_f(cgnsInd, cgnsBase, nZone, i, & - real(rotCenter, cgnsPerType), real(rotAngles, cgnsPerType), real(tlation, cgnsPerType), ierr) + real(rotCenter, cgnsPerType), real(rotAngles, cgnsPerType), & + real(tlation, cgnsPerType), ierr) if (ierr == CG_OK) then call readPeriodicSubface1to1(cgnsInd, cgnsBase, nZone, i, & cgnsDoms(nZone)%conn1to1(i)%connectName, & @@ -2127,7 +2128,8 @@ subroutine readBocos(cgnsInd, cgnsBase, nZone, & if (cgnsDoms(nZone)%bocoInfo(i)%BCType == bcNull) then write (errorMessage, strings) "Zone ", trim(cgnsDoms(nZone)%zoneName), & ", boundary face ", trim(cgnsDoms(nZone)%bocoInfo(i)%bocoName), & - ": Unknown user-defined boundary condition ", trim(cgnsDoms(nZone)%bocoInfo(i)%userDefinedName) + ": Unknown user-defined boundary condition ", & + trim(cgnsDoms(nZone)%bocoInfo(i)%userDefinedName) if (myID == 0) call terminate("readBocos", errorMessage) call mpi_barrier(ADflow_comm_world, ierr) end if @@ -2520,7 +2522,8 @@ subroutine readBCDataArrays(nArr, arr, DirNeu) print "(a)", "# Warning" print strings, "# Zone ", trim(cgnsDoms(nZone)%zoneName), & ", boundary subface ", trim(cgnsDoms(nZone)%bocoInfo(i)%bocoName) - print strings, "# BC data set ", trim(arr(k)%arrayName), ": No units specified, assuming SI units" + print strings, "# BC data set ", & + trim(arr(k)%arrayName), ": No units specified, assuming SI units" print "(a)", "#" end if end if @@ -2785,7 +2788,8 @@ subroutine readPeriodicSubface(cgnsInd, cgnsBase, zone, conn, & ! Check if this is a periodic boundary. call cg_conn_periodic_read_f(cgnsInd, cgnsBase, zone, conn, & - real(rotCenter, cgnsPerType), real(rotAngles, cgnsPerType), real(tlation, cgnsPerType), ierr) + real(rotCenter, cgnsPerType), real(rotAngles, cgnsPerType), & + real(tlation, cgnsPerType), ierr) testPeriodic: if (ierr == CG_OK) then @@ -2926,7 +2930,8 @@ subroutine readPeriodicSubface1to1(cgnsInd, cgnsBase, zone, conn, & ! Check if this is a periodic boundary. call cg_1to1_periodic_read_f(cgnsInd, cgnsBase, zone, conn, & - real(rotCenter, cgnsPerType), real(rotAngles, cgnsPerType), real(tlation, cgnsPerType), ierr) + real(rotCenter, cgnsPerType), real(rotAngles, cgnsPerType), & + real(tlation, cgnsPerType), ierr) testPeriodic: if (ierr == CG_OK) then @@ -2964,7 +2969,8 @@ subroutine readPeriodicSubface1to1(cgnsInd, cgnsBase, zone, conn, & print "(a)", "#" print "(a)", "# Warning" - print strings, "# Zone ", trim(cgnsDoms(zone)%zonename), ", 1to1 connectivity ", trim(connectName), & + print strings, "# Zone ", & + trim(cgnsDoms(zone)%zonename), ", 1to1 connectivity ", trim(connectName), & ": No unit specified for periodic angles, assuming radians." print "(a)", "#" diff --git a/src/preprocessing/preprocessingAPI.F90 b/src/preprocessing/preprocessingAPI.F90 index 00221e62f..e3090a8a0 100644 --- a/src/preprocessing/preprocessingAPI.F90 +++ b/src/preprocessing/preprocessingAPI.F90 @@ -1011,7 +1011,8 @@ subroutine checkSymmetry(level) print "(a)", "# Warning" print stringSpace, "# Symmetry boundary face", trim(cgnsDoms(i)%bocoInfo(j)%bocoName), & "of zone", trim(cgnsDoms(i)%zonename), "is not planar." - write (*, stringSci5) "# Maximum deviation from the mean normal: ", real(fact), " degrees" + write (*, stringSci5) "# Maximum deviation from the mean normal: ", & + real(fact), " degrees" print "(a)", "#" end if @@ -1364,7 +1365,8 @@ subroutine setSurfaceFamilyInfo write (*, "(2(A, I4), *(A))") "CGNS Block ", i, ", boundary condition ", j, ", of type ", & trim(BCTypeName(cgnsDoms(i)%bocoInfo(j)%BCTypeCGNS)), & " does not have a family. Based on the boundary condition type,", & - " a name of: '", trim(defaultFamName(cgnsDoms(i)%bocoInfo(j)%BCTypeCGNS)), "' will be used." + " a name of: '", trim(defaultFamName(cgnsDoms(i)%bocoInfo(j)%BCTypeCGNS)), & + "' will be used." end if cgnsDoms(i)%bocoInfo(j)%wallBCName = trim(defaultFamName(cgnsDoms(i)%bocoInfo(j)%BCTypeCGNS)) end if diff --git a/src/solver/fluxes.F90 b/src/solver/fluxes.F90 index b64b7458e..533a969bb 100644 --- a/src/solver/fluxes.F90 +++ b/src/solver/fluxes.F90 @@ -1229,7 +1229,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 @@ -1238,7 +1239,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 @@ -1247,7 +1249,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 @@ -1256,7 +1259,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 @@ -1306,7 +1310,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 @@ -1315,7 +1320,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 @@ -1324,7 +1330,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 @@ -1333,7 +1340,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 @@ -1382,7 +1390,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 @@ -1391,7 +1400,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 @@ -1400,7 +1410,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 @@ -1409,7 +1420,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 @@ -4025,7 +4037,8 @@ subroutine inviscidDissFluxScalarApprox ! of the face. dss2 = abs((shockSensor(i + 2, j, k) - two * shockSensor(i + 1, j, k) + shockSensor(i, j, k)) & - / (shockSensor(i + 2, j, k) + two * shockSensor(i + 1, j, k) + shockSensor(i, j, k) + sslim)) + / (shockSensor(i + 2, j, k) + two * shockSensor(i + 1, j, k) + & + shockSensor(i, j, k) + sslim)) ! Compute the dissipation coefficients for this face. @@ -4110,7 +4123,8 @@ subroutine inviscidDissFluxScalarApprox ! of the face. dss2 = abs((shockSensor(i, j + 2, k) - two * shockSensor(i, j + 1, k) + shockSensor(i, j, k)) & - / (shockSensor(i, j + 2, k) + two * shockSensor(i, j + 1, k) + shockSensor(i, j, k) + sslim)) + / (shockSensor(i, j + 2, k) + two * shockSensor(i, j + 1, k) + & + shockSensor(i, j, k) + sslim)) ! Compute the dissipation coefficients for this face. @@ -4190,7 +4204,8 @@ subroutine inviscidDissFluxScalarApprox ! of the face. dss2 = abs((shockSensor(i, j, k + 2) - two * shockSensor(i, j, k + 1) + shockSensor(i, j, k)) & - / (shockSensor(i, j, k + 2) + two * shockSensor(i, j, k + 1) + shockSensor(i, j, k) + sslim)) + / (shockSensor(i, j, k + 2) + two * shockSensor(i, j, k + 1) + & + shockSensor(i, j, k) + sslim)) ! Compute the dissipation coefficients for this face. @@ -4434,7 +4449,8 @@ subroutine inviscidDissFluxMatrixApprox ! of the face. dp2 = abs((shockSensor(i + 2, j, k) - two * shockSensor(i + 1, j, k) + shockSensor(i, j, k)) & - / (omega * (shockSensor(i + 2, j, k) + two * shockSensor(i + 1, j, k) + shockSensor(i, j, k)) & + / (omega * (shockSensor(i + 2, j, k) + & + two * shockSensor(i + 1, j, k) + shockSensor(i, j, k)) & + oneMinOmega * (abs(shockSensor(i + 2, j, k) - shockSensor(i + 1, j, k)) & + abs(shockSensor(i + 1, j, k) - shockSensor(i, j, k))) + plim)) @@ -4612,7 +4628,8 @@ subroutine inviscidDissFluxMatrixApprox ! of the face. dp2 = abs((shockSensor(i, j + 2, k) - two * shockSensor(i, j + 1, k) + shockSensor(i, j, k)) & - / (omega * (shockSensor(i, j + 2, k) + two * shockSensor(i, j + 1, k) + shockSensor(i, j, k)) & + / (omega * (shockSensor(i, j + 2, k) + & + two * shockSensor(i, j + 1, k) + shockSensor(i, j, k)) & + oneMinOmega * (abs(shockSensor(i, j + 2, k) - shockSensor(i, j + 1, k)) & + abs(shockSensor(i, j + 1, k) - shockSensor(i, j, k))) + plim)) @@ -4790,7 +4807,8 @@ subroutine inviscidDissFluxMatrixApprox ! of the face. dp2 = abs((shockSensor(i, j, k + 2) - two * shockSensor(i, j, k + 1) + shockSensor(i, j, k)) & - / (omega * (shockSensor(i, j, k + 2) + two * shockSensor(i, j, k + 1) + shockSensor(i, j, k)) & + / (omega * (shockSensor(i, j, k + 2) + & + two * shockSensor(i, j, k + 1) + shockSensor(i, j, k)) & + oneMinOmega * (abs(shockSensor(i, j, k + 2) - shockSensor(i, j, k + 1)) & + abs(shockSensor(i, j, k + 1) - shockSensor(i, j, k))) + plim)) diff --git a/src/solver/residuals.F90 b/src/solver/residuals.F90 index e12c01109..83e5abf8c 100644 --- a/src/solver/residuals.F90 +++ b/src/solver/residuals.F90 @@ -259,7 +259,8 @@ subroutine residual_block A55 = one * ((-(resM**2)) / 2) B11 = A11 * (gamma(i, j, k) - 1) * q / 2 + A12 * (-velXrho) & - / w(i, j, k, irho) + A13 * (-velYrho) / w(i, j, k, irho) + A14 * (-velZrho) / w(i, j, k, irho) & + / w(i, j, k, irho) + A13 * (-velYrho) / w(i, j, k, irho) + & + A14 * (-velZrho) / w(i, j, k, irho) & + A15 * (((gamma(i, j, k) - 1) * q / 2) - SoS**2) B12 = A11 * (1 - gamma(i, j, k)) * velXrho + A12 * 1 / w(i, j, k, irho) & + A15 * (1 - gamma(i, j, k)) * velXrho @@ -281,7 +282,8 @@ subroutine residual_block B25 = A21 * (gamma(i, j, k) - 1) + A25 * (gamma(i, j, k) - 1) B31 = A31 * (gamma(i, j, k) - 1) * q / 2 + A32 * (-velXrho) & - / w(i, j, k, irho) + A33 * (-velYrho) / w(i, j, k, irho) + A34 * (-velZrho) / w(i, j, k, irho) & + / w(i, j, k, irho) + A33 * (-velYrho) / w(i, j, k, irho) + & + A34 * (-velZrho) / w(i, j, k, irho) & + A35 * (((gamma(i, j, k) - 1) * q / 2) - SoS**2) B32 = A31 * (1 - gamma(i, j, k)) * velXrho + A32 & / w(i, j, k, irho) + A35 * (1 - gamma(i, j, k)) * velXrho @@ -1094,7 +1096,8 @@ subroutine computedwDADI viscTermI, viscTermJ, viscTermK real(kind=realType), dimension(:, :, :), pointer :: qq_i, qq_j, qq_k, & - cc_i, cc_j, cc_k, spectral_i, spectral_j, spectral_k, dual_dt + cc_i, cc_j, cc_k, spectral_i, & + spectral_j, spectral_k, dual_dt real(kind=realType), dimension(ie, 5) :: bbi, cci, ddi, ffi real(kind=realType), dimension(je, 5) :: bbj, ccj, ddj, ffj real(kind=realType), dimension(ke, 5) :: bbk, cck, ddk, ffk @@ -1380,7 +1383,9 @@ subroutine computedwDADI do j = 2, jl bbj(j, n) = bbj(j, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realType), zero) ddj(j, n) = ddj(j, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realType), zero) - ccj(j,n)=one+ccj(j,n)*dual_dt(i,j,k)*max(real(iblank(i,j,k),realType),zero)+spectral_i(i,j,k)+spectral_k(i,j,k) + ccj(j, n) = one + ccj(j, n) * dual_dt(i, j, k) * & + max(real(iblank(i, j, k), realType), zero) + & + spectral_i(i, j, k) + spectral_k(i, j, k) ffj(j, n) = dw(i, j, k, n) end do end do @@ -1514,7 +1519,9 @@ subroutine computedwDADI do i = 2, il bbi(i, n) = bbi(i, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realType), zero) ddi(i, n) = ddi(i, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realType), zero) - cci(i,n)=one+cci(i,n)*dual_dt(i,j,k)*max(real(iblank(i,j,k),realType),zero)+spectral_j(i,j,k)+spectral_k(i,j,k) + cci(i, n) = one + cci(i, n) * dual_dt(i, j, k) * & + max(real(iblank(i, j, k), realType), zero) + & + spectral_j(i, j, k) + spectral_k(i, j, k) ffi(i, n) = dw(i, j, k, n) end do end do @@ -1648,7 +1655,9 @@ subroutine computedwDADI do k = 2, kl bbk(k, n) = bbk(k, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realType), zero) ddk(k, n) = ddk(k, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realType), zero) - cck(k,n)=one+cck(k,n)*dual_dt(i,j,k)*max(real(iblank(i,j,k),realType),zero)+spectral_i(i,j,k)+spectral_j(i,j,k) + cck(k, n) = one + cck(k, n) * dual_dt(i, j, k) * & + max(real(iblank(i, j, k), realType), zero) + & + spectral_i(i, j, k) + spectral_j(i, j, k) ffk(k, n) = dw(i, j, k, n) end do end do diff --git a/src/solver/solverUtils.F90 b/src/solver/solverUtils.F90 index c5ac50c6a..2136d930e 100644 --- a/src/solver/solverUtils.F90 +++ b/src/solver/solverUtils.F90 @@ -446,20 +446,32 @@ subroutine gridVelocitiesFineLevel_TS_block(nn, sps) do j = 1, je_l do i = 1, ie_l - x_vc = eighth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k - 1, 1) + flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 1) & - + flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 1) + flowDoms(nn, 1, mm)%x(i, j, k - 1, 1) & - + flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 1) + flowDoms(nn, 1, mm)%x(i, j - 1, k, 1) & - + flowDoms(nn, 1, mm)%x(i - 1, j, k, 1) + flowDoms(nn, 1, mm)%x(i, j, k, 1)) - - y_vc = eighth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k - 1, 2) + flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 2) & - + flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 2) + flowDoms(nn, 1, mm)%x(i, j, k - 1, 2) & - + flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 2) + flowDoms(nn, 1, mm)%x(i, j - 1, k, 2) & - + flowDoms(nn, 1, mm)%x(i - 1, j, k, 2) + flowDoms(nn, 1, mm)%x(i, j, k, 2)) - - z_vc = eighth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k - 1, 3) + flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 3) & - + flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 3) + flowDoms(nn, 1, mm)%x(i, j, k - 1, 3) & - + flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 3) + flowDoms(nn, 1, mm)%x(i, j - 1, k, 3) & - + flowDoms(nn, 1, mm)%x(i - 1, j, k, 3) + flowDoms(nn, 1, mm)%x(i, j, k, 3)) + x_vc = eighth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k - 1, 1) + & + flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 1) & + + flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 1) + & + flowDoms(nn, 1, mm)%x(i, j, k - 1, 1) & + + flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 1) + & + flowDoms(nn, 1, mm)%x(i, j - 1, k, 1) & + + flowDoms(nn, 1, mm)%x(i - 1, j, k, 1) + & + flowDoms(nn, 1, mm)%x(i, j, k, 1)) + + y_vc = eighth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k - 1, 2) + & + flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 2) & + + flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 2) + & + flowDoms(nn, 1, mm)%x(i, j, k - 1, 2) & + + flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 2) + & + flowDoms(nn, 1, mm)%x(i, j - 1, k, 2) & + + flowDoms(nn, 1, mm)%x(i - 1, j, k, 2) + & + flowDoms(nn, 1, mm)%x(i, j, k, 2)) + + z_vc = eighth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k - 1, 3) + & + flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 3) & + + flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 3) + & + flowDoms(nn, 1, mm)%x(i, j, k - 1, 3) & + + flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 3) + & + flowDoms(nn, 1, mm)%x(i, j - 1, k, 3) & + + flowDoms(nn, 1, mm)%x(i - 1, j, k, 3) + & + flowDoms(nn, 1, mm)%x(i, j, k, 3)) s(i, j, k, 1) = s(i, j, k, 1) + dscalar(1, sps, mm) * x_vc s(i, j, k, 2) = s(i, j, k, 2) + dscalar(1, sps, mm) * y_vc @@ -533,12 +545,18 @@ subroutine gridVelocitiesFineLevel_TS_block(nn, sps) do j = 1, je_l do i = 0, ie_l - x_fc = fourth * (flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 1) + flowDoms(nn, 1, mm)%x(i, j, k, 1) & - + flowDoms(nn, 1, mm)%x(i, j - 1, k, 1) + flowDoms(nn, 1, mm)%x(i, j, k - 1, 1)) - y_fc = fourth * (flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 2) + flowDoms(nn, 1, mm)%x(i, j, k, 2) & - + flowDoms(nn, 1, mm)%x(i, j - 1, k, 2) + flowDoms(nn, 1, mm)%x(i, j, k - 1, 2)) - z_fc = fourth * (flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 3) + flowDoms(nn, 1, mm)%x(i, j, k, 3) & - + flowDoms(nn, 1, mm)%x(i, j - 1, k, 3) + flowDoms(nn, 1, mm)%x(i, j, k - 1, 3)) + x_fc = fourth * (flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 1) + & + flowDoms(nn, 1, mm)%x(i, j, k, 1) & + + flowDoms(nn, 1, mm)%x(i, j - 1, k, 1) + & + flowDoms(nn, 1, mm)%x(i, j, k - 1, 1)) + y_fc = fourth * (flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 2) + & + flowDoms(nn, 1, mm)%x(i, j, k, 2) & + + flowDoms(nn, 1, mm)%x(i, j - 1, k, 2) + & + flowDoms(nn, 1, mm)%x(i, j, k - 1, 2)) + z_fc = fourth * (flowDoms(nn, 1, mm)%x(i, j - 1, k - 1, 3) + & + flowDoms(nn, 1, mm)%x(i, j, k, 3) & + + flowDoms(nn, 1, mm)%x(i, j - 1, k, 3) + & + flowDoms(nn, 1, mm)%x(i, j, k - 1, 3)) sFaceI(i, j, k) = sFaceI(i, j, k) & + dscalar(1, sps, mm) * x_fc * sI(i, j, k, 1) & @@ -554,12 +572,18 @@ subroutine gridVelocitiesFineLevel_TS_block(nn, sps) do j = 0, je_l do i = 1, ie_l - x_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 1) + flowDoms(nn, 1, mm)%x(i, j, k, 1) & - + flowDoms(nn, 1, mm)%x(i - 1, j, k, 1) + flowDoms(nn, 1, mm)%x(i, j, k - 1, 1)) - y_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 2) + flowDoms(nn, 1, mm)%x(i, j, k, 2) & - + flowDoms(nn, 1, mm)%x(i - 1, j, k, 2) + flowDoms(nn, 1, mm)%x(i, j, k - 1, 2)) - z_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 3) + flowDoms(nn, 1, mm)%x(i, j, k, 3) & - + flowDoms(nn, 1, mm)%x(i - 1, j, k, 3) + flowDoms(nn, 1, mm)%x(i, j, k - 1, 3)) + x_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 1) + & + flowDoms(nn, 1, mm)%x(i, j, k, 1) & + + flowDoms(nn, 1, mm)%x(i - 1, j, k, 1) + & + flowDoms(nn, 1, mm)%x(i, j, k - 1, 1)) + y_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 2) + & + flowDoms(nn, 1, mm)%x(i, j, k, 2) & + + flowDoms(nn, 1, mm)%x(i - 1, j, k, 2) + & + flowDoms(nn, 1, mm)%x(i, j, k - 1, 2)) + z_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j, k - 1, 3) + & + flowDoms(nn, 1, mm)%x(i, j, k, 3) & + + flowDoms(nn, 1, mm)%x(i - 1, j, k, 3) + & + flowDoms(nn, 1, mm)%x(i, j, k - 1, 3)) sFaceJ(i, j, k) = sFaceJ(i, j, k) & + dscalar(1, sps, mm) * x_fc * sJ(i, j, k, 1) & @@ -575,12 +599,18 @@ subroutine gridVelocitiesFineLevel_TS_block(nn, sps) do j = 1, je_l do i = 1, ie_l - x_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 1) + flowDoms(nn, 1, mm)%x(i, j, k, 1) & - + flowDoms(nn, 1, mm)%x(i, j - 1, k, 1) + flowDoms(nn, 1, mm)%x(i - 1, j, k, 1)) - y_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 2) + flowDoms(nn, 1, mm)%x(i, j, k, 2) & - + flowDoms(nn, 1, mm)%x(i, j - 1, k, 2) + flowDoms(nn, 1, mm)%x(i - 1, j, k, 2)) - z_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 3) + flowDoms(nn, 1, mm)%x(i, j, k, 3) & - + flowDoms(nn, 1, mm)%x(i, j - 1, k, 3) + flowDoms(nn, 1, mm)%x(i - 1, j, k, 3)) + x_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 1) + & + flowDoms(nn, 1, mm)%x(i, j, k, 1) & + + flowDoms(nn, 1, mm)%x(i, j - 1, k, 1) + & + flowDoms(nn, 1, mm)%x(i - 1, j, k, 1)) + y_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 2) + & + flowDoms(nn, 1, mm)%x(i, j, k, 2) & + + flowDoms(nn, 1, mm)%x(i, j - 1, k, 2) + & + flowDoms(nn, 1, mm)%x(i - 1, j, k, 2)) + z_fc = fourth * (flowDoms(nn, 1, mm)%x(i - 1, j - 1, k, 3) + & + flowDoms(nn, 1, mm)%x(i, j, k, 3) & + + flowDoms(nn, 1, mm)%x(i, j - 1, k, 3) + & + flowDoms(nn, 1, mm)%x(i - 1, j, k, 3)) sFaceK(i, j, k) = sFaceK(i, j, k) & + dscalar(1, sps, mm) * x_fc * sK(i, j, k, 1) & @@ -1066,14 +1096,21 @@ subroutine gridVelocitiesFineLevel_block(useOldCoor, t, sps, nn) ! Determine the coordinates of the face center, ! which are stored in xc. - xc(1) = fourth*(flowDoms(nn, groundLevel, sps)%x(i,j-1,k-1,1) + flowDoms(nn, groundLevel, sps)%x(i,j,k-1,1) & - + flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 1) + flowDoms(nn, groundLevel, sps)%x(i, j, k, 1)) - xc(2) = fourth*(flowDoms(nn, groundLevel, sps)%x(i,j-1,k-1,2) + flowDoms(nn, groundLevel, sps)%x(i,j,k-1,2) & - + flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 2) + flowDoms(nn, groundLevel, sps)%x(i, j, k, 2)) - xc(3) = fourth*(flowDoms(nn, groundLevel, sps)%x(i,j-1,k-1,3) + flowDoms(nn, groundLevel, sps)%x(i,j,k-1,3) & - + flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 3) + flowDoms(nn, groundLevel, sps)%x(i, j, k, 3)) - - call cellFaceVelocities(xc, rotCenter, rotRate, velxGrid, velyGrid, velzGrid, derivRotationMatrix, sc) + xc(1) = fourth * (flowDoms(nn, groundLevel, sps)%x(i, j - 1, k - 1, 1) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k - 1, 1) & + + flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 1) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k, 1)) + xc(2) = fourth * (flowDoms(nn, groundLevel, sps)%x(i, j - 1, k - 1, 2) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k - 1, 2) & + + flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 2) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k, 2)) + xc(3) = fourth * (flowDoms(nn, groundLevel, sps)%x(i, j - 1, k - 1, 3) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k - 1, 3) & + + flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 3) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k, 3)) + + call cellFaceVelocities(xc, rotCenter, rotRate, & + velxGrid, velyGrid, velzGrid, derivRotationMatrix, sc) ! Store the dot product of grid velocity sc and ! the normal ss in sFace. @@ -1092,14 +1129,21 @@ subroutine gridVelocitiesFineLevel_block(useOldCoor, t, sps, nn) ! Determine the coordinates of the face center, ! which are stored in xc. - xc(1) = fourth * (flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 1) + flowDoms(nn, groundLevel, sps)%x(i, j, k - 1, 1) & - + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k - 1, 1) + flowDoms(nn, groundLevel, sps)%x(i, j, k, 1)) - xc(2) = fourth * (flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 2) + flowDoms(nn, groundLevel, sps)%x(i, j, k - 1, 2) & - + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k - 1, 2) + flowDoms(nn, groundLevel, sps)%x(i, j, k, 2)) - xc(3) = fourth * (flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 3) + flowDoms(nn, groundLevel, sps)%x(i, j, k - 1, 3) & - + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k - 1, 3) + flowDoms(nn, groundLevel, sps)%x(i, j, k, 3)) - - call cellFaceVelocities(xc, rotCenter, rotRate, velxGrid, velyGrid, velzGrid, derivRotationMatrix, sc) + xc(1) = fourth * (flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 1) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k - 1, 1) & + + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k - 1, 1) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k, 1)) + xc(2) = fourth * (flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 2) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k - 1, 2) & + + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k - 1, 2) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k, 2)) + xc(3) = fourth * (flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 3) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k - 1, 3) & + + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k - 1, 3) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k, 3)) + + call cellFaceVelocities(xc, rotCenter, rotRate, & + velxGrid, velyGrid, velzGrid, derivRotationMatrix, sc) ! Store the dot product of grid velocity sc and ! the normal ss in sFace. @@ -1118,14 +1162,21 @@ subroutine gridVelocitiesFineLevel_block(useOldCoor, t, sps, nn) ! Determine the coordinates of the face center, ! which are stored in xc. - xc(1) = fourth * (flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 1) + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 1) & - + flowDoms(nn, groundLevel, sps)%x(i - 1, j - 1, k, 1) + flowDoms(nn, groundLevel, sps)%x(i, j, k, 1)) - xc(2) = fourth * (flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 2) + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 2) & - + flowDoms(nn, groundLevel, sps)%x(i - 1, j - 1, k, 2) + flowDoms(nn, groundLevel, sps)%x(i, j, k, 2)) - xc(3) = fourth * (flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 3) + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 3) & - + flowDoms(nn, groundLevel, sps)%x(i - 1, j - 1, k, 3) + flowDoms(nn, groundLevel, sps)%x(i, j, k, 3)) - - call cellFaceVelocities(xc, rotCenter, rotRate, velxGrid, velyGrid, velzGrid, derivRotationMatrix, sc) + xc(1) = fourth * (flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 1) + & + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 1) & + + flowDoms(nn, groundLevel, sps)%x(i - 1, j - 1, k, 1) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k, 1)) + xc(2) = fourth * (flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 2) + & + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 2) & + + flowDoms(nn, groundLevel, sps)%x(i - 1, j - 1, k, 2) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k, 2)) + xc(3) = fourth * (flowDoms(nn, groundLevel, sps)%x(i, j - 1, k, 3) + & + flowDoms(nn, groundLevel, sps)%x(i - 1, j, k, 3) & + + flowDoms(nn, groundLevel, sps)%x(i - 1, j - 1, k, 3) + & + flowDoms(nn, groundLevel, sps)%x(i, j, k, 3)) + + call cellFaceVelocities(xc, rotCenter, rotRate, & + velxGrid, velyGrid, velzGrid, derivRotationMatrix, sc) ! Store the dot product of grid velocity sc and ! the normal ss in sFace. diff --git a/src/solver/surfaceIntegrations.F90 b/src/solver/surfaceIntegrations.F90 index d39af2d5f..bde08441d 100644 --- a/src/solver/surfaceIntegrations.F90 +++ b/src/solver/surfaceIntegrations.F90 @@ -23,7 +23,8 @@ subroutine getCostFunctions(globalVals, funcValues) ! Working real(kind=realType) :: fact, factMoment, ovrNTS real(kind=realType), dimension(3, nTimeIntervalsSpectral) :: force, forceP, forceV, forceM, & - moment, cForce, cForceP, cForceV, cForceM, cMoment, coFx, coFy, coFz + moment, cForce, cForceP, cForceV, cForceM, & + cMoment, coFx, coFy, coFz real(kind=realType), dimension(3) :: VcoordRef, VFreestreamRef real(kind=realType) :: mAvgPtot, mAvgTtot, mAvgRho, mAvgPs, mFlow, mAvgMn, mAvga, & mAvgVx, mAvgVy, mAvgVz, gArea, mAvgVi @@ -161,8 +162,10 @@ subroutine getCostFunctions(globalVals, funcValues) funcValues(costFuncAxisMoment) = funcValues(costFuncAxisMoment) + ovrNTS * globalVals(iAxisMoment, sps) funcValues(costFuncSepSensorAvgX) = funcValues(costFuncSepSensorAvgX) + ovrNTS * globalVals(iSepAvg, sps) - funcValues(costFuncSepSensorAvgY) = funcValues(costFuncSepSensorAvgY) + ovrNTS * globalVals(iSepAvg + 1, sps) - funcValues(costFuncSepSensorAvgZ) = funcValues(costFuncSepSensorAvgZ) + ovrNTS * globalVals(iSepAvg + 2, sps) + funcValues(costFuncSepSensorAvgY) = funcValues(costFuncSepSensorAvgY) + & + ovrNTS * globalVals(iSepAvg + 1, sps) + funcValues(costFuncSepSensorAvgZ) = funcValues(costFuncSepSensorAvgZ) + & + ovrNTS * globalVals(iSepAvg + 2, sps) funcValues(costFuncArea) = funcValues(costFuncArea) + ovrNTS * globalVals(iArea, sps) funcValues(costFuncFlowPower) = funcValues(costFuncFlowPower) + ovrNTS * globalVals(iPower, sps) @@ -206,8 +209,10 @@ subroutine getCostFunctions(globalVals, funcValues) gArea = globalVals(iArea, sps) if (gArea /= zero) then ! area averaged pressure - funcValues(costFuncAAvgPTot) = funcValues(costFuncAAvgPTot) + ovrNTS * globalVals(iAreaPTot, sps) / gArea - funcValues(costFuncAAvgPs) = funcValues(costFuncAAvgPs) + ovrNTS * globalVals(iAreaPs, sps) / gArea + funcValues(costFuncAAvgPTot) = funcValues(costFuncAAvgPTot) + & + ovrNTS * globalVals(iAreaPTot, sps) / gArea + funcValues(costFuncAAvgPs) = funcValues(costFuncAAvgPs) + & + ovrNTS * globalVals(iAreaPs, sps) / gArea end if funcValues(costFuncMdot) = funcValues(costFuncMdot) + ovrNTS * mFlow diff --git a/src/solver/userSurfaceIntegrations.F90 b/src/solver/userSurfaceIntegrations.F90 index af72b0070..760cca989 100644 --- a/src/solver/userSurfaceIntegrations.F90 +++ b/src/solver/userSurfaceIntegrations.F90 @@ -183,7 +183,8 @@ subroutine integrateUserSurfaces_d(localValues, localValuesd, famList, sps) flowDoms(nn, 1, sps)%realCommVars(iVy + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%w(:, :, :, iVy) flowDoms(nn, 1, sps)%realCommVars(iVz + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%w(:, :, :, iVz) flowDoms(nn, 1, sps)%realCommVars(iZippFlowP + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%P(:, :, :) - flowDoms(nn, 1, sps)%realCommVars(iZippFlowGamma + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%gamma(:, :, :) + flowDoms(nn, 1, sps)%realCommVars(iZippFlowGamma + nZippFlowComm)%var => & + flowDomsd(nn, 1, sps)%gamma(:, :, :) ! flowDoms(nn, 1, sps)%realCommVars(iZippFlowSface+nZippFlowComm)%var => Not Implemented flowDoms(nn, 1, sps)%realCommVars(iZippFlowX + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%x(:, :, :, 1) @@ -266,7 +267,8 @@ subroutine integrateUserSurfaces_d(localValues, localValuesd, famList, sps) fams = surf%famID ! Perform the actual integration - call flowIntegrationZipper_d(surf%isInflow, surf%conn, fams, vars, varsd, localValues, localValuesd, & + call flowIntegrationZipper_d(surf%isInflow, surf%conn, fams, vars, & + varsd, localValues, localValuesd, & famList, sps, ptValid) deallocate (ptValid, vars, varsd, fams) end if @@ -333,7 +335,8 @@ subroutine integrateUserSurfaces_b(localValues, localValuesd, famList, sps) flowDoms(nn, 1, sps)%realCommVars(iVy + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%w(:, :, :, iVy) flowDoms(nn, 1, sps)%realCommVars(iVz + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%w(:, :, :, iVz) flowDoms(nn, 1, sps)%realCommVars(iZippFlowP + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%P(:, :, :) - flowDoms(nn, 1, sps)%realCommVars(iZippFlowGamma + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%gamma(:, :, :) + flowDoms(nn, 1, sps)%realCommVars(iZippFlowGamma + nZippFlowComm)%var => & + flowDomsd(nn, 1, sps)%gamma(:, :, :) ! flowDoms(nn, 1, sps)%realCommVars(iZippFlowSface+nZippFlowComm)%var => Not Implemented flowDoms(nn, 1, sps)%realCommVars(iZippFlowX + nZippFlowComm)%var => flowDomsd(nn, 1, sps)%x(:, :, :, 1) @@ -413,7 +416,8 @@ subroutine integrateUserSurfaces_b(localValues, localValuesd, famList, sps) fams = surf%famID ! Perform the actual (reverse) integration - call flowIntegrationZipper_b(surf%isInflow, surf%conn, fams, vars, varsd, localValues, localValuesd, & + call flowIntegrationZipper_b(surf%isInflow, surf%conn, fams, vars, & + varsd, localValues, localValuesd, & famList, sps, ptValid) ! Accumulate into the receive buffers @@ -765,7 +769,8 @@ subroutine performInterpolation(pts, oBlocks, useDual, comm) ! Call the standard tree search call containmentTreeSearchSinglePoint(oBlock%ADT, xc, intInfo, uvw, & - oBlock%qualDonor, nInterpol, BB, frontLeaves, frontLeavesNew, failed) + oBlock%qualDonor, nInterpol, BB, & + frontLeaves, frontLeavesNew, failed) ! Make sure this point is not garbage. if (intInfo(1) >= 0) then @@ -1375,28 +1380,36 @@ subroutine commUserIntegrationSurfaceVars_b(recvBuffer, recvBufferd, varStart, v ! Accumulate back onto the derivative variables. flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1, k1) = & - flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1, k1) + weight(1) * sendBufferd(jj) + flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1, k1) + & + weight(1) * sendBufferd(jj) flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1, k1) = & - flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1, k1) + weight(2) * sendBufferd(jj) + flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1, k1) + & + weight(2) * sendBufferd(jj) flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1 + 1, k1) = & - flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1 + 1, k1) + weight(3) * sendBufferd(jj) + flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1 + 1, k1) + & + weight(3) * sendBufferd(jj) flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1 + 1, k1) = & - flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1 + 1, k1) + weight(4) * sendBufferd(jj) + flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1 + 1, k1) + & + weight(4) * sendBufferd(jj) flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1, k1 + 1) = & - flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1, k1 + 1) + weight(5) * sendBufferd(jj) + flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1, k1 + 1) + & + weight(5) * sendBufferd(jj) flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1, k1 + 1) = & - flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1, k1 + 1) + weight(6) * sendBufferd(jj) + flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1, k1 + 1) + & + weight(6) * sendBufferd(jj) flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1 + 1, k1 + 1) = & - flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1 + 1, k1 + 1) + weight(7) * sendBufferd(jj) + flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1, j1 + 1, k1 + 1) + & + weight(7) * sendBufferd(jj) flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1 + 1, k1 + 1) = & - flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1 + 1, k1 + 1) + weight(8) * sendBufferd(jj) + flowDoms(d1, 1, 1)%realCommVars(k + nZippFlowComm)%var(i1 + 1, j1 + 1, k1 + 1) + & + weight(8) * sendBufferd(jj) end if end do end do donorLoop diff --git a/src/utils/haloExchange.F90 b/src/utils/haloExchange.F90 index 275260815..11d9c2665 100644 --- a/src/utils/haloExchange.F90 +++ b/src/utils/haloExchange.F90 @@ -3448,7 +3448,8 @@ subroutine flowIntegrationZipperComm_b(isInflow, vars, varsd, sps) localPtr(ii) = zero end if case (iZippFlowX, iZippFlowY, iZippFlowZ) - xxd(i + 1, j + 1, iVar - iZippFlowX + 1) = xxd(i + 1, j + 1, iVar - iZippFlowX + 1) + localPtr(ii) + xxd(i + 1, j + 1, iVar - iZippFlowX + 1) = & + xxd(i + 1, j + 1, iVar - iZippFlowX + 1) + localPtr(ii) end select end do end do @@ -3734,7 +3735,8 @@ subroutine wallIntegrationZipperComm_b(vars, varsd, sps) case (iZippWallX, iZippWallY, iZippWallZ) ! The +1 is due to pointer offset - xxd(i + 1, j + 1, iVar - iZippWallX + 1) = xxd(i + 1, j + 1, iVar - iZippWallX + 1) + localPtr(ii) + xxd(i + 1, j + 1, iVar - iZippWallX + 1) = & + xxd(i + 1, j + 1, iVar - iZippWallX + 1) + localPtr(ii) end select end do diff --git a/src/utils/utils.F90 b/src/utils/utils.F90 index dd44fd2ce..1652107c9 100644 --- a/src/utils/utils.F90 +++ b/src/utils/utils.F90 @@ -1195,13 +1195,21 @@ subroutine computeRootBendingMoment(cf, cm, bendingMoment) bendingMoment = zero if (liftIndex == 2) then !z out wing sum momentx,momentz - elasticMomentx = cm(1) + cf(2)*(pointRefEC(3)-pointRef(3))/lengthref-cf(3)*(pointRefEC(2)-pointRef(2))/lengthref - elasticMomentz = cm(3) - cf(2)*(pointRefEC(1)-pointref(1))/lengthref+cf(1)*(pointRefEC(2)-pointRef(2))/lengthref + elasticMomentx = cm(1) + cf(2) * (pointRefEC(3) - & + pointRef(3)) / lengthref - cf(3) * & + (pointRefEC(2) - pointRef(2)) / lengthref + elasticMomentz = cm(3) - cf(2) * (pointRefEC(1) - & + pointref(1)) / lengthref + cf(1) * & + (pointRefEC(2) - pointRef(2)) / lengthref bendingMoment = sqrt(elasticMomentx**2 + elasticMomentz**2) elseif (liftIndex == 3) then !y out wing sum momentx,momenty - elasticMomentx = cm(1) + cf(3)*(pointrefEC(2)-pointRef(2))/lengthref+cf(3)*(pointrefEC(3)-pointref(3))/lengthref - elasticMomenty = cm(2) + cf(3)*(pointRefEC(1)-pointRef(1))/lengthref+cf(1)*(pointrefEC(3)-pointRef(3))/lengthref + elasticMomentx = cm(1) + cf(3) * (pointrefEC(2) - & + pointRef(2)) / lengthref + & + cf(3) * (pointrefEC(3) - pointref(3)) / lengthref + elasticMomenty = cm(2) + cf(3) * (pointRefEC(1) - & + pointRef(1)) / lengthref + & + cf(1) * (pointrefEC(3) - pointRef(3)) / lengthref bendingMoment = sqrt(elasticMomentx**2 + elasticMomenty**2) end if @@ -1388,7 +1396,8 @@ subroutine computeTSDerivatives(force, moment, coef0, dcdalpha, & !now compute dCl/dalpha do i = 1, 8 - call computeLeastSquaresRegression(BaseCoef(:, i), intervalAlpha, nTimeIntervalsSpectral, dcdAlpha(i), coef0(i)) + call computeLeastSquaresRegression(BaseCoef(:, i), & + intervalAlpha, nTimeIntervalsSpectral, dcdAlpha(i), coef0(i)) end do ! now subtract off estimated cl,cmz and use remainder to compute @@ -1401,7 +1410,9 @@ subroutine computeTSDerivatives(force, moment, coef0, dcdalpha, & !now compute dCi/dalphadot do i = 1, 8 - call computeLeastSquaresRegression(ResBaseCoef(:,i),intervalAlphadot,nTimeIntervalsSpectral,dcdalphadot(i),Coef0dot(i)) + call computeLeastSquaresRegression(ResBaseCoef(:, i), & + intervalAlphadot, nTimeIntervalsSpectral, & + dcdalphadot(i), Coef0dot(i)) end do a = sqrt(gammaInf * pInfDim / rhoInfDim) diff --git a/src/wallDistance/wallDistance.F90 b/src/wallDistance/wallDistance.F90 index 2017fb5da..965f4b3d3 100644 --- a/src/wallDistance/wallDistance.F90 +++ b/src/wallDistance/wallDistance.F90 @@ -1860,7 +1860,8 @@ subroutine determineWallAssociation(level, sps) coor(4) = large call minDistancetreeSearchSinglePoint(walls(c)%ADT, coor, & - intInfo2, uvw2, dummy, 0, BB, frontLeaves, frontLeavesNew) + intInfo2, uvw2, dummy, 0, BB, & + frontLeaves, frontLeavesNew) cellID2 = intInfo2(3) if (uvw2(4) < nearWallDist**2) then diff --git a/tests/reg_tests/test_functionals.py b/tests/reg_tests/test_functionals.py index 990030700..4631ee626 100644 --- a/tests/reg_tests/test_functionals.py +++ b/tests/reg_tests/test_functionals.py @@ -36,6 +36,7 @@ "ref_file": "funcs_euler_scalar_jst_tut_wing.json", "aero_prob": copy.deepcopy(ap_tutorial_wing), "N_PROCS": 1, + "dot_prod_tol": 1e-10, }, # scalar JST { @@ -51,6 +52,7 @@ }, "ref_file": "funcs_euler_scalar_jst_tut_wing.json", "aero_prob": ap_tutorial_wing, + "dot_prod_tol": 1e-10, }, # Matrix JST { @@ -70,6 +72,7 @@ }, "ref_file": "funcs_euler_matrix_jst_tut_wing.json", "aero_prob": ap_tutorial_wing, + "dot_prod_tol": 1e-10, }, # Upwind { @@ -88,6 +91,7 @@ }, "ref_file": "funcs_euler_upwind_tut_wing.json", "aero_prob": ap_tutorial_wing, + "dot_prod_tol": 1e-10, }, # Tutorial wing random block order { @@ -103,6 +107,7 @@ }, "ref_file": "funcs_euler_scalar_jst_rand_tut_wing.json", "aero_prob": ap_tutorial_wing, + "dot_prod_tol": 1e-10, }, # Tutorial wing laminar { @@ -124,6 +129,7 @@ }, "ref_file": "funcs_laminar_tut_wing.json", "aero_prob": ap_tutorial_wing_laminar, + "dot_prod_tol": 1e-10, }, # Tutorial wing RANS { @@ -151,6 +157,7 @@ }, "ref_file": "funcs_rans_tut_wing.json", "aero_prob": ap_tutorial_wing, + "dot_prod_tol": 2e-10, }, # Tutorial wing random RANS { @@ -178,6 +185,7 @@ }, "ref_file": "funcs_rans_rand_tut_wing.json", "aero_prob": ap_tutorial_wing, + "dot_prod_tol": 1e-10, }, # CRM WBT { @@ -201,6 +209,7 @@ }, "ref_file": "funcs_euler_scalar_jst_CRM_WBT.json", "aero_prob": ap_CRM, + "dot_prod_tol": 1e-10, }, ] ) @@ -266,7 +275,7 @@ def test_jac_vec_prod_bwd(self): utils.assert_bwd_mode_allclose(self.handler, self.CFDSolver, self.ap, tol=1e-10) def test_dot_products(self): - utils.assert_dot_products_allclose(self.handler, self.CFDSolver, tol=1e-10) + utils.assert_dot_products_allclose(self.handler, self.CFDSolver, tol=self.dot_prod_tol) if __name__ == "__main__":