Writing objective function sensitivities into OpenFOAM cell format when running in parallel

[uttamcadambi07]

Hi,

I am faced with a bit of an issue when it comes to writing out the objective function sensitivities.

I understand that the file "totalDerivHist.txt" writes out the total derivative of the objective function with respect to the design variable in a list format.

But I want this in the OpenFOAM cell ordered format so that I can visualize it on Paraview. I did this option in the runScript.py

"writeSensMap": ["beta"]

And surely, I was getting the cell ordered dictionary "sens_FI_beta".

But if you see the file, there is no connectivity information when I decompose it into say 4 cores (see file sens_FI_beta.txt ) you will notice that the dictionary terminates with boundary{} but no information on the patch.

I went into the code for pyDAFoam.py and found this

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
        """

        self.cdRoot()

        workingDir = self.rootDir
        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":
                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" % len(sensList) / 3)
                fSens.write("(\n")
                counterI = 0
                for i in range(len(sensList) / 3):
                    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()

So if you see, this file does the writing of the sensitivity map and there is nothing about the decomposing and reconstruction of the boundary fields.

There are three solutions for this

1. Try to visualize the case on Paraview WITHOUT reconstructing (I believe there is an option for this but I will investigate).
2. Use the totalDerivHist.txt file and try to re-order it to match the cell ordering.
3. Modify the source code to include boundary patch information.

Please let me know what you think and how I can proceed to fix this issue.

[uttamcadambi07]

For now, I do this procedure:

1. I load the designVariables.txt file into python
2. I load the betaFieldInversion dictionary file into python
3. For a given value in designVariables.txt, I check where that value is by comparing it with betaFieldInversion
4. I load the totaDerivHist.txt file into python
5. I use the indices from 3. and reorder the values of the sensitivities
6. I write it out into sens_FI_beta (which is now not in openfoam format but the header and footer can be easily added in gedit or any text editor like vim)

Here is the python script if anyone is interested - source_code.txt

Of course you have to make changes to this

time = '0.00000049'
field_type = 'betaFieldInversion'
file_DV = 'designVariableHist.txt'
file_sens = 'totalDerivHist.txt'
file_beta_ref = time + '/' + field_type

and also make sure that the python script is in the case directory.

[friedenhe]

You can also check the discussion from this thread: #102