I am debugging my Fortran code for my UMAT subroutine, and I noticed that some of my parameters keep changing throughout a subroutine which is strange to me. For example I have the following subroutine :
SUBROUTINE fabric_tensor_Mandel(vm1,vm2,vm3,v_m11,v_m22,v_m33, + sigma_0_plus,sigma_0_minus,vkisi_0,Tau_0,p,q,ro,AA,A,F_F, + vMATRIX_F,NTENS) C INCLUDE 'ABA_PARAM.INC' C PARAMETER (M=3) DIMENSION AA(NTENS,NTENS),A(NTENS),F_F(NTENS,NTENS),F(M,M), + sigma_plus(M),sigma_minus(M),Tau(M,M),vkisi(M,M),v_m(M), + vMM1(M,M),vMM2(M,M),vMM3(M,M),vMM1_v(NTENS),vMM2_v(NTENS), + vMM3_v(NTENS),vM11_dyad(NTENS,NTENS),vM22_dyad(NTENS,NTENS), + vM33_dyad(NTENS,NTENS),vM12_dyad(NTENS,NTENS),vMATRIX_F(M,M), + vM13_dyad(NTENS,NTENS),vM23_dyad(NTENS,NTENS), + vM21_dyad(NTENS,NTENS),vM31_dyad(NTENS,NTENS), + vM32_dyad(NTENS,NTENS),vM12_Sym(NTENS,NTENS), + vM13_Sym(NTENS,NTENS),vM23_Sym(NTENS,NTENS),FFA(NTENS), + vM21_Sym(NTENS,NTENS),vM31_Sym(NTENS,NTENS),FF_inv(NTENS,NTENS), + vM32_Sym(NTENS,NTENS),FF_inv_F(NTENS),F_V(NTENS), + FF_MATRIX(NTENS,NTENS) C C m=[m1;m2;m3]; v_m=(/vm1,vm2,vm3/) C Mi = mi (dyad) mi (eq.38) C MMi=mii*transpose(mii); CALL DYADICPROD(v_m11,v_m11,vMM1,M) CALL DYADICPROD(v_m22,v_m22,vMM2,M) CALL DYADICPROD(v_m33,v_m33,vMM3,M) C %convert to voigt notation C [MM1_V]= convert_matrix_vector(MM1); CALL TENS2VEC(vMM1,M,vMM1_v,NTENS) CALL TENS2VEC(vMM2,M,vMM2_v,NTENS) CALL TENS2VEC(vMM3,M,vMM3_v,NTENS) C (eq. 15-18 Barakuie) DO i=1,3 DO j=1,3 sigma_plus(i) = sigma_0_plus*(ro**p)*((v_m(i))(2*q)) sigma_minus(i) = sigma_0_minus*(ro**p)*((v_m(i))(2*q)) Tau(i,j) = Tau_0*(ro**p)*((v_m(i))(q))*((v_m(j))(q)) vkisi(i,j)=vkisi_0*((v_m(i))(2*q))/((v_m(j))(2*q)) END DO END DO C Construt F (eq.37) DO i=1,3 DO j=1,3 F(i,j)=((1.d0/sigma_plus(1))-(1.d0/sigma_minus(1)))*vMM1(i,j)+ 1 ((1.d0/sigma_plus(2))-(1.d0/sigma_minus(2)))*vMM2(i,j)+ 2 ((1.d0/sigma_plus(3))-(1.d0/sigma_minus(3)))*vMM3(i,j) END DO END DO C Convert to voigt notation C [F_V]= convert_matrix_vector(F); CALL TENS2VEC(F,M,F_V,NTENS) C Construct FF (eq. 36) C Construct the required dyadic product C [M11_dyad]= MM1_V*transpose(MM1_V); CALL DYADICPROD(vMM1_v,vMM1_v,vM11_dyad,NTENS) CALL DYADICPROD(vMM2_v,vMM2_v,vM22_dyad,NTENS) CALL DYADICPROD(vMM3_v,vMM3_v,vM33_dyad,NTENS) CALL DYADICPROD(vMM1_v,vMM2_v,vM12_dyad,NTENS) CALL DYADICPROD(vMM1_v,vMM3_v,vM13_dyad,NTENS) CALL DYADICPROD(vMM2_v,vMM1_v,vM21_dyad,NTENS) CALL DYADICPROD(vMM2_v,vMM3_v,vM23_dyad,NTENS) CALL DYADICPROD(vMM3_v,vMM1_v,vM31_dyad,NTENS) CALL DYADICPROD(vMM3_v,vMM2_v,vM32_dyad,NTENS) C Construct the required symmetric product CALL Symmetric_Zysset_Mandel(vMM1_v,vMM2_v,vM12_Sym) CALL Symmetric_Zysset_Mandel(vMM1_v,vMM3_v,vM13_Sym) CALL Symmetric_Zysset_Mandel(vMM2_v,vMM1_v,vM21_Sym) CALL Symmetric_Zysset_Mandel(vMM2_v,vMM3_v,vM23_Sym) CALL Symmetric_Zysset_Mandel(vMM3_v,vMM1_v,vM31_Sym) CALL Symmetric_Zysset_Mandel(vMM3_v,vMM2_v,vM32_Sym) c Construct fourth-order tensor FF DO i=1,NTENS DO j=1,NTENS F_F(i,j)=((1.d0/(sigma_plus(1)*sigma_minus(1)))*vM11_dyad(i,j))+ 1 ((1.d0/(sigma_plus(2)*sigma_minus(2)))*vM22_dyad(i,j))+ 2 ((1.d0/(sigma_plus(3)*sigma_minus(3)))*vM33_dyad(i,j))- 3 ((1.d0/(sigma_plus(1)*sigma_minus(1)))*(vkisi(1,2)*vM12_dyad(i,j)+ 4 vkisi(1,3)*vM13_dyad(i,j)))- 5 ((1.d0/(sigma_plus(2)*sigma_minus(2)))*(vkisi(2,1)*vM21_dyad(i,j)+ 6 vkisi(2,3)*vM23_dyad(i,j)))- 7 ((1.d0/(sigma_plus(3)*sigma_minus(3)))*(vkisi(3,1)*vM31_dyad(i,j)+ 8 vkisi(3,2)*vM32_dyad(i,j)))+ 9 ((0.5d0/(Tau(1,2)**2.d0))*vM12_Sym(i,j)+ + (0.5d0/(Tau(1,3)**2.d0))*vM13_Sym(i,j)+ + (0.5d0/(Tau(2,1)**2.d0))*vM21_Sym(i,j) + + ((0.5d0/(Tau(2,3)**2.d0))*vM23_Sym(i,j)+ + (0.5d0/(Tau(3,1)**2.d0))*vM31_Sym(i,j)+ + (0.5d0/(Tau(3,2)**2.d0))*vM32_Sym(i,j))) END DO END DO DO i=1,NTENS DO j=1,NTENS FF_MATRIX(i,j)=F_F(i,j) END DO END DO C FF_inv=inv(FF); %compute inverse of 4th order tensor FF CALL INVERT(FF_MATRIX,FF_inv,IPIV,6,X,FSMAL,IOUT) C FF_inv_F=inv(FF)*F_V; CALL KMLT1(FF_inv,F_V,FF_inv_F,NTENS) C A = -0.5*FF_inv_F; %Construct A (3,3) by eq.35 DO K=1,NTENS A(K)=-0.5*FF_inv_F(K) END DO C FFA= FF*A; %double cntraction of 4th-order FF and 2nd order A to be used in eq. 34 CALL KMLT1(F_F,A,FFA,NTENS) C AFFA = transpose(A)*FFA; %(scalar) %Construct A:FF:A (scalar) in the denominator of eq. 34 CALL DOTPROD(A,FFA,AFFA,NTENS) C AA = FF./(1+AFFA); %Construct 4th-order tensor AA DO i=1,NTENS DO j=1,NTENS AA(i,j) =(1/(1+AFFA))*F_F(i,j); END DO END DO RETURN END
but, for example, the matrix F_F changes after the following line:
DO i=1,NTENS DO j=1,NTENS AA(i,j) =(1/(1+AFFA))*F_F(i,j); END DO END DO
Can someone help me why does this happen and what should I do to fix it? Thanks, Faezeh
I don't know if this is what you're experiencing, but I hope you do know that you can inadvertently change the value of zero or 1 or any other fundamental constant in FORTRAN simply by calling a subroutine and passing it a constant. All arguments in FORTRAN are passed as pointers. In C and VBA, for example, you can pass arguments as values or pointers. More recent FORTRAN compilers allow qualifiers that block changes to specific arguments but even that requires you to enforce this distinction in the code. Any arguments you don't want changed, search for "variable=?". If it's on the left side of an = then you're in trouble. You may pass variables through many levels of subroutines where they might be inadvertently changed. This most often happens when you're using code that somebody else wrote.
Also, be sure to double-check if you are reading the parameters, vectors etc. correctly - i.e. if their length (dimension) and reading order in the problematic subroutine is the same as in the master subroutine.
Thanks, Vid Merljak for your answer. What do you mean by reading order in the problematic subroutine should be the same as in the master subroutine?
The size of individual parameters and also the dimension of arrays. For example, 4-byte (32-bit) integer, 4-byte single precision floating point, 8-byte double precision floating point, etc. Also know that 3.14159 in FORTRAN is a 4-byte single precision floating point and 3.14159D0 is an 8-byte double precision floating point. Some compilers check for argument size agreement and some do not, while for others it depends on the compiling options. Some constants are automatically "promoted" (elevated from a single to double precision floating point, while others are not). It also depends on the compiler and the target process.
Dudley J Benton, exactly! A parameter - say, the first one - can be passed with name 'foo' and read (in the subroutine) as 'bar'. If for example 'foo' has DIM(5) and 'bar' has DIM(6,2) every other parameter will be read from an incorrect position in the working memory. Leading to ... disaster.
Also, Faezeh Iranmanesh, are you using 'implicit none'? Strongly recommended. See for example:
- https://stackoverflow.com/questions/24337413/where-to-put-implicit-none-in-fortran#24337734
- https://www.tutorialspoint.com/fortran/fortran_arrays.htm
My boss' boss was on the witness stand in a civil suit over a power plant that failed. My boss had written in his notes that the manufacturer's data was FUBAR. The attorney demanded, "Explain to the court what FUBAR means." He glared at my boss and began, "It's a military term like SNAFU," which only made it worse.
- Badges
- Science topic
- Similar topics
- Mathematics
- Algebra
- Matrix