diff --git a/idwarp/AxisymmetricMesh.py b/idwarp/AxisymmetricMesh.py new file mode 100644 index 0000000..aa4e079 --- /dev/null +++ b/idwarp/AxisymmetricMesh.py @@ -0,0 +1,774 @@ +import os +import shutil +import warnings +from dataclasses import dataclass + +import numpy as np +from baseclasses import BaseSolver +from baseclasses.utils import Error +from mpi4py import MPI +from numpy.typing import ArrayLike + +from .MExt import MExt + + +@dataclass(eq=True, unsafe_hash=True, frozen=True) +class MeshPatch: + idxStart: int + idxEnd: int + point: ArrayLike + normal: ArrayLike + nodeSize: int + cellSize: int + patchName: str + + +class AxisymmetricMesh(BaseSolver): + def __init__(self, options, comm=None, debug=False): + """ + Create the AxisymmetricMesh object. + + Parameters + ---------- + options : dictionary + A dictionary containing the options for the mesh movement + strategies. + + comm : MPI_INTRA_COMM + MPI communication (as obtained from mpi4py) on which to create + the AxisymmetricMesh object. If not provided, MPI_COMM_WORLD is used. + + debug : bool + Flag specifying if the MExt import is automatically deleted. + This needs to be true ONLY when a symbolic debugger is used. + """ + name = "IDWarp" + category = "Volume mesh warping" + + if comm is None: + comm = MPI.COMM_WORLD + + # Get the default options + defOpts = self._getDefaultOptions() + + # Initialize the base class + super().__init__(name, category, defaultOptions=defOpts, options=options, comm=comm) + + self.printOptions() + + # Check if warp was already set (complex inheritance issue) + if not hasattr(self, "warp"): + curDir = os.path.basename(os.path.dirname(os.path.realpath(__file__))) + self.warp = MExt("libidwarp", curDir, debug=debug)._module + + # Initialize PETSc + self.warp.initpetsc(self.comm.py2f()) + + # Set the dtype to 'd' + self.dtype = "d" + + # Set the fortran option values + self._setMeshOptions() + + # Initialize class attributes + self.warpInitialized = False + self.fileType = self.getOption("fileType") + fileName = self.getOption("gridFile") + + if not self.fileType == "CGNS": + raise Error("Only CGNS file type is supported for axisymmetric mesh warping") + + self.warp.readcgns(fileName) + + if not self.warp.cgnsgrid.cgnsstructured: + raise Error("Only structured CGNS files are supported for axisymmetric mesh warping") + + # Check the boco types on the root proc + if self.comm.rank == 0: + if 2 not in self.warp.cgnsgrid.bocoTypes and 17 not in self.warp.cgnsgrid.bocoTypes: + raise Error("No 'wedge' nor 'symmPolar' boundary condition found in the CGNS file") + + @staticmethod + def _getDefaultOptions(): + defOpts = { + # General options + "gridFile": [(str, type(None)), None], + "fileType": [str, ["CGNS"]], + "specifiedSurfaces": [(list, type(None)), None], + "aExp": [float, 3.0], + "bExp": [float, 5.0], + "LdefFact": [float, 1.0], + "alpha": [float, 0.25], + "errTol": [float, 0.0005], + "evalMode": [str, ["fast", "exact"]], + "symmTol": [float, 1e-6], + "useRotations": [bool, True], + "zeroCornerRotations": [bool, True], + "cornerAngle": [float, 30.0], + "restartFile": [(str, type(None)), None], + "bucketSize": [int, 8], + # Axisymmetric only options + "axiSymmetric": [bool, True], + "axiAngle": [float, 0.0], + "axiAxis": [list, [1.0, 0.0, 0.0]], + "axiPlane": [(list, type(None)), None], + "axiMirrorFamily": [(str, type(None)), None], + } + return defOpts + + def getSurfaceCoordinates(self): + """Returns all defined surface coordinates on this processor + + Returns + ------- + coords : numpy array size (N,3) + Specified surface coordinates residing on this + processor. This may be empty array, size (0,3) + """ + self._setInternalSurface() + coords = np.zeros((self.nSurf, 3), self.dtype) + self.warp.getsurfacecoordinates(np.ravel(coords)) + + return coords + + def setSurfaceCoordinates(self, coordinates): + """Sets all surface coordinates on this processor + + Parameters + ---------- + coordinates : numpy array, size(N, 3) + The coordinate to set. This MUST be exactly the same size as + the array obtained from getSurfaceCoordinates() + """ + if len(coordinates) != self.nSurf: + raise Error( + "Incorrect length of coordinates supplied to " + "setSurfaceCoordinates on proc %d. Expected " + "aray of length %d, received an array of length " + "%d." % (self.comm.rank, self.nSurf, len(coordinates)) + ) + + self.warp.setsurfacecoordinates(np.ravel(coordinates)) + + def setExternalMeshIndices(self, ind): + """ + Set the indicies defining the transformation of an external + solver grid to the original CGNS grid. This is required to use + USMesh functions that involve the word "Solver" and warpDeriv. + The indices must be zero-based. + + Parameters + ---------- + ind : numpy integer array + The list of indicies this processor needs from the common mesh file + """ + + self.warp.setexternalmeshindices(ind) + + def setSurfaceDefinition(self, pts, conn, faceSizes, cgnsBlockID=None): + """This is the master function that determines the definition of the + surface to be used for the mesh movement. This surface may be + supplied from an external solver (such as SUmb) or it may be + generated by IDWarp internally. + + Parameters + ---------- + pts : array, size (M, 3) + Nodes on this processor + conn : int array, size (sum(faceSizes)) + Connectivity of the nodes on this processor + faceSizes : int Array size (N) + Treat the conn array as a flat list with the faceSizes giving + connectivity offset for each element. + cgnsBlockID : dummy argument. + This argument is not used at all. It is here just to have the + same API as the IDWarpMulti class. + """ + conn = np.array(conn).flatten() + + # Get the restart file option + restartFile = "" if self.getOption("restartFile") is None else self.getOption("restartFile") + + self._setSymmetryConditionsAxisymm() + + allPts = np.vstack(self.comm.allgather(pts)) + allFaceSizes = np.hstack(self.comm.allgather(faceSizes)) + + # Be careful with conn, need to increment by the offset from the points + ptSizes = self.comm.allgather(len(pts)) + offsets = np.zeros(len(ptSizes), "intc") + offsets[1:] = np.cumsum(ptSizes)[:-1] + allConn = np.hstack(self.comm.allgather(conn + offsets[self.comm.rank])) + + # Next point reduce all nodes: + uniquePts, link, nUnique = self.warp.pointreduce(allPts.T, 1e-12) + uniquePts = uniquePts.T[0:nUnique] + + # Now sort them by z, y and x to ensure that the tree is + # always created the same way. It is currently sensitivte to + # the ordering of the original nodes. Ie, if you reorder the + # nodes you get a differnet tree which gives slightly + # differnet warp values/derivatives. + + ind = np.lexsort(uniquePts.T) + uniquePts = uniquePts[ind] + + # Update the link array such that it points back to the original list. + inv = np.zeros_like(ind, "intc") + rng = np.arange(len(ind)).astype("intc") + + inv[ind[rng]] = rng + link = inv[link - 1] + 1 + self.warp.initializewarping(np.ravel(pts), uniquePts.T, link, allFaceSizes, allConn, restartFile) + self.nSurf = len(pts) + self.warpInitialized = True + + # Clean-up patch data no longer necessary. + self.warp.deallocatepatchio() + + def setSurfaceDefinitionFromFile(self, surfFile): + """Set the surface definition for the warping from a + multiblock PLOT3D surface file + + Parameters + ---------- + surfFile : filename of multiblock PLOT3D surface file. + """ + + # Read the PLOT3D surface file + self.warp.readplot3dsurface(surfFile) + + if self.comm.rank == 0: + # Extract out the data we want from the module. Note that + # becuase we are in python, convert conn to 0 based. + pts = self.warp.plot3dsurface.pts.T.copy() + conn = self.warp.plot3dsurface.conn.T.copy() - 1 + faceSizes = 4 * np.ones(len(conn)) + else: + pts = np.zeros((0, 3)) + conn = np.zeros((0, 4)) + faceSizes = 4 * np.ones(0) + + # Run the regular setSurfaceDefinition routine. Note that only + # the root proc has non-zero length arrays. That's ok. + self.setSurfaceDefinition(pts, conn, faceSizes) + + # Cleanup the fortran memory we used to read the surface + self.warp.plot3dsurface.pts = None + self.warp.plot3dsurface.conn = None + + def setSurfaceFromFile(self, surfFile): + """Update the internal surface surface coordinates using an + external PLOT3D surface file. This can be used in an analogous + way to setSurfaceDefinitionFromFile. The 'sense' of the file + must be same as the file used with + setSurfaceDefinitionFromFile. That means, the same number of + blocks, with the same sizes, in the same order. + + Parameters + ---------- + surfFile: filename of multiblock PLOT3D surface file' + """ + + # Read the PLOT3D surface file + self.warp.readplot3dsurface(surfFile) + + if self.comm.rank == 0: + # Extract out the data we want from the module + pts = self.warp.plot3dsurface.pts.T.copy() + else: + pts = np.zeros((0, 3)) + + # Do the regular update + self.setSurfaceCoordinates(pts) + + # Cleanup the memory we used to read the surface + self.warp.plot3dsurface.pts = None + self.warp.plot3dsurface.conn = None + + def getSolverGrid(self): + """Return the current grid in the order specified by + setExternalMeshIndices(). This is the main routine for + returning the defomed mesh to the external CFD solver. + + Returns + ------- + solverGrid, numpy array, real: The resulting grid. + The output is returned in flatted 1D coordinate + format. The len of the array is 3*len(indices) as + set by setExternalMeshIndices() + """ + solverGrid = np.zeros(self.warp.griddata.solvermeshdof, self.dtype) + warpGrid = self.getWarpGrid() + self.warp.warp_to_solver_grid(warpGrid, solverGrid) + + return solverGrid + + def getWarpGrid(self): + """ + Return the current grid. This funtion is typically unused. See + getSolverGrid for the more useful interface functionality. + + Returns + ------- + warpGrid, numpy array, real: The resulting grid. + The output is returned in flatted 1D coordinate + format. + """ + return self.warp.getvolumecoordinates(self.warp.griddata.warpmeshdof) + + def getCommonGrid(self): + """Return the grid in the original ordering. This is required for the + OpenFOAM tecplot writer since the connectivity is only known + in this ordering. + """ + return self.warp.getcommonvolumecoordinates(self.warp.griddata.commonmeshdof) + + def getdXs(self): + """Return the current values in dXs. This is the result from a + mesh-warp derivative computation. + + Returns + ------- + dXs : numpy array + The specific components of dXs. size(N,3). This the same + size as the array obtained with getSurfaceCoordiantes(). N may + be zero if this processor does not have any surface coordinates. + """ + dXs = np.zeros((self.nSurf, 3), self.dtype) + self.warp.getdxs(np.ravel(dXs)) + + return dXs + + def warpMesh(self): + """ + Run the applicable mesh warping strategy. + + This will update the volume coordinates to match surface + coordinates set with setSurfaceCoordinates() + """ + # Initialize internal surface if one isn't already set. + self._setInternalSurface() + + # Set Options + self._setMeshOptions() + + # Now run mesh warp command + self.warp.warpmesh() + + def warpDeriv(self, dXv, solverVec=True): + """Compute the warping derivative (dXv/dXs^T)*Vec (where vec is the + dXv argument to this function. + + This is the main routine to compute the mesh warping + derivative. + + Parameters + ---------- + dXv : numpy array + Vector of size external solver_grid. This is typically + obtained from the external solver's dRdx^T * psi + calculation. + + solverVec : logical + Flag to indicate that the dXv vector is in the solver + ordering and must be converted to the warp ordering + first. This is the usual approach and thus defaults to + True. + + Returns + ------- + None. The resulting calculation is available from the getdXs() + function. + + """ + if solverVec: + dXvWarp = np.zeros(self.warp.griddata.warpmeshdof, self.dtype) + self.warp.solver_to_warp_grid(dXv, dXvWarp) + else: + dXvWarp = dXv + self.warp.warpderiv(dXvWarp) + + def warpDerivFwd(self, dXs, solverVec=True): + """Compute the forward mode warping derivative + + This routine is not used for "regular" optimization; it is + used for matrix-free type optimization. dXs is assumed to be + the the peturbation on all the surface nodes. + + Parameters + ---------- + dXs : array, size Nsx3 + This is the forward mode peturbation seed. Same size as the + surface mesh from getSurfaceCoordinates(). + solverVec : bool + Whether or not to convert to the solver ordering. + + Returns + ------- + dXv : array + The peturbation on the volume meshes. It may be + in warp ordering or solver ordering depending on + the solverVec flag. + """ + dXvWarp = np.zeros(self.warp.griddata.warpmeshdof, dtype=self.dtype) + self.warpMesh() + self.warp.warpderivfwd(np.ravel(dXs), dXvWarp) + + if solverVec: + dXv = np.zeros(self.warp.griddata.solvermeshdof, dtype=self.dtype) + self.warp.warp_to_solver_grid(dXvWarp, dXv) + return dXv + else: + return dXvWarp + + def verifyWarpDeriv(self, dXv=None, solverVec=True, dofStart=0, dofEnd=10, h=1e-6, randomSeed=314): + """Run an internal verification of the solid warping + derivatives""" + + if dXv is None: + np.random.seed(randomSeed) # 'Random' seed to ensure runs are same + dXvWarp = np.random.random(self.warp.griddata.warpmeshdof) # , dtype=self.dtype) + else: + if solverVec: + dXvWarp = np.zeros(self.warp.griddata.warpmeshdof, self.dtype) + self.warp.solver_to_warp_grid(dXv, dXvWarp) + else: + dXvWarp = dXv + + self.warp.verifywarpderiv(dXvWarp, dofStart, dofEnd, h) + + # ========================================================================== + # Output Functionality + # ========================================================================== + def writeGrid(self, fileName=None): + """ + Write the current grid to the correct format + + Parameters + ---------- + fileName : str or None + + Filename for grid. Should end in .cgns for CGNS files. For + PLOT3D whatever you want. It is not optional for + CGNS/PLOT3D. It is not required for OpenFOAM meshes. This + call will update the 'points' file. + """ + if self.fileType == "CGNS": + if fileName is None: + raise Error("fileName is not optional for writeGrid with gridType of CGNS") + + # Copy the default and then write + if self.comm.rank == 0: + shutil.copy(self.getOption("gridFile"), fileName) + self.warp.writecgns(fileName) + + # ========================================================================= + # Internal Private Functions + # ========================================================================= + def _setInternalSurface(self): + """This function is used by default if setSurfaceDefinition() is not + set BEFORE an operation is requested that requires this + information. + """ + + if self.warpInitialized: + return + + if self.comm.rank == 0: + warnings.warn( + "Using internally generated IDWarp surfaces. If " + "this mesh object is to be used with an " + "external solver, ensure the mesh object is " + "passed to the solver immediatedly after it is created. " + "The external solver must then call " + "'setExternalMeshIndices()' and 'setSurfaceDefinition()' " + "routines.", + stacklevel=2, + ) + + conn = [] + pts = [] + faceSizes = [] + + if self.comm.rank == 0: + # Do the necessary fortran preprocessing + self.warp.processstructuredpatches() + + fullConn = self.warp.cgnsgrid.surfaceconn - 1 + fullPts = self.warp.cgnsgrid.surfacepoints + + nPatch = self.warp.cgnsgrid.getnpatch() + fullPatchNames = [] + for i in range(nPatch): + fullPatchNames.append(self.warp.cgnsgrid.getsurf(i + 1).strip().lower()) + # We now have all surfaces belonging to + # boundary conditions. We need to decide which + # ones to use depending on what the user has + # told us. + surfaceFamilies = set() + if self.getOption("specifiedSurfaces") is None: + # Use all wall surfaces: + for i in range(len(self.warp.cgnsgrid.surfaceiswall)): + if self.warp.cgnsgrid.surfaceiswall[i]: + surfaceFamilies.add(fullPatchNames[i].lower()) + else: + # The user has supplied a list of surface families + for name in self.getOption("specifiedSurfaces"): + surfaceFamilies.add(name.lower()) + + usedFams = set() + if self.warp.cgnsgrid.cgnsstructured: + if self.comm.rank == 0: + # Pull out data and convert to 0-based ordering + fullPatchSizes = self.warp.cgnsgrid.surfacesizes.T + + # Now we loop back over the "full" versions of + # thing and just take the ones that correspond + # to the families we are using. + + curNodeIndex = 0 + curCellIndex = 0 + curOffset = 0 + for i in range(len(fullPatchNames)): + curNodeSize = fullPatchSizes[i][0] * fullPatchSizes[i][1] + curCellSize = (fullPatchSizes[i][0] - 1) * (fullPatchSizes[i][1] - 1) + + if fullPatchNames[i] in surfaceFamilies: + # Keep track of the families we've actually used + usedFams.add(fullPatchNames[i]) + pts.extend(fullPts[curNodeIndex : curNodeIndex + curNodeSize * 3]) + conn.extend(fullConn[curCellIndex : curCellIndex + curCellSize * 4] - curOffset) + else: + # If we skipped, we increment the offset + curOffset += curNodeSize + + curNodeIndex += curNodeSize * 3 + curCellIndex += curCellSize * 4 + + # end for (root proc) + + # Run the common surface definition routine + pts = np.array(pts).reshape((len(pts) // 3, 3)) + faceSizes = 4 * np.ones(len(conn) // 4, "intc") + self.setSurfaceDefinition(pts=pts, conn=np.array(conn, "intc"), faceSizes=faceSizes) + + # Check if all supplied family names were actually + # used. The user probably wants to know if a family + # name was specified incorrectly. + if self.comm.rank == 0: + if usedFams < surfaceFamilies: + missing = list(surfaceFamilies - usedFams) + warnings.warn( + "Not all specified surface families that " + "were given were found the CGNS file. " + "The families not found are %s." % (repr(missing)), + stacklevel=2, + ) + if len(usedFams) == 0: + raise Error( + "No surfaces were selected. Check the names " + "given in the 'specifiedSurface' option. The " + "list of families is %s." % (repr(list(fullPatchNames))) + ) + + def _getFullPatchData(self): + self.warp.processstructuredpatches() + + fullConn = self.warp.cgnsgrid.surfaceconn - 1 + fullPts = self.warp.cgnsgrid.surfacepoints + + nPatch = self.warp.cgnsgrid.getnpatch() + fullPatchNames = [self.warp.cgnsgrid.getsurf(i + 1).strip().lower() for i in range(nPatch)] + fullPatchNames = [val.decode("utf-8") if isinstance(val, bytes) else val for val in fullPatchNames] + + return fullPatchNames, fullConn, fullPts + + def _getSymmetryFamilies(self, fullPatchNames): + symmFamilies = set() + symmFamilies.update([fullPatchNames[i] for i, isSymm in enumerate(self.warp.cgnsgrid.surfaceissymm) if isSymm]) + + # Decode family names that are bytes to convert them to strings + symmFamilies = {val.decode("utf-8") if isinstance(val, bytes) else val for val in symmFamilies} + + return symmFamilies + + def _getSymmetryPlanesFromGeometry(self): + fullPatchNames, fullConn, fullPts = self._getFullPatchData() + symmFamilies = self._getSymmetryFamilies(fullPatchNames) + symmPatches = [] + usedFams = set() + + fullPatchSizes = self.warp.cgnsgrid.surfacesizes.T + + curNodeIndex = 0 + curCellIndex = 0 + curOffset = 0 + + for i, patchName in enumerate(fullPatchNames): + # Get the node and cell sizes + curNodeSize = fullPatchSizes[i][0] * fullPatchSizes[i][1] + curCellSize = (fullPatchSizes[i][0] - 1) * (fullPatchSizes[i][1] - 1) + + if patchName in symmFamilies: + # Keep track of the families we've actually used + usedFams.add(patchName) + + # Determine the average normal for this patch + conn1 = fullConn[curCellIndex : curCellIndex + curCellSize * 4] + conn2 = fullConn[curCellIndex] + conn = conn1 - conn2 + 1 + + idxStart = curNodeIndex + idxEnd = idxStart + curNodeSize * 3 + pts = fullPts[curNodeIndex : curNodeIndex + curNodeSize * 3] + avgNorm = self.warp.averagenormal(pts, conn, 4 * np.ones(curCellSize, "intc")) + + # Decode patchName if it's a byte to convert it to a string + patchName = patchName.decode("utf-8") if isinstance(patchName, bytes) else patchName + patch = MeshPatch(idxStart, idxEnd, pts[:3], avgNorm, curNodeSize, curCellSize, patchName) + symmPatches.append(patch) + else: + curOffset += curNodeSize + + curNodeIndex += curNodeSize * 3 + curCellIndex += curCellSize * 4 + + self._checkUsedFamilies(usedFams, symmFamilies) + + return symmPatches, usedFams + + def _checkUsedFamilies(self, usedFams, symmFamilies): + if usedFams < symmFamilies: + missing = list(symmFamilies - usedFams) + warnings.warn( + "Not all specified symm families that " + "were given were found the CGNS file. " + "The families not found are %s." % (repr(missing)), + stacklevel=2, + ) + + def _printSymmetryPlanes(self, pts, normals): + print("+-------------------- Symmetry Planes -------------------+") + print("| Point Normal |") + for i in range(len(pts)): + print( + "| (%7.3f %7.3f %7.3f) (%7.3f %7.3f %7.3f) |" + % ( + np.real(pts[i][0]), + np.real(pts[i][1]), + np.real(pts[i][2]), + np.real(normals[i][0]), + np.real(normals[i][1]), + np.real(normals[i][2]), + ) + ) + print("+--------------------------------------------------------+") + + def _setSymmetryConditionsAxisymm(self): + # Initialize these variables on all procs + pts = [] + normals = [] + mirrorFamName = None + if self.comm.rank == 0: + mirrorFamName = self.getOption("axiMirrorFamily") + + if mirrorFamName is None: + raise Error("Must specify an 'axiMirrorFamily' for axisymmetric mesh warping.") + + # We need the information for all symm planes for axisymmetric warping + symmPatches, usedFams = self._getSymmetryPlanesFromGeometry() + + # Should only have two symmetry planes for an axisymm mesh + if len(usedFams) != 2: + raise Error("Could not detect 2 unique symmetry planes for an axisymmetric mesh.") + + mirrorPatches = [patch for patch in symmPatches if patch.patchName == mirrorFamName] + rotPatches = [patch for patch in symmPatches if patch.patchName != mirrorFamName] + + # Get the first two patches to represent the symmetry planes + plane1 = mirrorPatches[0] + plane2 = rotPatches[0] + + pts.append(plane1.point) + normals.append(plane1.normal) + + self._printSymmetryPlanes([plane1.point, plane2.point], [plane1.normal, plane2.normal]) + + # Broadcast the mirror family name to all procs + mirrorFamName = self.comm.bcast(mirrorFamName) + self.warp.griddata.setmirrorfamily(mirrorFamName) + + # We only want to set the mirror plane as the symmetry plane in the + # fortran layer. This "tricks" the mesh warping algorithm so that + # the mirroring will get equal weights to preserve the symmetry plane + # warping. + pts = self.comm.bcast(np.array(pts)) + normals = self.comm.bcast(np.array(normals)) + axiPlane = self.getOption("axiPlane") + if axiPlane is None: + raise Error( + "Must supply a point and normal vector using the " + "'axiPlane' option to perform axisymmetric mesh warping." + ) + + # Set a "phantom" plane that preserves zero-warping along the degenerate line BC + pts = np.vstack((pts, axiPlane[0])) + normals = np.vstack((normals, axiPlane[1])) + + self.warp.setsymmetryplanes(pts.T, normals.T) + + # ========================================================================= + # Utility functions + # ========================================================================= + def _setMeshOptions(self): + """Private function to set the options currently in + self.options to the corresponding place in Fortran""" + + self.warp.gridinput.ldeffact = self.getOption("LdefFact") + self.warp.gridinput.alpha = self.getOption("alpha") + self.warp.gridinput.aexp = self.getOption("aExp") + self.warp.gridinput.bexp = self.getOption("bExp") + self.warp.gridinput.symmtol = self.getOption("symmTol") + self.warp.gridinput.userotations = self.getOption("useRotations") + self.warp.gridinput.zerocornerrotations = self.getOption("zeroCornerRotations") + self.warp.gridinput.cornerangle = self.getOption("cornerAngle") + self.warp.gridinput.errtol = self.getOption("errTol") + self.warp.kd_tree.bucket_size = self.getOption("bucketSize") + self.warp.griddata.axisymm = self.getOption("axiSymmetric") + if self.getOption("evalMode") == "fast": + self.warp.gridinput.evalmode = self.warp.gridinput.eval_fast + elif self.getOption("evalMode") == "exact": + self.warp.gridinput.evalmode = self.warp.gridinput.eval_exact + + if self.getOption("axiSymmetric"): + if self.getOption("axiAngle") != 0.0: + angle = self.getOption("axiAngle") + self.warp.griddata.axisymmangle = np.deg2rad(angle) # Convert to radians + else: + raise Error("Must supply an 'axiSymmAngle' for axisymmetric mesh warping.") + + if any([i != 0.0 for i in self.getOption("axiAxis")]): + self.warp.griddata.axisymmaxis = np.array(self.getOption("axiAxis"), self.dtype) + else: + raise Error("Must supply an 'axiAxis' for axisymmetric mesh warping.") + + def _processFortranStringArray(self, strArray): + """Getting arrays of strings out of Fortran can be kinda nasty. This + takes the array and returns a nice python list of strings""" + shp = strArray.shape + arr = strArray.reshape((shp[1], shp[0]), order="F") + tmp = [] + + for i in range(arr.shape[1]): + tmp.append("") + for j in range(arr.shape[0]): + tmp[-1] += arr[j, i].decode() + tmp[-1] = tmp[-1].strip().lower() + + return tmp + + def __del__(self): + """Release all the mesh warping memory. This should be called + automatically when the object is garbage collected.""" + self.warp.releasememory() diff --git a/idwarp/AxisymmetricMesh_C.py b/idwarp/AxisymmetricMesh_C.py new file mode 100644 index 0000000..bcb7a30 --- /dev/null +++ b/idwarp/AxisymmetricMesh_C.py @@ -0,0 +1,22 @@ +import os + +from . import MExt +from .AxisymmetricMesh import AxisymmetricMesh + + +class AxisymmetricMesh_C(AxisymmetricMesh): + """ + Represents a (Complex) Axisymmetric mesh + """ + + def __init__(self, *args, **kwargs): + """Initialize the object.""" + + debug = False + if "debug" in kwargs: + debug = kwargs["debug"] + + curDir = os.path.basename(os.path.dirname(os.path.realpath(__file__))) + self.warp = MExt.MExt("libidwarp_cs", curDir, debug=debug)._module + AxisymmetricMesh.__init__(self, *args, **kwargs) + self.dtype = "D" diff --git a/idwarp/__init__.py b/idwarp/__init__.py index 67fd5a4..d12f92c 100644 --- a/idwarp/__init__.py +++ b/idwarp/__init__.py @@ -1,8 +1,10 @@ __version__ = "2.6.2" +from .AxisymmetricMesh import AxisymmetricMesh +from .AxisymmetricMesh_C import AxisymmetricMesh_C from .UnstructuredMesh import USMesh from .UnstructuredMesh_C import USMesh_C from .MultiUnstructuredMesh import MultiUSMesh from .MultiUnstructuredMesh_C import MultiUSMesh_C -__all__ = ["USMesh", "MultiUSMesh", "USMesh_C", "MultiUSMesh_C"] +__all__ = ["AxisymmetricMesh", "USMesh", "MultiUSMesh", "AxisymmetricMesh_C", "USMesh_C", "MultiUSMesh_C"] diff --git a/src/IO/readStructuredCGNS.F90 b/src/IO/readStructuredCGNS.F90 index 75b650b..ab5402e 100644 --- a/src/IO/readStructuredCGNS.F90 +++ b/src/IO/readStructuredCGNS.F90 @@ -28,8 +28,9 @@ subroutine readStructuredCGNS(cg) #endif integer(kind=intType), dimension(:), allocatable :: surfaceNodes, localSurfaceNodes integer(kind=intType), dimension(:, :), allocatable :: sizes - integer(kind=intType) :: nSurf, nodeCount, nConn, ni, nj, nx, ny, nz + integer(kind=intType) :: nSurf, nodeCount, nConn, ni, nj, nx, ny, nz, bocoIdx integer(kind=intType) :: status(MPI_STATUS_SIZE) + integer(kind=intType) :: nTotalBocos logical :: lowFace ! Only do reading on root proc: @@ -66,6 +67,8 @@ subroutine readStructuredCGNS(cg) nSurf = 0 nConn = 0 + nTotalBocos = 0 + bocoIdx = 0 zoneLoop1: do iZone = 1, nZones call cg_zone_read_f(cg, base, iZone, zoneName, dims, ierr) if (ierr .eq. CG_ERROR) call cg_error_exit_f @@ -77,6 +80,7 @@ subroutine readStructuredCGNS(cg) call cg_nbocos_f(cg, base, iZone, nbocos, ierr) if (ierr .eq. ERROR) call cg_error_exit_f allocate (zones(iZone)%bocos(nbocos)) + nTotalBocos = nTotalBocos + nbocos bocoLoop1: do boco = 1, nbocos @@ -137,6 +141,9 @@ subroutine readStructuredCGNS(cg) allocate (surfacePoints(3 * nSurf)) allocate (surfaceConn(4 * nConn)) + ! Allocate the flattened BC types storage + allocate (bocoTypes(nTotalBocos)) + ! Reset the counters here nSurf = 0 offSet = 0 @@ -191,7 +198,9 @@ subroutine readStructuredCGNS(cg) end do bocoLoop2: do boco = 1, size(zones(iZone)%bocos) + bocoIdx = bocoIdx + 1 bocoType = zones(iZone)%bocos(boco)%type + bocoTypes(bocoIdx) = bocoType pts = zones(iZone)%bocos(boco)%ptRange ! Flag the surface nodes @@ -284,7 +293,7 @@ subroutine readStructuredCGNS(cg) if (ierr .eq. CG_ERROR) call cg_error_exit_f deallocate (sizes) else - ! Allocate these to zero so we can just blindly dealloc later + !Allocate these to zero so we can just blindly dealloc later allocate (surfacePoints(0), surfaceConn(0)) end if rootProc @@ -292,6 +301,73 @@ subroutine readStructuredCGNS(cg) call MPI_bcast(nNodes, 1, MPI_INTEGER, 0, warp_comm_world, ierr) call EChk(ierr, __FILE__, __LINE__) + ! For axisymm cases we need to communicate some of the zone data to every proc + axisymmetric: if (axisymm) then + ! Total number of zones + call MPI_bcast(nZones, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Allocate the zones + if (myid /= 0) then + allocate (zones(nZones)) + end if + + ! Communicate the zone data + do i = 1, nZones + ! Communicate the zone data + call MPI_bcast(zones(i)%il, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%jl, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%kl, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%name, maxCGNSNameLen, MPI_CHARACTER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%nVertices, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%nElements, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Get the number of bocos for this zone from the root proc + if (myID == 0) nbocos = size(zones(i)%bocos) + + ! Communicate the number of bocos + call MPI_bcast(nbocos, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + if (myid /= 0) then + allocate (zones(i)%bocos(nbocos)) + + ! Nullify section pointers since we won't have any for structured mesh + nullify (zones(i)%sections) + end if + + do j = 1, nbocos + ! Communicate all boco data + call MPI_bcast(zones(i)%bocos(j)%name, maxCGNSNameLen, MPI_CHARACTER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%bocos(j)%type, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%bocos(j)%ptRange, 6, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%bocos(j)%family, maxCGNSNameLen, MPI_CHARACTER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%bocos(j)%nBCElem, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + call MPI_bcast(zones(i)%bocos(j)%nBCNodes, 1, MPI_INTEGER, 0, warp_comm_world, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Allocate empty data so we don't get errors trying to deallocate later + if (myid /= 0) then + allocate (zones(i)%bocos(j)%BCElements(0)) + nullify (zones(i)%bocos(j)%elemPtr) + nullify (zones(i)%bocos(j)%elemConn) + allocate (zones(i)%bocos(j)%elemNodes(0, 0)) + end if + end do + end do + end if axisymmetric + rootProc2: if (myid == 0) then istart = 0 + 1 diff --git a/src/f2py/libidwarp.pyf b/src/f2py/libidwarp.pyf index c8bb084..5a477da 100644 --- a/src/f2py/libidwarp.pyf +++ b/src/f2py/libidwarp.pyf @@ -165,12 +165,20 @@ python module libidwarp ! in use constants integer(kind=inttype) :: warpmeshdof integer(kind=inttype) :: solvermeshdof - integer(Kind=inttype) :: commonmeshdof + integer(kind=inttype) :: commonmeshdof + logical, optional :: axisymm=.false. + real(kind=realtype) :: axisymmangle + real(kind=realtype) dimension(3) :: axisymmaxis + + subroutine setmirrorfamily(famname) ! in :idwarp:../modules/gridData.f90 + character*(*) intent(in) :: famname + end subroutine setmirrorfamily end module griddata module cgnsgrid real(kind=realtype), dimension(:), allocatable :: surfacepoints integer(kind=inttype), dimension(:), allocatable :: surfaceconn + integer(kind=inttype), dimension(:), allocatable :: bocoTypes logical, dimension(:), allocatable :: surfaceiswall logical, dimension(:), allocatable :: surfaceissymm integer(kind=inttype), dimension(:,:), allocatable :: surfacesizes @@ -204,7 +212,6 @@ python module libidwarp ! in module kd_tree integer(kind=inttype) :: bucket_size end module kd_tree - end interface #ifdef USE_COMPLEX end python module libidwarp_cs diff --git a/src/modules/cgnsGrid.F90 b/src/modules/cgnsGrid.F90 index 085e9e2..3409690 100644 --- a/src/modules/cgnsGrid.F90 +++ b/src/modules/cgnsGrid.F90 @@ -45,6 +45,9 @@ module cgnsGrid ! List of the zones type(zoneDataType), dimension(:), allocatable :: zones + ! All BC types (Need this for axisymmetric meshes) + integer(kind=intType), dimension(:), allocatable :: bocoTypes + ! Deduced information of the surfaces real(kind=realType), dimension(:), allocatable :: surfacePoints logical, dimension(:), allocatable :: surfaceIsWall diff --git a/src/modules/gridData.F90 b/src/modules/gridData.F90 index a4299e0..025d98e 100644 --- a/src/modules/gridData.F90 +++ b/src/modules/gridData.F90 @@ -39,7 +39,7 @@ module gridData #ifndef USE_TAPENADE real(kind=realType), pointer, dimension(:) :: XsPtrb, XsPtrd - real(kind=realType), pointer, dimension(:) :: XvPtrb, XVPtrd + real(kind=realType), pointer, dimension(:) :: XvPtrb, XvPtrd #endif ! Sizes of the three different mesh sizes: @@ -60,6 +60,18 @@ module gridData integer(kind=intType) :: nLoop real(kind=realType), dimension(:, :), allocatable :: symmPts, symmNormals + ! Axisymmetric Information + logical :: axiSymm = .false. + real(kind=realType) :: axiSymmAngle + real(kind=realType) :: axiSymmAxis(3) + character(len=maxStringLen) :: mirrorFamily + +contains + + subroutine setMirrorFamily(famName) + character(len=*), intent(in) :: famName + mirrorFamily = famName + end subroutine setMirrorFamily end module gridData module plot3dSurface diff --git a/src/modules/kd_tree.F90 b/src/modules/kd_tree.F90 index f41da44..a531385 100644 --- a/src/modules/kd_tree.F90 +++ b/src/modules/kd_tree.F90 @@ -1637,7 +1637,7 @@ subroutine computeNodalProperties(tp, initialPoint) else ! Now get the rotation Matrix if (useRotations .and. .not. tp%isCorner(i)) then - call getRotationMatrix3d(tp%normals0(:, i), tp%normals(:, i), tp%Mi(:, :, i)) + call getRotationMatrixVecVec(tp%normals0(:, i), tp%normals(:, i), tp%Mi(:, :, i)) else ! Set identity tp%Mi(:, :, i) = zero diff --git a/src/utils/releaseMemory.F90 b/src/utils/releaseMemory.F90 index 8a9f174..4bcc414 100644 --- a/src/utils/releaseMemory.F90 +++ b/src/utils/releaseMemory.F90 @@ -106,6 +106,11 @@ subroutine releaseMemory ! Delete zone array itself deallocate (zones) + ! Delete the boco types array + if (allocated(bocoTypes)) then + deallocate (bocoTypes) + end if + ! Delete the surface-stuff in the cgns grid if necessary if (allocated(surfacePoints)) then deallocate (surfacePoints) diff --git a/src/utils/vectorUtils.f90 b/src/utils/vectorUtils.f90 index 9c1abfa..6444fe3 100644 --- a/src/utils/vectorUtils.f90 +++ b/src/utils/vectorUtils.f90 @@ -28,7 +28,7 @@ subroutine getMag(V, mag) end subroutine getMag -subroutine getRotationMatrix3d(v1, v2, Mi) +subroutine getRotationMatrixVecVec(v1, v2, Mi) use precision use constants implicit none @@ -100,4 +100,54 @@ subroutine getRotationMatrix3d(v1, v2, Mi) Mi = Mi + sin(angle) * A + (one - cos(angle)) * C -end subroutine getRotationMatrix3d +end subroutine getRotationMatrixVecVec + +subroutine getRotationMatrixAngleAxis(angle, axis, Mi) + use precision + use constants + implicit none + + ! Input variables + real(kind=realType), intent(in) :: axis(3), angle + + ! Output variables + real(kind=realType), intent(out) :: Mi(3, 3) + + ! Working variables + real(kind=realType), dimension(3, 3) :: A, C + real(kind=realType) :: overNorm, axisNorm(3), axisMag + + call getMag(axis, axisMag) + overNorm = 1 / axisMag + axisNorm = axis * overNorm + + ! Compute the rotation matrix + A(1, 1) = zero + A(1, 2) = -axisNorm(3) + A(1, 3) = axisNorm(2) + A(2, 1) = axisNorm(3) + A(2, 2) = zero + A(2, 3) = -axisNorm(1) + A(3, 1) = -axisNorm(2) + A(3, 2) = axisNorm(1) + A(3, 3) = zero + + !C = A*A + C(1, 1) = A(1, 1) * A(1, 1) + A(1, 2) * A(2, 1) + A(1, 3) * A(3, 1) + C(1, 2) = A(1, 1) * A(1, 2) + A(1, 2) * A(2, 2) + A(1, 3) * A(3, 2) + C(1, 3) = A(1, 1) * A(1, 3) + A(1, 2) * A(2, 3) + A(1, 3) * A(3, 3) + C(2, 1) = A(2, 1) * A(1, 1) + A(2, 2) * A(2, 1) + A(2, 3) * A(3, 1) + C(2, 2) = A(2, 1) * A(1, 2) + A(2, 2) * A(2, 2) + A(2, 3) * A(3, 2) + C(2, 3) = A(2, 1) * A(1, 3) + A(2, 2) * A(2, 3) + A(2, 3) * A(3, 3) + C(3, 1) = A(3, 1) * A(1, 1) + A(3, 2) * A(2, 1) + A(3, 3) * A(3, 1) + C(3, 2) = A(3, 1) * A(1, 2) + A(3, 2) * A(2, 2) + A(3, 3) * A(3, 2) + C(3, 3) = A(3, 1) * A(1, 3) + A(3, 2) * A(2, 3) + A(3, 3) * A(3, 3) + + ! Rodrigues formula for the rotation matrix + Mi = zero + Mi(1, 1) = one + Mi(2, 2) = one + Mi(3, 3) = one + + Mi = Mi + sin(angle) * A + (one - cos(angle)) * C +end subroutine getRotationMatrixAngleAxis diff --git a/src/warp/Makefile b/src/warp/Makefile index 15fab6b..ef48152 100644 --- a/src/warp/Makefile +++ b/src/warp/Makefile @@ -11,6 +11,7 @@ include ${RULES_FILE} vpath %.o $(OBJDIR) FILES_TO_COMPLEXIFY = \ + axisymm.F90\ getdXs.F90\ getElementProps.F90\ getSurfaceCoordinates.F90\ diff --git a/src/warp/axisymm.F90 b/src/warp/axisymm.F90 new file mode 100644 index 0000000..a890b7e --- /dev/null +++ b/src/warp/axisymm.F90 @@ -0,0 +1,476 @@ +subroutine copyRotateVolumeCoordinates() + use precision + use gridData + use cgnsGrid + implicit none + + ! Working variables + integer(kind=intType) :: i, j, k, ii, jj, istart, iend, ierr, iiLocal + integer(kind=intType) :: iZone, nZones, iBoco, ptRange(3, 2), rotDir, iPlane + integer(kind=intType) :: iNode, dir1, dir2, rotOffset, nNodes, mirrorIdx + integer(kind=intType), dimension(:), allocatable :: insertIdxs + real(kind=realType) :: Mi(3, 3), coords(3), rotCoords(3) + real(kind=realType), dimension(:), allocatable :: insertCoords + real(kind=realType), dimension(:), pointer :: coordsLocal + + nZones = size(zones) + + ! Get the rotation matrix + call getRotationMatrixAngleAxis(axiSymmAngle, axiSymmAxis, Mi) + + ! Get the range of the volume nodes on this proc + call VecGetOwnershipRange(Xv, istart, iend, ierr) + istart = istart + 1 ! Convert to 1-based indexing + call EChk(ierr, __FILE__, __LINE__) + + ! Get the local array of coordinates + call VecGetArrayF90(Xv, coordsLocal, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ii = 0 + jj = 0 + ! Loop over the zones in the grid + zoneLoop: do iZone = 1, nZones + rotDir = 0 ! Reset the rotation direction + + ! Loop over the BCs to see if there is a mirror plane BC + bocoLoop: do iBoco = 1, size(zones(iZone)%bocos) + ! Check if the BCs for this zone contain the mirror family + if (zones(iZone)%bocos(iBoco)%family .eq. mirrorFamily) then + ! Determine the rotation direction + ! This is the direction in which the zone is one cell wide + ! The point range for the BC will be equal in this direction + ! signifying that the zone is one cell wide in this direction + ptRange = zones(iZone)%bocos(iBoco)%ptRange + if (ptRange(1, 1) == ptRange(1, 2)) then + mirrorIdx = ptRange(1, 1) + rotDir = 1 ! i-direction + else if (ptRange(2, 1) == ptRange(2, 2)) then + mirrorIdx = ptRange(2, 1) + rotDir = 2 ! j-direction + else if (ptRange(3, 1) == ptRange(3, 2)) then + mirrorIdx = ptRange(3, 1) + rotDir = 3 ! k-direction + end if + end if + end do bocoLoop + + if (rotDir == 0) then + ! No mirror plane BC found, skip this zone + cycle zoneLoop + end if + + ! Determine the plane directions to loop over and set the rotation + ! plane offset (These come from the loop structure in readStructuredCGNS.F90) + if (rotDir == 1) then + ! The zone is one cell wide in the i-direction (final loop in readStructuredCGNS.F90) + dir1 = zones(izone)%jl + dir2 = zones(izone)%kl + rotOffset = 1 + else if (rotDir == 2) then + ! The zone is one cell wide in the j-direction (second loop in readStructuredCGNS.F90) + dir1 = zones(iZone)%il + dir2 = zones(iZone)%kl + rotOffset = zones(iZone)%il + else if (rotDir == 3) then + ! The zone is one cell wide in the k-direction (first loop in readStructuredCGNS.F90) + dir1 = zones(izone)%il + dir2 = zones(izone)%jl + rotOffset = zones(iZone)%il * zones(iZone)%jl + end if + + ! Number of nodes on a single plane for this zone + nNodes = dir1 * dir2 + + ! Loop over the nodes on both planes in this zone (mirror and rotated) + nodeLoop: do iNode = 1, 2 * nNodes + ii = ii + 1 ! Global index + + if (3 * ii - 2 .lt. istart .or. 3 * ii .gt. iend) then + ! ii is not in the local ownership range so skip + cycle nodeLoop + end if + + ! Determine if the node is on the first or second plane (1 or 0) + iPlane = mod((iNode - 1) / rotOffset + 1, 2) + if (mirrorIdx .eq. 1 .and. iPlane .eq. 1) then + ! Mirror plane is the first plane and we are on the mirror plane + jj = ii + rotOffset + else if (mirrorIdx .eq. 2 .and. iPlane .eq. 0) then + ! Mirror plane is the second plane and we are on the mirror plane + jj = ii - rotOffset + else + ! This node is on the rot plane so skip + cycle nodeLoop + end if + + ! Adjust ii to be in the local indexing + iiLocal = 3 * ii - istart + 1 + + ! Copy the x, y, z coordinates + coords(1) = coordsLocal(iiLocal - 2) ! x + coords(2) = coordsLocal(iiLocal - 1) ! y + coords(3) = coordsLocal(iiLocal) ! z + + ! Apply the rotation + rotCoords = matmul(Mi, coords) + + ! Add the rotated coords to the array of coords to insert + ! into the global vector + if (.not. allocated(insertCoords)) then + allocate (insertCoords(3)) + insertCoords = rotCoords + else + insertCoords = [insertCoords, rotCoords] + end if + + ! Add the index to the array of indices to insert + ! into the global vector + if (.not. allocated(insertIdxs)) then + allocate (insertIdxs(3)) + insertIdxs = (/3 * jj - 2, 3 * jj - 1, 3 * jj/) + else + insertIdxs = [insertIdxs, (/3 * jj - 2, 3 * jj - 1, 3 * jj/)] + end if + + end do nodeLoop + + end do zoneLoop + + ! Insert the rotated coordinates into the global vector + ! We only do this collective communication once at the end if insertCoords is allocated + if (allocated(insertCoords)) then + call VecSetValues(Xv, size(insertCoords), insertIdxs - 1, insertCoords, INSERT_VALUES, ierr) + call EChk(ierr, __FILE__, __LINE__) + end if + + ! Need to call VecAssemblyBegin/End + call VecAssemblyBegin(Xv, ierr) + call EChk(ierr, __FILE__, __LINE__) + call VecAssemblyEnd(Xv, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Restore the array to deallocate the coordsLocal pointer + call VecRestoreArrayF90(Xv, coordsLocal, ierr) + call EChk(ierr, __FILE__, __LINE__) + +end subroutine copyRotateVolumeCoordinates + +subroutine copyRotateVolumeCoordinates_d() + use precision + use gridData + use cgnsGrid + implicit none + + ! Working variables + integer(kind=intType) :: i, j, k, ii, jj, istart, iend, ierr, iiLocal + integer(kind=intType) :: iZone, nZones, iBoco, ptRange(3, 2), rotDir, iPlane + integer(kind=intType) :: iNode, dir1, dir2, rotOffset, nNodes, mirrorIdx + integer(kind=intType), dimension(:), allocatable :: insertIdxs + real(kind=realType) :: Mi(3, 3), coords_d(3), rotCoords_d(3) + real(kind=realType), dimension(:), allocatable :: insertCoords_d + real(kind=realType), dimension(:), pointer :: coordsLocal_d + + nZones = size(zones) + + ! Get the rotation matrix + call getRotationMatrixAngleAxis(axiSymmAngle, axiSymmAxis, Mi) + + ! Get the range of the volume nodes on this proc + call VecGetOwnershipRange(dXv, istart, iend, ierr) + istart = istart + 1 ! Convert to 1-based indexing + call EChk(ierr, __FILE__, __LINE__) + + ! Get the local array of coordinates + call VecGetArrayF90(dXv, coordsLocal_d, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ii = 0 + jj = 0 + ! Loop over the zones in the grid + zoneLoop: do iZone = 1, nZones + rotDir = 0 ! Reset the rotation direction + + ! Loop over the BCs to see if there is a mirror plane BC + bocoLoop: do iBoco = 1, size(zones(iZone)%bocos) + ! Check if the BCs for this zone contain the mirror family + if (zones(iZone)%bocos(iBoco)%family .eq. mirrorFamily) then + ! Determine the rotation direction + ! This is the direction in which the zone is one cell wide + ! The point range for the BC will be equal in this direction + ! signifying that the zone is one cell wide in this direction + ptRange = zones(iZone)%bocos(iBoco)%ptRange + if (ptRange(1, 1) == ptRange(1, 2)) then + mirrorIdx = ptRange(1, 1) + rotDir = 1 ! i-direction + else if (ptRange(2, 1) == ptRange(2, 2)) then + mirrorIdx = ptRange(2, 1) + rotDir = 2 ! j-direction + else if (ptRange(3, 1) == ptRange(3, 2)) then + mirrorIdx = ptRange(3, 1) + rotDir = 3 ! k-direction + end if + end if + end do bocoLoop + + if (rotDir == 0) then + ! No mirror plane BC found, skip this zone + cycle zoneLoop + end if + + ! Determine the plane directions to loop over and set the rotation + ! plane offset (These come from the loop structure in readStructuredCGNS.F90) + if (rotDir == 1) then + ! The zone is one cell wide in the i-direction (final loop in readStructuredCGNS.F90) + dir1 = zones(izone)%jl + dir2 = zones(izone)%kl + rotOffset = 1 + else if (rotDir == 2) then + ! The zone is one cell wide in the j-direction (second loop in readStructuredCGNS.F90) + dir1 = zones(iZone)%il + dir2 = zones(iZone)%kl + rotOffset = zones(iZone)%il + else if (rotDir == 3) then + ! The zone is one cell wide in the k-direction (first loop in readStructuredCGNS.F90) + dir1 = zones(izone)%il + dir2 = zones(izone)%jl + rotOffset = zones(iZone)%il * zones(iZone)%jl + end if + + ! Number of nodes on a single plane for this zone + nNodes = dir1 * dir2 + + ! Loop over the nodes on both planes in this zone (mirror and rotated) + nodeLoop: do iNode = 1, 2 * nNodes + ii = ii + 1 ! Global index + + if (3 * ii - 2 .lt. istart .or. 3 * ii .gt. iend) then + ! ii is not in the local ownership range so skip + cycle nodeLoop + end if + + ! Determine if the node is on the mirror plane or the rotated plane (1 or 0) + iPlane = mod((iNode - 1) / rotOffset + 1, 2) + if (mirrorIdx .eq. 1 .and. iPlane .eq. 1) then + ! Mirror plane is the first plane and we are on the mirror plane + jj = ii + rotOffset + else if (mirrorIdx .eq. 2 .and. iPlane .eq. 0) then + ! Mirror plane is the second plane and we are on the mirror plane + jj = ii - rotOffset + else + ! This node is on the rot plane so skip + cycle nodeLoop + end if + + ! Adjust ii to be in the local indexing + iiLocal = 3 * ii - istart + 1 + + ! Copy the x, y, z coordinates + coords_d(1) = coordsLocal_d(iiLocal - 2) ! x + coords_d(2) = coordsLocal_d(iiLocal - 1) ! y + coords_d(3) = coordsLocal_d(iiLocal) ! z + + ! Apply the rotation + rotCoords_d = matmul(Mi, coords_d) + + ! Add the rotated coords to the array of coords to insert + ! into the global vector + if (.not. allocated(insertCoords_d)) then + allocate (insertCoords_d(3)) + insertCoords_d = rotCoords_d + else + insertCoords_d = [insertCoords_d, rotCoords_d] + end if + + ! Add the index to the array of indices to insert + ! into the global vector + if (.not. allocated(insertIdxs)) then + allocate (insertIdxs(3)) + insertIdxs = (/3 * jj - 2, 3 * jj - 1, 3 * jj/) + else + insertIdxs = [insertIdxs, (/3 * jj - 2, 3 * jj - 1, 3 * jj/)] + end if + + end do nodeLoop + + end do zoneLoop + + ! Insert the rotated coordinates into the global vector + ! We only do this collective communication once at the end if insertCoords is allocated + if (allocated(insertCoords_d)) then + call VecSetValues(dXv, size(insertCoords_d), insertIdxs - 1, insertCoords_d, INSERT_VALUES, ierr) + call EChk(ierr, __FILE__, __LINE__) + end if + + ! Need to call VecAssemblyBegin/End + call VecAssemblyBegin(dXv, ierr) + call EChk(ierr, __FILE__, __LINE__) + call VecAssemblyEnd(dXv, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Restore the array to deallocate the coordsLocal pointer + call VecRestoreArrayF90(dXv, coordsLocal_d, ierr) + call EChk(ierr, __FILE__, __LINE__) + +end subroutine copyRotateVolumeCoordinates_d + +subroutine copyRotateVolumeCoordinates_b() + use precision + use gridData + use cgnsGrid + implicit none + + ! Working variables + integer(kind=intType) :: i, j, k, ii, jj, istart, iend, ierr, iiLocal + integer(kind=intType) :: iZone, nZones, iBoco, ptRange(3, 2), rotDir, iPlane + integer(kind=intType) :: iNode, dir1, dir2, rotOffset, nNodes, mirrorIdx + integer(kind=intType), dimension(:), allocatable :: insertIdxs + real(kind=realType) :: Mi(3, 3), coords_b(3), rotCoords_b(3) + real(kind=realType), dimension(:), allocatable :: insertCoords_b + real(kind=realType), dimension(:), pointer :: coordsLocal_b + + nZones = size(zones) + + ! Get the rotation matrix for the reverse rotation (-angle) + call getRotationMatrixAngleAxis(-axiSymmAngle, axiSymmAxis, Mi) + + ! Get the range of the volume nodes on this proc + call VecGetOwnershipRange(dXv, istart, iend, ierr) + istart = istart + 1 ! Convert to 1-based indexing + call EChk(ierr, __FILE__, __LINE__) + + ! Get the local array of coordinates + call VecGetArrayF90(dXv, coordsLocal_b, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ii = 0 + jj = 0 + ! Loop over the zones in the grid + zoneLoop: do iZone = 1, nZones + rotDir = 0 ! Reset the rotation direction + + ! Loop over the BCs to see if there is a mirror plane BC + bocoLoop: do iBoco = 1, size(zones(iZone)%bocos) + ! Check if the BCs for this zone contain the mirror family + if (zones(iZone)%bocos(iBoco)%family .eq. mirrorFamily) then + ! Determine the rotation direction + ! This is the direction in which the zone is one cell wide + ! The point range for the BC will be equal in this direction + ! signifying that the zone is one cell wide in this direction + ptRange = zones(iZone)%bocos(iBoco)%ptRange + if (ptRange(1, 1) == ptRange(1, 2)) then + mirrorIdx = ptRange(1, 1) + rotDir = 1 ! i-direction + else if (ptRange(2, 1) == ptRange(2, 2)) then + mirrorIdx = ptRange(2, 1) + rotDir = 2 ! j-direction + else if (ptRange(3, 1) == ptRange(3, 2)) then + mirrorIdx = ptRange(3, 1) + rotDir = 3 ! k-direction + end if + end if + end do bocoLoop + + if (rotDir == 0) then + ! No mirror plane BC found, skip this zone + cycle zoneLoop + end if + + ! Determine the plane directions to loop over and set the rotation + ! plane offset (These come from the loop structure in readStructuredCGNS.F90) + if (rotDir == 1) then + ! The zone is one cell wide in the i-direction (final loop in readStructuredCGNS.F90) + dir1 = zones(izone)%jl + dir2 = zones(izone)%kl + rotOffset = 1 + else if (rotDir == 2) then + ! The zone is one cell wide in the j-direction (second loop in readStructuredCGNS.F90) + dir1 = zones(iZone)%il + dir2 = zones(iZone)%kl + rotOffset = zones(iZone)%il + else if (rotDir == 3) then + ! The zone is one cell wide in the k-direction (first loop in readStructuredCGNS.F90) + dir1 = zones(izone)%il + dir2 = zones(izone)%jl + rotOffset = zones(iZone)%il * zones(iZone)%jl + end if + + ! Number of nodes on a single plane for this zone + nNodes = dir1 * dir2 + + ! Loop over the nodes on both planes in this zone (mirror and rotated) + nodeLoop: do iNode = 1, 2 * nNodes + ii = ii + 1 ! Global index + + if (3 * ii - 2 .lt. istart .or. 3 * ii .gt. iend) then + ! ii is not in the local ownership range so skip + cycle nodeLoop + end if + + ! Determine if the node is on plane 1 or 2 (1 or 0) + iPlane = mod((iNode - 1) / rotOffset + 1, 2) + if (mirrorIdx .eq. 2 .and. iPlane .eq. 1) then + ! On a rotated plane that is plane 1, offset forward + jj = ii + rotOffset + else if (mirrorIdx .eq. 1 .and. iPlane .eq. 0) then + ! On a rotated plane this is plane 0, offset backward + jj = ii - rotOffset + else + ! This node is on the rot plane so skip + cycle nodeLoop + end if + + ! Adjust ii to be in the local indexing + iiLocal = 3 * ii - istart + 1 + + ! Copy the x, y, z coordinates + coords_b(1) = coordsLocal_b(iiLocal - 2) ! x + coords_b(2) = coordsLocal_b(iiLocal - 1) ! y + coords_b(3) = coordsLocal_b(iiLocal) ! z + + ! Set the local coords on the rotated plane to zero after the copy + coordsLocal_b(iiLocal - 2:iiLocal) = zero + + ! Apply the rotation + rotCoords_b = matmul(Mi, coords_b) + + ! Add the rotated coords to the array of coords to insert + ! into the global vector + if (.not. allocated(insertCoords_b)) then + allocate (insertCoords_b(3)) + insertCoords_b = rotCoords_b + else + insertCoords_b = [insertCoords_b, rotCoords_b] + end if + + ! Add the index to the array of indices to insert + ! into the global vector + if (.not. allocated(insertIdxs)) then + allocate (insertIdxs(3)) + insertIdxs = (/3 * jj - 2, 3 * jj - 1, 3 * jj/) + else + insertIdxs = [insertIdxs, (/3 * jj - 2, 3 * jj - 1, 3 * jj/)] + end if + + end do nodeLoop + + end do zoneLoop + + ! Add the rotated coordinates into the global vector + ! We only do this collective communication once at the end if insertCoords is allocated + if (allocated(insertCoords_b)) then + call VecSetValues(dXv, size(insertCoords_b), insertIdxs - 1, insertCoords_b, ADD_VALUES, ierr) + call EChk(ierr, __FILE__, __LINE__) + end if + + ! Need to call VecAssemblyBegin/End + call VecAssemblyBegin(dXv, ierr) + call EChk(ierr, __FILE__, __LINE__) + call VecAssemblyEnd(dXv, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Restore the array to deallocate the coordsLocal pointer + call VecRestoreArrayF90(dXv, coordsLocal_b, ierr) + call EChk(ierr, __FILE__, __LINE__) + +end subroutine copyRotateVolumeCoordinates_b diff --git a/src/warp/warpDeriv.F90 b/src/warp/warpDeriv.F90 index 469f44e..cc97bec 100644 --- a/src/warp/warpDeriv.F90 +++ b/src/warp/warpDeriv.F90 @@ -48,6 +48,25 @@ subroutine warpDeriv(dXv_f, ndof_warp) ! dXv_f is the actual reverse seed so copy: XvPtrb(:) = dXv_f + if (axiSymm) then + call VecGetOwnershipRange(dXv, istart, iend, ierr) + call EChk(ierr, __FILE__, __LINE__) + + call VecSetValues(dXv, ndof_warp, (/(j, j=istart, iend)/), dXv_f, INSERT_VALUES, ierr) + call EChk(ierr, __FILE__, __LINE__) + + call VecAssemblyBegin(dXv, ierr) + call EChk(ierr, __FILE__, __LINE__) + call VecAssemblyEnd(dXv, ierr) + call EChk(ierr, __FILE__, __LINE__) + + call copyRotateVolumeCoordinates_b() + + ! Copy the rotated seeds into XvPtrb + call VecGetArrayF90(dXv, XvPtrb, ierr) + call EChk(ierr, __FILE__, __LINE__) + end if + ! Zero the output surface derivative vector dXs call vecZeroEntries(dXs, ierr) call EChk(ierr, __FILE__, __LINE__) @@ -75,7 +94,7 @@ subroutine warpDeriv(dXv_f, ndof_warp) ! Numerator derivative oden = one / denominator(j) - numb = xvptrb(3 * j - 2:3 * j) * oden + numb = XvPtrb(3 * j - 2:3 * j) * oden call getMirrorNumerator_b(numb, kk) if (evalMode == EVAL_EXACT) then @@ -113,9 +132,12 @@ subroutine warpDeriv(dXv_f, ndof_warp) deallocate (dxssummed) ! Deallocate the extra space - deallocate (XvPtrb, XsPtrb) + deallocate (XsPtrb) ! Restore all the arrays + call VecRestoreArrayF90(dXv, XvPtrb, ierr) + call EChk(ierr, __FILE__, __LINE__) + call VecRestoreArrayF90(XsLocal, XsPtr, ierr) call EChk(ierr, __FILE__, __LINE__) diff --git a/src/warp/warpDerivFwd.F90 b/src/warp/warpDerivFwd.F90 index fc58025..91da3f7 100644 --- a/src/warp/warpDerivFwd.F90 +++ b/src/warp/warpDerivFwd.F90 @@ -11,7 +11,7 @@ subroutine warpDerivFwd(Xsdot, cDof, Xvdot, meshDOF) real(kind=realType), intent(in) :: xsdot(cdof) ! Output Arguments - real(kind=realType), intent(inout), dimension(meshDOF) :: XvDot + real(kind=realType), intent(out), dimension(meshDOF) :: XvDot #ifndef USE_COMPLEX ! Working Data integer(kind=intType) :: i, j, kk, istart, iend, ierr, nVol @@ -30,6 +30,10 @@ subroutine warpDerivFwd(Xsdot, cDof, Xvdot, meshDOF) call EChk(ierr, __FILE__, __LINE__) end do + ! Create the volume seed vector as a copy of the volume nodes + call VecGetArrayF90(dXv, XvPtrd, ierr) + call EChk(ierr, __FILE__, __LINE__) + ! While we only set local values, we STILL have to call ! VecAssemblyBegin/End call VecAssemblyBegin(dXs, ierr) @@ -58,7 +62,7 @@ subroutine warpDerivFwd(Xsdot, cDof, Xvdot, meshDOF) call allocDerivValues(mytrees(1)%tp) call zeroDeriv(mytrees(1)%tp%root) - nVol = size(XvPtr) / 3 + nVol = size(XvPtrd) / 3 call computeNodalProperties_d(mytrees(1)%tp, .False.) mytrees(1)%tp%Ldef = mytrees(1)%tp%Ldef0 * LdefFact mytrees(1)%tp%alphaToBexp = alpha**bExp @@ -93,12 +97,15 @@ subroutine warpDerivFwd(Xsdot, cDof, Xvdot, meshDOF) updateLoop: do j = 1, nVol oden = one / denominator(j) - Xvdot(3 * j - 2) = numerator(1, j) * oden - Xvdot(3 * j - 1) = numerator(2, j) * oden - Xvdot(3 * j) = numerator(3, j) * oden + XvPtrd(3 * j - 2) = numerator(1, j) * oden + XvPtrd(3 * j - 1) = numerator(2, j) * oden + XvPtrd(3 * j) = numerator(3, j) * oden end do updateLoop ! Restore all the arrays + call VecRestoreArrayF90(dXv, XvPtrd, ierr) + call EChk(ierr, __FILE__, __LINE__) + call VecRestoreArrayF90(XsLocal, XsPtr, ierr) call EChk(ierr, __FILE__, __LINE__) @@ -107,5 +114,22 @@ subroutine warpDerivFwd(Xsdot, cDof, Xvdot, meshDOF) call VecRestoreArrayF90(Xv0, Xv0Ptr, ierr) call EChk(ierr, __FILE__, __LINE__) + + ! Copy and rotate the seeds if we have axisymmetry + if (axiSymm) then + call copyRotateVolumeCoordinates_d() + end if + + call VecGetArrayF90(dXv, XvPtrd, ierr) + call EChk(ierr, __FILE__, __LINE__) + + ! Copy values from XvPtr to XvDot + do i = 1, meshDOF + XvDot(i) = XvPtrd(i) + end do + + call VecRestoreArrayF90(dXv, XvPtrd, ierr) + call EChk(ierr, __FILE__, __LINE__) + #endif end subroutine warpDerivFwd diff --git a/src/warp/warpMesh.F90 b/src/warp/warpMesh.F90 index 7be909d..623d7cd 100644 --- a/src/warp/warpMesh.F90 +++ b/src/warp/warpMesh.F90 @@ -79,5 +79,10 @@ subroutine warpMesh() call VecRestoreArrayF90(Xv, XvPtr, ierr) call EChk(ierr, __FILE__, __LINE__) + ! If we have an axisymmetric mesh we need to rotate and copy + if (axiSymm) then + call copyRotateVolumeCoordinates() + end if + end subroutine warpMesh