topic Re: Abaqus VUMAT error during compilation in Intel® Fortran Compiler

Abaqus VUMAT error during compilation

Domv — Wed, 17 Jan 2024 16:45:54 GMT

I need help to solve the following problem, i understood there is an array already allocated, but i don't know how find it, can someone help me?

forrtl: severe (151): allocatable array is already allocated
Image PC Routine Line Source
libifcoremd.dll 00007FF97DE2CC0E Unknown Unknown Unknown
explicitU.dll 00007FF94DCE6F49 Unknown Unknown Unknown
explicitU.dll 00007FF94DCE2E1E Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92CFE2347 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92D08438D Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92D0ABD6A Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92C6F2874 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92C6CA58C Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92C5BB0E8 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E16B832 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E16033C Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E1A63AB Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E1A8F2F Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E19BA86 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E19E314 Unknown Unknown Unknown
ABQSMAPacModule.d 00007FF92E1B08F9 Unknown Unknown Unknown
package.exe 00007FF78275FD09 Unknown Unknown Unknown
package.exe 00007FF78275DD86 Unknown Unknown Unknown
package.exe 00007FF782761210 Unknown Unknown Unknown
KERNEL32.DLL 00007FF9B0AD244D Unknown Unknown Unknown
ntdll.dll 00007FF9B272DFB8 Unknown Unknown Unknown

@Steve_Lionel

Re: Abaqus VUMAT error during compilation

andrew_4619 — Wed, 17 Jan 2024 17:43:20 GMT

well the error message explains itself. Why not post your VUMAT routine and suggestions on how to fix it might be possible then

Re: Abaqus VUMAT error during compilation

Domv — Wed, 17 Jan 2024 17:52:29 GMT

parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
3 four=4.d0,
2 third=one/three, half= .5d0, twoThirds=two/three,
3 threeHalfs = 1.5d0, fourth=one/four, six=6.d0)

character*80 cmname
C
do i = 1,nblock

C Read material properties

a_a=props(1)

a_b=a_a

b_a=props(2)

b_b=b_a

d_a=props(3)

d_b=d_a

pmin_a=props(4)

pmin_b=pmin_a

c_a=props(5)

c_b=c_a

s1=props(6)

s2=props(7)

s3=props(8)

sT_a=props(9)

sT_b=sT_a

sC_a=props(10)

sC_b=sC_a

c ----------

c Deformation gradiEnt before the increment

F_old(1,1)=defgradOld(i,1)
F_old(1,2)=defgradOld(i,4)
F_old(1,3)=defgradOld(i,9)
F_old(2,1)=defgradOld(i,7)
F_old(2,2)=defgradOld(i,2)
F_old(2,3)=defgradOld(i,5)
F_old(3,1)=defgradOld(i,6)
F_old(3,2)=defgradOld(i,8)
F_old(3,3)=defgradOld(i,3)

c Deformation gradiant after the increment

F_new(1,1)=defgradNew(i,1)
F_new(1,2)=defgradNew(i,4)
F_new(1,3)=defgradNew(i,9)
F_new(2,1)=defgradNew(i,7)
F_new(2,2)=defgradNew(i,2)
F_new(2,3)=defgradNew(i,5)
F_new(3,1)=defgradNew(i,6)
F_new(3,2)=defgradNew(i,8)
F_new(3,3)=defgradNew(i,3)

c Stretch tensor after the increment

xU(1,1)=stretchNew(i,1)
xU(2,2)=stretchNew(i,2)
xU(3,3)=stretchNew(i,3)
xU(1,2)=stretchNew(i,4)
xU(2,3)=stretchNew(i,5)
xU(3,1)=stretchNew(i,6)
xU(2,1)=xU(1,2)
xU(3,2)=xU(2,3)
xU(1,3)=xU(3,1)

c Inverse of Stretch

CALL CalDet(xU,xJU)
CALL CalInv(xU,xJU,U_1)

c We need to calculate the rotational matrix: rot=FU^-1
c
c Rotation Matrix

rot=matmul(F_new,U_1)

c Right cauchy stress tensor

C_new=matmul(TRANSPOSE(F_new),F_new)

c Inizial fiber direction

a0_old(1,1)=one
a0_old(1,2)=zero
a0_old(1,3)=zero

b0_old(1,1)=zero
b0_old(1,2)=one
b0_old(1,3)=zero

a0=matmul(a0_old,F_old)
b0=matmul(b0_old,F_old)

c Calcolo invariante I4

aI4=sum(matmul(a0,C_new)*a0)

c Calcolo invariante I6

aI6=sum(matmul(a0,C_new)*b0)

c Calcolo invariante I8

aI8=sum(matmul(b0,C_new)*b0)

c Calcolo invariante I10

aI10=acos(sum(a0*b0))-acos(aI6/sqrt(aI4*aI8))

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

WTa_bar=a_a*(SQRT(aI4)-1)**three+b_a*(SQRT(aI4)-one)**two

WTb_bar=a_b*(SQRT(aI8)-one)**three+b_b*(SQRT(aI8)-one)**two

kTC_a=(one-pmin_a)*exp(c_a*(SQRT(aI8)-one))+pmin_a

kTC_b=(one-pmin_b)*exp(c_b*(SQRT(aI4)-one))+pmin_b

WTa=WTa_bar*kTC_a

WCa=d_a*(SQRT(aI4)-one)**two

WTb=WTb_bar*kTC_b

WCb=d_b*(SQRT(aI8)-one)**two

WS=half*s1*aI10**two+third*s2*aI10**two+fourth*s3*aI10**four

WSTa=third*aI10**three*sT_a*(aI4-one)**two

WSTb=third*aI10**three*sT_b*(aI8-one)**two

kSCa=exp(sC_a*(aI4-one)**three)

kSCb=exp(sC_b*(aI8-one)**three)

C Definizione vettore vectW

vectW(1,1)=WTa
vectW(1,2)=WTb
vectW(1,3)=WS
vectW(1,4)=WSTa
vectW(1,5)=WSTb
vectW(1,6)=WCb
vectW(1,7)=WS*kSCb
vectW(1,8)=WSTa*kSCb
vectW(1,9)=WCa
vectW(1,10)=WS*kSCa
vectW(1,11)=WSTb*kSCa
vectW(1,12)=WS*kSCa*kSCb

c Definizione vettori vectI4_full:

C a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b

vectI4_full(1,1)=(a_a*(threehalfs*SQRT(aI4)-
3 three+threehalfs*one/SQRT(aI4))+b_a*(one-
1 one/SQRT(aI4)))*((one-pmin_a)*exp(c_a*(SQRT(aI8)-one))-one)

vectI4_full(1,2)=(a_b*(SQRT(aI8)-1)**three+
1 b_b*(SQRT(aI8)-one)**two)*(one-
3 pmin_b)*(c_b*(half*one/SQRT(aI4)))*exp(c_b*(SQRT(aI4)-one))

vectI4_full(1,3)=(twothirds*aI10**three*sT_a*(aI4-one))*(one+kSCb)

vectI4_full(1,4)=(one-one/SQRT(aI4))*d_a

vectI4_full(1,5)=(WS+WSTb+WS*kSCb)*(three*aI4**one-six*aI4+
1 three)*exp(sC_a*(aI4**three-three*aI4**two+three*aI4-one))

c Definizione vettori vectI8_full:

C a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b

vectI8_full(1,1)=(a_a*(SQRT(aI4)-one)**three+b_a*(SQRT(aI4)-
1 one)**two)*(one-pmin_a)*(half*c_a/SQRT(aI8))*exp(c_b*(SQRT(aI8)-
2 one))

vectI8_full(1,2)=(a_b*(threehalfs*SQRT(aI8)-
2 three+three/(two*SQRT(aI8)))+b_b*(one-
1 one/SQRT(aI8)))*((one-pmin_b)*exp(c_b*(SQRT(aI4)-one)+pmin_b))

vectI8_full(1,3)=(twothirds*aI10**three*sT_b*(aI8-one))*(one+kSCa)

vectI8_full(1,4)=d_b*(one-one/SQRT(aI8))

vectI8_full(1,5)=(WS+WSTa+WS*kSCa)*((three*aI8**two-six*aI8+
2 three)*exp(sC_b*(aI8**three-three*aI8**two+three*aI4-one)))

c Definizione vettori vectI10_full:

C a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b

vectI10_full(1,1)=(s1*aI10+s2*aI10**two+s3*aI10**three)*(one+kSCb+
1 kSCa+kSCa*kSCb)

vectI10_full(1,2)=(one+kSCb)*(aI10**two*sT_a*(SQRT(aI4)-one)**two)

vectI10_full(1,3)=(one+kSCa)*(aI10**two*sT_b*(SQRT(aI8)-one)**two)

C -----------------------
C IFCICLE W_tot=vectW(1,12)*vect_Wtot(1,12)
c
c vect_Wtot(,,,,,,,,,,,)
c W_tot=matmul(vectW(1,12)*vect_Wtot(1,12))
c
c vectI4=(,,,,)
c vectI8=(,,,,)
c vectI10=(,,)

c DWDI4=matmul(vectI4_full,transpose(vectI4))
c DWDI8=matmul(vectI8_full,transpose(vectI8))
c DWDI10=matmul(vectI10_full,transpose(vectI10))

C IF aI4.GT.one.AND.aI8.GT.one

IF(aI4.GT.one) THEN
IF(aI8.GT.one) THEN

c vect_Wtot=(0,1,1,1,1,0,0,0,0,0,0,0)
vect_Wtot(1,1)=zero
vect_Wtot(1,2)=one
vect_Wtot(1,3)=one
vect_Wtot(1,4)=one
vect_Wtot(1,5)=one
vect_Wtot(1,6)=zero
vect_Wtot(1,7)=zero
vect_Wtot(1,8)=zero
vect_Wtot(1,9)=zero
vect_Wtot(1,10)=zero
vect_Wtot(1,11)=zero
vect_Wtot(1,12)=zero

W_tot=sum(vectW*vect_Wtot)

c vectI4=(1,1,1,0,0)

vectI4(1,1)=one
vectI4(1,1)=one
vectI4(1,1)=one
vectI4(1,1)=zero
vectI4(1,1)=zero

c vectI8=(1,1,1,0,0)

vectI8(1,1)=one
vectI8(1,1)=one
vectI8(1,1)=one
vectI8(1,1)=zero
vectI8(1,1)=zero

c vectI10=(1,1,1)

vectI10(1,1)=one
vectI10(1,1)=one
vectI10(1,1)=one

DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8
DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10

END IF
END IF

IF(aI4.GT.one) THEN
IF(aI8.LT.one.OR.aI8.EQ.one) THEN

c vect_Wtot(1,0,0,0,0,1,1,1,0,0,0,0)

vect_Wtot(1,1)=one
vect_Wtot(1,2)=zero
vect_Wtot(1,3)=zero
vect_Wtot(1,4)=zero
vect_Wtot(1,5)=zero
vect_Wtot(1,6)=one
vect_Wtot(1,7)=one
vect_Wtot(1,8)=one
vect_Wtot(1,9)=zero
vect_Wtot(1,10)=zero
vect_Wtot(1,11)=zero
vect_Wtot(1,12)=zero

W_tot=sum(vectW*vect_Wtot)

c vectI4=(1,0,0,1,1)

vectI4(1,1)=one
vectI4(1,1)=zero
vectI4(1,1)=zero
vectI4(1,1)=one
vectI4(1,1)=one

c vectI8=(1,0,0,1,1)

vectI8(1,1)=one
vectI8(1,1)=zero
vectI8(1,1)=zero
vectI8(1,1)=one
vectI8(1,1)=one

c vectI10=(1,1,0)

vectI10(1,1)=one
vectI10(1,1)=one
vectI10(1,1)=zero

DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8
DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10

END IF
END IF

C IF aI4 LT and EQ one and one and aI8 GT one

IF(aI4.LT.one.OR.aI4.EQ.one) THEN

IF(aI8.GT.one) THEN

c vect_Wtot(0,1,0,0,0,0,0,0,1,1,1,0)

vect_Wtot(1,1)=zero
vect_Wtot(1,2)=one
vect_Wtot(1,3)=zero
vect_Wtot(1,4)=zero
vect_Wtot(1,5)=zero
vect_Wtot(1,6)=zero
vect_Wtot(1,7)=zero
vect_Wtot(1,8)=zero
vect_Wtot(1,9)=one
vect_Wtot(1,10)=one
vect_Wtot(1,11)=one
vect_Wtot(1,12)=zero

W_tot=sum(vectW*vect_Wtot)

c vectI4=(0,1,1,0,0)

vectI4(1,1)=zero
vectI4(1,1)=one
vectI4(1,1)=one
vectI4(1,1)=zero
vectI4(1,1)=zero

c vectI8=(0,1,0,1,1)

vectI8(1,1)=zero
vectI8(1,1)=one
vectI8(1,1)=zero
vectI8(1,1)=one
vectI8(1,1)=one

c vectI10=(1,0,1)

vectI10(1,1)=one
vectI10(1,1)=zero
vectI10(1,1)=one

DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8
DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10

END IF
END IF

c IF aI4 LT one and EQ one and aI8 LT one and EQ one

IF(aI4.LT.one.OR.aI4.EQ.one) THEN
IF(aI8.LT.one.OR.aI8.EQ.one) THEN

c vect_Wtot=(0,0,0,0,0,1,0,0,1,0,0,1)

vect_Wtot(1,1)=zero
vect_Wtot(1,2)=zero
vect_Wtot(1,3)=zero
vect_Wtot(1,4)=zero
vect_Wtot(1,5)=zero
vect_Wtot(1,6)=one
vect_Wtot(1,7)=zero
vect_Wtot(1,8)=zero
vect_Wtot(1,9)=one
vect_Wtot(1,10)=zero
vect_Wtot(1,11)=zero
vect_Wtot(1,12)=one

W_tot=sum(vectW*vect_Wtot)

c vectI4=(0,0,0,1,1)

vectI4(1,1)=zero
vectI4(1,1)=zero
vectI4(1,1)=zero
vectI4(1,1)=one
vectI4(1,1)=one

c vectI8=(0,0,0,1,1)

vectI8(1,1)=zero
vectI8(1,1)=zero
vectI8(1,1)=zero
vectI8(1,1)=one
vectI8(1,1)=one

c vectI10=(1,0,0)

vectI10(1,1)=one
vectI10(1,1)=zero
vectI10(1,1)=zero

DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4
DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8
DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10

END IF
END IF

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

c tensor products

CALL KronProd(a0,a0,a0xa0)

CALL KronProd(b0,b0,b0xb0)

CALL KronProd(a0,b0,a0xb0)

CALL KronProd(b0,b0,b0xa0)

aS1=a0xa0*DW(1,1)+b0xb0*DW(1,2)
aS2=DW(1,3)*(one/SQRT(aI4*aI8-aI6**two))
aS3=aI6*(a0xa0/aI4+b0xb0/aI8)

S=two*aS1+aS2*(a0xb0 + b0xa0 - aS3)

sigma_co=one/aJ*matmul(F_new,matmul(S,transpose(F_new)))

sigma_new=matmul(matmul(transpose(rot),sigma_co),rot)

stressNew(i,1)=sigma_new(1,1)
stressNew(i,2)=sigma_new(2,2)
stressNew(i,3)=sigma_new(3,3)
stressNew(i,4)=sigma_new(1,2)
stressNew(i,5)=sigma_new(2,3)
stressNew(i,6)=sigma_new(3,1)

c stateNew(i,1)=W_tot/rho

c deallocate(C)

return
end do

end

c calcolo del determinante

subroutine CalDet(A,aJ)

implicit none
REAL*8 A(3,3), aJ

aJ = A(1,1) * A(2,2) * A(3,3) -
1 A(1,2) * A(2,1) * A(3,3) +
2 + A(1,2) * A(2,3) * A(3,1) +
1 A(1,3) * A(3,2) * A(2,1) -
2 A(1,3) * A(3,1) * A(2,2) -
3 A(2,3) * A(3,2) * A(1,1)

return
end

c calcolo della matrice inversa

subroutine CalInv(A,aJ,B)
include 'vaba_param.inc'
dimension A(3,3),B(3,3)

B(1,1)=A(2,2)*A(3,3)-A(2,3)*A(3,2)
B(1,2)=A(1,3)*A(3,2)-A(1,2)*A(3,3)
B(1,3)=A(1,2)*A(2,3)-A(1,3)*A(2,2)
B(2,1)=A(2,3)*A(3,1)-A(2,1)*A(3,3)
B(2,2)=A(1,1)*A(3,3)-A(1,3)*A(3,1)
B(2,3)=A(1,3)*A(2,1)-A(1,1)*A(2,3)
B(3,1)=A(2,1)*A(3,2)-A(2,2)*A(3,1)
B(3,2)=A(1,2)*A(3,1)-A(1,1)*A(3,2)
B(3,3)=A(1,1)*A(2,2)-A(1,2)*A(2,1)

B=B/aJ

return
end

c calcolo prodotto tensioriale

function KronProd(A,B) result(C)
IMPLICIT NONE
real, dimension (:,:), intent(in) :: A, B
real, dimension (:,:), allocatable :: C
integer :: i = 0, j = 0, k = 0, l = 0
integer :: m = 0, n = 0, p = 0, q = 0

allocate(C(size(A,1)*size(B,1),size(A,2)*size(B,2)))
C = 0

do i = 1,size(A,1)
do j = 1,size(A,2)
n=(i-1)*size(B,1) + 1
m=n+size(B,1) - 1
p=(j-1)*size(B,2) + 1
q=p+size(B,2) - 1
C(n:m,p:q) = A(i,j)*B
enddo
enddo

end function KronProd

Re: Abaqus VUMAT error during compilation

Domv — Wed, 17 Jan 2024 18:45:31 GMT

subroutine vumat( C Read only (unmodifiable)variables - 1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, 2 stepTime, totalTime, dt, cmname, coordMp, charLength, 3 props, density, strainInc, relSpinInc, 4 tempOld, stretchOld, defgradOld, fieldOld, 5 stressOld, stateOld, enerInternOld, enerInelasOld, 6 tempNew, stretchNew, defgradNew, fieldNew, C Write only (modifiable) variables - 7 stressNew, stateNew, enerInternNew, enerInelasNew ) C include 'vaba_param.inc' C dimension props(nprops), density(nblock), coordMp(nblock,*), 1 charLength(nblock), strainInc(nblock,ndir+nshr), 2 relSpinInc(nblock,nshr), tempOld(nblock), 3 stretchOld(nblock,ndir+nshr), 4 defgradOld(nblock,ndir+nshr+nshr), 5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr), 6 stateOld(nblock,nstatev), enerInternOld(nblock), 7 enerInelasOld(nblock), tempNew(nblock), 8 stretchNew(nblock,ndir+nshr), 8 defgradNew(nblock,ndir+nshr+nshr), 9 fieldNew(nblock,nfieldv), 1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), 2 enerInternNew(nblock), enerInelasNew(nblock) C dimension F_new(3,3), xU(3,3), U_1(3,3), rot(3,3), 3 C_new(3,3), a0(1,3), b0(1,3), a0_old(1,3), b0_old(1,3), 4 F_old(3,3), DWDI4(1,5), DWDI8(1,5), DWDI10(1,5), 5 vectI4(1,5),vectI8(1,5),vectI10(1,5), 6 vectI4_full(1,5),vectI8_full(1,5),vectI10_full(1,5), 7 S(3,3), sigma_co(3,3), sigma_new(3,3), vectW(1,12), 8 vect_Wtot(1,12), D(3,3), a0xa0(3,3), 3 b0xb0(3,3), a0xb0(3,3), b0xa0(3,3), aS1(3,3), 4 aS3(3,3), DW(1,3) c real*8 a_a, a_b, b_a, b_b, d_a, d_b, pmin_a, pmin_b, cc + c_a, c_b, s1, s2, s3, sT_a, sT_b, sC_a, sC_b, rho, xJU, cc + aI4, aI6, aI8, aI10, WTa_bar, WTb_bar, kTC_a, kTC_b, c + WTa, WCa, WTb, WCb, WS, WSTa, WSTb, kSCa, kSCb, W_tot parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0, 3 four=4.d0, 2 third=one/three, half= .5d0, twoThirds=two/three, 3 threeHalfs = 1.5d0, fourth=one/four, six=6.d0) character*80 cmname C do i = 1,nblock C Read material properties a_a=props(1) a_b=a_a b_a=props(2) b_b=b_a d_a=props(3) d_b=d_a pmin_a=props(4) pmin_b=pmin_a c_a=props(5) c_b=c_a s1=props(6) s2=props(7) s3=props(8) sT_a=props(9) sT_b=sT_a sC_a=props(10) sC_b=sC_a c ---------- c Deformation gradiEnt before the increment F_old(1,1)=defgradOld(i,1) F_old(1,2)=defgradOld(i,4) F_old(1,3)=defgradOld(i,9) F_old(2,1)=defgradOld(i,7) F_old(2,2)=defgradOld(i,2) F_old(2,3)=defgradOld(i,5) F_old(3,1)=defgradOld(i,6) F_old(3,2)=defgradOld(i,8) F_old(3,3)=defgradOld(i,3) c Deformation gradiant after the increment F_new(1,1)=defgradNew(i,1) F_new(1,2)=defgradNew(i,4) F_new(1,3)=defgradNew(i,9) F_new(2,1)=defgradNew(i,7) F_new(2,2)=defgradNew(i,2) F_new(2,3)=defgradNew(i,5) F_new(3,1)=defgradNew(i,6) F_new(3,2)=defgradNew(i,8) F_new(3,3)=defgradNew(i,3) c Stretch tensor after the increment xU(1,1)=stretchNew(i,1) xU(2,2)=stretchNew(i,2) xU(3,3)=stretchNew(i,3) xU(1,2)=stretchNew(i,4) xU(2,3)=stretchNew(i,5) xU(3,1)=stretchNew(i,6) xU(2,1)=xU(1,2) xU(3,2)=xU(2,3) xU(1,3)=xU(3,1) c Inverse of Stretch CALL CalDet(xU,xJU) CALL CalInv(xU,xJU,U_1) c We need to calculate the rotational matrix: rot=FU^-1 c c Rotation Matrix rot=matmul(F_new,U_1) c Right cauchy stress tensor C_new=matmul(TRANSPOSE(F_new),F_new) c Inizial fiber direction a0_old(1,1)=one a0_old(1,2)=zero a0_old(1,3)=zero b0_old(1,1)=zero b0_old(1,2)=one b0_old(1,3)=zero a0=matmul(a0_old,F_old) b0=matmul(b0_old,F_old) c Calcolo invariante I4 aI4=sum(matmul(a0,C_new)*a0) c Calcolo invariante I6 aI6=sum(matmul(a0,C_new)*b0) c Calcolo invariante I8 aI8=sum(matmul(b0,C_new)*b0) c Calcolo invariante I10 aI10=acos(sum(a0*b0))-acos(aI6/sqrt(aI4*aI8)) C ----------------------- WTa_bar=a_a*(SQRT(aI4)-1)**three+b_a*(SQRT(aI4)-one)**two WTb_bar=a_b*(SQRT(aI8)-one)**three+b_b*(SQRT(aI8)-one)**two kTC_a=(one-pmin_a)*exp(c_a*(SQRT(aI8)-one))+pmin_a kTC_b=(one-pmin_b)*exp(c_b*(SQRT(aI4)-one))+pmin_b WTa=WTa_bar*kTC_a WCa=d_a*(SQRT(aI4)-one)**two WTb=WTb_bar*kTC_b WCb=d_b*(SQRT(aI8)-one)**two WS=half*s1*aI10**two+third*s2*aI10**two+fourth*s3*aI10**four WSTa=third*aI10**three*sT_a*(aI4-one)**two WSTb=third*aI10**three*sT_b*(aI8-one)**two kSCa=exp(sC_a*(aI4-one)**three) kSCb=exp(sC_b*(aI8-one)**three) C Definizione vettore vectW vectW(1,1)=WTa vectW(1,2)=WTb vectW(1,3)=WS vectW(1,4)=WSTa vectW(1,5)=WSTb vectW(1,6)=WCb vectW(1,7)=WS*kSCb vectW(1,8)=WSTa*kSCb vectW(1,9)=WCa vectW(1,10)=WS*kSCa vectW(1,11)=WSTb*kSCa vectW(1,12)=WS*kSCa*kSCb c Definizione vettori vectI4_full: C a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b vectI4_full(1,1)=(a_a*(threehalfs*SQRT(aI4)- 3 three+threehalfs*one/SQRT(aI4))+b_a*(one- 1 one/SQRT(aI4)))*((one-pmin_a)*exp(c_a*(SQRT(aI8)-one))-one) vectI4_full(1,2)=(a_b*(SQRT(aI8)-1)**three+ 1 b_b*(SQRT(aI8)-one)**two)*(one- 3 pmin_b)*(c_b*(half*one/SQRT(aI4)))*exp(c_b*(SQRT(aI4)-one)) vectI4_full(1,3)=(twothirds*aI10**three*sT_a*(aI4-one))*(one+kSCb) vectI4_full(1,4)=(one-one/SQRT(aI4))*d_a vectI4_full(1,5)=(WS+WSTb+WS*kSCb)*(three*aI4**one-six*aI4+ 1 three)*exp(sC_a*(aI4**three-three*aI4**two+three*aI4-one)) c Definizione vettori vectI8_full: C a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b vectI8_full(1,1)=(a_a*(SQRT(aI4)-one)**three+b_a*(SQRT(aI4)- 1 one)**two)*(one-pmin_a)*(half*c_a/SQRT(aI8))*exp(c_b*(SQRT(aI8)- 2 one)) vectI8_full(1,2)=(a_b*(threehalfs*SQRT(aI8)- 2 three+three/(two*SQRT(aI8)))+b_b*(one- 1 one/SQRT(aI8)))*((one-pmin_b)*exp(c_b*(SQRT(aI4)-one)+pmin_b)) vectI8_full(1,3)=(twothirds*aI10**three*sT_b*(aI8-one))*(one+kSCa) vectI8_full(1,4)=d_b*(one-one/SQRT(aI8)) vectI8_full(1,5)=(WS+WSTa+WS*kSCa)*((three*aI8**two-six*aI8+ 2 three)*exp(sC_b*(aI8**three-three*aI8**two+three*aI4-one))) c Definizione vettori vectI10_full: C a_a a_b b_a b_b d_a d_b pmin_a pmin_b c_a c_b s1 s2 s3 sT_a sT_b sC_a sC_b vectI10_full(1,1)=(s1*aI10+s2*aI10**two+s3*aI10**three)*(one+kSCb+ 1 kSCa+kSCa*kSCb) vectI10_full(1,2)=(one+kSCb)*(aI10**two*sT_a*(SQRT(aI4)-one)**two) vectI10_full(1,3)=(one+kSCa)*(aI10**two*sT_b*(SQRT(aI8)-one)**two) C ----------------------- C IFCICLE W_tot=vectW(1,12)*vect_Wtot(1,12) c c vect_Wtot(,,,,,,,,,,,) c W_tot=matmul(vectW(1,12)*vect_Wtot(1,12)) c c vectI4=(,,,,) c vectI8=(,,,,) c vectI10=(,,) c DWDI4=matmul(vectI4_full,transpose(vectI4)) c DWDI8=matmul(vectI8_full,transpose(vectI8)) c DWDI10=matmul(vectI10_full,transpose(vectI10)) C IF aI4.GT.one.AND.aI8.GT.one IF(aI4.GT.one) THEN IF(aI8.GT.one) THEN c vect_Wtot=(0,1,1,1,1,0,0,0,0,0,0,0) vect_Wtot(1,1)=zero vect_Wtot(1,2)=one vect_Wtot(1,3)=one vect_Wtot(1,4)=one vect_Wtot(1,5)=one vect_Wtot(1,6)=zero vect_Wtot(1,7)=zero vect_Wtot(1,8)=zero vect_Wtot(1,9)=zero vect_Wtot(1,10)=zero vect_Wtot(1,11)=zero vect_Wtot(1,12)=zero W_tot=sum(vectW*vect_Wtot) c vectI4=(1,1,1,0,0) vectI4(1,1)=one vectI4(1,1)=one vectI4(1,1)=one vectI4(1,1)=zero vectI4(1,1)=zero c vectI8=(1,1,1,0,0) vectI8(1,1)=one vectI8(1,1)=one vectI8(1,1)=one vectI8(1,1)=zero vectI8(1,1)=zero c vectI10=(1,1,1) vectI10(1,1)=one vectI10(1,1)=one vectI10(1,1)=one DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4 DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8 DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10 END IF END IF IF(aI4.GT.one) THEN IF(aI8.LT.one.OR.aI8.EQ.one) THEN c vect_Wtot(1,0,0,0,0,1,1,1,0,0,0,0) vect_Wtot(1,1)=one vect_Wtot(1,2)=zero vect_Wtot(1,3)=zero vect_Wtot(1,4)=zero vect_Wtot(1,5)=zero vect_Wtot(1,6)=one vect_Wtot(1,7)=one vect_Wtot(1,8)=one vect_Wtot(1,9)=zero vect_Wtot(1,10)=zero vect_Wtot(1,11)=zero vect_Wtot(1,12)=zero W_tot=sum(vectW*vect_Wtot) c vectI4=(1,0,0,1,1) vectI4(1,1)=one vectI4(1,1)=zero vectI4(1,1)=zero vectI4(1,1)=one vectI4(1,1)=one c vectI8=(1,0,0,1,1) vectI8(1,1)=one vectI8(1,1)=zero vectI8(1,1)=zero vectI8(1,1)=one vectI8(1,1)=one c vectI10=(1,1,0) vectI10(1,1)=one vectI10(1,1)=one vectI10(1,1)=zero DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4 DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8 DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10 END IF END IF C IF aI4 LT and EQ one and one and aI8 GT one IF(aI4.LT.one.OR.aI4.EQ.one) THEN IF(aI8.GT.one) THEN c vect_Wtot(0,1,0,0,0,0,0,0,1,1,1,0) vect_Wtot(1,1)=zero vect_Wtot(1,2)=one vect_Wtot(1,3)=zero vect_Wtot(1,4)=zero vect_Wtot(1,5)=zero vect_Wtot(1,6)=zero vect_Wtot(1,7)=zero vect_Wtot(1,8)=zero vect_Wtot(1,9)=one vect_Wtot(1,10)=one vect_Wtot(1,11)=one vect_Wtot(1,12)=zero W_tot=sum(vectW*vect_Wtot) c vectI4=(0,1,1,0,0) vectI4(1,1)=zero vectI4(1,1)=one vectI4(1,1)=one vectI4(1,1)=zero vectI4(1,1)=zero c vectI8=(0,1,0,1,1) vectI8(1,1)=zero vectI8(1,1)=one vectI8(1,1)=zero vectI8(1,1)=one vectI8(1,1)=one c vectI10=(1,0,1) vectI10(1,1)=one vectI10(1,1)=zero vectI10(1,1)=one DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4 DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8 DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10 END IF END IF c IF aI4 LT one and EQ one and aI8 LT one and EQ one IF(aI4.LT.one.OR.aI4.EQ.one) THEN IF(aI8.LT.one.OR.aI8.EQ.one) THEN c vect_Wtot=(0,0,0,0,0,1,0,0,1,0,0,1) vect_Wtot(1,1)=zero vect_Wtot(1,2)=zero vect_Wtot(1,3)=zero vect_Wtot(1,4)=zero vect_Wtot(1,5)=zero vect_Wtot(1,6)=one vect_Wtot(1,7)=zero vect_Wtot(1,8)=zero vect_Wtot(1,9)=one vect_Wtot(1,10)=zero vect_Wtot(1,11)=zero vect_Wtot(1,12)=one W_tot=sum(vectW*vect_Wtot) c vectI4=(0,0,0,1,1) vectI4(1,1)=zero vectI4(1,1)=zero vectI4(1,1)=zero vectI4(1,1)=one vectI4(1,1)=one c vectI8=(0,0,0,1,1) vectI8(1,1)=zero vectI8(1,1)=zero vectI8(1,1)=zero vectI8(1,1)=one vectI8(1,1)=one c vectI10=(1,0,0) vectI10(1,1)=one vectI10(1,1)=zero vectI10(1,1)=zero DW(1,1)=sum(vectI4_full*vectI4) ! DWDI4 DW(1,2)=sum(vectI8_full*vectI8) ! DWDI8 DW(1,3)=sum(vectI10_full*vectI10) ! DWDI10 END IF END IF C ----------------------- c tensor products CALL KronProd(a0,a0,a0xa0) CALL KronProd(b0,b0,b0xb0) CALL KronProd(a0,b0,a0xb0) CALL KronProd(b0,b0,b0xa0) aS1=a0xa0*DW(1,1)+b0xb0*DW(1,2) aS2=DW(1,3)*(one/SQRT(aI4*aI8-aI6**two)) aS3=aI6*(a0xa0/aI4+b0xb0/aI8) S=two*aS1+aS2*(a0xb0 + b0xa0 - aS3) sigma_co=one/aJ*matmul(F_new,matmul(S,transpose(F_new))) sigma_new=matmul(matmul(transpose(rot),sigma_co),rot) stressNew(i,1)=sigma_new(1,1) stressNew(i,2)=sigma_new(2,2) stressNew(i,3)=sigma_new(3,3) stressNew(i,4)=sigma_new(1,2) stressNew(i,5)=sigma_new(2,3) stressNew(i,6)=sigma_new(3,1) c stateNew(i,1)=W_tot/rho c deallocate(C) return end do end c calcolo del determinante subroutine CalDet(A,aJ) implicit none REAL*8 A(3,3), aJ aJ = A(1,1) * A(2,2) * A(3,3) - 1 A(1,2) * A(2,1) * A(3,3) + 2 + A(1,2) * A(2,3) * A(3,1) + 1 A(1,3) * A(3,2) * A(2,1) - 2 A(1,3) * A(3,1) * A(2,2) - 3 A(2,3) * A(3,2) * A(1,1) return end c calcolo della matrice inversa subroutine CalInv(A,aJ,B) include 'vaba_param.inc' dimension A(3,3),B(3,3) B(1,1)=A(2,2)*A(3,3)-A(2,3)*A(3,2) B(1,2)=A(1,3)*A(3,2)-A(1,2)*A(3,3) B(1,3)=A(1,2)*A(2,3)-A(1,3)*A(2,2) B(2,1)=A(2,3)*A(3,1)-A(2,1)*A(3,3) B(2,2)=A(1,1)*A(3,3)-A(1,3)*A(3,1) B(2,3)=A(1,3)*A(2,1)-A(1,1)*A(2,3) B(3,1)=A(2,1)*A(3,2)-A(2,2)*A(3,1) B(3,2)=A(1,2)*A(3,1)-A(1,1)*A(3,2) B(3,3)=A(1,1)*A(2,2)-A(1,2)*A(2,1) B=B/aJ return end c calcolo prodotto tensioriale function KronProd(A,B) result(C) IMPLICIT NONE real, dimension (:,:), intent(in) :: A, B real, dimension (:,:), allocatable :: C integer :: i = 0, j = 0, k = 0, l = 0 integer :: m = 0, n = 0, p = 0, q = 0 allocate(C(size(A,1)*size(B,1),size(A,2)*size(B,2))) C = 0 do i = 1,size(A,1) do j = 1,size(A,2) n=(i-1)*size(B,1) + 1 m=n+size(B,1) - 1 p=(j-1)*size(B,2) + 1 q=p+size(B,2) - 1 C(n:m,p:q) = A(i,j)*B enddo enddo end function KronProd

Re: Abaqus VUMAT error during compilation

andrew_4619 — Wed, 17 Jan 2024 22:59:08 GMT

Well it appears that the message "forrtl: severe (151): allocatable array is already allocated" is not generated by the code you pasted, though it maybe your code sets things that cause other parts of the ABAQUS solver to crash. I see many bugs... For example

c vectI4=(1,1,1,0,0) vectI4(1,1)=one vectI4(1,1)=one vectI4(1,1)=one vectI4(1,1)=zero vectI4(1,1)=zero c vectI8=(1,1,1,0,0) vectI8(1,1)=one vectI8(1,1)=one vectI8(1,1)=one vectI8(1,1)=zero vectI8(1,1)=zero c vectI10=(1,1,1) vectI10(1,1)=one vectI10(1,1)=one vectI10(1,1)=one

you are setting vectI4(1,1) five times as this array is declared (1,5) you really want:

c vectI4=(1,1,1,0,0) vectI4(1,1)=one vectI4(1,2)=one vectI4(1,3)=one vectI4(1,4)=zero vectI4(1,5)=zero ....

but is would make more sense to have

vectI4(1,:) = [ one, one, one, zero, zero ] vectI8(1,:) = [ one, one, one, zero, zero ] vectI10(1,1:3) = one c note vectI10(1,4) and (1,5) are not initialised is that intentional

Fortran has syntax for whole array and array slice operations.

Re: Abaqus VUMAT error during compilation

Domv — Wed, 17 Jan 2024 23:21:14 GMT

I really don't know how i made this mistake, i've already correct and still receiving the same error, with abaqus message:

ERROR in job messaging system: Error in connection to analysis.

Error in job Job-1: Abaqus/Explicit Packager exited with an error - Please see the status file for possible error messages if the file exists.
Job Job-1 aborted due to errors.

Re: Abaqus VUMAT error during compilation

andrew_4619 — Wed, 17 Jan 2024 23:40:31 GMT

There are I suspect more errors in the code I only had a quick look at it. You can get some Fortran help here but not much help with ABAQUS. Did you start by having a working sample VUMAT that ran OK? That should be the first test.