Skip to content
Snippets Groups Projects

[Refractor] - Numerical Tangent as an exteral function call

Merged [Refractor] - Numerical Tangent as an exteral function call
fftmad_memory into main
Merged Noels Ludovic requested to merge fftmad_memory into main
Viewing commit f8cf6170
Show latest version
1 file
+ 172
− 48
Compare changes
  • Side-by-side
  • Inline
MaterialModels/FiniteStrain/Finite_VEVP/umat.f
+ 172
− 48
@@ -33,7 +33,7 @@ C !--------------------------------------------------------------
C ! UMAT SUBROUTINE
C !--------------------------------------------------------------
SUBROUTINE umat_vevp(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,
3 DPRED,CMNAME,
@@ -42,24 +42,33 @@ C !--------------------------------------------------------------
6 CELENT,DFGRD0,DFGRD1,NOEL,NPT,
7 LAYER,KSPT,KSTEP,KINC)
C INCLUDE 'ABA_PARAM.INC'
INCLUDE 'ABA_PARAM.INC'
CHARACTER*80 CMNAME
REAL*8 STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),DFGRD0(3,3), DFGRD1(3,3)
3 PROPS(NPROPS), COORDS(3),DROT(3,3),DFGRD0(3,3),
4 DFGRD1(3,3)
REAL*8 SSE, SPD, SCD, RPL, DTIME, TEMP, CELENT, DRPLDT
REAL*8 PNEWDT, DTEMP, R, dR, ddR
IF (PROPS(1) .GT. 2.D0) THEN
call VEVP(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,PROPS,
call VEVP(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,PROPS,
& NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1,
& PNEWDT)
return
end
return
ELSE
call NEOHOOK(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,PROPS,
& NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1)
return
END IF
end
C !--------------------------------------------------------------
@@ -83,7 +92,7 @@ C !--------------------------------------------------------------
REAL(prec), INTENT(IN) :: DFGRD0(3,3), DFGRD1(3,3)
REAL(prec), INTENT(INOUT) :: statev(nstatv)
REAL(prec), INTENT(OUT) :: stress(ntens)
REAL(prec), INTENT(OUT) :: ddsdde(ntens,ntens), PNEWDT
REAL(prec), INTENT(OUT) :: ddsdde(ntens,ntens), PNEWDT
!Declaration of local variables
INTEGER :: ii, jj, O6, order
@@ -101,7 +110,8 @@ C !--------------------------------------------------------------
! VE trial variables
REAL(prec) :: F_ve_tr(3, 3), C_ve_tr(3, 3), D_ve_tr(3, 3)
REAL(prec) :: E_ve_tr(3, 3), dlogC_ve_tr(3,3,3,3)
REAL(prec) :: ddlogC_ve_tr(3,3,3,3,3,3), AA_tr(8,3,3)
!REAL(prec) :: ddlogC_ve_tr(3,3,3,3,3,3)
REAL(prec) :: AA_tr(8,3,3)
REAL(prec) :: BB_tr(8),GGe, KKe, kappa_tr(3, 3), S_tr(3, 3)
REAL(prec) :: temp1(3,3), tau_tr(3, 3)
! Variables for trial yield
@@ -111,7 +121,7 @@ C !--------------------------------------------------------------
REAL(prec) :: F_hat(6, 3, 3), J, F_ve_tr_hat(6, 3, 3)
REAL(prec) :: C_ve_tr_hat(6, 3, 3), D_ve_tr_hat(6,3,3)
REAL(prec) :: E_ve_tr_hat(6,3,3), dlogC_ve_tr_hat(6,3,3,3,3)
REAL(prec) :: ddlogC_ve_tr_hat(6,3,3,3,3,3,3)
!REAL(prec) :: ddlogC_ve_tr_hat(6,3,3,3,3,3,3)
REAL(prec) :: AA_tr_hat(6,8,3,3)
REAL(prec) :: BB_tr_hat(6, 8), GGe_hat(6), KKe_hat(6)
REAL(prec) :: kappa_tr_hat(6, 3, 3), S_tr_hat(6, 3, 3)
@@ -124,7 +134,8 @@ C !--------------------------------------------------------------
REAL(prec) :: exp_GQ(3, 3), F_vp(3,3), F_vp_inv(3,3)
! VE corrected variables
REAL(prec) :: F_ve(3, 3), C_ve(3, 3), D_ve(3, 3), E_ve(3, 3)
REAL(prec) :: dlogC_ve(3,3,3,3), ddlogC_ve(3,3,3,3,3,3)
REAL(prec) :: dlogC_ve(3,3,3,3)
!REAL(prec) :: ddlogC_ve(3,3,3,3,3,3)
REAL(prec) :: AA(8,3,3), BB(8), kappa(3, 3), b(3, 3)
REAL(prec) :: S(3, 3), temp2(3, 3), tau(3, 3)
! VP variables perturbated
@@ -143,14 +154,15 @@ C !--------------------------------------------------------------
REAL(prec) :: F_ve_hat(6,3,3), C_ve_hat(6,3,3)
REAL(prec) :: D_ve_hat(6,3,3), E_ve_hat(6,3,3)
REAL(prec) :: dlogC_ve_hat(6,3,3,3,3)
REAL(prec) :: ddlogC_ve_hat(6,3,3,3,3,3,3), AA_hat(6,8,3,3)
REAL(prec) :: AA_hat(6,8,3,3)
!REAL(prec) :: ddlogC_ve_hat(6,3,3,3,3,3,3)
REAL(prec) :: BB_hat(6,8), kappa_hat(6,3,3), S_hat(6,3,3)
REAL(prec) :: temp2_hat(6,3,3), tau_hat(6,3,3)
REAL(prec) :: tau_hat_v(6,6)
! Tollerances
REAL(prec), PARAMETER :: TOLL=0.001D0, TOLL_G=0.999D-11
INTEGER, PARAMETER :: MAX_i=500
REAL(prec), PARAMETER :: TOLL=0.0001D0, TOLL_G=0.999D-6
INTEGER, PARAMETER :: MAX_i=100
!Define 2nd order identity tensor
data I_mat(1,:) /1.D0, 0.D0, 0.D0/
@@ -161,11 +173,11 @@ C !--------------------------------------------------------------
! for the first time step (initial condition) in FFTMAD.
! Can be commented out for ABAQUS, If you comment out, then
! in abaqus use the inp file to set initial values for F_vp.
C IF ((KINC .EQ. 1) .AND. (KSTEP .EQ. 1)) THEN
C STATEV(1) = 1.D0
C STATEV(2) = 1.D0
C STATEV(3) = 1.D0
C END IF
IF ((KINC .EQ. 1) .AND. (KSTEP .EQ. 1)) THEN
STATEV(1) = 1.D0
STATEV(2) = 1.D0
STATEV(3) = 1.D0
END IF
! State variables at previous time step
CALL voit2mat(STATEV(1:9), F_vp_n(:,:))
@@ -255,7 +267,7 @@ C END IF
! Approximate VE log Strain (power series) at trial state
D_ve_tr(:,:) = C_ve_tr(:, :) - I_mat(:,:)
CALL approx_log(D_ve_tr(:,:), order, E_ve_tr(:, :),
1 dlogC_ve_tr(:,:,:,:), ddlogC_ve_tr(:,:,:,:,:,:))
1 dlogC_ve_tr(:,:,:,:))
E_ve_tr(:, :) = 0.5D0 * E_ve_tr(:, :)
! Calculate the corotated kirchoff stress at trial state along
@@ -296,8 +308,7 @@ C END IF
! Approximate VE log Strain at trial state
D_ve_tr_hat(O6, :,:) = C_ve_tr_hat(O6, :, :) - I_mat(:,:)
CALL approx_log(D_ve_tr_hat(O6, :,:), order,
1 E_ve_tr_hat(O6, :, :), dlogC_ve_tr_hat(O6,:,:,:,:),
2 ddlogC_ve_tr_hat(O6,:,:,:,:,:,:))
1 E_ve_tr_hat(O6, :, :), dlogC_ve_tr_hat(O6,:,:,:,:))
E_ve_tr_hat(O6,:, :) = 0.5D0 * E_ve_tr_hat(O6, :, :)
! Calculate the corotated kirchoff stress at trial state
@@ -375,7 +386,7 @@ C END IF
DO ii = 1, 6
DO jj = 1, 6
DDSDDE(ii, jj) = (1.D0 / (J * TOLL))
1 * (tau_tr_hat_v(jj, ii) - J*STRESS(ii))
1 * (tau_tr_hat_v(jj, ii) - J*STRESS(ii))
END DO
END DO
@@ -400,14 +411,14 @@ C END IF
GQ(:,:) = GAMMA * Q(:,:)
CALL approx_exp(GQ(:,:), order, exp_GQ(:,:))
CALL mat2dot(exp_GQ(:,:), F_vp_n(:,:), F_vp(:,:))
CALL matInv(F_vp(:,:), F_vp_inv(:,:))
CALL matInverse(F_vp(:,:), F_vp_inv(:,:))
CALL mat2dot(DFGRD1(:,:), F_vp_inv(:,:), F_ve(:, :))
CALL mat2Tdot(F_ve(:,:), F_ve(:,:), C_ve(:,:))
! Approximate VE log Strain (power series) at the corrected state
D_ve(:,:) = C_ve(:, :) - I_mat(:,:)
CALL approx_log(D_ve(:,:), order, E_ve(:, :),
1 dlogC_ve(:,:,:,:), ddlogC_ve(:,:,:,:,:,:))
1 dlogC_ve(:,:,:,:))
E_ve(:, :) = 0.5D0 * E_ve(:, :)
! Calculate the corotated kirchoff stress at corrected state along
@@ -513,7 +524,7 @@ C END IF
CALL approx_exp(GQ_hat(O6,:,:), order, exp_GQ_hat(O6,:,:))
CALL mat2dot(exp_GQ_hat(O6, :,:), F_vp_n(:,:),
1 F_vp_hat(O6,:,:))
CALL matInv(F_vp_hat(O6,:,:), F_vp_inv_hat(O6,:,:))
CALL matInverse(F_vp_hat(O6,:,:), F_vp_inv_hat(O6,:,:))
CALL mat2dot(F_hat(O6, :,:), F_vp_inv_hat(O6,:,:),
1 F_ve_hat(O6, :, :))
CALL mat2Tdot(F_ve_hat(O6,:,:), F_ve_hat(O6,:,:),
@@ -523,8 +534,7 @@ C END IF
! Approximate VE log Strain (power series) at the purturbated state
D_ve_hat(O6,:,:) = C_ve_hat(O6,:, :) - I_mat(:,:)
CALL approx_log(D_ve_hat(O6,:,:), order, E_ve_hat(O6,:, :),
1 dlogC_ve_hat(O6,:,:,:,:),
2 ddlogC_ve_hat(O6,:,:,:,:,:,:))
1 dlogC_ve_hat(O6,:,:,:,:))
E_ve_hat(O6,:, :) = 0.5D0 * E_ve_hat(O6,:, :)
! Calculate the corotated kirchoff stress at purturbed state along
@@ -570,6 +580,119 @@ C END IF
RETURN
END SUBROUTINE VEVP
C !--------------------------------------------------------------
C ! UMAT SUBROUTINE FOR ISOTROPIC NEO HOOKIAN MATERIAL MODEL
C !--------------------------------------------------------------
SUBROUTINE NEOHOOK(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,
1 PROPS,NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1)
IMPLICIT NONE
INTEGER, PARAMETER :: double=kind(1.d0)
INTEGER, PARAMETER :: prec=double
INTEGER, INTENT(IN) :: ntens,nprops,nstatv
INTEGER, INTENT(IN) :: kinc,kstep,noel
REAL(prec), INTENT(IN) :: stran(ntens), dstran(ntens)
REAL(prec), INTENT(IN) :: props(nprops), dtime
REAL(prec), INTENT(IN) :: DFGRD0(3,3), DFGRD1(3,3)
REAL(prec), INTENT(INOUT) :: statev(nstatv)
REAL(prec), INTENT(OUT) :: stress(ntens)
REAL(prec), INTENT(OUT) :: ddsdde(ntens,ntens)
!List of internal variables
INTEGER :: ii, jj, mm, K1, ll, m2v(3, 3)
REAL(prec) :: I(6), B(6), B_bar(6), dev_B_bar(6)
REAL(prec) :: F_inv(3, 3)
REAL(prec) :: J, I_bar, fac1, fac2, fac3, fac4, cE, cnu
!Decleration of constants
REAL(prec), PARAMETER :: ZERO=0.D0, ONE=1.D0, TWO=2.D0
REAL(prec), PARAMETER :: THREE=3.D0, FOUR= 4.D0, SIX=6.D0
REAL(prec), PARAMETER :: NINE=9.D0
!Declare matrix to voit mapping
data m2v/ 1, 4, 5,
1 4, 2, 6,
2 5, 6, 3/
!2nd Order Identity
data I/ ONE, ONE, ONE, ZERO, ZERO, ZERO/
!Get material properties
cE=PROPS(1)
cnu=PROPS(2)
!Calculate determinant of deformation gradient
CALL determinant(DFGRD1(:,:), J)
!Calculate inverse of deformation gradient
CALL matInverse(DFGRD1(:,:), F_inv(:, :))
!Calculate finger tensor B
B(:) = ZERO
DO K1 = 1, 3
DO ii = 1, 3
DO jj = 1, 3
IF (ii .EQ. jj) THEN
B(m2v(ii, jj)) = B(m2v(ii, jj)) + DFGRD1(ii, K1)
1 * DFGRD1(jj, K1)
ELSE
B(m2v(ii, jj)) = B(m2v(ii, jj)) + 0.5d0
1 * DFGRD1(ii, K1) * DFGRD1(jj, K1)
END IF
END DO
END DO
END DO
B_bar(:) = B(:) * J ** (-TWO / THREE)
I_bar = B_bar(1) + B_bar(2) + B_bar(3)
dev_B_bar(:) = B_bar(:) - (ONE / THREE) * I_bar * I(:)
!Return Cauchy stress to Abaqus
STRESS(:) = 2 * cE * (J - ONE) * J * I(:)
1 + 2 * cnu * dev_B_bar(:)
!Algorithmic constants for the Tangent
fac1 = -(FOUR / THREE) * cnu * J ** (-ONE)
fac2 = -(TWO * cE * (J - ONE) * J
1 - (TWO/THREE) * cnu * I_bar) * J ** (-ONE)
fac3 = (TWO * cE * (J - ONE) * J
1 + (FOUR / NINE) * cnu * I_bar + TWO * cE * J ** TWO)
2 * J ** (-ONE)
fac4 = TWO * cnu * J ** (-ONE)
!Tangent for ABAQUS
DDSDDE(:,:) = ZERO
DO ii = 1, 3
DO jj = 1, 3
DO ll = 1, 3
DO mm = 1, 3
DDSDDE(m2v(ii, jj), m2v(ll, mm)) =
1 DDSDDE(m2v(ii, jj), m2v(ll, mm))
2 + fac1 * (B_bar(m2v(ii, jj)) * I(m2v(ll, mm))
3 + I(m2v(ii, jj)) * B_bar(m2v(ll, mm)))
4 + (fac4 - TWO * cnu * J**(-ONE))
5 * B_bar(m2v(ii, ll)) * I(m2v(jj, mm))
6 + TWO * cnu * (dev_B_bar(m2v(ii, ll))
7 * I(m2v(jj, mm))
8 + I(m2v(ii, ll)) * dev_B_bar(m2v(jj, mm)))
9 + fac2 * I(m2v(ii, mm)) * I(m2v(ll, jj))
1 + fac3 * I(m2v(ii, jj)) * I(m2v(ll, mm))
2 + (fac2 + FOUR * cE * (J - ONE))
3 * I(m2v(ii, ll)) * I(m2v(jj, mm))
END DO
END DO
END DO
END DO
RETURN
END SUBROUTINE NEOHOOK
!--------------------------------------------------------------
! Helper function to convert a voit array to matrix
!--------------------------------------------------------------
@@ -643,7 +766,7 @@ C END IF
!--------------------------------------------------------------
! Helper function to compute inverse of 3x3 matrix
!--------------------------------------------------------------
SUBROUTINE matInv(matrix, inv)
SUBROUTINE matInverse(matrix, inv)
IMPLICIT NONE
INTEGER, PARAMETER :: double=kind(1.d0)
@@ -668,7 +791,7 @@ C END IF
inv(:, :) = inv(:, :)/ J
RETURN
END SUBROUTINE matInv
END SUBROUTINE matInverse
!------------------------------------------------------------
@@ -843,7 +966,7 @@ C END IF
! Compute inverse of int variable F_vp
CALL matInv(F_vp(:, :), F_vp_inv(:, :))
CALL matInverse(F_vp(:, :), F_vp_inv(:, :))
! Compute F_ve
CALL mat2dot(F(:, :), F_vp_inv(:, :), F_ve(:, :))
@@ -929,7 +1052,7 @@ C END IF
! Power series approximation
!--------------------------------------------------------------
SUBROUTINE approx_log(AA, order, logAA, dlogAA, ddlogAA)
SUBROUTINE approx_log(AA, order, logAA, dlogAA)!, ddlogAA)
IMPLICIT NONE
INTEGER, PARAMETER :: double=kind(1.d0)
INTEGER, PARAMETER :: prec=double
@@ -937,13 +1060,14 @@ C END IF
INTEGER, INTENT(IN) :: order
REAL(prec), INTENT(IN) :: AA(3, 3)
REAL(prec), INTENT(OUT) :: logAA(3, 3), dlogAA(3, 3, 3, 3)
REAL(prec), INTENT(OUT) :: ddlogAA(3, 3, 3, 3, 3, 3)
!REAL(prec), INTENT(OUT) :: ddlogAA(3, 3, 3, 3, 3, 3)
! List of internal variables
INTEGER :: ii, jj, nn, mm, ll, rr, ss, pp, qq
INTEGER :: ii, jj, nn, mm, ll, rr, ss
REAL(prec) :: coeffs(order + 1), I_mat(3, 3)
REAL(prec) :: FF(3, 3), dFF(3, 3, 3, 3)
REAL(prec) :: ddFF(3, 3, 3, 3, 3, 3), II_mat(3, 3, 3, 3)
!REAL(prec) :: ddFF(3, 3, 3, 3, 3, 3)
REAL(prec) :: II_mat(3, 3, 3, 3)
REAL(prec) :: temp1(3, 3, 3, 3), temp2(3, 3, 3, 3)
! Define 2nd order identity tensor
@@ -990,12 +1114,12 @@ C END IF
END DO
dlogAA(:, :, :, :) = 0.D0
ddlogAA(:, :, :, :, :, :) = 0.D0
!ddlogAA(:, :, :, :, :, :) = 0.D0
DO ii = 1, order
FF(:, :) = 0.D0
dFF(:, :, :, :) = 0.D0
ddFF(:, :, :, :, :, :) = 0.D0
!ddFF(:, :, :, :, :, :) = 0.D0
nn = order + 1 - ii
@@ -1012,14 +1136,14 @@ C END IF
1 + dlogAA(jj, ll, rr, ss) * AA(ll, mm)
2 + logAA(jj, ll) * II_mat(ll, mm, rr, ss)
DO pp = 1, 3
DO qq = 1, 3
ddFF(jj, mm, rr, ss, pp, qq) = ddFF(jj, mm, rr, ss, pp, qq)
1 + ddlogAA(jj, ll, rr, ss, pp, qq) * AA(ll, mm)
2 + dlogAA(jj, ll, rr, ss) * II_mat(ll, mm, pp, qq)
3 + dlogAA(jj, ll, pp, qq) * II_mat(ll, mm, rr, ss)
END DO
END DO
! DO pp = 1, 3
! DO qq = 1, 3
! ddFF(jj, mm, rr, ss, pp, qq) = ddFF(jj, mm, rr, ss, pp, qq)
! 1 + ddlogAA(jj, ll, rr, ss, pp, qq) * AA(ll, mm)
! 2 + dlogAA(jj, ll, rr, ss) * II_mat(ll, mm, pp, qq)
! 3 + dlogAA(jj, ll, pp, qq) * II_mat(ll, mm, rr, ss)
! END DO
! END DO
END DO
END DO
@@ -1030,7 +1154,7 @@ C END IF
END DO
logAA(:, :) = FF(:, :)
dlogAA(:,:,:,:) = dFF(:,:,:,:)
ddlogAA(:,:,:,:,:,:) = ddFF(:, :, :, :, :, :)
!ddlogAA(:,:,:,:,:,:) = ddFF(:, :, :, :, :, :)
END DO
END SUBROUTINE approx_log
@@ -1399,8 +1523,8 @@ C END IF
gma_i = gma + Dgma
! Update hardening variables on gma_i
CALL hardn(PROPS, NPROPS, gma_0, gma_i, sigma_c, sigma_t,
1 HHc, HHt, HHb)
CALL hardn(PROPS, NPROPS, gma_0, gma_i, sigma_c,
1 sigma_t, HHc, HHt, HHb)
! Update Drucker Pragger coefficients along with ecc factor m
CALL DPcoeff(alpha, sigma_c, sigma_t, m, a0, a1, a2)
@@ -1421,7 +1545,7 @@ C END IF
gma = gma_i
! Debuging statements for N.R. Comment out for big sims
! !Debuging statements for N.R. Comment out for big sims
! WRITE(*,*) "Iter : ", iter_G
! WRITE(*,*) "GAMMA : ", GAMMA
! WRITE(*,*) "RESI : ", ABS(F)