inigrav__eos_8F_source.html

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/.

!||====================================================================

!||    inigrav_eos        ../starter/source/initial_conditions/inigrav/inigrav_eos.F

!||--- called by ------------------------------------------------------

!||    inigrav_load       ../starter/source/initial_conditions/inigrav/inigrav_load.F

!||--- calls      -----------------------------------------------------

!||--- uses       -----------------------------------------------------

!||====================================================================


      SUBROUTINE inigrav_eos(NELG, NEL, NG, MATID, IPM, GRAV0, RHO0, DEPTH, PM, BUFMAT,

     .                       ELBUF_TAB,PSURF,LIST,PRES,IMAT,MLW,NPF,TF,NUMMAT,MAT_PARAM)

C-----------------------------------------------

C   M o d u l e s

C-----------------------------------------------

      USE elbufdef_mod

      USE matparam_def_mod

      USE eosmain_mod , ONLY : eosmain

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      "param_c.inc"

#include      "com01_c.inc"

#include      "tabsiz_c.inc"

C-----------------------------------------------

C   D u m m y   A r g u m e n t s

C-----------------------------------------------

      INTEGER, INTENT(IN)    :: NUMMAT

      INTEGER, INTENT(IN)    :: NEL, NG, MATID, IPM(NPROPMI, *), LIST(NEL),NELG, IMAT, MLW

      my_real, INTENT(IN)    :: grav0, depth(*), pm(npropm, *), bufmat(*),psurf

      my_real, INTENT(INOUT) :: pres(nel)

      TYPE(elbuf_struct_), DIMENSION(NGROUP), TARGET, INTENT(IN) :: ELBUF_TAB

      my_real, INTENT(IN) :: rho0

      INTEGER,INTENT(IN)::NPF(SNPC)

      my_real,INTENT(IN)::tf(stf)

      TYPE(matparam_struct_) ,DIMENSION(NUMMAT), INTENT(IN) :: MAT_PARAM

C-----------------------------------------------

C   L o c a l   V a r i a b l e s

C-----------------------------------------------

      INTEGER :: I,K

      INTEGER :: EOSTYPE, ITER, MAX_ITER, MAT(NEL), REMAINING_ELTS, NBMAT, NVAREOS, NVARTMP_EOS

      my_real :: bid(nel) !bfrac not used int his context

      my_real :: tol, incr

      my_real :: func(nel), dfunc(nel), error(nel), rhozero(nel), rho(nel),

     .           mu(nel), mu2(nel), espe(nel), dvol(nel), df(nel), vol(nel), psh(nel),

     .           dpdrho(nel), dpdm(nel), dpde(nel), theta(nel), ecold(nel), p0(nel),

     .           grun(nel),pref(nel), off(nel), eint(nel), sig(nel,6),rho_old,v0(nel)


      LOGICAL :: CONT, CONVERGED(NEL)

      TYPE(g_bufel_), POINTER :: GBUF

      TYPE(L_BUFEL_), POINTER :: LBUF

      TYPE(buf_eos_), POINTER :: EBUF

C-----------------------------------------------

C   S o u r c e   L i n e s

C-----------------------------------------------


C LIST IS SUBGROUP TO TREAT : ONLY ELEM WITH RELEVANT PARTS ARE KEPT

C NEL IS ISEZ OF LIST

C NELG IS SIZE OF ORIGINAL GROUP : needed to shift indexes in GBUF%SIG & MBUF%VAR

C IMAT : relevant for multifluid eos. imat=0 => monomaterial law


      gbuf => elbuf_tab(ng)%GBUF

      NULLIFY(lbuf)

      NULLIFY(ebuf)

      ebuf =>  elbuf_tab(ng)%BUFLY(1)%EOS(1, 1, 1)

      nvareos = elbuf_tab(ng)%BUFLY(1)%NVAR_EOS

      nvartmp_eos = elbuf_tab(ng)%BUFLY(1)%NVARTMP_EOS

      IF(imat>0)THEN

        lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)

        ebuf =>  elbuf_tab(ng)%BUFLY(imat)%EOS(1, 1, 1)

        nvareos = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS

        nvartmp_eos = elbuf_tab(ng)%BUFLY(imat)%NVARTMP_EOS

      ENDIF


      eostype          = mat_param(matid)%IEOS

      eostype          = ipm(4, matid)    ! still used temporarily in this section

      iter             = 0

      max_iter         = 30

      cont             =.true.

      error(1:nel)     = ep30

      tol              = em4

      converged(1:nel) =.false.

      dvol(1:nel)      = zero

      theta(1:nel)     = zero


      !===================================================================!

      !                     BUFFER INITIALIZATION                         !

      !===================================================================!

      ! local array instead of pointers since they may be reshaped with a

      ! subgroup (1:NEL) c (1:NELG) depending on element PART (possibly several PARTS in a given group)

      DO k = 1, nel

           rhozero(k) = pm(1,matid)

           psh(k)     = pm(88, matid)

           mat(k)     = matid

           pref(k)    = psurf - rho0 * grav0 * depth(k)   !RHO in : RHO0_mixture (gravity model applied to the mixture not to the single submaterial)

           vol(k)     = one   !  Real volume value not useful here, just an indicator : ONE = compute, ZERO = do not compute

           mu(k)      = zero

           sig(k,1)   = -pm(31,matid)

           sig(k,2)   = -pm(31,matid)

           sig(k,3)   = -pm(31,matid)

           sig(k,4)   = zero

           sig(k,5)   = zero

           sig(k,6)   = zero

      ENDDO


      IF(imat==0)THEN   ! WORK IN GBUF

        DO k = 1, nel

            i         = list(k)

           rho(k)     = gbuf%RHO(k)

           off(k)     = gbuf%OFF(i)

           eint(k)    = gbuf%EINT(i)

           v0(k)      = gbuf%VOL(i)

        ENDDO

      ELSE ! THEN IMAT>0    WORK IN LBUF

        DO k = 1, nel

           i          = list(k)

           rho(k)     = pm(1,matid) !submaterial

           off(k)     = lbuf%OFF(i)

           eint(k)    = lbuf%EINT(i)

           v0(k)      = lbuf%VOL(i)

        ENDDO

      ENDIF

      !EINT=Eint/V in this context with law 151

      !===================================================================!

      !                     ITERATIONS                                    !

      !===================================================================!

      remaining_elts = nel                      !Number of elements for which convergence is not achieved

      DO WHILE(cont .AND. iter < max_iter)   !Start nonlinear solver loop

         iter = iter + 1

         DO k = 1, nel

            mu(k) = rho(k) / rhozero(k) - one

            df(k) = rhozero(k) / rho(k)

            eint(k) = eint(k) / df(k) !with IFLAG==2 we need ESPE=Df*EINT

         ENDDO

         CALL eosmain(2     , nel        , eostype  , pm     , off  , eint,

     .                rho   , rhozero    , mu       , mu2    , espe ,

     .                dvol  , df         , vol      , mat    , psh  ,

     .                pres  , dpdm       , dpde     , theta  ,

     .                bufmat, sig        , mu       , mlw    ,

     .                npf   , tf         , ebuf%VAR , nvareos, mat_param(matid),

     .                bid  , nvartmp_eos, ebuf%VARTMP)

          DO k=1,nel

            eint(k)=eint(k)*df(k) ! set back to Eint/V

          ENDDO

         !===================================================================!

         !                     CHECK CONVERGENCY                             !

         !===================================================================!

         !loop over retained element in the current group I \in LIST(1:NEL)

         DO k = 1, nel

            IF (.NOT. converged(k)) THEN   ! No need to take the square root, only the squared sound speed is needed here

               dpdrho(k) = dpdm(k) / rhozero(k)  !DPDM is total derivative (variation along isentropic path)

               !GRUN(K) = DPDE(K) / (ONE + MU(K))

               func(k) = pres(k)-pref(k)

               dfunc(k) = dpdrho(k) !- grun(k) * pres(k) / rho(k)

               incr = - func(k) / dfunc(k)

               rho_old = rho(k)

               rho(k) = rho(k) + incr

               error(k)= max( abs(incr)/rhozero(k) , abs(func(k)) )


               incr =   - (pres(k)+psh(k)+half*dpdrho(k)*incr)*(rhozero(k)/rho(k)-rhozero(k)/rho_old)

               eint(k) = eint(k) + incr   !EINT is in fact Eint/V with law 151 here

               IF(dpde(k)==zero)eint(k)=zero


               IF (error(k) < tol .OR. iter == max_iter) THEN

                  converged(k) = .true.

                  remaining_elts = remaining_elts - 1

                  vol(k) = zero

                  i = list(k)

                  df(k) = rhozero(i) / rho(k)

                  IF(imat==0)THEN

                    gbuf%RHO(i) = rho(k)

                    gbuf%SIG(i           ) = - pres(k)

                    gbuf%SIG(i +     nelg) = - pres(k)

                    gbuf%SIG(i + 2 * nelg) = - pres(k)

                    gbuf%EINT(i) = eint(k)/df(k)     !Eint/V with law 151 here     => EINT*DF == Eint/V0

                  ELSE

                    lbuf%RHO(i) = rho(k)

                    lbuf%SIG(i           ) = - pres(k)

                    lbuf%SIG(i +     nelg) = - pres(k)

                    lbuf%SIG(i + 2 * nelg) = - pres(k)

                    lbuf%EINT(i) = eint(k)/df(k)     !Eint/V with law 151 here

                  ENDIF

               ENDIF

            ENDIF

         ENDDO

         cont = (remaining_elts /= 0)

      ENDDO                     ! WHILE


      END SUBROUTINE inigrav_eos

my_real
#define my_real
Definition cppsort.cpp:32

inigrav_eos
subroutine inigrav_eos(nelg, nel, ng, matid, ipm, grav0, rho0, depth, pm, bufmat, elbuf_tab, psurf, list, pres, imat, mlw, npf, tf, nummat, mat_param)
Definition inigrav_eos.F:32

max
#define max(a, b)
Definition macros.h:21

eosmain_mod
Definition eosmain.F:35

eosmain_mod::eosmain
subroutine eosmain(iflag, nel, eostyp, pm, off, eint, rho, rho0, mu, mu2, espe, dvol, df, vnew, mat, psh, pnew, dpdm, dpde, theta, bufmat, sig, mu_bak, mlw, npf, tf, vareos, nvareos, mat_param, bfrac, nvartmp, vartmp)
Definition eosmain.F:102