Fixed a bug for vector sens map and using vector field variables as D…

…V. (#513)
mdolab · Nov 1, 2023 · 767274c · 767274c
1 parent 2e981bd
commit 767274c
Show file tree

Hide file tree

Showing 4 changed files with 33 additions and 9 deletions.
diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py
@@ -1814,13 +1814,14 @@ def writeFieldSensitivityMap(self, objFuncName, designVarName, solutionTime, fie
                 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" % len(sensList) / 3)
+                fSens.write("%d\n" % size)
                 fSens.write("(\n")
                 counterI = 0
-                for i in range(len(sensList) / 3):
+                for i in range(size):
                     fSens.write("(")
                     for j in range(3):
                         fSens.write("%g " % sensList[counterI])

diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C
@@ -8048,11 +8048,11 @@ void DASolver::setFieldValue4GlobalCellI(
 
     if (meshPtr_->thisDb().foundObject<volVectorField>(fieldName))
     {
-        if (daIndexPtr_->globalCellVectorNumbering.isLocal(globalCellI))
+        if (daIndexPtr_->globalCellNumbering.isLocal(globalCellI))
         {
             volVectorField& field =
                 const_cast<volVectorField&>(meshPtr_->thisDb().lookupObject<volVectorField>(fieldName));
-            label localCellI = daIndexPtr_->globalCellVectorNumbering.toLocal(globalCellI);
+            label localCellI = daIndexPtr_->globalCellNumbering.toLocal(globalCellI);
             field[localCellI][compI] = val;
         }
     }

diff --git a/tests/refs/DAFoam_Test_DASimpleFoamFieldRef.txt b/tests/refs/DAFoam_Test_DASimpleFoamFieldRef.txt
@@ -1,13 +1,15 @@
 Dictionary Key: FI
-@value     5.67958608369146 1e-08 1e-10
+@value    5.679586083422453 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
 Dictionary Key: FI
 Dictionary Key: alpha
 @value                    0 0.0001 1e-06
 Dictionary Key: alphaPorosity
-@value    329.0823053537178 0.0001 1e-06
+@value     329.081743624962 0.0001 1e-06
 Dictionary Key: beta
-@value    2.497448257761139 0.0001 1e-06
+@value    2.497448257892173 0.0001 1e-06
+Dictionary Key: fvSource
+@value    50.51907226688485 0.0001 1e-06
 Dictionary Key: fail
 @value                    0 0.0001 1e-06
diff --git a/tests/runTests_DASimpleFoamField.py b/tests/runTests_DASimpleFoamField.py
@@ -26,7 +26,7 @@
 if gcomm.rank == 0:
     os.system("rm -rf 0 processor*")
     os.system("cp -r 0.incompressible 0")
-    # rename reference U field 
+    # rename reference U field
     os.system("cp 0/varRefFieldInversion 0/UData")
 
 replace_text_in_file("0/varRefFieldInversion", "    object      varRefFieldInversion;", "    object      UData;")
@@ -45,7 +45,7 @@
     "designSurfaces": ["wing"],
     "primalMinResTol": 1e-12,
     "writeJacobians": ["all"],
-    "writeSensMap": ["betaSA", "alphaPorosity"],
+    "writeSensMap": ["betaSA", "alphaPorosity", "fvSource"],
     "primalBC": {
         "U0": {"variable": "U", "patches": ["inout"], "value": [U0, 0.0, 0.0]},
         "p0": {"variable": "p", "patches": ["inout"], "value": [p0]},
@@ -107,6 +107,7 @@
         "beta": {"designVarType": "Field", "fieldName": "betaFieldInversion", "fieldType": "scalar"},
         "alphaPorosity": {"designVarType": "Field", "fieldName": "alphaPorosity", "fieldType": "scalar"},
         "alpha": {"designVarType": "AOA", "patches": ["inout"], "flowAxis": "x", "normalAxis": "y"},
+        "fvSource": {"designVarType": "Field", "fieldName": "fvSource", "fieldType": "vector"},
     },
 }
 
@@ -132,16 +133,29 @@ def alpha(val, geo):
     DASolver.setOption("primalBC", {"U0": {"variable": "U", "patches": ["inout"], "value": inletU}})
     DASolver.updateDAOption()
 
+
 def betaFieldInversion(val, geo):
     for idxI, v in enumerate(val):
         DASolver.setFieldValue4GlobalCellI(b"betaFieldInversion", v, idxI)
         DASolver.updateBoundaryConditions(b"betaFieldInversion", b"scalar")
 
+
 def alphaPorosity(val, geo):
     for idxI, v in enumerate(val):
         DASolver.setFieldValue4GlobalCellI(b"alphaPorosity", v, idxI)
         DASolver.updateBoundaryConditions(b"alphaPorosity", b"scalar")
 
+
+def fvSource(val, geo):
+    nVals = int(len(val) // 3)
+    for cellI in range(nVals):
+        for i in range(3):
+            idxI = cellI * 3 + i
+            v = val[idxI]
+            DASolver.setFieldValue4GlobalCellI(b"fvSource", float(v), cellI, i)
+            DASolver.updateBoundaryConditions(b"fvSource", b"vector")
+
+
 # select points
 DVGeo.addGlobalDV("alpha", [alpha0], alpha, lower=-10.0, upper=10.0, scale=1.0)
 nCells = 4032
@@ -151,6 +165,9 @@ def alphaPorosity(val, geo):
 alphaPorosity0 = np.zeros(nCells, dtype="d")
 DVGeo.addGlobalDV("alphaPorosity", value=alphaPorosity0, func=alphaPorosity, lower=0, upper=100.0, scale=1.0)
 
+fvSource0 = np.zeros(nCells * 3, dtype="d")
+DVGeo.addGlobalDV("fvSource", value=fvSource0, func=fvSource, lower=-10, upper=10.0, scale=1.0)
+
 # DAFoam
 DASolver = PYDAFOAM(options=aeroOptions, comm=gcomm)
 DASolver.setDVGeo(DVGeo)
@@ -195,6 +212,10 @@ def alphaPorosity(val, geo):
     funcsSens["FI"]["alphaPorosity"] = np.zeros(1, "d")
     funcsSens["FI"]["alphaPorosity"][0] = np.linalg.norm(alphaPorositySens)
 
+    fvSourceSens = funcsSens["FI"]["fvSource"]
+    funcsSens["FI"]["fvSource"] = np.zeros(1, "d")
+    funcsSens["FI"]["fvSource"][0] = np.linalg.norm(fvSourceSens)
+
     # Do not consider alpha deriv, just assign it to 0
     funcsSens["FI"]["alpha"] = 0