Issue #18: Mass scaling engine for superelements

src/vpmCommon/supElTypeModule.f90

@@ -144,6 +144,8 @@ module SupElTypeModule
      real(dp) :: dmpScl(2)  !< Current and previous scaling of mDmp0 and kDmp0
      integer  :: dmpSclIdx  !< Index to structural damping scaling function
 
+     real(dp) :: massScl(2) !< Current and previous mass scaling factors
+     integer  :: massSclIdx !< Index to mass scaling function
      real(dp) :: stifScl(2) !< Current and previous stiffness scaling factors
      integer  :: stifSclIdx !< Index to stiffness scaling function
 
@@ -356,6 +358,7 @@ subroutine WriteSupElType (sup,io,complexity)
     write(io,3) 'mDmpFactor    =', sup%mDmpFactor
     write(io,3) 'kDmpFactor    =', sup%kDmpFactor
     write(io,*) 'dmpScale(id)  =', sup%dmpSclIdx
+    write(io,*) 'massScale(id) =', sup%massSclIdx
     write(io,*) 'stifScale(id) =', sup%stifSclIdx
     write(io,*) 'savePos       =', sup%savePos
     write(io,*) 'saveVar       =', sup%saveVar
@@ -666,8 +669,10 @@ subroutine NullifySupEl (sup)
     sup%mDmp0 = 0.0_dp
     sup%kDmp0 = 0.0_dp
     sup%dmpScl = 1.0_dp
+    sup%massScl = 1.0_dp
     sup%stifScl = 1.0_dp
     sup%dmpSclIdx = 0
+    sup%massSclIdx = 0
     sup%stifSclIdx = 0
     sup%savePos = .false.
     sup%saveVar = .false.
@@ -1273,6 +1278,7 @@ module SupElNamelistModule
   integer :: recoveryFlag          !< Flag for enabling stress recovery
   integer :: BC(maxGenDofs_p)      !< Boundary conditions for component modes
   integer :: strDmpEngineId        !< Base ID for structural damping function
+  integer :: massEngineId          !< Base ID for mass scaling function
   integer :: stiffEngineId         !< Base ID for stiffness scaling function
 
   integer :: savePos    !< Flag indicating whether position should be saved
@@ -1307,7 +1313,7 @@ module SupElNamelistModule
        &            refTriad2Id, offset2, &
        &            refTriad3Id, offset3, &
        &            stressStiffFlag, massCorrFlag, recoveryFlag, &
-       &            stiffEngineId, stiffScale, massScale, &
+       &            stiffEngineId, massEngineId, stiffScale, massScale, &
        &            dragParams, slamParams, &
        &            strDmpEngineId, alpha1, alpha2, alpha3, alpha4, BC, &
        &            savePos, saveVar
@@ -1343,7 +1349,8 @@ subroutine read_SUP_EL (infp,stat)
     stressStiffFlag=-1; massCorrFlag=-1; recoveryFlag=-1
     stiffScale=1.0_dp; massScale=1.0_dp; dragParams=0.0_dp; slamParams=0.0_dp
     alpha1=0.0_dp; alpha2=0.0_dp; alpha3=0.0_dp; alpha4=0.0_dp
-    stiffEngineId=0; strDmpEngineId=0; BC=1; savePos=0; saveVar=0
+    massEngineId=0; stiffEngineId=0; strDmpEngineId=0
+    BC=1; savePos=0; saveVar=0
 
     read(infp,nml=SUP_EL,iostat=stat)
 

src/vpmSolver/addInSysModule.f90

@@ -40,7 +40,7 @@ module AddInSysModule
   !> @param Rhs Right-hand-side force vector
   !>
   !> @details The Newton matrix is a linear combination of the mass-, damping-
-  !> and stiffness matrices, Nmat = scaleM*M  + scaleC*C + scaleK*K.
+  !> and stiffness matrices, Nmat = scaleM*M + scaleC*C + scaleK*K.
   !> If @a Rhs is provided, contributions due to prescribed motions,
   !> if any, will be added.
   !>
@@ -142,6 +142,7 @@ subroutine BuildNewtonMat (Nmat, scaleM, scaleC, scaleK, mech, sam, iter, &
     do i = 1, size(mech%sups)
 
        if ( reComputeSupMat .or. mech%sups(i)%stressStiffFlag(1) == 1 .or. &
+            abs(mech%sups(i)%massScl(1)-mech%sups(i)%massScl(2)) > eps_p .or. &
             abs(mech%sups(i)%stifScl(1)-mech%sups(i)%stifScl(2)) > eps_p ) then
 
           call BuildSupNewtonMat (newTanStiff, scaleM, scaleC, scaleK, &
@@ -353,7 +354,8 @@ subroutine BuildMassMat (Mmat, mech, sam, ierr)
 
     do i = 1, size(mech%sups)
        if (associated(mech%sups(i)%Mmat)) then
-          call addInSupMat (mech%sups(i)%Mmat, Mmat, mech%sups(i), sam, ierr)
+          call addInSupMat (mech%sups(i)%Mmat, Mmat, mech%sups(i), sam, ierr, &
+               &            scale=mech%sups(i)%massScl(1))
           if (ierr < 0) goto 900
        end if
        if (associated(mech%sups(i)%Mamat)) then

src/vpmSolver/initiateSupElTypeModule.f90

@@ -551,7 +551,7 @@ subroutine ReadSupEls (infp,env,triads,sups,err)
           end if
        end if
 
-       scaledStrDmp = scaledStructDamp .or. stiffEngineId > 0
+       scaledStrDmp = scaledStructDamp .or. stiffEngineId>0 .or. massEngineId>0
        if (scaledStrDmp .or. .not.haveDamping(0)) then
 
           !! Apply stiffness and/or mass scaling factors
@@ -610,6 +610,19 @@ subroutine ReadSupEls (infp,env,triads,sups,err)
           sups(idIn)%genDOFs%alpha2 = alpha4(1:numGenDOFs)
        end if
 
+       if (massEngineId > 0) then
+
+          !! Time-dependent mass scaling
+          if (haveDamping(2)) then
+             err = err - 1
+             call ReportInputError('SUP_EL',idIn,sups(idIn)%id, &
+                  'Stiffness-proportional damping can not be combined with '// &
+                  'time-dependent mass scaling.')
+          else
+             sups(idIn)%massSclIdx = massEngineId
+          end if
+
+       end if
        if (stiffEngineId > 0) then
 
           !! Time-dependent stiffness scaling
@@ -903,6 +916,14 @@ subroutine InitiateSupEls2 (sam,sys,sups,masses,engines,err)
                   &                'Non-existing engine referred')
           end if
        end if
+       if (sups(idIn)%massSclIdx > 0) then
+          if (.not. associated(GetPtrToId(engines,sups(idIn)%massSclIdx, &
+               &                          index=sups(idIn)%massSclIdx))) then
+             lerr = lerr - 1
+             call ReportInputError('SUP_EL',idIn,sups(idIn)%id, &
+                  &                'Non-existing engine referred')
+          end if
+       end if
        if (sups(idIn)%stifSclIdx > 0) then
           if (.not. associated(GetPtrToId(engines,sups(idIn)%stifSclIdx, &
                &                          index=sups(idIn)%stifSclIdx))) then

src/vpmSolver/supElRoutinesModule.f90

@@ -298,7 +298,8 @@ subroutine UpdateSupEls (sups,supLoads,env,beta,gamma,time,timeStep, &
 
        if (associated(sups(i)%Mmat)) then
           !! Calculate the inertia forces in the superelement, Fi = M*uldd
-          call scaledMatmul (m,n,1.0_dp,sups(i)%Mmat,sups(i)%uldd,sups(i)%Fi)
+          call scaledMatmul (m,n,sups(i)%massScl(1), &
+               &             sups(i)%Mmat,sups(i)%uldd,sups(i)%Fi)
        else
           call DCOPY (m,0.0_dp,0,sups(i)%Fi(1),1)
        end if
@@ -335,7 +336,7 @@ subroutine UpdateSupEls (sups,supLoads,env,beta,gamma,time,timeStep, &
 
        if (associated(sups(i)%fg)) then
           !! Calculate the gravitational forces which is an external force
-          call scaledMatmul (m,3,sups(i)%stifScl(1),sups(i)%fg, &
+          call scaledMatmul (m,3,sups(i)%massScl(1),sups(i)%fg, &
                &             matmul(env%gravity,sups(i)%supTr(:,1:3)),sups(i)%Q)
        else
           call DCOPY (m,0.0_dp,0,sups(i)%Q(1),1)
@@ -698,7 +699,8 @@ end subroutine UpdateSupElLoad
   !> @param[out] ierr Error flag
   !>
   !> @details Effective with time-dependent structural damping coefficients.
-  !> This subroutine also updates the time-dependent stiffness scaling factor.
+  !> This subroutine also updates the time-dependent stiffness- and mass
+  !> scaling factors.
   !>
   !> @callgraph @callergraph
   !>
@@ -739,6 +741,10 @@ subroutine UpdateSupElDamping (sups,engs,alpha,ierr)
        else
           newDmp = .false.
        end if
+       if (sups(i)%massSclIdx > 0) then
+          sups(i)%massScl(2) = sups(i)%massScl(1)
+          sups(i)%massScl(1) = EngineValue(engs(sups(i)%massSclIdx),ierr)
+       end if
        if (sups(i)%stifSclIdx > 0) then
           sups(i)%stifScl(2) = sups(i)%stifScl(1)
           sups(i)%stifScl(1) = EngineValue(engs(sups(i)%stifSclIdx),ierr)
@@ -1271,7 +1277,7 @@ subroutine BuildSupNewtonMat (newTangent,scaleM,scaleC,scaleK,sup,ierr)
 
        !! Structural mass matrix contribution
        call DCOPY (n,sup%Mmat(1,1),1,sup%Nmat(1,1),1)
-       call DSCAL (n,scaleM,sup%Nmat(1,1),1)
+       call DSCAL (n,scaleM*sup%massScl(1),sup%Nmat(1,1),1)
 
     else
        call DCOPY (n,0.0_dp,0,sup%Nmat(1,1),1)