multi_compute_dt.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    multi_compute_dt   ../engine/source/multifluid/multi_compute_dt.F
!||--- called by ------------------------------------------------------
!||    alemain            ../engine/source/ale/alemain.F
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod       ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    initbuf_mod        ../engine/share/resol/initbuf.F
!||    multi_fvm_mod      ../common_source/modules/ale/multi_fvm_mod.f90
!||====================================================================
      SUBROUTINE multi_compute_dt(DT2T, ELBUF_TAB, IPARG, ITASK, IXS, IXQ, IXTG,
     .     PM, IPM, MULTI_FVM, WGRID, XGRID,
     .     NELTST, ITYPTST)
C-----------------------------------------------
C     M o d u l e s
C-----------------------------------------------
      USE initbuf_mod      
      USE elbufdef_mod
      USE multi_fvm_mod
C-----------------------------------------------
C     I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C     C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "param_c.inc"
#include      "task_c.inc"
#include      "mvsiz_p.inc"
#include      "scr18_c.inc"         
#include      "units_c.inc"        
#include      "com08_c.inc"           
#include      "comlock.inc"        
C-----------------------------------------------
C     D u m m y   A r g u m e n t s
C-----------------------------------------------
      my_real, INTENT(OUT) :: dt2t
      TYPE(elbuf_struct_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
      INTEGER, INTENT(IN) :: IPARG(NPARG, *)
      INTEGER, INTENT(IN) :: ITASK ! SMP TASK
      INTEGER, INTENT(IN) :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
      INTEGER, INTENT(IN) :: IPM(NPROPMI, *)
      my_real, INTENT(IN) :: pm(npropm, *)
      TYPE(multi_fvm_struct), INTENT(INOUT) :: MULTI_FVM
      my_real, INTENT(IN) :: wgrid(*), xgrid(3, *)
      INTEGER, INTENT(OUT) :: ITYPTST, NELTST
C-----------------------------------------------
C     L o c a l   V a r i a b l e s
C-----------------------------------------------
      TYPE(g_bufel_), POINTER :: GBUF
      INTEGER :: NG, NEL, II, JJ, KFACE, I, J, NB_FACE, NFT, ITY
      INTEGER :: IPLA
      INTEGER :: NBMAT, IMAT
      INTEGER :: NODE1, NODE2, NODE3, NODE4, 
     .     node5, node6, node7, node8
      my_real :: w1(3), w2(3), w3(3), w4(3), 
     .     w5(3), w6(3), w7(3), w8(3)
      my_real :: x1(3), x2(3), x3(3), x4(3), 
     .     x5(3), x6(3), x7(3), x8(3)
      my_real :: wfac(1:3), surf
      my_real :: lambdaii, lambdaf, normuii, normujj
      my_real :: fii(5), fjj(5), normal_vel, normal_vel2, vii(5), vjj(5), vel2
      my_real :: dtel(mvsiz), nx, ny, nz
      INTEGER :: ISOLNOD, MATLAW
      LOGICAL :: l_FOUND_LOWER
 
C     Time step
      dt2t = zero
      ityptst = 0
      neltst = 0
 
      DO ng = itask + 1, ngroup, nthread
         matlaw = iparg(1, ng)         
         IF (matlaw == 151) THEN
            nel = iparg(2, ng)
            nft = iparg(3, ng)
            ity = iparg(5, ng)
            isolnod = iparg(28, ng)
            gbuf => elbuf_tab(ng)%GBUF
C     DELTAX is to be kept
            gbuf%DELTAX(1:nel) = zero
C     Number of faces in an element
            nb_face = 6 
            IF (ity == 2) THEN
               nb_face = 4
            ELSEIF (ity == 7) THEN
               nb_face = 3
            ENDIF
            
C     Flux computation
            dtel(1:nel) = zero
            DO ii = 1, nel
               i = ii + nft               
C     Face KFACE
               DO kface = 1, nb_face              
                  nx = multi_fvm%FACE_DATA%NORMAL(1, kface, i)
                  ny = multi_fvm%FACE_DATA%NORMAL(2, kface, i)
                  nz = multi_fvm%FACE_DATA%NORMAL(3, kface, i)
                  wfac(1:3) = multi_fvm%FACE_DATA%WFAC(1:3, kface, i)
                  surf = multi_fvm%FACE_DATA%SURF(kface, i)
C     Time step
                  normal_vel2 = (multi_fvm%VEL(1, i) - wfac(1)) * nx + 
     .                 (multi_fvm%VEL(2, i) - wfac(2)) * ny + 
     .                 (multi_fvm%VEL(3, i) - wfac(3)) * nz
                  dtel(ii) = max(dtel(ii), 
     .                 surf / gbuf%VOL(ii) * (multi_fvm%SOUND_SPEED(i) + abs(normal_vel2)) / dtfac1(102))
                  gbuf%DELTAX(ii) = max(gbuf%DELTAX(ii), surf / gbuf%VOL(ii))
               ENDDO            !KFACE
               gbuf%DELTAX(ii) = one / gbuf%DELTAX(ii)
            ENDDO ! ii = 1, nel
C----------------------
C     Globalize time step for this group
C----------------------
            l_found_lower=.false.
            DO ii = 1, nel
               IF(dtel(ii)>zero)gbuf%DT(ii) = one/dtel(ii)
               IF (dtel(ii) > dt2t) THEN
                  l_found_lower=.true.
                  dt2t = dtel(ii)
                  ityptst = ity
                  IF (multi_fvm%SYM == 0) THEN
                     neltst = ixs(nixs, ii + nft)
                  ELSE
                     IF (ity == 2) THEN
C     QUADS
                        neltst = ixq(nixq, ii + nft)
                     ELSEIF (ity == 7) THEN
C     TRIANGLES
                        neltst = ixtg(nixtg, ii + nft)
                     ENDIF
                  ENDIF
               ENDIF
            ENDDO
            
            !CHECK IF LOWER THAN DTMIN
            IF(l_found_lower .AND. dt2t/=zero)THEN
              IF(one/dt2t<dtmin1(102))THEN
                tstop = tt
#include "lockon.inc"
                WRITE(iout,*) ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR ALE/EULER CELL',neltst
                WRITE(istdo,*)' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR ALE/EULER CELL',neltst
#include "lockoff.inc"          
              ENDIF  
            ENDIF           
            
         ENDIF
      ENDDO  ! NG = ITASK + 1, NGROUP, NTHREAD
C----------------------
C     Global time step
C----------------------
      IF (dt2t > zero) THEN
         dt2t = one / dt2t
      ELSE
         dt2t = ep30
      ENDIF
C----------------------
C     Boundary fluxes
C----------------------
      END SUBROUTINE multi_compute_dt