Skip to content

Commit

Permalink
Re-implemented the sensitivity map and changed the renameSolution for…
Browse files Browse the repository at this point in the history 
…mat (#566)

* Added new sensMap funcs.

* Made the renameSolution consistent with OF's naming convention.
  • Loading branch information
@friedenhe
friedenhe authored Jan 19, 2024
1 parent 458148b commit aa18fb4
Show file tree
Hide file tree
Showing 5 changed files with 318 additions and 222 deletions.
244 changes: 22 additions & 222 deletions dafoam/pyDAFoam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1556,7 +1556,7 @@ def getInternalDVDict(self):
internalDVDict = {}
for dvName in self.internalDV:
internalDVDict[dvName] = self.internalDV[dvName]["value"]

return internalDVDict

def addInternalDV(self, dvName, dvInit, dvFunc, lower, upper, scale):
Expand Down Expand Up @@ -1872,224 +1872,6 @@ def _writeOFCaseFiles(self):

return

def writeFieldSensitivityMap(self, objFuncName, designVarName, solutionTime, fieldType, sensVec):
"""
Save the field sensitivity dObjFunc/dDesignVar map to disk.
Parameters
----------
objFuncName : str
Name of the objective function
designVarName : str
Name of the design variable
solutionTime : float
The solution time where the sensitivity will be save
fieldType : str
The type of the field, either scalar or vector
sensVec : petsc vec
The Petsc vector that contains the sensitivity
"""

workingDir = os.getcwd()
if self.parallel:
sensDir = "processor%d/%.8f/" % (self.rank, solutionTime)
else:
sensDir = "%.8f/" % solutionTime

sensDir = os.path.join(workingDir, sensDir)

sensList = []
Istart, Iend = sensVec.getOwnershipRange()
for idxI in range(Istart, Iend):
sensList.append(sensVec[idxI])

# write sens
if not os.path.isfile(os.path.join(sensDir, "sens_%s_%s" % (objFuncName, designVarName))):
fSens = open(os.path.join(sensDir, "sens_%s_%s" % (objFuncName, designVarName)), "w")
if fieldType == "scalar":
self._writeOpenFoamHeader(fSens, "volScalarField", sensDir, "sens_%s_%s" % (objFuncName, designVarName))
fSens.write("dimensions [0 0 0 0 0 0 0];\n")
fSens.write("internalField nonuniform List<scalar>\n")
fSens.write("%d\n" % len(sensList))
fSens.write("(\n")
for i in range(len(sensList)):
fSens.write("%g\n" % sensList[i])
fSens.write(")\n")
fSens.write(";\n")
elif fieldType == "vector":
size = len(sensList) // 3
self._writeOpenFoamHeader(fSens, "volVectorField", sensDir, "sens_%s_%s" % (objFuncName, designVarName))
fSens.write("dimensions [0 0 0 0 0 0 0];\n")
fSens.write("internalField nonuniform List<vector>\n")
fSens.write("%d\n" % size)
fSens.write("(\n")
counterI = 0
for i in range(size):
fSens.write("(")
for j in range(3):
fSens.write("%g " % sensList[counterI])
counterI = counterI + 1
fSens.write(")\n")
fSens.write(")\n")
fSens.write(";\n")
else:
raise Error("fieldType %s not valid! Options are: scalar or vector" % fieldType)

fSens.write("boundaryField\n")
fSens.write("{\n")
fSens.write(' "(.*)"\n')
fSens.write(" {\n")
fSens.write(" type zeroGradient;\n")
fSens.write(" }\n")
fSens.write("}\n")
fSens.close()

def writeSurfaceSensitivityMap(self, objFuncName, designVarName, solutionTime):
"""
Save the sensitivity dObjFunc/dXs map to disk. where Xs is the wall surface mesh coordinate
Parameters
----------
objFuncName : str
Name of the objective function
designVarName : str
Name of the design variable
solutionTime : float
The solution time where the sensitivity will be save
"""

dFdXs = self.mesh.getdXs()
dFdXs = self.mapVector(dFdXs, self.allWallsGroup, self.allWallsGroup)

pts = self.getSurfaceCoordinates(self.allWallsGroup)
conn, faceSizes = self.getSurfaceConnectivity(self.allWallsGroup)
conn = np.array(conn).flatten()

workingDir = os.getcwd()
if self.parallel:
meshDir = "processor%d/%.11f/polyMesh/" % (self.rank, solutionTime)
sensDir = "processor%d/%.11f/" % (self.rank, solutionTime)
else:
meshDir = "%.11f/polyMesh/" % solutionTime
sensDir = "%.11f/" % solutionTime

meshDir = os.path.join(workingDir, meshDir)
sensDir = os.path.join(workingDir, sensDir)

if not os.path.isdir(sensDir):
try:
os.mkdir(sensDir)
except Exception:
raise Error("Can not make a directory at %s" % sensDir)
if not os.path.isdir(meshDir):
try:
os.mkdir(meshDir)
except Exception:
raise Error("Can not make a directory at %s" % meshDir)

# write points
if not os.path.isfile(os.path.join(meshDir, "points")):
fPoints = open(os.path.join(meshDir, "points"), "w")
self._writeOpenFoamHeader(fPoints, "dictionary", meshDir, "points")
fPoints.write("%d\n" % len(pts))
fPoints.write("(\n")
for i in range(len(pts)):
fPoints.write("(%g %g %g)\n" % (float(pts[i][0]), float(pts[i][1]), float(pts[i][2])))
fPoints.write(")\n")
fPoints.close()

# write faces
if not os.path.isfile(os.path.join(meshDir, "faces")):
fFaces = open(os.path.join(meshDir, "faces"), "w")
self._writeOpenFoamHeader(fFaces, "dictionary", meshDir, "faces")
counterI = 0
fFaces.write("%d\n" % len(faceSizes))
fFaces.write("(\n")
for i in range(len(faceSizes)):
fFaces.write("%d(" % faceSizes[i])
for j in range(faceSizes[i]):
fFaces.write(" %d " % conn[counterI])
counterI += 1
fFaces.write(")\n")
fFaces.write(")\n")
fFaces.close()

# write owner
if not os.path.isfile(os.path.join(meshDir, "owner")):
fOwner = open(os.path.join(meshDir, "owner"), "w")
self._writeOpenFoamHeader(fOwner, "dictionary", meshDir, "owner")
fOwner.write("%d\n" % len(faceSizes))
fOwner.write("(\n")
for i in range(len(faceSizes)):
fOwner.write("0\n")
fOwner.write(")\n")
fOwner.close()

# write neighbour
if not os.path.isfile(os.path.join(meshDir, "neighbour")):
fNeighbour = open(os.path.join(meshDir, "neighbour"), "w")
self._writeOpenFoamHeader(fNeighbour, "dictionary", meshDir, "neighbour")
fNeighbour.write("%d\n" % len(faceSizes))
fNeighbour.write("(\n")
for i in range(len(faceSizes)):
fNeighbour.write("0\n")
fNeighbour.write(")\n")
fNeighbour.close()

# write boundary
if not os.path.isfile(os.path.join(meshDir, "boundary")):
fBoundary = open(os.path.join(meshDir, "boundary"), "w")
self._writeOpenFoamHeader(fBoundary, "dictionary", meshDir, "boundary")
fBoundary.write("1\n")
fBoundary.write("(\n")
fBoundary.write(" allWalls\n")
fBoundary.write(" {\n")
fBoundary.write(" type wall;\n")
fBoundary.write(" nFaces %d;\n" % len(faceSizes))
fBoundary.write(" startFace 0;\n")
fBoundary.write(" }\n")
fBoundary.write(")\n")
fBoundary.close()

# write sens
if not os.path.isfile(os.path.join(sensDir, "sens_%s_%s" % (objFuncName, designVarName))):
fSens = open(os.path.join(sensDir, "sens_%s_%s" % (objFuncName, designVarName)), "w")
self._writeOpenFoamHeader(fSens, "volVectorField", sensDir, "sens_%s_%s" % (objFuncName, designVarName))
fSens.write("dimensions [0 0 0 0 0 0 0];\n")
fSens.write("internalField uniform (0 0 0);\n")

counterI = 0
fSens.write("boundaryField\n")
fSens.write("{\n")
fSens.write(" allWalls\n")
fSens.write(" {\n")
fSens.write(" type wall;\n")
fSens.write(" value nonuniform List<vector>\n")
fSens.write("%d\n" % len(faceSizes))
fSens.write("(\n")
counterI = 0
for i in range(len(faceSizes)):
sensXMean = 0.0
sensYMean = 0.0
sensZMean = 0.0
for j in range(faceSizes[i]):
idxI = conn[counterI]
sensXMean += dFdXs[idxI][0]
sensYMean += dFdXs[idxI][1]
sensZMean += dFdXs[idxI][2]
counterI += 1
sensXMean /= faceSizes[i]
sensYMean /= faceSizes[i]
sensZMean /= faceSizes[i]
fSens.write("(%f %f %f)\n" % (sensXMean, sensYMean, sensZMean))
fSens.write(")\n")
fSens.write(";\n")
fSens.write(" }\n")
fSens.write("}\n")
fSens.close()

def writePetscVecMat(self, name, vecMat, mode="Binary"):
"""
Write Petsc vectors or matrices
Expand Down Expand Up @@ -2234,6 +2016,17 @@ def calcTotalDerivsFFD(self, objFuncName, designVarName, dFScaling=1.0, accumula

if self.DVGeo is not None and self.DVGeo.getNDV() > 0:
dFdFFD = self.mapdXvTodFFD(totalDerivXv)
# check if we need to write the sens map
if designVarName in self.getOption("writeSensMap"):
dFdXs = self.mesh.getdXs()
dFdXs = self.mapVector(dFdXs, self.allWallsGroup, self.designSurfacesGroup)
Xs = self.getSurfaceCoordinates(self.allWallsGroup)
dFdXsFlatten = dFdXs.flatten()
XsFlatten = Xs.flatten()
size = len(dFdXsFlatten)
timeName = float(self.nSolveAdjoints) / 1e8
name = "sens_" + objFuncName + "_" + designVarName
self.solver.writeSensMapSurface(name, dFdXsFlatten, XsFlatten, size, timeName)
# assign the total derivative to self.adjTotalDeriv
if self.adjTotalDeriv[objFuncName][designVarName] is None:
self.adjTotalDeriv[objFuncName][designVarName] = np.zeros(nDVs, self.dtype)
Expand Down Expand Up @@ -2285,6 +2078,13 @@ def calcTotalDerivsField(self, objFuncName, designVarName, fieldType, dFScaling=
totalDeriv.scale(-1.0)
totalDeriv.axpy(1.0, dFdField)

# check if we need to save the sensitivity maps
if designVarName in self.getOption("writeSensMap"):
timeName = float(self.nSolveAdjoints) / 1e8
dFdFieldArray = self.vec2Array(totalDeriv)
name = "sens_" + objFuncName + "_" + designVarName
self.solver.writeSensMapField(name, dFdFieldArray, fieldType, timeName)

# assign the total derivative to self.adjTotalDeriv
if self.adjTotalDeriv[objFuncName][designVarName] is None:
self.adjTotalDeriv[objFuncName][designVarName] = np.zeros(nDVs, self.dtype)
Expand Down Expand Up @@ -2377,7 +2177,7 @@ def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accum
if designVarName in xDV:
nDVs = len(xDV[designVarName])
parameters = xDV[designVarName].copy(order="C")

if designVarName in self.internalDV:
nDVs = len(self.internalDV[designVarName]["init"])
parameters = self.internalDV[designVarName]["value"]
Expand Down Expand Up @@ -3422,7 +3222,7 @@ def deletePrevPrimalSolTime(self):
def renameSolution(self, solIndex):
"""
Rename the primal solution folder to specific format for post-processing. The renamed time has the
format like 0.00000001, 0.00000002, etc. One can load these intermediate shapes and fields and
format like 1e-8, 2e-8, etc. One can load these intermediate shapes and fields and
plot them in paraview.
The way it is implemented is that we sort the solution folder and consider the largest time folder
as the solution folder and rename it
Expand Down Expand Up @@ -3456,7 +3256,7 @@ def renameSolution(self, solIndex):
renamed = False
return solutionTime, renamed

distTime = "%.8f" % (solIndex / 1e8)
distTime = "%g" % (solIndex / 1e8)

src = os.path.join(checkPath, solutionTime)
dst = os.path.join(checkPath, distTime)
Expand Down
Loading

0 comments on commit aa18fb4

Please sign in to comment.