Shell Element DT is too bigger #113
Unanswered
Myungjun-god
asked this question in
Q&A
Replies: 0 comments
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Uh oh!
There was an error while loading. Please reload this page.
I am currently planning to conduct a simple tension-followed-by-compression test using a Shell Element model.
The issue is that the DT value comes out high, and the process simply stops after a single cycle. I need assistance in improving this aspect.
Below is my code : Two rad file (Starter, Engine) & Two fortran code
plz help me!
#RADIOSS STARTER
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BEGIN
Job1_VM2D_2x2_Shell_Bauschinger
2022 0
Mg mm s
Mg mm s
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 1. CONTROL CARDS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/TITLE
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 2. MATERIALS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/MAT/USER01/1
VM2D_SwiftVoce_Al
RHO_I
E NU
SWIFTK SWIFTE0 SWIFTN
VOCESAT VOCE0 VOCEC ALPHA
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 3. NODES (2x2 mesh, quarter-symmetry model: origin = symmetry intersection):
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/NODE
1 0 0 0
2 1.0 0 0
3 2.0 0 0
4 0 1.0 0
5 1.0 1.0 0
6 2.0 1.0 0
7 0 2.0 0
8 1.0 2.0 0
9 2.0 2.0 0
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 4. PARTS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PART/1
Part-1
prop_ID mat_ID subset_ID Thick
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/SHELL/1
shell_ID node1 node2 node3 node4
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 5. BOUNDARY CONDITIONS (Quarter symmetry):
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
xsymm (X=0 plane, nodes 1,4,7): Abaqus DOF(1,5,6)=0 -> Tx,Ry,Rz fixed -> Trarot = 100 011
/BCS/1
xsymm
Trarot skew_ID grnod_ID
100 011 0 11
/GRNOD/NODE/11
xsymm_nodes
1
4
7
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
ysymm (Y=0 plane, nodes 1,2,3): Abaqus DOF(2,4,6)=0 -> Ty,Rx,Rz fixed -> Trarot = 010 101
/BCS/2
ysymm
Trarot skew_ID grnod_ID
010 101 0 12
/GRNOD/NODE/12
ysymm_nodes
1
2
3
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
Out-of-plane stabilization: fix Tz on all nodes
/BCS/3
out_of_plane
Trarot skew_ID grnod_ID
001 000 0 13
/GRNOD/NODE/13
all_nodes
1
2
3
4
5
6
7
8
9
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 6. GEOMETRICAL SETS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PROP/SHELL/1
Section-1
Ishell Ismstr Ish3n Idrill P_thick_fail
hm hf hr dm dn
N Istrain Thick Ashear Ithick Iplas Ipos
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 7. FUNCTIONS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/FUNCT/1
Loading_DispY_TensionThenReturn
X Y
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 8. IMPOSED DISPLACEMENT (free edge, Y direction, nodes 7,8,9):
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/GRNOD/NODE/21
loading_nodes
7
8
9
/IMPDISP/1
loading_dispY
#fct_IDT Dir Skew_ID sens_ID grnd_ID icoor
1 Y 0 0 21
Ascale_x Fscale_Y Tstart Tstop
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 9. TIME HISTORIES:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/TH/SHEL/1
TH_Shell
var1 var2 var3 var4 var5 var6 var7 var8 var9 var10
DEF
elem_ID skew_ID elem_name
/TH/SHEL/2
TH_Shell2
var1 var2 var3 var4 var5 var6 var7 var8 var9 var10
DEF
elem_ID skew_ID elem_name
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/TH/NODE/3
TH_Loading_Node
var1 var2 var3 var4 var5 var6 var7 var8 var9 var10
DEF
NODid Iskew NODname
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/END
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/ANIM/DT
0 0.01
/ANIM/ELEM/VONM
/ANIM/ELEM/EPSP
/ANIM/SHELL/TENS/STRESS/ALL
/ANIM/SHELL/TENS/STRAIN/ALL
/PRINT/-100/55
/RUN/Job1_VM2D_2x2_Shell_Bauschinger/1/
2.0
/STOP
Emax Mmax Nmax NTH NANIM NERR_POSIT
0 0 0 1 1 0
/TFILE/0
dT_HIS
0.000100
/VERS/2022
/DT/NODA/CST/0
0.9 5e-07
#-------------------------------------
UNCOMMENT LINES BELOW FOR H3D OUTPUT
#-------------------------------------
#/H3D/NODA/VEL
#/H3D/SHELL/TENS/STRESS/NPT=ALL
#/H3D/SHELL/TENS/STRAIN/NPT=LOWER
#/H3D/SHELL/TENS/STRAIN/NPT=UPPER
#/H3D/ELEM/EPSP/NPT=UPPER
#/H3D/ELEM/EPSP/NPT=LOWER
#/H3D/DT
#0.000000000000000 0.001000000000000
lecmuser.f
. NUVAR,IFUNC,MAXFUNC,NFUNC,STIFINT,
. USERBUF)
C-----------------------------------------------------------------------
C OpenRadioss starter user material reader for Abaqus VUMAT:
C VM_VUMAT_2D.for
C
C Abaqus PROPS / OpenRadioss UPARAM order:
C 1 E
C 2 NU
C 3 SWIFTK
C 4 SWIFTE0
C 5 SWIFTN
C 6 VOCESAT
C 7 VOCE0
C 8 VOCEC
C 9 ALPHA
C
C State variables:
C UVAR(:,1) = equivalent plastic strain
C UVAR(:,2) = Newton iteration count
C-----------------------------------------------------------------------
USE LAW_USER
IMPLICIT NONE
INTEGER IIN,IOUT,MAXUPARAM,NUPARAM,NUVAR,MAXFUNC,NFUNC
INTEGER IFUNC(MAXFUNC)
DOUBLE PRECISION UPARAM(MAXUPARAM),STIFINT
TYPE(ULAWBUF) :: USERBUF
C
DOUBLE PRECISION E,NU,SWIFTK,SWIFTE0,SWIFTN
DOUBLE PRECISION VOCESAT,VOCE0,VOCEC,ALPHA
C
READ(IIN,'(2F20.0)') E,NU
READ(IIN,'(3F20.0)') SWIFTK,SWIFTE0,SWIFTN
READ(IIN,'(4F20.0)') VOCESAT,VOCE0,VOCEC,ALPHA
C
NUPARAM = 9
NUVAR = 2
NFUNC = 0
C
UPARAM(1) = E
UPARAM(2) = NU
UPARAM(3) = SWIFTK
UPARAM(4) = SWIFTE0
UPARAM(5) = SWIFTN
UPARAM(6) = VOCESAT
UPARAM(7) = VOCE0
UPARAM(8) = VOCEC
UPARAM(9) = ALPHA
C
STIFINT = E
C
WRITE(IOUT,1000)
1000 FORMAT(
& 5X,' VM 2D plane-stress Swift-Voce user material',/,
& 5X,' Converted from Abaqus VUMAT: VM_VUMAT_2D.for',/,
& 5X,' NUPARAM=9, NUVAR=2',//)
C
RETURN
END
1 NEL ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,
2 NPF ,NGL ,TF ,TIME ,TIMESTEP ,
3 UPARAM ,RHO0 ,AREA ,EINT ,SHF ,
4 SOUNDSP,VISCMAX,PLA ,UVAR ,OFF ,
5 SIGY ,USERBUF)
C-----------------------------------------------------------------------
C OpenRadioss SHELL user material for Abaqus VUMAT:
C VM_VUMAT_2D.for
C
C Use this file for /MAT/USER01 assigned to shell elements.
C Shell user laws in the SDK use the C-suffixed routine name:
C LUSER01C + LAW_USERSH + ULAWCINTBUF
C
C UPARAM order follows Abaqus PROPS(1:9):
C E, NU, SWIFTK, SWIFTE0, SWIFTN, VOCESAT, VOCE0, VOCEC, ALPHA
C
C UVAR usage:
C UVAR(:,1) = equivalent plastic strain
C UVAR(:,2) = Newton iteration count
C-----------------------------------------------------------------------
USE LAW_USERSH
IMPLICIT NONE
C
INTEGER NEL,NUPARAM,NUVAR,NFUNC
INTEGER IFUNC(NFUNC),NPF(),NGL(NEL)
DOUBLE PRECISION TF(),TIME
DOUBLE PRECISION TIMESTEP(NEL),UPARAM(NUPARAM),RHO0(NEL)
DOUBLE PRECISION AREA(NEL),EINT(2,NEL),SHF(NEL)
DOUBLE PRECISION SOUNDSP(NEL),VISCMAX(NEL),PLA(NEL)
DOUBLE PRECISION UVAR(NEL,NUVAR),OFF(NEL),SIGY(NEL)
TYPE(ULAWCINTBUF) :: USERBUF
C
INTEGER I,J,KNRITER
DOUBLE PRECISION ZERO,ONE,TWO,THREE,TOL
PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,THREE=3.0D0)
PARAMETER (TOL=1.0D-8)
C
DOUBLE PRECISION EMOD,ENU,SWIFTK,SWIFTE0,SWIFTN
DOUBLE PRECISION VOCESAT,VOCE0,VOCEC,ALPHA
DOUBLE PRECISION EMAT(4,4),V(4,4),TEMPVAL
DOUBLE PRECISION DSTRAN(4),STRESST(4),STRESS(4)
DOUBLE PRECISION FLOW(4),DFLOW(4,4),DSTRESS(4)
DOUBLE PRECISION TEMPVEC(4),RESIVEC(5),FRESIVEC(4)
DOUBLE PRECISION EYE4(4,4),Q(4,4),QINV(4,4),ADYAA(4,4)
DOUBLE PRECISION STRAINP,STRESSY,APRIME,STRESSEFF
DOUBLE PRECISION DLAM,DDLAM,YF,RESIDUAL,RCRITERION
DOUBLE PRECISION CSOUND
C
DOUBLE PRECISION
. SIGOXX(MVSIZ),SIGOYY(MVSIZ),SIGOXY(MVSIZ),
. DEPSXX(MVSIZ),DEPSYY(MVSIZ),DEPSXY(MVSIZ),
. SIGNXX(MVSIZ),SIGNYY(MVSIZ),SIGNXY(MVSIZ),
. DPLA(MVSIZ),ETSE(MVSIZ),THKN(MVSIZ)
C
IF (NUPARAM .LT. 9) THEN
WRITE(,) 'ERROR LUSER01C_VM2D: NUPARAM must be >= 9'
RETURN
ENDIF
C
EMOD = UPARAM(1)
ENU = UPARAM(2)
SWIFTK = UPARAM(3)
SWIFTE0 = UPARAM(4)
SWIFTN = UPARAM(5)
VOCESAT = UPARAM(6)
VOCE0 = UPARAM(7)
VOCEC = UPARAM(8)
ALPHA = UPARAM(9)
C
C Copy input state from OpenRadioss shell user-law buffer.
SIGOXX(1:NEL) = USERBUF%SIGOXX(1:NEL)
SIGOYY(1:NEL) = USERBUF%SIGOYY(1:NEL)
SIGOXY(1:NEL) = USERBUF%SIGOXY(1:NEL)
DEPSXX(1:NEL) = USERBUF%DEPSXX(1:NEL)
DEPSYY(1:NEL) = USERBUF%DEPSYY(1:NEL)
DEPSXY(1:NEL) = USERBUF%DEPSXY(1:NEL)
SIGNXX(1:NEL) = USERBUF%SIGNXX(1:NEL)
SIGNYY(1:NEL) = USERBUF%SIGNYY(1:NEL)
SIGNXY(1:NEL) = USERBUF%SIGNXY(1:NEL)
DPLA(1:NEL) = USERBUF%DPLA(1:NEL)
ETSE(1:NEL) = USERBUF%ETSE(1:NEL)
THKN(1:NEL) = USERBUF%THKN(1:NEL)
C
C Plane-stress elastic stiffness, matching the Abaqus VUMAT.
EMAT = ZERO
TEMPVAL = EMOD/(ONE-ENUENU)
EMAT(1,1) = TEMPVAL
EMAT(1,2) = TEMPVALENU
EMAT(2,1) = TEMPVALENU
EMAT(2,2) = TEMPVAL
EMAT(4,4) = TEMPVAL((ONE-ENU)/TWO)
C
C Sound speed for stability time step and hourglass forces.
C Required by Radioss doc. c = sqrt(E/((1-nu^2)rho0))
C
C Plane-stress Von Mises matrix. Component 3 is unused sigma33.
V = ZERO
V(1,1) = ONE
V(1,2) = -ONE/TWO
V(2,1) = -ONE/TWO
V(2,2) = ONE
V(4,4) = THREE
C
DO I = 1,NEL
IF (RHO0(I) .GT. ZERO) THEN
CSOUND = DSQRT(EMOD/((ONE-ENUENU)*RHO0(I)))
ELSE
CSOUND = ZERO
ENDIF
SOUNDSP(I) = CSOUND
VISCMAX(I) = ZERO
C
DSTRAN(1) = DEPSXX(I)
DSTRAN(2) = DEPSYY(I)
DSTRAN(3) = ZERO
DSTRAN(4) = DEPSXY(I)
C
IF (TIME .LT. 5.0D-4) THEN
WRITE(,) 'DBG TIME=',TIME,' I=',I,
1 ' DEPSXX=',DEPSXX(I),' DEPSYY=',DEPSYY(I),
2 ' SIGOXX=',SIGOXX(I),' SIGOYY=',SIGOYY(I)
ENDIF
C
STRESS(1) = SIGOXX(I)
STRESS(2) = SIGOYY(I)
STRESS(3) = ZERO
STRESS(4) = SIGOXY(I)
C
STRESST = STRESS + MATMUL(EMAT,DSTRAN)
CALL SVHARDEN(STRAINP,SWIFTK,SWIFTE0,SWIFTN,
1 VOCESAT,VOCE0,VOCEC,ALPHA,STRESSY,APRIME)
CALL VMSTRESS2D(STRESST,V,STRESSEFF)
STRESS = STRESST
YF = STRESSEFF - STRESSY
C
IF (YF .LE. TOL) THEN
SIGY(I) = STRESSY
SIGNXX(I) = STRESST(1)
SIGNYY(I) = STRESST(2)
SIGNXY(I) = STRESST(4)
DPLA(I) = ZERO
ETSE(I) = ONE
PLA(I) = STRAINP
IF (NUVAR .GE. 1) UVAR(I,1) = STRAINP
IF (NUVAR .GE. 2) UVAR(I,2) = ZERO
CYCLE
ENDIF
C
KNRITER = 1
DLAM = ZERO
RESIDUAL = 1000.0D0
RCRITERION = 1.0D-2
C
DO WHILE (RESIDUAL .GT. RCRITERION)
IF (KNRITER .GT. 200) THEN
WRITE(,) 'WARNING LUSER01C_VM2D: NR not converged'
EXIT
ENDIF
KNRITER = KNRITER + 1
C
CALL VMSTRESS2D(STRESS,V,STRESSEFF)
IF (STRESSEFF .LT. 1.0D-30) STRESSEFF = 1.0D-30
FLOW = MATMUL(V,STRESS)/STRESSEFF
CALL DYADIC(FLOW,FLOW,4,ADYAA)
DFLOW = (V-ADYAA)/STRESSEFF
C
CALL SVHARDEN(STRAINP,SWIFTK,SWIFTE0,SWIFTN,
1 VOCESAT,VOCE0,VOCEC,ALPHA,STRESSY,APRIME)
C
TEMPVEC = MATMUL(EMAT,FLOW)
DO J = 1,4
FRESIVEC(J) = STRESS(J)-STRESST(J)+DLAMTEMPVEC(J)
ENDDO
YF = STRESSEFF - STRESSY
C
CALL MAKEEYE(4,EYE4)
Q = EYE4 + DLAMMATMUL(EMAT,DFLOW)
CALL DOINV(Q,4,QINV)
C
DDLAM = (YF-DOT_PRODUCT(FLOW,MATMUL(QINV,FRESIVEC)))
1 /(DOT_PRODUCT(FLOW,MATMUL(QINV,
2 MATMUL(EMAT,FLOW)))+APRIME)
DSTRESS = -MATMUL(QINV,FRESIVEC
1 + DDLAMMATMUL(EMAT,FLOW))
C
DLAM = DLAM + DDLAM
STRESS = STRESS + DSTRESS
STRAINP = STRAINP + DDLAM
C
TEMPVEC = MATMUL(EMAT,FLOW)
DO J = 1,4
RESIVEC(J) = STRESS(J)-STRESST(J)+DLAMTEMPVEC(J)
ENDDO
CALL VMSTRESS2D(STRESS,V,STRESSEFF)
CALL SVHARDEN(STRAINP,SWIFTK,SWIFTE0,SWIFTN,
1 VOCESAT,VOCE0,VOCEC,ALPHA,STRESSY,APRIME)
RESIVEC(5) = STRESSEFF - STRESSY
CALL VECNORM(RESIVEC,5,2,RESIDUAL)
ENDDO
C
SIGY(I) = STRESSY
SIGNXX(I) = STRESS(1)
SIGNYY(I) = STRESS(2)
SIGNXY(I) = STRESS(4)
DPLA(I) = DLAM
ETSE(I) = APRIME/(APRIME+EMOD)
PLA(I) = STRAINP
C
C Thickness update for plastic incompressibility:
C dEPS_zz_plastic = -(dEPS_xx_plastic + dEPS_yy_plastic)
C Plastic strain increments from associative flow:
C dEPS_pl = DLAM * dF/dSIG where F = von Mises yield
C For plane stress vM: dEPS33_pl = -(dEPS11_pl + dEPS22_pl)
C Use FLOW directions computed during the last NR iteration.
THKN(I) = THKN(I) * (ONE - DLAM*(FLOW(1)+FLOW(2)))
C
IF (NUVAR .GE. 1) UVAR(I,1) = STRAINP
IF (NUVAR .GE. 2) UVAR(I,2) = DBLE(KNRITER)
ENDDO
C
C Copy results back to the OpenRadioss shell buffer.
USERBUF%SIGNXX(1:NEL) = SIGNXX(1:NEL)
USERBUF%SIGNYY(1:NEL) = SIGNYY(1:NEL)
USERBUF%SIGNXY(1:NEL) = SIGNXY(1:NEL)
USERBUF%DPLA(1:NEL) = DPLA(1:NEL)
USERBUF%ETSE(1:NEL) = ETSE(1:NEL)
USERBUF%THKN(1:NEL) = THKN(1:NEL)
C
RETURN
END
C=======================================================================
SUBROUTINE SVHARDEN(EP,SK,SE0,SN,VSAT,V0,VC,AA,SY,HAPRIME)
IMPLICIT NONE
DOUBLE PRECISION EP,SK,SE0,SN,VSAT,V0,VC,AA,SY,HAPRIME
DOUBLE PRECISION SWIFT_VAL,VOCE_VAL,SWIFT_D,VOCE_D
DOUBLE PRECISION ONE
PARAMETER (ONE=1.0D0)
C
SWIFT_VAL = SK*((SE0+EP)SN)
VOCE_VAL = VSAT-(VSAT-V0)DEXP(-VCEP)
SY = AASWIFT_VAL + (ONE-AA)VOCE_VAL
C
SWIFT_D = SNSK((SE0+EP)(SN-ONE))
VOCE_D = VC*(VSAT-V0)DEXP(-VCEP)
HAPRIME = AA*SWIFT_D + (ONE-AA)*VOCE_D
RETURN
END
C=======================================================================
SUBROUTINE VMSTRESS2D(STRESS,V,STRESSVM)
IMPLICIT NONE
DOUBLE PRECISION STRESS(4),V(4,4),STRESSVM
DOUBLE PRECISION TEMPVEC(4),TEMPVAL
C
TEMPVEC = MATMUL(V,STRESS)
TEMPVAL = DOT_PRODUCT(STRESS,TEMPVEC)
IF (TEMPVAL .LT. 0.0D0) TEMPVAL = 0.0D0
STRESSVM = DSQRT(TEMPVAL)
RETURN
END
C=======================================================================
SUBROUTINE MAKEEYE(N,EYE)
IMPLICIT NONE
INTEGER N,I,J
DOUBLE PRECISION EYE(N,N)
C
DO I = 1,N
DO J = 1,N
EYE(I,J) = 0.0D0
ENDDO
EYE(I,I) = 1.0D0
ENDDO
RETURN
END
C=======================================================================
SUBROUTINE DOINV(A,N,AINV)
IMPLICIT NONE
INTEGER N,I,J,K,PIVOTROW
DOUBLE PRECISION A(N,N),AINV(N,N),L(N,N),U(N,N)
DOUBLE PRECISION P(N,N),LINV(N,N),UINV(N,N)
DOUBLE PRECISION PIVOT,TEMPVAL,MAXVAL
C
CALL MAKEEYE(N,L)
CALL MAKEEYE(N,P)
U = A
LINV = 0.0D0
UINV = 0.0D0
C
DO J = 1,N
MAXVAL = ABS(U(J,J))
PIVOTROW = J
DO I = J+1,N
IF (ABS(U(I,J)) .GT. MAXVAL) THEN
MAXVAL = ABS(U(I,J))
PIVOTROW = I
ENDIF
ENDDO
C
IF (PIVOTROW .NE. J) THEN
DO K = 1,N
TEMPVAL = U(J,K)
U(J,K) = U(PIVOTROW,K)
U(PIVOTROW,K) = TEMPVAL
TEMPVAL = P(J,K)
P(J,K) = P(PIVOTROW,K)
P(PIVOTROW,K) = TEMPVAL
ENDDO
DO K = 1,J-1
TEMPVAL = L(J,K)
L(J,K) = L(PIVOTROW,K)
L(PIVOTROW,K) = TEMPVAL
ENDDO
ENDIF
C
PIVOT = U(J,J)
IF (ABS(PIVOT) .LT. 1.0D-30) PIVOT = 1.0D-30
DO I = J+1,N
L(I,J) = U(I,J)/PIVOT
DO K = J,N
U(I,K) = U(I,K)-L(I,J)*U(J,K)
ENDDO
ENDDO
ENDDO
C
DO I = 1,N
LINV(I,I) = 1.0D0
DO J = 1,I-1
TEMPVAL = 0.0D0
DO K = J,I-1
TEMPVAL = TEMPVAL + L(I,K)*LINV(K,J)
ENDDO
LINV(I,J) = -TEMPVAL
ENDDO
ENDDO
C
DO I = N,1,-1
PIVOT = U(I,I)
IF (ABS(PIVOT) .LT. 1.0D-30) PIVOT = 1.0D-30
UINV(I,I) = 1.0D0/PIVOT
DO J = I+1,N
TEMPVAL = 0.0D0
DO K = I+1,J
TEMPVAL = TEMPVAL + U(I,K)*UINV(K,J)
ENDDO
UINV(I,J) = -TEMPVAL/PIVOT
ENDDO
ENDDO
C
AINV = MATMUL(MATMUL(UINV,LINV),P)
RETURN
END
C=======================================================================
SUBROUTINE VECNORM(VEC,N,P,NORM)
IMPLICIT NONE
INTEGER N,P,I
DOUBLE PRECISION VEC(N),NORM,SUMV
C
SUMV = 0.0D0
IF (P .EQ. 2) THEN
DO I = 1,N
SUMV = SUMV + VEC(I)*VEC(I)
ENDDO
NORM = DSQRT(SUMV)
ELSE
DO I = 1,N
SUMV = SUMV + ABS(VEC(I))
ENDDO
NORM = SUMV
ENDIF
RETURN
END
C=======================================================================
SUBROUTINE DYADIC(VECA,VECB,N,MAT)
IMPLICIT NONE
INTEGER N,I,J
DOUBLE PRECISION VECA(N),VECB(N),MAT(N,N)
C
DO I = 1,N
DO J = 1,N
MAT(I,J) = VECA(I)*VECB(J)
ENDDO
ENDDO
RETURN
END
All reactions