Ad actuator (#98)

* Added a script to scale plot3d.

* Changed plot3d2tecplot.py

* Changed the design variable hist output to JSON format.

* Added the AD calculation for ACT derivs.

* Fixed the act deriv in parallel.

* Updated the version.

* Fixed the ACT deriv in parallel. Updated the test.

* Fixed the parallel BC derivs.

* The AOA deriv is working in parallel.

* Updated tests.

dafoam/scripts/dafoam_plot3d2tecplot.py

@@ -2,7 +2,7 @@
 """
 This script converts a Plot3D file to a Tecplot file
 
-Usage: python dafoam-plot3d2tecplot.py plot3d_file_input.xyz tecplot_file_output.dat
+Usage: python dafoam_plot3d2tecplot.py plot3d_file_input.xyz tecplot_file_output.dat
 """
 
 import sys

dafoam/scripts/dafoam_plot3dscale.py

@@ -0,0 +1,35 @@
+#!/usr/bin/env python
+"""
+This script scale the coordinates in a plot3D file
+
+Usage: python dafoam_plot3dscale.py plot3d_file_input.xyz plot3d_file_output.xyz 2 2 2
+This will scale the x, y, and z coordinates by a factor of 2
+"""
+
+import sys
+from pygeo import *
+
+inputFileName = sys.argv[1]
+outputFileName = sys.argv[2]
+scaleX = float(sys.argv[3])
+scaleY = float(sys.argv[4])
+scaleZ = float(sys.argv[5])
+
+print(
+    "Scaling %s to %s with scaleX: %g scaleY: %g scaleZ: %g ...."
+    % (inputFileName, outputFileName, scaleX, scaleY, scaleZ)
+)
+
+ffd = pyBlock("plot3d", fileName=inputFileName, FFD=True)
+
+for ivol in range(ffd.nVol):
+    vals = ffd.vols[ivol].coef
+    vals[:, :, :, 0] = vals[:, :, :, 0] * scaleX
+    vals[:, :, :, 1] = vals[:, :, :, 1] * scaleY
+    vals[:, :, :, 2] = vals[:, :, :, 2] * scaleZ
+ffd.writePlot3dCoef(outputFileName)
+
+print(
+    "Scaling %s to %s with scaleX: %g scaleY: %g scaleZ: %g Done!"
+    % (inputFileName, outputFileName, scaleX, scaleY, scaleZ)
+)

setup.py

@@ -31,6 +31,7 @@
         "dafoam/scripts/dafoam_matgetvalues.py",
         "dafoam/scripts/dafoam_vecgetvalues.py",
         "dafoam/scripts/dafoam_plot3d2tecplot.py",
+        "dafoam/scripts/dafoam_plot3dscale.py",
     ],
     install_requires=[
         "numpy>=1.16.4",

src/adjoint/DAFvSource/DAFvSource.C

@@ -110,6 +110,52 @@ bool DAFvSource::writeData(Ostream& os) const
     // do nothing
     return true;
 }
+
+void DAFvSource::syncDAOptionToActuatorDVs()
+{
+    /*
+    Description:
+        Synchronize the values in DAOption and actuatorDiskDVs_. 
+        We need to synchronize the values defined in fvSource from DAOption to actuatorDiskDVs_
+        NOTE: we need to call this function whenever we change the actuator design variables 
+        during optimization. This is needed because we need to use actuatorDiskDVs_ in AD 
+    */
+
+    // now we need to initialize actuatorDiskDVs_
+    dictionary fvSourceSubDict = daOption_.getAllOptions().subDict("fvSource");
+    word diskName0 = fvSourceSubDict.toc()[0];
+    word source0 = fvSourceSubDict.subDict(diskName0).getWord("source");
+
+    if (source0 == "cylinderAnnulusSmooth")
+    {
+        forAll(fvSourceSubDict.toc(), idxI)
+        {
+            word diskName = fvSourceSubDict.toc()[idxI];
+
+            // sub dictionary with all parameters for this disk
+            dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
+
+            // now read in all parameters for this actuator disk
+            scalarList centerList;
+            diskSubDict.readEntry<scalarList>("center", centerList);
+
+            // we have 9 design variables for each disk
+            scalarList dvList(9);
+            dvList[0] = centerList[0];
+            dvList[1] = centerList[1];
+            dvList[2] = centerList[2];
+            dvList[3] = diskSubDict.getScalar("outerRadius");
+            dvList[4] = diskSubDict.getScalar("innerRadius");
+            dvList[5] = diskSubDict.getScalar("scale");
+            dvList[6] = diskSubDict.getScalar("POD");
+            dvList[7] = diskSubDict.getScalar("expM");
+            dvList[8] = diskSubDict.getScalar("expN");
+
+            // set actuatorDiskDVs_
+            actuatorDiskDVs_.set(diskName, dvList);
+        }
+    }
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/adjoint/DAFvSource/DAFvSource.H

@@ -32,8 +32,7 @@ namespace Foam
 \*---------------------------------------------------------------------------*/
 
 class DAFvSource
-:
-    public regIOobject
+    : public regIOobject
 {
 
 private:
@@ -44,10 +43,9 @@ private:
     void operator=(const DAFvSource&);
 
 protected:
-
     /// model name
     const word& modelType_;
-    
+
     /// fvMesh
     const fvMesh& mesh_;
 
@@ -60,6 +58,9 @@ protected:
     /// DAIndex object
     const DAIndex& daIndex_;
 
+    /// the list of design variables for all the actuator disks
+    HashTable<List<scalar>> actuatorDiskDVs_;
+
 public:
     /// Runtime type information
     TypeName("DAFvSource");
@@ -104,6 +105,26 @@ public:
     /// compute the FvSource term
     virtual void calcFvSource(volVectorField& fvSource);
 
+    /// set a new value to the actuator disk design variable
+    void setActuatorDVs(
+        const word diskName,
+        const label dvI,
+        const scalar val)
+    {
+        actuatorDiskDVs_[diskName][dvI] = val;
+    }
+
+    /// get the value from the actuator disk design variable
+    scalar getActuatorDVs(
+        const word diskName,
+        const label dvI)
+    {
+        return actuatorDiskDVs_[diskName][dvI];
+    }
+
+    /// synchronize the values in DAOption and actuatorDiskDVs_
+    void syncDAOptionToActuatorDVs();
+
     /// virtual function for regIOobject
     bool writeData(Ostream& os) const;
 };

src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C

@@ -27,6 +27,12 @@ DAFvSourceActuatorDisk::DAFvSourceActuatorDisk(
     this->calcFvSourceCellIndices(fvSourceCellIndices_);
 
     printInterval_ = daOption.getOption<label>("printInterval");
+
+    // now we need to initialize actuatorDiskDVs_ by synchronizing the values
+    // defined in fvSource from DAOption to actuatorDiskDVs_
+    // NOTE: we need to call this function whenever we change the actuator
+    // design variables during optimization
+    this->syncDAOptionToActuatorDVs();
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -213,21 +219,20 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource)
             word diskName = fvSourceSubDict.toc()[idxI];
             dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
 
-            scalarList centerList;
             scalarList direction;
-            diskSubDict.readEntry<scalarList>("center", centerList);
             diskSubDict.readEntry<scalarList>("direction", direction);
-            vector center = {centerList[0], centerList[1], centerList[2]};
             vector dirNorm = {direction[0], direction[1], direction[2]};
             dirNorm = dirNorm / mag(dirNorm);
-            scalar outerRadius = diskSubDict.getScalar("outerRadius");
-            scalar innerRadius = diskSubDict.getScalar("innerRadius");
+            vector center = {
+                actuatorDiskDVs_[diskName][0], actuatorDiskDVs_[diskName][1], actuatorDiskDVs_[diskName][2]};
+            scalar outerRadius = actuatorDiskDVs_[diskName][3];
+            scalar innerRadius = actuatorDiskDVs_[diskName][4];
             word rotDir = diskSubDict.getWord("rotDir");
-            scalar scale = diskSubDict.getScalar("scale");
-            scalar POD = diskSubDict.getScalar("POD");
+            scalar scale = actuatorDiskDVs_[diskName][5];
+            scalar POD = actuatorDiskDVs_[diskName][6];
             scalar eps = diskSubDict.getScalar("eps");
-            scalar expM = diskSubDict.getScalar("expM");
-            scalar expN = diskSubDict.getScalar("expN");
+            scalar expM = actuatorDiskDVs_[diskName][7];
+            scalar expN = actuatorDiskDVs_[diskName][8];
             scalar rStarMin = diskSubDict.lookupOrDefault<scalar>("rStarMin", 0.02);
             scalar rStarMax = diskSubDict.lookupOrDefault<scalar>("rStarMax", 0.98);
             scalar epsR = diskSubDict.lookupOrDefault<scalar>("epsR", 0.02);
@@ -295,22 +300,22 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource)
                 if (rStar < rStarMin)
                 {
                     scalar dR2 = (rStar - rStarMin) * (rStar - rStarMin);
-                    scalar fR = fRMin * exp(-dR2 / epsR / epsR) ;
-                    fAxial = fR * exp(-dA2 / eps / eps) ;
+                    scalar fR = fRMin * exp(-dR2 / epsR / epsR);
+                    fAxial = fR * exp(-dA2 / eps / eps);
                     fCirc = fAxial * POD / constant::mathematical::pi / rStarMin;
                 }
                 else if (rStar >= rStarMin && rStar <= rStarMax)
                 {
                     scalar fR = pow(rStar, expM) * pow(1.0 - rStar, expN) * scale;
-                    fAxial = fR * exp(-dA2 / eps / eps) ;
+                    fAxial = fR * exp(-dA2 / eps / eps);
                     // we use Hoekstra's method to calculate the fCirc based on fAxial
                     fCirc = fAxial * POD / constant::mathematical::pi / rPrime;
                 }
                 else
                 {
                     scalar dR2 = (rStar - rStarMax) * (rStar - rStarMax);
-                    scalar fR = fRMax * exp(-dR2 / epsR / epsR) ;
-                    fAxial = fR * exp(-dA2 / eps / eps) ;
+                    scalar fR = fRMax * exp(-dR2 / epsR / epsR);
+                    fAxial = fR * exp(-dA2 / eps / eps);
                     fCirc = fAxial * POD / constant::mathematical::pi / rStarMax;
                 }
 

src/adjoint/DAObjFunc/DAObjFunc.C

@@ -287,7 +287,6 @@ scalar DAObjFunc::masterFunction(
     return objFuncValue;
 }
 
-/// calcluate the value of objective function
 scalar DAObjFunc::getObjFuncValue()
 {
     /*

src/adjoint/DAObjFunc/DAObjFunc.H

@@ -174,7 +174,7 @@ public:
         scalarList& objFuncCellValues,
         scalar& objFuncValue) = 0;
 
-    /// calcluate the value of objective function
+    /// calculate the value of objective function
     scalar getObjFuncValue();
 
     /// the master function to compute objective function given the state and point vectors

src/adjoint/DAObjFunc/DAObjFuncForce.C

@@ -228,12 +228,12 @@ void DAObjFuncForce::updateForceDir(vector& forceDir)
         {
             FatalErrorIn("") << "boundaryType: " << U.boundaryField()[patchI].type()
                              << " not supported!"
-                             << "Avaiable options are: fixedValue, inletOutlet"
+                             << "Available options are: fixedValue, inletOutlet"
                              << abort(FatalError);
         }
     }
 
-    // need to reduce the sum of force across all processors, this is becasue some of
+    // need to reduce the sum of force across all processors, this is because some of
     // the processor might not own the inoutRefPatchName_ so their flowDir will be -1e16, but
     // when calling the following reduce function, they will get the correct flowDir
     // computed by other processors

src/adjoint/DAPartDeriv/DAPartDerivdRdACTD.C

@@ -130,6 +130,9 @@ void DAPartDerivdRdACTD::calcPartDerivMat(
     scalar expM = diskModelSubDict.getScalar("expM");
     scalar expN = diskModelSubDict.getScalar("expN");
 
+    DAFvSource& fvSource = const_cast<DAFvSource&>(
+        mesh_.thisDb().lookupObject<DAFvSource>("DAFvSource"));
+
     for (label i = 0; i < nActDVs_; i++)
     {
 
@@ -189,6 +192,9 @@ void DAPartDerivdRdACTD::calcPartDerivMat(
             diskModelSubDict.set("expN", expN);
         }
 
+        // we need to synchronize the DAOption to actuatorDVs
+        fvSource.syncDAOptionToActuatorDVs();
+
         // Info<<daOption_.getAllOptions().subDict("fvSource")<<endl;
 
         // compute residual
@@ -250,13 +256,15 @@ void DAPartDerivdRdACTD::calcPartDerivMat(
             diskModelSubDict.set("expN", expN);
         }
 
+        // we need to synchronize the DAOption to actuatorDVs
+        fvSource.syncDAOptionToActuatorDVs();
+
         // Info<<daOption_.getAllOptions().subDict("fvSource")<<endl;
 
         // compute residual partial using finite-difference
         VecAXPY(resVec, -1.0, resVecRef);
         VecScale(resVec, rDeltaValue);
 
-
         // assign resVec to jacMat
         PetscInt Istart, Iend;
         VecGetOwnershipRange(resVec, &Istart, &Iend);

src/adjoint/DAResidual/DAResidualPisoFoam.H

@@ -15,7 +15,6 @@
 #include "addToRunTimeSelectionTable.H"
 #include "pisoControl.H"
 #include "adjustPhi.H"
-#include "DAFvSource.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C

@@ -146,7 +146,7 @@ label DARhoSimpleFoam::solvePrimal(
         if (printToScreen)
         {
             daTurbulenceModelPtr_->printYPlus();
-            
+
             this->printAllObjFuncs();
 
             Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
@@ -155,11 +155,10 @@ label DARhoSimpleFoam::solvePrimal(
         }
 
         runTime.write();
-
     }
 
     this->writeAssociatedFields();
-    
+
     this->calcPrimalResidualStatistics("print");
 
     // primal converged, assign the OpenFoam fields to the state vec wVec

src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.H

@@ -35,7 +35,6 @@ class DARhoSimpleFoam
 {
 
 protected:
-
     /// simple pointer
     autoPtr<simpleControl> simplePtr_;
 
@@ -66,7 +65,7 @@ protected:
     /// DASource pointer
     autoPtr<DAFvSource> daFvSourcePtr_;
 
-    /// fvSource term 
+    /// fvSource term
     autoPtr<volVectorField> fvSourcePtr_;
 
     /// fvSource term for the energy equation
@@ -105,7 +104,6 @@ public:
     virtual label solvePrimal(
         const Vec xvVec,
         Vec wVec);
-
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //