Bugfix for rotating frame sensitivities + add Tapenade checks (#121)

* first changes made * Updating the fork before jess-install * Only output from Rank0 * bug fix for COMPLEX makefile * Merged mdolab/adflow into default * Fwd() routine using master_d close * Fwd() routine using master_d close * Fwd() routine using master_d close * evalFunctionsSensFwd() now outputs in desires O.dict fmt * dMxdMach works. Still problems with FFD DVs * This is the working settings for dMx/dMach. The Fwd()-call n pyADflow.py is still hardcoded. This needs to be improved. * inviscidCentralFlux(), add to Makefile_tapenade * The dMx/dMach still works - but dMxdFFDpts have not been fixed all together. However, it seems from the older commits that these files we add in this commit are crucial. There are a few additions that should be located, but then it should (hopefully) work. Most likely, the error that remains to be fixed now is that some functions already in use are not differentiated correctly (just as in the previous case with inviscidcentralflux()). Hopefully a fix is coming up shortly. * Forward mode works for FFD and for Children. Bits and pieces have been hardcoded. Note; this commit will not work without the proper DVGeometry.py file from PyGeo... * Merged mdolab/adflow into default * normalvel_d seems to work. Still need slip() and grid(). * Now both normalvel() and gridvel() seems to work. We (hopefully) only need slipvel() and then we need to move them all to adjointextra_d.f90 * Now also slipvel() seems to work. Will now move them to adjointExtra.F90 * FWD() seems to give ok gradient, BWD is still off * cleaning * forward works, master_b close * Forward still good, Reverse still off. All _b files should however have been added. Weird * still trying to fix master_b. Master_d is still fine.... * master_d still works when removing changes made to bcroutines. Now the possible files containing the error for master_b has been narrowed down. * Master_b WORKS! for wind turbines * notes on how to do dot-tests added * added rot frame test - manual ref, to be retrained * updated ref for rot test * updated rot test input folder * finished pyADflow linting - new functs * missed linting * gridvelocitiesfine_block fixed * modified the make file * reset a deleted function that needds to be aded * fix function slipVelocitiesFineLevel_block * normalvelocities_block fixed * minor edits to remove inconstency * remove inconsistency of bcroutines * fix the inconsistency of residual * fix inconsistency of flowutils * more * add missing nes * fix indexing * comments cleanup in fortran code * added back addabsslice * moved Mads dev tips to documentation * removed unused functions from pyAdflow * turning on the complex rotating test * manual fix to the ref file - float to array * fixed aeroproblem definition * should have checked the whole file * added sens wrt mach * added mach to dv list * removed unused fast_b files and updated gitignore * renamed variable t to time * removed gammaInf from Makefile_tapenade * retrained tests (span func change) and added metadata * adding missing function geometry parameters for AD * add missing function derivative for initres and timestep_block, now FAD matches FD for pRpq for a rotating problem * pRpx achieves 1e-5 accuracy * fix missing reverse aded function * clean not used function * fixed API for totalSensitivityProd * re-trained rotating adjoint test * cleaned up debug lines * updated test - will it work? * tol adjustment - it's happening * added tapenade tests * added the shell script for tapenade tests * removed non-diferentiated TS routine * addressing Sabet's comments * cleaned up blockette.f90 * bump patch version * cleaned up BCroutines * removed low-speed preconditioner edits * fixed comments in outputMod * removing outdated documentation * retrained tests * forgot to update one value * replaced repeated code with single subroutine * differentiated new subroutine * removed commented code * removed unnecessary options from test + retrained metadata * updated options capitalization for all adjoint tests * attempting tigher tolerances for complex tests * Revert "attempting tigher tolerances for complex tests" This reverts commit c2e554f. * fixed autoEdit file to not convert unused files to fast_b * removed unused variable * typo * moved initres_block * moved initres_block again to fix diff * fixed spacing * fixed spacing forreal * moved initres_block back up because it is actually differentiated now * using MachCoeff for Cp normalization - should be equivalent to previous two approaches, I checked the cp curves again! * removed useless square root * Makefile tapenade: cleaned fluxes&viscousFlux * Makefile tapenade: cleaned flowUtils%computeSpeedOfSoundSquared * Makefile tapenade: cleaned various solverutils routines * Makefile tapenade: fixed solverutils routines * Makefile tapenade: removed timeref from solverutils routines * Makefile tapenade: removed pinf and rhoinf from solverutils routines * cleaned up comments in Cp calculation - outputMod.F90 * version bump * removed extra call to timeStep_block_d Co-authored-by: Sicheng He <[email protected]> Co-authored-by: Sicheng He <[email protected]>
mdolab · May 17, 2022 · bb0b139 · bb0b139
1 parent b3878fb
commit bb0b139
Show file tree

Hide file tree

Showing 37 changed files with 11,679 additions and 3,415 deletions.
diff --git a/.github/azure-pipelines.yaml b/.github/azure-pipelines.yaml
@@ -19,3 +19,4 @@ extends:
     GCC_CONFIG: config/defaults/config.LINUX_GFORTRAN.mk
     INTEL_CONFIG: config/defaults/config.LINUX_INTEL.mk
     COVERAGE: true
+    TAPENADE: true
diff --git a/.github/build_tapenade.sh b/.github/build_tapenade.sh
@@ -0,0 +1,5 @@
+#!/bin/bash
+set -e
+
+cd src/adjoint
+make -f Makefile_tapenade ad_forward ad_reverse ad_reverse_fast
diff --git a/adflow/__init__.py b/adflow/__init__.py
@@ -1,4 +1,4 @@
-__version__ = "2.7.2"
+__version__ = "2.7.3"
 
 from mpi4py import MPI
 

diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py
@@ -567,7 +567,7 @@ def addSlices(self, direction, positions, sliceType="relative", groupName=None):
 
         for i in range(len(positions)):
             # It is important to ensure each slice get a unique
-            # name...so we will number sequentially from pythhon
+            # name...so we will number sequentially from python
             j = self.nSlice + i + 1
             if sliceType == "relative":
                 sliceName = "Slice_%4.4d %s Para Init %s=%7.3f" % (j, groupTag, direction, positions[i])
@@ -3975,9 +3975,9 @@ def computeJacobianVectorProductFwd(
         # already existing (and possibly nonzero) xsdot and xvdot
         if xDvDot is not None or xSDot is not None:
             if xDvDot is not None and self.DVGeo is not None:
-                xsdot += self.DVGeo.totalSensitivityProd(
-                    xDvDot, self.curAP.ptSetName, self.comm, config=self.curAP.name
-                ).reshape(xsdot.shape)
+                xsdot += self.DVGeo.totalSensitivityProd(xDvDot, self.curAP.ptSetName, config=self.curAP.name).reshape(
+                    xsdot.shape
+                )
             if self.mesh is not None:
                 xsdot = self.mapVector(xsdot, self.meshFamilyGroup, self.designFamilyGroup, includeZipper=False)
                 xvdot += self.mesh.warpDerivFwd(xsdot)

diff --git a/doc/adjoint.rst b/doc/adjoint.rst
@@ -95,5 +95,4 @@ The user can give a custom function that uses the intrinsic ADflow functions. AD
 
 The boundary condition definition and sensitivities are stored under ``AeroProblem`` to facilitate the use of different boundary conditions for multipoint cases.
 
-The reverse master (``master_b``) routine is wrapped by a thin wrapper which is the ``computeMatrixFreeProductBwd`` subroutine defined in ``adjointAPI.F90``. This wrapper allocates memory etc. before the derivatives are calculated. 
-
+The reverse master (``master_b``) routine is wrapped by a thin wrapper which is the ``computeMatrixFreeProductBwd`` subroutine defined in ``adjointAPI.F90``. This wrapper allocates memory etc. before the derivatives are calculated.
diff --git a/src/NKSolver/blockette.F90 b/src/NKSolver/blockette.F90
@@ -226,6 +226,7 @@ subroutine blocketteRes(useDissApprox, useViscApprox, useUpdateIntermed, useFlow
              call applyAllTurbBCthisblock(.True.)
           end if
           call applyAllBC_block(.True.)
+
        end do
     end do
 

diff --git a/src/adjoint/Makefile_tapenade b/src/adjoint/Makefile_tapenade
@@ -103,7 +103,7 @@ adjointExtra%volume_block(x) > \
 \
 adjointExtra%metric_block(x) > \
                          (x, si, sj, sk) \
-wallDistance%updateWallDistancesQuickly(x, xSurf) > \
+wallDistance%updateWallDistancesQuickly(x, xSurf, d2wall) > \
 	                               (x, xSurf, d2wall) \
 \
 adjointExtra%boundaryNormals(si, sj, sk) > \
@@ -132,14 +132,14 @@ BCRoutines%bcSymmPolar1stHalo(xx, ww1, ww2, pp1, pp2, rlv1, rlv2, rev1, rev2) >
 BCRoutines%bcSymmPolar2ndHalo(xx, ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3) > \
                              (xx, ww0, ww3, pp0, pp3, rlv0, rlv3, rev0, rev3) \
 \
-BCRoutines%bcNSWallAdiabatic(ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2) > \
-                            (ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2) \
+BCRoutines%bcNSWallAdiabatic(ww0, ww1, ww2, pp0, pp1, pp2, pp3, rlv0, rlv1, rlv2, rev0, rev1, rev2, bcData%uSlip) > \
+                            (ww0, ww1, ww2, pp0, pp1, pp2, pp3, rlv0, rlv1, rlv2, rev0, rev1, rev2, bcData%uSlip) \
 \
-BCRoutines%bcNSWallIsothermal(ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, rGas, BCData%TNS_Wall) > \
-                             (ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, rGas, BCData%TNS_Wall) \
+BCRoutines%bcNSWallIsothermal(ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, rGas, BCData%TNS_Wall, bcData%uSlip) > \
+                             (ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, rGas, BCData%TNS_Wall, bcData%uSlip) \
 \
-BCRoutines%bcFarField(ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, pInfCorr, wInf, bcData%norm) > \
-                     (ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, pInfCorr, wInf, bcData%norm) \
+BCRoutines%bcFarField(ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, pInfCorr, wInf, bcData%norm, bcData%rface) > \
+                     (ww0, ww1, ww2, pp0, pp1, pp2, rlv0, rlv1, rlv2, rev0, rev1, rev2, pInfCorr, wInf, bcData%norm, bcData%rface) \
 \
 BCRoutines%bcEulerWall(ww0, ww1, ww2, pp0, pp1, pp2, pp3, rlv0, rlv1, rlv2, rev0, rev1, rev2, bcData%norm, ssi, ssj, ssk, ss) > \
                       (ww0, ww1, ww2, pp0, pp1, pp2, pp3, rlv0, rlv1, rlv2, rev0, rev1, rev2, bcData%norm, ssi, ssj, ssk, ss) \
@@ -156,14 +156,14 @@ turbBCRoutines%applyAllTurbBCThisBlock(rev, w, bvtj1, bvtj2, bvtk1, bvtk2, bvti1
 turbBCRoutines%bcTurbTreatment(w, rlv, d2wall, winf) > \
                               (w, rlv, d2wall, winf, bvtj1, bvtj2, bvtk1, bvtk2, bvti1, bvti2) \
 \
-solverUtils%timeStep_block(w, pInfCorr, rhoInf, si, sj, sk, p) > \
-              (w, pInfCorr, rhoInf, si, sj, sk, p, radi, radj, radk) \
+solverUtils%timeStep_block(w, pInfCorr, rhoInf, si, sj, sk, sFaceI, sFaceJ, sFaceK, p, radi, radj, radk, dtl, rlv, rev, vol) > \
+              (w, pInfCorr, rhoInf, si, sj, sk, sFaceI, sFaceJ, sFaceK, p, radi, radj, radk, dtl, rlv, rev, vol) \
 \
 sa%saSource(w, rlv, vol, si, sj, sk, timeRef, d2wall) > \
            (w, rlv, vol, si, sj, sk, timeREf, scratch) \
 \
-turbutils%turbAdvection(w, vol, si, sj, sk, scratch) > \
-                       (w, vol, si, sj, sk, scratch) \
+turbutils%turbAdvection(w, vol, si, sj, sk, scratch, sfacei, sfacej, sfacek) > \
+                       (w, vol, si, sj, sk, scratch, sfacei, sfacej, sfacek) \
 \
 sa%saViscous(w, vol, si, sj, sk, rlv, scratch) > \
             (w, vol, si, sj, sk, rlv, scratch) \
@@ -174,34 +174,37 @@ sa%saResScale(scratch, dw) > \
 residuals%sourceTerms_block(w, pref, uref, plocal, dw, actuatorRegions%force, actuatorRegions%heat) > \
                            (w, pref, uref, plocal, dw, actuatorRegions%force, actuatorRegions%heat) \
 \
-fluxes%inviscidCentralFlux(w, timeRef, vol, si, sj, sk, p, dw) > \
-                          (w, timeRef, vol, si, sj, sk, p, dw) \
+residuals%initres_block(dw, fw, flowDoms%w, flowDoms%vol, dscalar, dvector) > \
+                       (dw, fw, flowDoms%w, flowDoms%vol, dscalar, dvector) \
+\
+fluxes%inviscidCentralFlux(w, timeRef, vol, si, sj, sk, p, dw, sfacei, sfacej, sfacek) > \
+                          (w, timeRef, vol, si, sj, sk, p, dw, sfacei, sfacej, sfacek) \
 \
 fluxes%inviscidDissFluxScalar(w, pInfCorr, rhoInf, p, radi, radj, radk, fw) > \
-                             (w, pInfCOrr, rhoInf, p, fw) \
+                             (w, pInfCOrr, rhoInf, p, radi, radj, radk, fw) \
 \
-fluxes%inviscidDissFluxMatrix(w, pInfCorr, si, sj, sk, p, fw) > \
-                             (w, pInfCorr, si, sj, sk, p, fw) \
+fluxes%inviscidDissFluxMatrix(w, pInfCorr, si, sj, sk, p, fw, sfacei, sfacej, sfacek) > \
+                             (w, pInfCorr, si, sj, sk, p, fw, sfacei, sfacej, sfacek) \
 \
-fluxes%inviscidDissFluxScalarApprox(w, p, radi, radj, radk, fw) > \
+fluxes%inviscidDissFluxScalarApprox(w, p, radi, radj, radk, fw, rhoinf, pinfcorr) > \
                                    (w, p, fw) \
 \
-fluxes%inviscidDissFluxMatrixApprox(w, p, fw) > \
-                                   (w, p, fw) \
+fluxes%inviscidDissFluxMatrixApprox(w, p, fw, pinfcorr, sfacei, sfacej, sfacek, si, sj, sk) > \
+                                   (w, p, fw, pinfcorr, sfacei, sfacej, sfacek, si, sj, sk) \
 \
-fluxes%inviscidUpwindFlux(w, si, sj, sk, p, fw) > \
-                         (w, si, sj, sk, p, fw) \
+fluxes%inviscidUpwindFlux(w, si, sj, sk, p, fw, sfacei, sfacej, sfacek) > \
+                         (w, si, sj, sk, p, fw, sfacei, sfacej, sfacek) \
 \
 flowUtils%computeSpeedOfSoundSquared(w, p) > \
                                     (w, p, aa) \
 \
-flowUtils%allNodalGradients(w, vol, si, sj, sk, aa) > \
+flowUtils%allNodalGradients(w, vol, si, sj, sk, aa, ux, uy, uz, vx, vy, vz, wx, wy, wz, qx, qy, qz) > \
                            (w, vol, si, sj, sk, aa, ux, uy, uz, vx, vy, vz, wx, wy, wz, qx, qy, qz) \
 \
-fluxes%viscousFlux(w, x, si, sj, sk, rlv, rev, aa, ux, uy, uz, vx, vy, vz, wx, wy, wz, qx, qy, qz, fw) > \
-                  (w, x, si, sj, sk, fw, viscsubface%tau) \
+fluxes%viscousFlux(w, x, si, sj, sk, rlv, rev, aa, ux, uy, uz, vx, vy, vz, wx, wy, wz, qx, qy, qz, fw, viscsubface%tau, viscSubface%q) > \
+                  (w, x, si, sj, sk, fw, viscsubface%tau, viscSubface%q) \
 \
-fluxes%viscousFluxApprox(w, x, rlv, rev, aa, fw) > \
+fluxes%viscousFluxApprox(w, x, rlv, rev, aa, fw, si, sj, sk) > \
                         (w, x, fw) \
 \
 adjointExtra%sumDwandFw(fw, dw) > (dw) \
@@ -228,6 +231,17 @@ oversetUtilities%newtonUpdate(xCen, blk, frac) > \
 oversetUtilities%fracToWeights(frac) > \
 			      (weights) \
 \
+solverutils%slipVelocitiesFineLevel_block(pinf, timeref, rhoinf, veldirfreestream, machgrid, bcdata%uslip, flowDoms%x) > \
+										 (veldirfreestream, bcdata%uslip, flowDoms%x) \
+\
+solverutils%normalvelocities_block(sfacei, sfacej, sfacek, si, sj, sk, bcdata%rface) > \
+			      				  (sfacei, sfacej, sfacek, si, sj, sk, bcdata%rface) \
+\
+solverutils%gridVelocitiesFineLevel_block(flowDoms%x, si, sj, sk, s, sfacei, sfacej, sfacek, velDirFreestream, pinf, timeref, rhoinf, machgrid) > \
+                                         (flowDoms%x, si, sj, sk, s, sfacei, sfacej, sfacek, velDirFreestream) \
+\
+solverutils%cellFaceVelocities(xc, rotCenter, rotRate, velxGrid, velyGrid, velzGrid, derivRotationMatrix, sc) > \
+							  (xc, velxGrid, velyGrid, velzGrid, derivRotationMatrix, sc) \
 "
 
 #applyLowSpeedPreconditioner(w, P, dw) > (w, dw, P) \
@@ -243,6 +257,9 @@ stateOnlyRoutines = "\
 residuals%sourceTerms_block(w, dw) > \
                            (w, dw) \
 \
+residuals%initres_block(dw, fw, flowDoms%w) > \
+                       (dw, fw, flowDoms%w) \
+\
 fluxes%inviscidCentralFlux(w, p, dw) > \
                           (w, p, dw) \
 \

diff --git a/src/adjoint/adjointUtils.F90 b/src/adjoint/adjointUtils.F90
@@ -796,6 +796,7 @@ subroutine allocDerivativeValues(level)
                flowDomsd(nn, level, sps)%qz(il,jl,kl), &
                flowDomsd(nn, level, sps)%rlv(0:ib,0:jb,0:kb), &
                flowDomsd(nn, level, sps)%rev(0:ib,0:jb,0:kb), &
+               flowDomsd(nn, level, sps)%dtl(1:ie,1:je,1:ke), &
                flowDomsd(nn, level, sps)%radI(1:ie,1:je,1:ke), &
                flowDomsd(nn, level, sps)%radJ(1:ie,1:je,1:ke), &
                flowDomsd(nn, level, sps)%radK(1:ie,1:je,1:ke), &
@@ -942,6 +943,8 @@ subroutine zeroADSeeds(nn, level, sps)
     flowDomsd(nn, level, sps)%rlv = zero
     flowDomsd(nn, level, sps)%rev = zero
 
+    flowDomsd(nn, level, sps)%dtl = zero
+
     flowDomsd(nn, level, sps)%radI = zero
     flowDomsd(nn, level, sps)%radJ = zero
     flowDomsd(nn, level, sps)%radK = zero
@@ -2234,6 +2237,11 @@ subroutine setDiffSizes
     ISIZE2OFDrfgamma = jb + 1
     ISIZE3OFDrfgamma = kb + 1
 
+    ! dtl
+    ISIZE1OFDrfdtl = ie
+    ISIZE2OFDrfdtl = je
+    ISIZE3OFDrfdtl = ke
+
     ! radI
     ISIZE1OFDrfradI = ie
     ISIZE2OFDrfradI = je

diff --git a/src/adjoint/autoEdit/autoEditReverseFast.py b/src/adjoint/autoEdit/autoEditReverseFast.py
@@ -35,30 +35,31 @@
 print("Directory of output source files :", DIR_MOD)
 
 useful_modules = [
+    "bcpointers_b",
     "bcroutines_b",
-    "turbbcroutines_b",
-    "utils_b",
     "flowutils_b",
-    "walldistance_b",
-    "bcpointers_b",
+    "fluxes_b",
     "initializeflow_b",
-    "turbutils_b",
+    "residuals_b",
     "sa_b",
-    "fluxes_b",
     "solverutils_b",
-    "residuals_b",
     "surfaceintegrations_b",
+    "turbbcroutines_b",
+    "turbutils_b",
+    "utils_b",
+    "walldistance_b",
 ]
 
+FILE_IGNORE = ["adjointextra_b.f90", "bcdata_b.f90", "oversetutilities_b.f90", "zipperintegrations_b.f90"]
+
 for f in os.listdir(DIR_ORI):
-    if f.endswith(EXT):
+    if f not in FILE_IGNORE and f.endswith(EXT):
         # open original file in read mode
         file_object_ori = open(os.path.join(DIR_ORI, f), "r")
         print("\nParsing input file", file_object_ori.name)
 
         # read to whole file to string and reposition the pointer
         # at the first byte for future reading
-        all_src = file_object_ori.read()
         file_object_ori.seek(0)
 
         # First we want to dertmine if it is a 'useful' module or a
@@ -94,13 +95,6 @@
             # Just deal with lower case string
             line = line.lower()
 
-            # # Replace modules
-            # m = patt_modules.match(line)
-            # if m:
-            #     if not addedModule:
-            #         file_object_mod.write('  use myPushPopLib\n')
-            #         addedModule = True
-
             # Replace _cb on calls
             if "_cb" in line:
                 line = line.replace("_cb", "")