From 1ab1f6bf317b7000d1dc0ee4fde0be9d46447ba0 Mon Sep 17 00:00:00 2001 From: Anil Yildirim Date: Wed, 29 Mar 2023 09:57:41 +0200 Subject: [PATCH 1/8] first implementation of custom subset of active child FFDs for family groups. derivatives not added yet --- adflow/pyADflow.py | 97 +++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 88 insertions(+), 9 deletions(-) diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 24635c330..57d133e1e 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -3011,8 +3011,6 @@ def setSurfaceCoordinates(self, coordinates, groupName=None): def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): """Set the supplied aeroProblem to be used in ADflow""" - ptSetName = "adflow_%s_coords" % aeroProblem.name - newAP = False # Tell the user if we are switching aeroProblems if self.curAP != aeroProblem: @@ -3026,9 +3024,24 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): aeroProblem.adflowData except AttributeError: aeroProblem.adflowData = adflowFlowCase() - aeroProblem.ptSetName = ptSetName aeroProblem.surfMesh = self.getSurfaceCoordinates(self.designFamilyGroup) + # dictionary that holds the ptsetname for each family + ptSetNames = {} + activeChildrenFam = self.getOption("activechilddvgeofamilies") + if activeChildrenFam is None: + # we dont have a custom child dvgeo mapping. the surface family will be only + # designFamilyGroup and we will have a single pointset + ptSetName = f"adflow_{self.designFamilyGroup}_{aeroProblem.name}_coords" + ptSetNames[self.designFamilyGroup] = ptSetName + else: + # we have a custom surface family to child dvgeo mapping. + for familyName in activeChildrenFam.keys(): + ptSetName = f"adflow_{familyName}_{aeroProblem.name}_coords" + ptSetNames[familyName] = ptSetName + + aeroProblem.ptSetNames = ptSetNames + if self.curAP is not None: # If we have already solved something and are now # switching, save what we need: @@ -3056,14 +3069,75 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): # Now check if we have an DVGeo object to deal with: if self.DVGeo is not None: - # DVGeo appeared and we have not embedded points! - if ptSetName not in self.DVGeo.points: - coords0 = self.mapVector(self.coords0, self.allFamilies, self.designFamilyGroup, includeZipper=False) - self.DVGeo.addPointSet(coords0, ptSetName, **self.pointSetKwargs) + # we have a DVGeo added. check if we have already added the points for this AP + activeChildrenFam = self.getOption("activechilddvgeofamilies") + if activeChildrenFam is None: + # we have a single pointset + ptSetName = aeroProblem.ptSetNames[self.designFamilyGroup] + + if ptSetName not in self.DVGeo.points: + coords0 = self.mapVector(self.coords0, self.allFamilies, self.designFamilyGroup, includeZipper=False) + self.DVGeo.addPointSet(coords0, ptSetName, **self.pointSetKwargs) + else: + # we have custom pointsets + for family in activeChildrenFam.keys(): + ptSetName = aeroProblem.ptSetNames[family] + + if ptSetName not in self.DVGeo.points: + # coords0 is now the subset of this family + coords0 = self.mapVector(self.coords0, self.allFamilies, family, includeZipper=False) + # this pointset is added with a custom activeChildrenFamily + self.DVGeo.addPointSet( + coords0, + ptSetName, + activeChildren=activeChildrenFam[family], + **self.pointSetKwargs + ) # Check if our point-set is up to date: - if not self.DVGeo.pointSetUpToDate(ptSetName) or aeroProblem.adflowData.disp is not None: - coords = self.DVGeo.update(ptSetName, config=aeroProblem.name) + updateSurface = False + + # check for disp first. doing this check here is not the most efficient, + # but it results in simpler code + if aeroProblem.adflowData.disp is not None: + updateSurface = True + + if activeChildrenFam is None: + # we have a single pointset + ptSetName = aeroProblem.ptSetNames[self.designFamilyGroup] + + if not self.DVGeo.pointSetUpToDate(ptSetName): + updateSurface = True + + coords = self.DVGeo.update(ptSetName, config=aeroProblem.name) + + else: + # we have custom pointsets + # first figure out if we want to update the pointset; check each family we are tracking + for family in activeChildrenFam.keys(): + # return immediately if we are flagged for a surfafce update + if updateSurface: + break + + # check if any of the pointsets are out of date + ptSetName = aeroProblem.ptSetNames[family] + if not self.DVGeo.pointSetUpToDate(ptSetName): + updateSurface = True + + # if we have the updateSurface flag True, compute the coords array + if updateSurface: + # get the current design surface family. we will overwrite this as we go through the families + coords = self.mapVector(self.coords0, self.allFamilies, self.designFamilyGroup, includeZipper=False) + + for family in activeChildrenFam.keys(): + ptSetName = aeroProblem.ptSetNames[family] + familyCoords = self.DVGeo.update(ptSetName, config=aeroProblem.name) + # map this to the coords vector + self.mapVector(familyCoords, family, self.designFamilyGroup, vec2=coords, includeZipper=False) + + if updateSurface: + # the coords array is computed above. if we have the update surface flag enabled, + # run the surface update # Potentially add a fixed set of displacements to it. if aeroProblem.adflowData.disp is not None: @@ -4239,6 +4313,7 @@ def computeJacobianVectorProductFwd( # For the geometric xDvDot perturbation we accumulate into the # already existing (and possibly nonzero) xsdot and xvdot if xDvDot is not None or xSDot is not None: + # TODO fix if xDvDot is not None and self.DVGeo is not None: xsdot += self.DVGeo.totalSensitivityProd(xDvDot, self.curAP.ptSetName, config=self.curAP.name).reshape( xsdot.shape @@ -4560,6 +4635,7 @@ def computeJacobianVectorProductBwd( if self.mesh is not None: # Include geometric # derivatives if mesh is # present + # TODO fix if self.DVGeo is not None and self.DVGeo.getNDV() > 0: xdvbar.update( self.DVGeo.totalSensitivity(xsbar, self.curAP.ptSetName, self.comm, config=self.curAP.name) @@ -5198,6 +5274,7 @@ def _getDefaultOptions(): "meshSurfaceFamily": [(str, type(None)), None], "designSurfaceFamily": [(str, type(None)), None], "closedSurfaceFamilies": [(list, type(None)), None], + "activeChildDVGeoFamilies": [(dict, type(None)), None], # Output Parameters "storeRindLayer": [bool, True], "outputDirectory": [str, "./"], @@ -5478,6 +5555,7 @@ def _getImmutableOptions(self): "partitiononly", "meshsurfacefamily", "designsurfacefamily", + "activechilddvgeofamilies", "closedsurfacefamilies", "zippersurfacefamily", "cutcallback", @@ -5859,6 +5937,7 @@ def _getSpecialOptionLists(self): "liftindex", "meshsurfacefamily", "designsurfacefamily", + "activechilddvgeofamilies", "closedsurfacefamilies", "zippersurfacefamily", "outputsurfacefamily", From 13b30e1628384043d722d8c37219e11f634e46da Mon Sep 17 00:00:00 2001 From: Anil Yildirim Date: Wed, 5 Apr 2023 00:13:51 +0300 Subject: [PATCH 2/8] added reverse mode sensitivities w/ active children check --- adflow/pyADflow.py | 38 ++++++++++++++++++++++++++++++++------ 1 file changed, 32 insertions(+), 6 deletions(-) diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 192686418..d321e3d3b 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -4658,13 +4658,39 @@ def computeJacobianVectorProductBwd( if xDvDeriv: xdvbar = {} if self.mesh is not None: # Include geometric - # derivatives if mesh is - # present - # TODO fix + # derivatives if mesh is present if self.DVGeo is not None and self.DVGeo.getNDV() > 0: - xdvbar.update( - self.DVGeo.totalSensitivity(xsbar, self.curAP.ptSetName, self.comm, config=self.curAP.name) - ) + # we have a DVGeo added. check if we have already added the points for this AP + activeChildrenFam = self.getOption("activechilddvgeofamilies") + + if activeChildrenFam is None: + # we have a single pointset + ptSetName = self.curAP.ptSetNames[self.designFamilyGroup] + + xdvbar.update( + self.DVGeo.totalSensitivity(xsbar, ptSetName, self.comm, config=self.curAP.name) + ) + + else: + # we have custom pointsets + # start with an empty dict, we set the entries in the first pass, + # and add in the following passes + dvgeoSens = {} + for family in activeChildrenFam.keys(): + ptSetName = self.curAP.ptSetNames[family] + + # get the mapped xsbar + xsbarFamily = self.mapVector(xsbar, self.allFamilies, family, includeZipper=False) + + familySens = self.DVGeo.totalSensitivity(xsbarFamily, ptSetName, self.comm, config=self.curAP.name) + + for key, val in familySens.items(): + # not the best way to do this but should work... + if key in dvgeoSens: + dvgeoSens[key] += val + else: + dvgeoSens[key] = val + else: if self.comm.rank == 0: ADFLOWWarning( From 9da19f2b3f5759ee4201050c321d400629d88022 Mon Sep 17 00:00:00 2001 From: Anil Yildirim Date: Wed, 5 Apr 2023 02:50:33 +0300 Subject: [PATCH 3/8] bug fix --- adflow/pyADflow.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index d321e3d3b..38383df1a 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -4680,7 +4680,7 @@ def computeJacobianVectorProductBwd( ptSetName = self.curAP.ptSetNames[family] # get the mapped xsbar - xsbarFamily = self.mapVector(xsbar, self.allFamilies, family, includeZipper=False) + xsbarFamily = self.mapVector(xsbar, self.designFamilyGroup, family, includeZipper=False) familySens = self.DVGeo.totalSensitivity(xsbarFamily, ptSetName, self.comm, config=self.curAP.name) From be7d63f0a46e6ad3b75d9fc5139883821174c576 Mon Sep 17 00:00:00 2001 From: Anil Yildirim Date: Tue, 11 Apr 2023 11:11:49 +0200 Subject: [PATCH 4/8] fixed the derivative bug --- adflow/pyADflow.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 38383df1a..499a74448 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -4675,7 +4675,6 @@ def computeJacobianVectorProductBwd( # we have custom pointsets # start with an empty dict, we set the entries in the first pass, # and add in the following passes - dvgeoSens = {} for family in activeChildrenFam.keys(): ptSetName = self.curAP.ptSetNames[family] @@ -4686,10 +4685,10 @@ def computeJacobianVectorProductBwd( for key, val in familySens.items(): # not the best way to do this but should work... - if key in dvgeoSens: - dvgeoSens[key] += val + if key in xdvbar: + xdvbar[key] += val else: - dvgeoSens[key] = val + xdvbar[key] = val else: if self.comm.rank == 0: From 5beaff34ec845f16305356d90af218ae8c062e56 Mon Sep 17 00:00:00 2001 From: Anil Yildirim Date: Thu, 20 Jul 2023 13:48:41 -0400 Subject: [PATCH 5/8] modified the implementation to be more general --- adflow/pyADflow.py | 31 ++++++++++++------------------- 1 file changed, 12 insertions(+), 19 deletions(-) diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index bc072a343..39b95bd7d 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -3302,15 +3302,14 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): # dictionary that holds the ptsetname for each family ptSetNames = {} - activeChildrenFam = self.getOption("activechilddvgeofamilies") - if activeChildrenFam is None: + if self.customPointSetFamilies is None: # we dont have a custom child dvgeo mapping. the surface family will be only # designFamilyGroup and we will have a single pointset ptSetName = f"adflow_{self.designFamilyGroup}_{aeroProblem.name}_coords" ptSetNames[self.designFamilyGroup] = ptSetName else: # we have a custom surface family to child dvgeo mapping. - for familyName in activeChildrenFam.keys(): + for familyName in self.customPointSetFamilies.keys(): ptSetName = f"adflow_{familyName}_{aeroProblem.name}_coords" ptSetNames[familyName] = ptSetName @@ -3344,8 +3343,7 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): # Now check if we have an DVGeo object to deal with: if self.DVGeo is not None: # we have a DVGeo added. check if we have already added the points for this AP - activeChildrenFam = self.getOption("activechilddvgeofamilies") - if activeChildrenFam is None: + if self.customPointSetFamilies is None: # we have a single pointset ptSetName = aeroProblem.ptSetNames[self.designFamilyGroup] @@ -3354,18 +3352,18 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): self.DVGeo.addPointSet(coords0, ptSetName, **self.pointSetKwargs) else: # we have custom pointsets - for family in activeChildrenFam.keys(): + for family, familyKwargs in self.customPointSetFamilies.items(): ptSetName = aeroProblem.ptSetNames[family] if ptSetName not in self.DVGeo.points: # coords0 is now the subset of this family coords0 = self.mapVector(self.coords0, self.allFamilies, family, includeZipper=False) - # this pointset is added with a custom activeChildrenFamily + # this pointset is added with a custom familyKwargs self.DVGeo.addPointSet( coords0, ptSetName, - activeChildren=activeChildrenFam[family], - **self.pointSetKwargs + **self.pointSetKwargs, + **familyKwargs, ) # Check if our point-set is up to date: @@ -3376,7 +3374,7 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): if aeroProblem.adflowData.disp is not None: updateSurface = True - if activeChildrenFam is None: + if self.customPointSetFamilies is None: # we have a single pointset ptSetName = aeroProblem.ptSetNames[self.designFamilyGroup] @@ -3388,7 +3386,7 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): else: # we have custom pointsets # first figure out if we want to update the pointset; check each family we are tracking - for family in activeChildrenFam.keys(): + for family in self.customPointSetFamilies.keys(): # return immediately if we are flagged for a surfafce update if updateSurface: break @@ -3403,7 +3401,7 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): # get the current design surface family. we will overwrite this as we go through the families coords = self.mapVector(self.coords0, self.allFamilies, self.designFamilyGroup, includeZipper=False) - for family in activeChildrenFam.keys(): + for family in self.customPointSetFamilies.keys(): ptSetName = aeroProblem.ptSetNames[family] familyCoords = self.DVGeo.update(ptSetName, config=aeroProblem.name) # map this to the coords vector @@ -4918,9 +4916,7 @@ def computeJacobianVectorProductBwd( # derivatives if mesh is present if self.DVGeo is not None and self.DVGeo.getNDV() > 0: # we have a DVGeo added. check if we have already added the points for this AP - activeChildrenFam = self.getOption("activechilddvgeofamilies") - - if activeChildrenFam is None: + if self.customPointSetFamilies is None: # we have a single pointset ptSetName = self.curAP.ptSetNames[self.designFamilyGroup] @@ -4932,7 +4928,7 @@ def computeJacobianVectorProductBwd( # we have custom pointsets # start with an empty dict, we set the entries in the first pass, # and add in the following passes - for family in activeChildrenFam.keys(): + for family in self.customPointSetFamilies.keys(): ptSetName = self.curAP.ptSetNames[family] # get the mapped xsbar @@ -5586,7 +5582,6 @@ def _getDefaultOptions(): "meshSurfaceFamily": [(str, type(None)), None], "designSurfaceFamily": [(str, type(None)), None], "closedSurfaceFamilies": [(list, type(None)), None], - "activeChildDVGeoFamilies": [(dict, type(None)), None], # Output Parameters "storeRindLayer": [bool, True], "outputDirectory": [str, "./"], @@ -5874,7 +5869,6 @@ def _getImmutableOptions(self): "partitiononly", "meshsurfacefamily", "designsurfacefamily", - "activechilddvgeofamilies", "closedsurfacefamilies", "zippersurfacefamily", "cutcallback", @@ -6266,7 +6260,6 @@ def _getSpecialOptionLists(self): "liftindex", "meshsurfacefamily", "designsurfacefamily", - "activechilddvgeofamilies", "closedsurfacefamilies", "zippersurfacefamily", "outputsurfacefamily", From 9f047f42fd9851ec7264ef32037fac40ab601100 Mon Sep 17 00:00:00 2001 From: Anil Yildirim Date: Sun, 26 Nov 2023 12:58:33 +0100 Subject: [PATCH 6/8] add the missing jacvec modifications --- adflow/pyADflow.py | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 4f8ccf161..39bbbface 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -3358,7 +3358,9 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): ptSetName = aeroProblem.ptSetNames[self.designFamilyGroup] if ptSetName not in self.DVGeo.points: - coords0 = self.mapVector(self.coords0, self.allFamilies, self.designFamilyGroup, includeZipper=False) + coords0 = self.mapVector( + self.coords0, self.allFamilies, self.designFamilyGroup, includeZipper=False + ) self.DVGeo.addPointSet(coords0, ptSetName, **self.pointSetKwargs) else: # we have custom pointsets @@ -4610,11 +4612,23 @@ def computeJacobianVectorProductFwd( # For the geometric xDvDot perturbation we accumulate into the # already existing (and possibly nonzero) xsdot and xvdot if xDvDot is not None or xSDot is not None: - # TODO fix if xDvDot is not None and self.DVGeo is not None: - xsdot += self.DVGeo.totalSensitivityProd(xDvDot, self.curAP.ptSetName, config=self.curAP.name).reshape( - xsdot.shape - ) + if self.customPointSetFamilies is None: + # no custom pointset families, just process the entire dvdot and add to xsdot + xsdot += self.DVGeo.totalSensitivityProd( + xDvDot, self.curAP.ptSetName, config=self.curAP.name + ).reshape(xsdot.shape) + else: + # custom pointsets. accumulate local xsdots in the complete vector + for family in self.customPointSetFamilies.keys(): + # process this family's xsdot + ptSetName = self.curAP.ptSetNames[family] + xsdot_family = self.DVGeo.totalSensitivityProd( + xDvDot, ptSetName, self.comm, config=self.curAP.name + ) + # map it to an empty surface vector and accumulate + xsdot += self.mapVector(xsdot_family, family, self.designFamilyGroup) + if self.mesh is not None: xsdot = self.mapVector(xsdot, self.meshFamilyGroup, self.designFamilyGroup, includeZipper=False) xvdot += self.mesh.warpDerivFwd(xsdot) @@ -4951,7 +4965,9 @@ def computeJacobianVectorProductBwd( # get the mapped xsbar xsbarFamily = self.mapVector(xsbar, self.designFamilyGroup, family, includeZipper=False) - familySens = self.DVGeo.totalSensitivity(xsbarFamily, ptSetName, self.comm, config=self.curAP.name) + familySens = self.DVGeo.totalSensitivity( + xsbarFamily, ptSetName, self.comm, config=self.curAP.name + ) for key, val in familySens.items(): # not the best way to do this but should work... From bb2a6cd067973cf3c9caaa2bf74e6efb57fab881 Mon Sep 17 00:00:00 2001 From: Anil Yildirim Date: Wed, 29 Nov 2023 09:55:19 +0100 Subject: [PATCH 7/8] fix fwd jacvec product --- adflow/pyADflow.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 39bbbface..3f8bbdb79 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -4616,7 +4616,7 @@ def computeJacobianVectorProductFwd( if self.customPointSetFamilies is None: # no custom pointset families, just process the entire dvdot and add to xsdot xsdot += self.DVGeo.totalSensitivityProd( - xDvDot, self.curAP.ptSetName, config=self.curAP.name + xDvDot, self.curAP.ptSetNames[self.designFamilyGroup], config=self.curAP.name ).reshape(xsdot.shape) else: # custom pointsets. accumulate local xsdots in the complete vector From fc136955ed5b455121898f93200037dc30dcc86a Mon Sep 17 00:00:00 2001 From: Anil Yildirim Date: Thu, 16 Jan 2025 15:51:14 -0800 Subject: [PATCH 8/8] add custom ptset kwargs for the blanking surfaces when they are used in the full update mode --- adflow/pyADflow.py | 3 ++- doc/options.yaml | 7 +++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index af5b859d9..2d1aabddd 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -3346,6 +3346,7 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): # used in the explicit hole cutting multiple times in different # configurations. surfFile = self.blankingSurfDict[surf]["surfFile"] + surfPointSetKwargs = self.blankingSurfDict[surf]["pointSetKwargs"] surfPtSetName = f"points_{surfFile}" if surfPtSetName not in self.DVGeo.points: # we need to add the pointset to dvgeo. do it in parallel @@ -3367,7 +3368,7 @@ def setAeroProblem(self, aeroProblem, releaseAdjointMemory=True): # we already communicated the points when loading the file, # so just add them to dvgeo now procPts = surfPts[disp[self.comm.rank] : disp[self.comm.rank + 1]] - self.DVGeo.addPointSet(procPts, surfPtSetName, **self.pointSetKwargs) + self.DVGeo.addPointSet(procPts, surfPtSetName, **self.pointSetKwargs, **surfPointSetKwargs) # Check if our point-set is up to date: updateSurface = False diff --git a/doc/options.yaml b/doc/options.yaml index 18c34416d..137fc862a 100644 --- a/doc/options.yaml +++ b/doc/options.yaml @@ -629,6 +629,13 @@ explicitSurfaceCallback: # We use this approach so that the users can have complete # control over the surface mesh nodes after they are loaded. "coordXfer": coordXfer, + + # Optional entry: pointSetKwargs, + # This is a dictionary that contains custom kwargs that will be + # used with the addPointSet call when this surface is added to + # DVGeo. This is done when the solver does full overset updates + # between iterations, where we also track the blanking surfaces + # with DVGeo. }, # Similar entries for the fuselage.