OM_DVGEOCOMP tests and forward mode derivatives

.gitignore

@@ -6,4 +6,6 @@ reg_tests/pygeo_reg.orig
 doc/_build
 *testflo_report.out
 input_files
-.isort.cfg
+.isort.cfg
+*.vscode
+tests/reg_tests/reports/

pygeo/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "1.14.0"
+__version__ = "1.15.0"
 
 from .pyNetwork import pyNetwork
 from .pyGeo import pyGeo

pygeo/geo_utils/norm.py

@@ -27,17 +27,17 @@ def cross_b(a, b, crossb):
     bb[1] = bb[1] + a[0] * crossb[2]
     ab[1] = ab[1] - b[0] * crossb[2]
     bb[0] = bb[0] - a[1] * crossb[2]
-    crossb[2] = 0.0
+
     ab[2] = ab[2] + b[0] * crossb[1]
     bb[0] = bb[0] + a[2] * crossb[1]
     ab[0] = ab[0] - b[2] * crossb[1]
     bb[2] = bb[2] - a[0] * crossb[1]
-    crossb[1] = 0.0
+
     ab[1] = ab[1] + b[2] * crossb[0]
     bb[2] = bb[2] + a[1] * crossb[0]
     ab[2] = ab[2] - b[1] * crossb[0]
     bb[1] = bb[1] - a[2] * crossb[0]
-    crossb[0] = 0.0
+
     return ab, bb
 
 
@@ -48,8 +48,8 @@ def dot_b(a, b, dotb):
     ab = np.zeros_like(a)
     bb = np.zeros_like(b)
 
-    ab = ab + b * dotb
-    bb = bb + a * dotb
+    ab = b * dotb
+    bb = a * dotb
 
     return ab, bb
 

pygeo/mphys/mphys_dvgeo.py

@@ -201,7 +201,12 @@ def nom_addPointSet(self, points, ptName, add_output=True, DVGeoName=None, **kwa
         DVGeo = self.nom_getDVGeo(DVGeoName=DVGeoName)
 
         # add the points to the dvgeo object
-        DVGeo.addPointSet(points.reshape(len(points) // 3, 3), ptName, **kwargs)
+        dMaxGlobal = None
+        if isinstance(DVGeo, DVGeometryESP):
+            # DVGeoESP can return a value to check the pointset distribution
+            dMaxGlobal = DVGeo.addPointSet(points.reshape(len(points) // 3, 3), ptName, **kwargs)
+        else:
+            DVGeo.addPointSet(points.reshape(len(points) // 3, 3), ptName, **kwargs)
         self.omPtSetList.append(ptName)
 
         if isinstance(DVGeo, DVGeometry):
@@ -215,6 +220,8 @@ def nom_addPointSet(self, points, ptName, add_output=True, DVGeoName=None, **kwa
             # add an output to the om component
             self.add_output(ptName, distributed=True, val=points.flatten())
 
+        return dMaxGlobal
+
     def nom_add_point_dict(self, point_dict):
         # add every pointset in the dict, and set the ptset name as the key
         for k, v in point_dict.items():
@@ -569,7 +576,32 @@ def nom_addVSPVariable(self, component, group, parm, isComposite=False, DVGeoNam
         if not isComposite:
             self.add_input(dvName, distributed=False, shape=1, val=val)
 
-    def nom_addESPVariable(self, desmptr_name, isComposite=False, DVGeoName=None, **kwargs):
+    def nom_addESPVariable(self, desmptr_name, rows=None, cols=None, dh=0.001, isComposite=False, DVGeoName=None):
+        """
+        Add an ESP design variables to the DVGeometryESP object
+        Wrapper for :meth:`addVariable <.DVGeometryESP.addVariable>`
+        Input parameters are identical to those in wrapped function unless otherwise specified
+
+        Parameters
+        ----------
+        desmptr_name : str
+            See :meth:`addVariable <.DVGeometryESP.addVariable>`
+        rows : list or None, optional
+            See :meth:`addVariable <.DVGeometryESP.addVariable>`
+        cols : list or None, optional
+            See :meth:`addVariable <.DVGeometryESP.addVariable>`
+        dh : float, optional
+            See :meth:`addVariable <.DVGeometryESP.addVariable>`
+        isComposite : bool, optional
+            Whether this DV is to be included in the composite DVs, by default False
+        DVGeoName : string, optional
+            The name of the DVGeo to add DVs to, necessary if there are multiple DVGeo objects
+
+        Raises
+        ------
+        RuntimeError
+            Raised if the underlying DVGeo parameterization is not ESP-based
+        """
         # if we have multiple DVGeos use the one specified by name
         DVGeo = self.nom_getDVGeo(DVGeoName=DVGeoName)
 
@@ -578,7 +610,7 @@ def nom_addESPVariable(self, desmptr_name, isComposite=False, DVGeoName=None, **
             raise RuntimeError(f"Only ESP-based DVGeo objects can use ESP DVs, not type: {type(DVGeo).__name__}")
 
         # actually add the DV to ESP
-        DVGeo.addVariable(desmptr_name, **kwargs)
+        DVGeo.addVariable(desmptr_name, rows=rows, cols=cols, dh=dh)
 
         # get the value
         val = DVGeo.DVs[desmptr_name].value.copy()
@@ -962,10 +994,12 @@ def nom_setConstraintSurface(
         self.DVCon.setSurface(surface, name=name, addToDVGeo=addToDVGeo, DVGeoName=DVGeoName, surfFormat=surfFormat)
 
     def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
-        # only do the computations when we have more than zero entries in d_inputs in the reverse mode
-        ni = len(list(d_inputs.keys()))
+        # only do the computations when we have more than zero entries in d_inputs
+        # in the reverse mode or d_outputs in the forward mode
+        doRev = mode == "rev" and len(list(d_inputs.keys())) > 0
+        doFwd = mode == "fwd" and len(list(d_outputs.keys())) > 0
 
-        if mode == "rev" and ni > 0:
+        if doFwd or doRev:
             # this flag will be set to True after every compute call.
             # if it is true, we assume the design has changed so we re-run the sensitivity update
             # there can be hundreds of calls to this routine due to thickness constraints,
@@ -977,15 +1011,24 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
                 # set the flag to False so we dont run the update again if this is called w/o a compute in between
                 self.update_jac = False
 
+            # Directly do Jacobian vector product with the derivatives from DVConstraints
             for constraintname in self.constraintfuncsens:
                 for dvname in self.constraintfuncsens[constraintname]:
-                    if dvname in d_inputs:
+                    if constraintname in d_outputs and dvname in d_inputs:
                         dcdx = self.constraintfuncsens[constraintname][dvname]
-                        dout = d_outputs[constraintname]
-                        jvtmp = np.dot(np.transpose(dcdx), dout)
-                        d_inputs[dvname] += jvtmp
+                        if doFwd:
+                            din = d_inputs[dvname]
+                            jvtmp = np.dot(dcdx, din)
+                            d_outputs[constraintname] += jvtmp
+                        elif doRev:
+                            dout = d_outputs[constraintname]
+                            jvtmp = np.dot(np.transpose(dcdx), dout)
+                            d_inputs[dvname] += jvtmp
 
             for _, DVGeo in self.DVGeos.items():
+                dvs = DVGeo.getVarNames()
+
+                # Collet point sets from all DVGeos if DVGeometryMulti object
                 if isinstance(DVGeo, DVGeometryMulti):
                     ptSetNames = []
                     for comp in DVGeo.DVGeoDict.keys():
@@ -997,34 +1040,38 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
 
                 for ptSetName in ptSetNames:
                     if ptSetName in self.omPtSetList:
-                        dout = d_outputs[ptSetName].reshape(len(d_outputs[ptSetName]) // 3, 3)
-
-                        # only do the calc. if d_output is not zero on ANY proc
-                        local_all_zeros = np.all(dout == 0)
+                        # Process the seeds
+                        if doFwd:
+                            # Collect the d_inputs associated with the current DVGeo
+                            seeds = {}
+                            for k in d_inputs:
+                                if k in dvs:
+                                    seeds[k] = d_inputs[k]
+                        elif doRev:
+                            seeds = d_outputs[ptSetName].reshape(len(d_outputs[ptSetName]) // 3, 3)
+
+                        # only do the calc. if seeds are not zero on ANY proc
+                        local_all_zeros = np.all(seeds == 0)
                         global_all_zeros = np.zeros(1, dtype=bool)
                         # we need to communicate for this check otherwise we may hang
                         self.comm.Allreduce([local_all_zeros, MPI.BOOL], [global_all_zeros, MPI.BOOL], MPI.LAND)
 
-                        # global_all_zeros is a numpy array of size 1
+                        # Compute the Jacobian vector product
                         if not global_all_zeros[0]:
-                            # TODO totalSensitivityTransProd is broken. does not work with zero surface nodes on a proc
-                            # xdot = DVGeo.totalSensitivityTransProd(dout, ptSetName)
-                            xdot = DVGeo.totalSensitivity(dout, ptSetName)
-
-                            # loop over dvs and accumulate
-                            xdotg = {}
-                            for k in xdot:
-                                # check if this dv is present
-                                if k in d_inputs:
-                                    # do the allreduce
-                                    # TODO reove the allreduce when this is fixed in openmdao
-                                    # reduce the result ourselves for now. ideally, openmdao will do the reduction itself when this is fixed. this is because the bcast is also done by openmdao (pyoptsparse, but regardless, it is not done here, so reduce should also not be done here)
-                                    xdotg[k] = self.comm.allreduce(xdot[k], op=MPI.SUM)
-
-                                    # accumulate in the dict
-                                    # TODO
-                                    # because we only do one point set at a time, we always want the 0th
-                                    # entry of this array since dvgeo always behaves like we are passing
-                                    # in multiple objective seeds with totalSensitivity. we can remove the [0]
-                                    # once we move back to totalSensitivityTransProd
-                                    d_inputs[k] += xdotg[k][0]
+                            if doFwd:
+                                d_outputs[ptSetName] += DVGeo.totalSensitivityProd(seeds, ptSetName)
+                            elif doRev:
+                                # TODO totalSensitivityTransProd is broken. does not work with zero surface nodes on a proc
+                                # xdot = DVGeo.totalSensitivityTransProd(dout, ptSetName)
+                                xdot = DVGeo.totalSensitivity(seeds, ptSetName)
+
+                                # loop over dvs and accumulate
+                                xdotg = {}
+                                for k in xdot:
+                                    # check if this dv is present
+                                    if k in d_inputs:
+                                        # do the allreduce
+                                        xdotg[k] = self.comm.allreduce(xdot[k], op=MPI.SUM)
+
+                                        # accumulate in the dict
+                                        d_inputs[k] += xdotg[k][0]

pygeo/parameterization/DVGeo.py

@@ -2498,7 +2498,7 @@ def totalSensitivityProd(self, vec, ptSetName, config=None):
         newvec = np.zeros(self.getNDV(), self.dtype)
 
         i = 0
-        missingVars = set()  # set of variables
+        missingVars = set(names)  # set of variables to find in this DVGeo or its children
         for vecKey in vec:
             # check if the seed DV is actually a design variable for the DVGeo object
             if vecKey not in names:
@@ -2522,8 +2522,6 @@ def totalSensitivityProd(self, vec, ptSetName, config=None):
                     dv = geoObj.DV_listLocal[key]
                     missingVars.discard(key)
                 else:
-                    # keep track of DVs which are in the full name list but not in this DVGeo object
-                    missingVars.add(key)
                     continue
 
                 if key in vec:

tests/reg_tests/commonUtils.py

@@ -283,3 +283,53 @@ def spanX(val, geo):
         C[i, 0] *= val
 
     geo.restoreCoef(C, axis_key)
+
+
+def getShapeFunc(lidx, direction=None):
+    """
+    Get shape dictionaries for use with shape function DVs. Common to DVGeometry and MPhys DVGeo tests.
+    Requires local index from DVGeo object and optionally a three-element numpy array for the
+    positive direction vector and returns shapes.
+    """
+    shape_1 = {}
+    shape_2 = {}
+
+    k_center = 2
+    i_center = 1
+    n_chord = lidx.shape[0]
+
+    d_up = np.array([0.0, 1.0, 0.0])
+    if direction is not None:
+        d_up = direction
+        d_up /= np.linalg.norm(d_up)
+
+    for kk in [-1, 0, 1]:
+        if kk == 0:
+            k_weight = 1.0
+        else:
+            k_weight = 0.5
+
+        for ii in range(n_chord):
+            # compute the chord weight. we want the shape to peak at i_center
+            if ii == i_center:
+                i_weight = 1.0
+            elif ii < i_center:
+                # we are ahead of the center point
+                i_weight = ii / i_center
+            else:
+                # we are behind the center point
+                i_weight = (n_chord - ii - 1) / (n_chord - i_center - 1)
+
+            # scale direction by the i and k weights
+            d_up_scaled = d_up * k_weight * i_weight
+
+            # get this point's global index and add to the dictionary with the direction vector.
+            gidx_up = lidx[ii, 1, kk + k_center]
+            gidx_down = lidx[ii, 0, kk + k_center]
+
+            shape_1[gidx_up] = d_up_scaled
+            # the lower face is perturbed with a separate dictionary
+            shape_2[gidx_down] = -d_up_scaled
+
+    shapes = [shape_1, shape_2]
+    return shapes

tests/reg_tests/test_DVGeometry.py

@@ -1401,51 +1401,7 @@ def test_shape_functions(self, train=False, refDeriv=False):
             # coef indices, which is what we need to set the shape
             lidx = DVGeo.getLocalIndex(0)
 
-            # create the dictionaries
-            shape_1 = {}
-            shape_2 = {}
-
-            k_center = 2
-            i_center = 1
-            n_chord = lidx.shape[0]
-
-            for kk in [-1, 0, 1]:
-                if kk == 0:
-                    k_weight = 1.0
-                else:
-                    k_weight = 0.5
-
-                for ii in range(n_chord):
-                    # compute the chord weight. we want the shape to peak at i_center
-                    if ii == i_center:
-                        i_weight = 1.0
-                    elif ii < i_center:
-                        # we are ahead of the center point
-                        i_weight = ii / i_center
-                    else:
-                        # we are behind the center point
-                        i_weight = (n_chord - ii - 1) / (n_chord - i_center - 1)
-
-                    # get the direction vectors with unit length
-                    dir_up = np.array([0.0, 1.0, 0.0])
-                    # dir down can also be defined as an upwards pointing vector. Then, the DV itself
-                    # getting a negative value means the surface would move down etc. For now, we define
-                    # the vector as its pointing down, so a positive DV value moves the surface down.
-                    dir_down = np.array([0.0, -1.0, 0.0])
-
-                    # scale them by the i and k weights
-                    dir_up *= k_weight * i_weight
-                    dir_down *= k_weight * i_weight
-
-                    # get this point's global index and add to the dictionary with the direction vector.
-                    gidx_up = lidx[ii, 1, kk + k_center]
-                    gidx_down = lidx[ii, 0, kk + k_center]
-
-                    shape_1[gidx_up] = dir_up
-                    # the lower face is perturbed with a separate dictionary
-                    shape_2[gidx_down] = dir_down
-
-            shapes = [shape_1, shape_2]
+            shapes = commonUtils.getShapeFunc(lidx)
             DVGeo.addShapeFunctionDV("shape_func", shapes)
 
             # test derivatives