Time spectral with flexible grid

adflow/pyADflow.py

@@ -2594,9 +2594,9 @@ def getIblankCheckSum(self, fileName=None):
     #   i.e. an Aerostructural solver
     # =========================================================================
 
-    def getSurfaceCoordinates(self, groupName=None, includeZipper=True):
+    def getSurfaceCoordinates(self, groupName=None, includeZipper=True, TS=0):
         # This is an alias for getSurfacePoints
-        return self.getSurfacePoints(groupName, includeZipper)
+        return self.getSurfacePoints(groupName, includeZipper, TS)
 
     def getPointSetName(self, apName):
         """
@@ -2940,6 +2940,19 @@ def _setAeroProblemData(self, aeroProblem, firstCall=False):
         # 4. Periodic Parameters --- These are not checked/verified
         # and come directly from aeroProblem. Make sure you specify
         # them there properly!!
+        if self.getOption("useExternalDynamicMesh"):
+            # The inputs here are mainly used for TS stability problem.
+            # For general aeroelastic problem, we rely on more general settings.
+            # Currently, polynomial motion is not supported for aeroelastic problem.
+            AP.degreePol = 0
+            AP.coefPol = [0.0]
+            #  no odd number of instance allowed!
+            AP.degreeFourier = (self.adflow.inputtimespectral.ntimeintervalsspectral - 1) // 2
+            AP.cosCoefFourier = [0.0]
+            AP.sinCoefFourier = [0.0]
+        # if use time spectral, we need one of the following
+        # to be true to get a correct time period.
+        # the number of time instance is set directly.
         if self.getOption("alphaMode"):
             self.adflow.inputmotion.degreepolalpha = int(AP.degreePol)
             self.adflow.inputmotion.coefpolalpha = AP.coefPol
@@ -2983,6 +2996,9 @@ def _setAeroProblemData(self, aeroProblem, firstCall=False):
             self.adflow.inputmotion.degreefouryrot = AP.degreeFourier
             self.adflow.inputmotion.coscoeffouryrot = AP.cosCoefFourier
             self.adflow.inputmotion.sincoeffouryrot = AP.sinCoefFourier
+        elif self.getOption("useExternalDynamicMesh"):
+            # if it is an aeroelastic case
+            self.adflow.inputtimespectral.omegafourier = AP.omegaFourier
 
         # Set any possible BC Data coming out of the aeroProblem
         nameArray, dataArray, groupArray, groupNames, empty = self._getBCDataFromAeroProblem(AP)
@@ -3604,7 +3620,10 @@ def updateGeometryInfo(self, warpMesh=True):
                 self.adflow.killsignals.routinefailed = False
                 self.adflow.killsignals.fatalFail = False
                 self.updateTime = time.time() - timeA
-                if newGrid is not None:
+                if newGrid is not None and not self.getOption("useExternalDynamicMesh"):
+                    # since updateGeometryInfo assumes one time slice
+                    # when using ts with externally defined mesh, the grid is provided externally;
+                    # updateGeometryInfo mainly serves to update the metrics and etc.
                     self.adflow.warping.setgrid(newGrid)
             # Update geometric data, depending on the type of simulation
             if self.getOption("equationMode").lower() == "unsteady":
@@ -4762,6 +4781,8 @@ def _getDefaultOptions():
             "windAxis": [bool, False],
             "alphaFollowing": [bool, True],
             "TSStability": [bool, False],
+            "useTSInterpolatedGridVelocity": [bool, False],
+            "useExternalDynamicMesh": [bool, False],
             # Convergence Parameters
             "L2Convergence": [float, 1e-8],
             "L2ConvergenceRel": [float, 1e-16],
@@ -5117,6 +5138,7 @@ def _getOptionMap(self):
             "windaxis": ["stab", "usewindaxis"],
             "alphafollowing": ["stab", "tsalphafollowing"],
             "tsstability": ["stab", "tsstability"],
+            "usetsinterpolatedgridvelocity": ["ts", "usetsinterpolatedgridvelocity"],
             # Convergence Parameters
             "l2convergence": ["iter", "l2conv"],
             "l2convergencerel": ["iter", "l2convrel"],
@@ -5272,6 +5294,7 @@ def _getSpecialOptionLists(self):
             "cutcallback",
             "infchangecorrection",
             "skipafterfailedadjoint",
+            "useexternaldynamicmesh",
         }
 
         # Deprecated options that may be in old scripts and should not be used.

doc/options.yaml

@@ -554,6 +554,14 @@ useALE:
   desc: >
     Use the arbitrary Lagrangian-Eulerian (ALE) formulation for unsteady mesh deformations.
 
+useExternalDynamicMesh:
+  desc: >
+    Use externally provided deformed mesh for the time spectral equation. Each time instance gets its own deformed grid. This can be used for both an aerodynamic case or an aeroelastic case.
+
+useTSInterpolatedGridVelocity:
+  desc: >
+    Use spectral differentiation to compute the grid velocity for the time spectral equation.
+
 useGridMotion:
   desc: >
     Whether or not a rigid body motion of the grid has been specified.

src/f2py/adflow.pyf

@@ -331,6 +331,18 @@ python module libadflow
            integer(kind=inttype), optional,intent(in),check(len(famlist)>=nfamlist),depend(famlist) :: nfamlist=len(famlist)
            integer(kind=inttype) intent(in) :: sps
          end subroutine getsurfaceperturbation
+
+         subroutine setgridforoneinstance(grid,sps) ! in warping.F90:warping
+           use blockpointers, only: jl,x,il,ndom,kl
+           use section, only: sections,nsections
+           use inputphysics, only: equationmode
+           use preprocessingapi, only: xhalo
+           use utils, only: setpointers
+           use constants
+           real(kind=realtype) dimension(:),intent(in) :: grid
+           integer intent(in) :: sps
+         end subroutine setgridforoneinstance
+
        end module warping
 
        module surfaceutils
@@ -1176,11 +1188,13 @@ python module libadflow
          real(kind=realtype) allocatable,dimension(:,:,:) :: dscalar
          real(kind=realtype) allocatable,dimension(:,:,:) :: dvector
          real(kind=realtype) :: dtunsteadyrestartspectral
+         logical :: usetsinterpolatedgridvelocity
          logical :: writeunsteadyrestartspectral
          integer(kind=inttype) :: nunsteadysolspectral
          logical :: writeunsteadyvolspectral
          logical :: writeunsteadysurfspectral
          real(kind=realtype) allocatable,dimension(:,:,:) :: rotmatrixspectral
+         real(kind=realtype) :: omegafourier
        end module inputtimespectral
        module inputoverset
          logical :: useoversetloadbalance

src/modules/inputParam.F90

@@ -667,6 +667,9 @@ module inputTimeSpectral
 
   real(kind=realType), dimension(:,:,:), allocatable :: &
        rotMatrixSpectral
+  logical :: useTSInterpolatedGridVelocity
+
+  real(kind=realType) :: omegaFourier
 
 end module inputTimeSpectral
 

src/partitioning/partitioning.F90

@@ -849,6 +849,7 @@ subroutine timePeriodSpectral
          omegaFourYRot, omegaFourZRot, degreeFourMach, degreeFourAlpha, &
          degreeFourBeta
     use inputPhysics, only : equationMOde, flowType
+    use inputTimeSpectral, only : omegaFourier
     use section, only : sections, nSections
     use utils, only : terminate
     implicit none
@@ -969,6 +970,20 @@ subroutine timePeriodSpectral
           endif
        endif
 
+       ! aeroelastic case
+       if(omegaFourier > 0) then 
+         tt = two*pi/omegaFourier
+
+          ! Check if a time period was already determined. If so, try
+          ! to determine a common time. Otherwise just copy the data.
+
+          if( timeDetermined ) then
+             timePeriod = commonTimeSpectral(timePeriod, tt)
+          else
+             timePeriod     = tt
+             timeDetermined = .true.
+          endif
+       end if
 
 !!$         ! Altitude.
 !!$

src/preprocessing/preprocessingAPI.F90

@@ -3851,6 +3851,7 @@ subroutine updateGridVelocitiesAllLevels
        endif
 
        call gridVelocitiesFineLevel(.false., t, mm)
+
        call gridVelocitiesCoarseLevels(mm)
        call normalVelocitiesAllLevels(mm)
 
@@ -3889,7 +3890,10 @@ subroutine updatePeriodicInfoAllLevels
 
     call timePeriodSpectral
     call timeRotMatricesSpectral
-    call fineGridSpectralCoor
+    ! solve for the new grid only for rigid rotation with analytical deformation case
+    if (.NOT. usetsinterpolatedgridvelocity) then 
+       call fineGridSpectralCoor
+    end if
     call timeSpectralMatrices