AD: UA: no linearization for UA states that are off (OpenFAST#1716)

* AD: UA extrapinterp and input check

* AD: UA using BDF2 time integration scheme

* AD: UA using smaller increments for linearization of continuous states

* AD: UA not linearizing UA states that are off

* AD: changing r-test pointer (first pass)

* RTest: more verbose outputs for lin comparison and only abort at the end

* AD: changing r-test pointer (new baselines)

* AD: reverting to ABM4, DBF2 seem to cause problems for linearization

* Update of r-test pointer

* Update of r-test pointer (igen)

* Update of r-test pointer (to dev)

modules/aerodyn/src/UnsteadyAero.f90

@@ -761,22 +761,16 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg )
    end if
 
    ! Compute derivative step size
-   ! ---  ADENV
-   !p%dx    = 0.5_R8Ki * D2R_D
-   !p%dx(4) = 0.0001_R8Ki
-   ! --- ADLEG
-   p%dx = 2.0_R8Ki * D2R_D 
-   p%dx(4) = 0.001 ! x4 is a number between 0 and 1, so we need this to be small
-
-
-
+   p%dx    = 0.5_R8Ki * D2R_D
+   p%dx(4) = 0.0001_R8Ki
+
    p%UA_off_forGood = .false. ! flag that determines if UA should be turned off for the whole simulation
    if (allocated(InitInp%UAOff_innerNode)) then
       do j=1,min(size(p%UA_off_forGood,2), size(InitInp%UAOff_innerNode)) !blade
          do i=1,min(InitInp%UAOff_innerNode(j),size(p%UA_off_forGood,1)) !node
 !            call WrScr( 'Warning: Turning off Unsteady Aerodynamics on inner node (node '//trim(num2lstr(i))//', blade '//trim(num2lstr(j))//')' )
             p%UA_off_forGood(i,j) = .true.
-            !p%lin_nx = p%lin_nx - UA_NumLinStates
+            p%lin_nx = p%lin_nx - UA_NumLinStates
          end do
       end do
    end if
@@ -787,7 +781,7 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg )
 !            call WrScr( 'Warning: Turning off Unsteady Aerodynamics on outer node (node '//trim(num2lstr(i))//', blade '//trim(num2lstr(j))//')' )
             if (.not. p%UA_off_forGood(i,j)) then
                p%UA_off_forGood(i,j) = .true.
-               !p%lin_nx = p%lin_nx - UA_NumLinStates
+               p%lin_nx = p%lin_nx - UA_NumLinStates
             end if
          end do
       end do
@@ -801,7 +795,7 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg )
             if (ErrStat2 > ErrID_None) then
                call WrScr( 'Warning: Turning off Unsteady Aerodynamics because '//trim(ErrMsg2)//' (node '//trim(num2lstr(i))//', blade '//trim(num2lstr(j))//')' )
                p%UA_off_forGood(i,j) = .true.
-               !p%lin_nx = p%lin_nx - UA_NumLinStates
+               p%lin_nx = p%lin_nx - UA_NumLinStates
             end if
          end if
 
@@ -818,7 +812,7 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg )
       n = 1
       do j=1,size(p%UA_off_forGood,2) !blade
          do i=1,size(p%UA_off_forGood,1) !node
-            !if (.not. p%UA_off_forGood(i,j)) then
+            if (.not. p%UA_off_forGood(i,j)) then
                do k=1,UA_NumLinStates
                   p%lin_xIndx(n,1) = i ! node
                   p%lin_xIndx(n,2) = j ! blade
@@ -830,7 +824,7 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg )
                   endif
                   n = n + 1
                end do
-            !end if
+            end if
          end do
       end do
    end if
@@ -2297,7 +2291,19 @@ subroutine UA_UpdateStates( i, j, t, n, u, uTimes, p, x, xd, OtherState, AFInfo,
          call HGM_Steady( i, j, u_interp, p, x%element(i,j), AFInfo, ErrStat2, ErrMsg2 )
       end if
 
+      ! get inputs at t+dt
+      CALL UA_Input_ExtrapInterp( u, utimes, u_interp_raw, t+p%dt, ErrStat2, ErrMsg2 )
+         CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
+         IF ( ErrStat >= AbortErrLev ) RETURN
+
+         ! make sure that u%u is not zero (this previously turned off UA for the entire simulation. 
+         ! Now, we keep it on, but we don't want the math to blow up when we divide by u%u)
+      call UA_fixInputs(u_interp_raw, u_interp, ErrStat2, ErrMsg2)
+         call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
+
+      ! update states to value at t+dt:
       call UA_ABM4( i, j, t, n, u, utimes, p, x, OtherState, AFInfo, m, ErrStat2, ErrMsg2 )
+      !call UA_BDF2( i, j, t, n, u_interp, p, x, OtherState, AFInfo, m, ErrStat2, ErrMsg2 )
          call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
 
       if (.not. p%ShedEffect) then

reg_tests/executeOpenfastLinearRegressionCase.py

@@ -182,128 +182,164 @@ def indent(msg, sindent='\t'):
     exitWithError("An expected local solution file does not exist:")
 
 ### test for regression (compare lin files only)
-try:
-    for i, f in enumerate(localOutFiles):
-        local_file = os.path.join(testBuildDirectory, f)
-        baseline_file = os.path.join(targetOutputDirectory, f)
+
+def compareLin(f):
+    Errors     = []
+    ElemErrors = []
+
+    local_file = os.path.join(testBuildDirectory, f)
+    baseline_file = os.path.join(targetOutputDirectory, f)
+    local_file2 = local_file # Ellude the long filename
+
+    # Set a prefix for all errors to identify where it comes from
+    basename = os.path.splitext(os.path.basename(local_file))[0]
+    ext2 = os.path.splitext(basename)[1][1:] +'.lin' # '.1' or '.AD' '.BD'
+    errPrefix = CasePrefix[:-1]+ext2+': '
+
+    if verbose:
+        print(errPrefix+'file_ref:', baseline_file)
+        print(errPrefix+'file_new:', local_file)
+
+    def newError(msg):
+        msg=errPrefix+msg
         if verbose:
-            print(CasePrefix+'ref:', baseline_file)
-            print(CasePrefix+'new:', local_file)
-
-        # verify both files have the same number of lines
-        local_file_line_count = file_line_count(local_file)
-        baseline_file_line_count = file_line_count(baseline_file)
-        if local_file_line_count != baseline_file_line_count:
-            Err="Local and baseline solutions have different line counts in"
-            Err+="\n\tFile1:{}".format(local_file)
-            Err+="\n\tFile2:{}\n\n".format(baseline_file)
-            raise Exception(Err)
-
-        # open both files
-        floc = FASTLinearizationFile(local_file)
-        fbas = FASTLinearizationFile(baseline_file)
-
-        # --- Test that they have the same variables
-        kloc = floc.keys()
-        kbas = fbas.keys()
-        try:
-            np.testing.assert_equal(kloc, kbas)
-        except Exception as e:
-            Err = 'Different keys in local linfile.\n'
-            Err+= '\tNew:{}\n'.format(kloc)
-            Err+= '\tRef:{}\n'.format(kbas)
-            Err+= '\tin linfile: {}.\n'.format(local_file)
-            raise Exception(Err)
-
-        # --- Compare 10 first frequencies and damping ratios in 'A' matrix
-        if 'A' in fbas.keys(): 
-            Abas = fbas['A']
-            Aloc = floc['A']
-            # Note: we could potentially reorder states like MBC does, but no need for freq/damping
-            _, zeta_bas, _, freq_bas = eigA(Abas, nq=None, nq1=None, sort=True)
-            _, zeta_loc, _, freq_loc = eigA(Aloc, nq=None, nq1=None, sort=True)
-
-            if len(freq_bas)==0:
-                # We use complex eigenvalues instead of frequencies/damping
-                # If this fails often, we should discard this test.
-                _, Lambda = eig(Abas, sort=False)
-                v_bas = np.diag(Lambda)
-                _, Lambda = eig(Aloc, sort=False)
-                v_loc = np.diag(Lambda)
-
-                if verbose:
-                    print(CasePrefix+'val_ref:', v_bas[:7])
-                    print(CasePrefix+'val_new:', v_loc[:7])
-                try:
-                    np.testing.assert_allclose(v_bas[:10], v_loc[:10], rtol=rtol_f, atol=atol_f)
-                except Exception as e:
-                    raise Exception('Failed to compare A-matrix frequencies\n\tLinfile: {}.\n\tException: {}'.format(local_file, indent(e.args[0])))
-            else:
-
-                #if verbose:
-                print(CasePrefix+'freq_ref:', np.around(freq_bas[:8]    ,5), '[Hz]')
-                print(CasePrefix+'freq_new:', np.around(freq_loc[:8]    ,5),'[Hz]')
-                print(CasePrefix+'damp_ref:', np.around(zeta_bas[:8]*100,5), '[%]')
-                print(CasePrefix+'damp_new:', np.around(zeta_loc[:8]*100,5), '[%]')
-
-                try:
-                    np.testing.assert_allclose(freq_loc[:10], freq_bas[:10], rtol=rtol_f, atol=atol_f)
-                except Exception as e:
-                    raise Exception('Failed to compare A-matrix frequencies\n\tLinfile: {}.\n\tException: {}'.format(local_file, indent(e.args[0])))
-
-
-                if caseName=='Ideal_Beam_Free_Free_Linear':
-                    # The free-free case is a bit weird, same frequencies but dampings are +/- a value
-                    zeta_loc=np.abs(zeta_loc)
-                    zeta_bas=np.abs(zeta_bas)
-
-
-                try:
-                    # Note: damping ratios in [%]
-                    np.testing.assert_allclose(zeta_loc[:10]*100, zeta_bas[:10]*100, rtol=rtol_d, atol=atol_d)
-                except Exception as e:
-                    raise Exception('Failed to compare A-matrix damping ratios\n\tLinfile: {}.\n\tException: {}'.format(local_file, indent(e.args[0])))
-
-
-
-        # --- Compare individual matrices/vectors
-        KEYS= ['A','B','C','D','dUdu','dUdy']
-        KEYS+=['x','y','u','xdot']
-        for k,v in fbas.items():
-            if k in KEYS and v is not None:
-                if verbose:
-                    print(CasePrefix+'key:', k)
-                # Arrays
-                Mloc=np.atleast_2d(floc[k])
-                Mbas=np.atleast_2d(fbas[k])
-
-                # --- Compare dimensions
-                try:
-                    np.testing.assert_equal(Mloc.shape, Mbas.shape)
-                except Exception as e:
-                    Err = 'Different dimensions for variable `{}`.\n'.format(k)
-                    Err+= '\tNew:{}\n'.format(Mloc.shape)
-                    Err+= '\tRef:{}\n'.format(Mbas.shape)
-                    Err+= '\tLinfile: {}.\n'.format(local_file)
-                    raise Exception(Err)
-
-
-                # We for loop below to get the first element that mismatch
-                # Otherwise, do: np.testing.assert_allclose(floc[k], fbas[k], rtol=rtol, atol=atol)
-                for i in range(Mbas.shape[0]):
-                    for j in range(Mbas.shape[1]):
-                        # Old method:
-                        #if not isclose(Mloc[i,j], Mbas[i,j], rtol=rtol, atol=atol):
-                        #    sElem = 'Element [{},{}], new : {}, baseline: {}'.format(i+1,j+1,Mloc[i,j], Mbas[i,j])
-                        #    raise Exception('Failed to compare variable `{}`, {} \n\tLinfile: {}.'.format(k, sElem, local_file)) #, e.args[0]))
-                        try:
-                            np.testing.assert_allclose(Mloc[i,j], Mbas[i,j], rtol=rtol, atol=atol)
-                        except Exception as e:
-                            sElem = 'Element [{},{}], new : {}, baseline: {}'.format(i+1,j+1,Mloc[i,j], Mbas[i,j])
-                            raise Exception('Failed to compare variable `{}`, {} \n\tLinfile: {}.\n\tException: {}'.format(k, sElem, local_file, indent(e.args[0])))
-
-except Exception as e:
-    exitWithError(e.args[0])
+            print(msg)
+        Errors.append(msg)
+
+
+
+
+    # verify both files have the same number of lines
+    local_file_line_count = file_line_count(local_file)
+    baseline_file_line_count = file_line_count(baseline_file)
+    if local_file_line_count != baseline_file_line_count:
+        Err="Local and baseline solutions have different line counts in"
+        Err+="\n\tFile1:{}".format(local_file)
+        Err+="\n\tFile2:{}\n\n".format(baseline_file)
+        newError(Err)
+        #raise Exception(Err)
+
+    # open both files
+    floc = FASTLinearizationFile(local_file)
+    fbas = FASTLinearizationFile(baseline_file)
+
+    # --- Test that they have the same variables
+    kloc = floc.keys()
+    kbas = fbas.keys()
+    try:
+        np.testing.assert_equal(kloc, kbas)
+    except Exception as e:
+        Err = 'Different keys in local linfile.\n'
+        Err+= '\tNew:{}\n'.format(kloc)
+        Err+= '\tRef:{}\n'.format(kbas)
+        Err+= '\tin linfile: {}.\n'.format(local_file2)
+        newError(Err)
+        #raise Exception(Err)
+
+    # --- Compare 10 first frequencies and damping ratios in 'A' matrix
+    if 'A' in fbas.keys(): 
+        Abas = fbas['A']
+        Aloc = floc['A']
+        # Note: we could potentially reorder states like MBC does, but no need for freq/damping
+        _, zeta_bas, _, freq_bas = eigA(Abas, nq=None, nq1=None, sort=True, fullEV=True)
+        _, zeta_loc, _, freq_loc = eigA(Aloc, nq=None, nq1=None, sort=True, fullEV=True)
+
+        if len(freq_bas)==0:
+            # We use complex eigenvalues instead of frequencies/damping
+            # If this fails often, we should discard this test.
+            _, Lambda = eig(Abas, sort=False)
+            v_bas = np.diag(Lambda)
+            _, Lambda = eig(Aloc, sort=False)
+            v_loc = np.diag(Lambda)
+
+            if verbose:
+                print(errPrefix+'val_ref:', v_bas[:7])
+                print(errPrefix+'val_new:', v_loc[:7])
+            try:
+                np.testing.assert_allclose(v_bas[:10], v_loc[:10], rtol=rtol_f, atol=atol_f)
+            except Exception as e:
+                Err='Failed to compare A-matrix frequencies\n\tLinfile: {}.\n\tException: {}'.format(local_file2, indent(e.args[0]))
+                newError(Err)
+        else:
+
+            #if verbose:
+            print(errPrefix+'freq_ref:', np.around(freq_bas[:8]    ,5), '[Hz]')
+            print(errPrefix+'freq_new:', np.around(freq_loc[:8]    ,5),'[Hz]')
+            print(errPrefix+'damp_ref:', np.around(zeta_bas[:8]*100,5), '[%]')
+            print(errPrefix+'damp_new:', np.around(zeta_loc[:8]*100,5), '[%]')
+
+            try:
+                np.testing.assert_allclose(freq_loc[:10], freq_bas[:10], rtol=rtol_f, atol=atol_f)
+            except Exception as e:
+                Err = 'Failed to compare A-matrix frequencies\n\tLinfile: {}.\n\tException: {}'.format(local_file2, indent(e.args[0]))
+                newError(Err)
+
+
+            if caseName=='Ideal_Beam_Free_Free_Linear':
+                # The free-free case is a bit weird, smae frequencies but dampings are +/- a value
+                zeta_loc=np.abs(zeta_loc)
+                zeta_bas=np.abs(zeta_bas)
+
+
+            try:
+                # Note: damping ratios in [%]
+                np.testing.assert_allclose(zeta_loc[:10]*100, zeta_bas[:10]*100, rtol=rtol_d, atol=atol_d)
+            except Exception as e:
+                Err = 'Failed to compare A-matrix damping ratios\n\tLinfile: {}.\n\tException: {}'.format(local_file2, indent(e.args[0]))
+                newError(Err)
+
+    # --- Compare individual matrices/vectors
+    KEYS= ['A','B','C','D','dUdu','dUdy']
+    KEYS+=['x','y','u','xdot']
+    for k,v in fbas.items():
+        if k in KEYS and v is not None:
+            if verbose:
+                print(errPrefix+'key:', k)
+            # Arrays
+            Mloc=np.atleast_2d(floc[k])
+            Mbas=np.atleast_2d(fbas[k])
+
+            # --- Compare dimensions
+            try:
+                np.testing.assert_equal(Mloc.shape, Mbas.shape)
+                dimEqual=True
+            except Exception as e:
+                Err = 'Different dimensions for variable `{}`.\n'.format(k)
+                Err+= '\tNew:{}\n'.format(Mloc.shape)
+                Err+= '\tRef:{}\n'.format(Mbas.shape)
+                Err+= '\tLinfile: {}.\n'.format(local_file2)
+                newError(Err)
+                dimEqual=False
+
+            if not dimEqual:
+                # We don't compare elements if shapes are different
+                continue
+
+            # We for loop below to get the first element that mismatch
+            # Otherwise, do: np.testing.assert_allclose(floc[k], fbas[k], rtol=rtol, atol=atol)
+            for i in range(Mbas.shape[0]):
+                for j in range(Mbas.shape[1]):
+                    # Old method:
+                    #if not isclose(Mloc[i,j], Mbas[i,j], rtol=rtol, atol=atol):
+                    #    sElem = 'Element [{},{}], new : {}, baseline: {}'.format(i+1,j+1,Mloc[i,j], Mbas[i,j])
+                    #    raise Exception('Failed to compare variable `{}`, {} \n\tLinfile: {}.'.format(k, sElem, local_file)) #, e.args[0]))
+                    try:
+                        np.testing.assert_allclose(Mloc[i,j], Mbas[i,j], rtol=rtol, atol=atol)
+                    except Exception as e:
+                        sElem = 'Element [{},{}], new : {}, baseline: {}'.format(i+1,j+1,Mloc[i,j], Mbas[i,j])
+                        Err=errPrefix+'Failed to compare variable `{}`, {} \n\tLinfile: {}.\n\tException: {}'.format(k, sElem, local_file2, indent(e.args[0]))
+                        ElemErrors.append(Err)
+    return Errors, ElemErrors
+
+Errors=[]
+for i, f in enumerate(localOutFiles):
+    ErrorsLoc, ElemErrorsLoc = compareLin(f)
+    Errors += ErrorsLoc
+    if len(ElemErrorsLoc)>0:
+        Errors += ElemErrorsLoc[:3] # Just a couple of them
+
+if len(Errors)>0:
+    exitWithError('See errors below: \n'+'\n'.join(Errors))
 
 # passing case
 sys.exit(0)

reg_tests/lib/eva.py

@@ -158,18 +158,18 @@ def eigA(A, nq=None, nq1=None, fullEV=False, normQ=None, sort=True):
 
     if m!=n:
         raise Exception('Matrix needs to be squared')
-    if nq is None:
-        if nq1 is None:
-            nq1=0
-        nq = int((n-nq1)/2)
-    else:
-        nq1 = n-2*nq
-    if n!=2*nq+nq1 or nq1<0:
-        raise Exception('Number of 1st and second order dofs should match the matrix shape (n= 2*nq + nq1')
     Q, Lambda = eig(A, sort=False)
     v = np.diag(Lambda)
 
     if not fullEV:
+        if nq is None:
+            if nq1 is None:
+                nq1=0
+            nq = int((n-nq1)/2)
+        else:
+            nq1 = n-2*nq
+        if n!=2*nq+nq1 or nq1<0:
+            raise Exception('Number of 1st and second order dofs should match the matrix shape (n= 2*nq + nq1')
         Q=np.delete(Q, slice(nq,2*nq), axis=0)
 
     # Selecting eigenvalues with positive imaginary part (frequency)