!DEC$ FREEFORM

INCLUDE 'mkl_service.f90'
INCLUDE 'Additional_routines.f'
INCLUDE 'Type_trainingparameters.f'
!Abaqus utility routines that have to be replaced
INCLUDE 'Abaqus_Routines.f'

PROGRAM UMAT_Driver
    ! UMAT_Driver is a program that mimics Abaqus in calling an arbitrary strain or
    ! stress test trajectory; it calls a material routine as Abaqus throught the UMAT
    ! interface and defines all variables accordingly; currently the thermal variables
    ! have no meaning; UEXTERNALDB must be present even if not needed
    
    USE ABQINTERFACE
    USE type_trainingparameters
    USE linear_interpolation
    USE DSTRAN_DFGRD1_exponential_map
    
    IMPLICIT NONE
    
    INTEGER,DIMENSION(:),ALLOCATABLE:: ipiv
    INTEGER:: info
    REAL(KIND=AbqRK):: STRESS_criterion
    INTEGER(KIND=4),DIMENSION(3):: simulation_start_time,simulation_end_time
    INTEGER:: i,k,l,m,o
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE::timesteps,timesteps_temp
    INTEGER,DIMENSION(:),ALLOCATABLE::increments,increments_temp
    LOGICAL:: STEP_finished
    CHARACTER(LEN=32):: arg
    INTEGER:: nCPUs !number of processors to parallize UMATs
    CHARACTER(LEN=256)::  OUTDIR
    INTEGER:: LENOUTDIR
    REAL(KIND=AbqRK),DIMENSION(3,3),PARAMETER:: Ident=reshape([1.0_AbqRK,0.0_AbqRK,0.0_AbqRK,&
                                                               0.0_AbqRK,1.0_AbqRK,0.0_AbqRK,&
                                                               0.0_AbqRK,0.0_AbqRK,1.0_AbqRK],[3,3])
    !---------------------variables from Abaqus interface----------------------
    INTEGER(KIND=AbqIK)::NPT,KINC,NOEL
    INTEGER(KIND=AbqIK)::KSPT=1,LAYER=1 !has currently no meaning
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: STRESS,STRESS_t,STRAN,DSTRAN,&
                                                STRESS_set,STRESS_rel
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: STATEV,STATEV_t
    REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE:: DDSDDE
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE::DDSDDT,DRPLDE
    REAL(KIND=AbqRK),DIMENSION(2):: TIME
    REAL(KIND=AbqRK),DIMENSION(1):: PREDEF,DPRED
    REAL(KIND=AbqRK),DIMENSION(3):: COORDS
    REAL(KIND=AbqRK),DIMENSION(3,3):: DFGRD0,DFGRD1,DROT
    INTEGER,DIMENSION(4):: JSTEP
    REAL(KIND=AbqRK):: DTIME,TEMP,DTEMP,SSE,SPD,SCD,CELENT,DRPLDT,PNEWDT,RPL
    !--------------------hard coded simulation parameters---------------------
    INTEGER:: NPT_max=27 !(maximal number of IP's of Abaqus Elements)
    INTEGER:: kinc_max=8
    REAL(KIND=AbqRK):: tol_stress_iter=0.0001_AbqRK ! tolerance in interation for
                                                    ! DSTRAN when STRESS is given
    
    WRITE(*,*) 'Running UMAT_Driver Version 0.1'
    
    !------------------get the simulation arguments-----------------------------
    nCPUs=1 !default values
    training%Jobname=''
    DO i=1,COMMAND_ARGUMENT_COUNT() !get arguments from command line
        CALL GET_COMMAND_ARGUMENT(i,arg)
        IF (arg(1:5)=='cpus=') THEN
            READ(arg(6:LEN_TRIM(arg)),*) nCPUs
        ELSE IF (arg(1:4)=='job=') THEN
            READ(arg(5:LEN_TRIM(arg)),*) training%Jobname
            DO k=1,LEN_TRIM(training%Jobname)
                IF (training%Jobname(k:k)=='.') THEN
                    IF (training%Jobname(k:k+3)=='.inp') THEN
                        training%Jobname=training%Jobname(1:k-1)
                        WRITE(*,*) '.inp has been removed from the Jobname.'
                    END IF
                END IF
            END DO
        ELSE
            STOP 'Error: unrecognized keyword. Allowed are: cpus, job'
        END IF
    END DO
    
    IF (training%Jobname=='') STOP 'Error: Jobname must be specified'
    !---------------------------------------------------------------------------
    
    !open files to write information to
    CALL GETOUTDIR(OUTDIR,LENOUTDIR)
    OPEN(unit=6,file=OUTDIR(1:LENOUTDIR)//'/'//TRIM(training%Jobname)//'.dat',action='WRITE',position='append',status='replace')
    OPEN(unit=7,file=OUTDIR(1:LENOUTDIR)//'/'//TRIM(training%Jobname)//'.msg',action='WRITE',position='append',status='replace')
    
    !write the current time to the Messagefile
    CALL ITIME(simulation_start_time)
    WRITE(7,fmt='(A29,I2,A1,I2,A1,I2)') 'Starting simulation at time: ',simulation_start_time(1),':',simulation_start_time(2),':',simulation_start_time(3)
    
    !at first get the simulation input parameters
    CALL training%read_data()
    
    !call UEXTERNALDB at the start of the Analysis
    TIME=0.0_AbqRK; DTIME=0.0_AbqRK;
    CALL UEXTERNALDB(0,0,TIME,DTIME,0,0)
    
    call omp_set_num_threads(nCPUs) !set up the parallel envionment
    
    !$omp parallel do &
    !$omp private(NPT,NOEL,KINC,LAYER,TIME,DTIME,STRESS,STRESS_t,STRAN,DSTRAN,STATEV,STATEV_t,&
      DDSDDE,DDSDDT,DRPLDE,PREDEF,DPRED,COORDS,DFGRD0,DFGRD1,DROT,TEMP,DTEMP,SSE,SPD,&
      SCD,CELENT,DRPLDT,PNEWDT,RPL,JSTEP,ipiv,info,STRESS_set,STRESS_rel,o,i,k,l,m,&
      STRESS_criterion,timesteps,timesteps_temp,increments,increments_temp,STEP_finished) &
    !$omp shared(training,tol_stress_iter,kinc_max,NPT_max)
    
    DO o=1,SIZE(training%Step) !go over all training steps in parallel
    
        IF (.NOT.(ALLOCATED(STRESS))) ALLOCATE(STRESS(training%NTENS),STRESS_t(training%NTENS),STRAN(training%NTENS),STRESS_set(training%NTENS),&
                                               STRESS_rel(training%NTENS),DSTRAN(training%NTENS),DDSDDE(training%NTENS,training%NTENS),&
                                               DDSDDT(training%NTENS),DRPLDE(training%NTENS),STATEV(training%NSTATV),STATEV_t(training%NSTATV))
        
        !initialization at the beginning of a timestep
        IF (ALLOCATED(timesteps)) DEALLOCATE(timesteps); ALLOCATE(timesteps(0))
        IF (ALLOCATED(increments)) DEALLOCATE(increments); ALLOCATE(increments(0))
        STEP_finished=.FALSE.
        NOEL=1+(o-1)/NPT_max
        NPT=MOD(o-1,NPT_max)+1
        JSTEP=0; JSTEP(3)=training%Step(o)%nlgeom
        KINC=1
        TIME(1)=0.0_AbqRK; TIME(2)=0.0_AbqRK
        DFGRD0=Ident; DFGRD1=Ident; DROT=Ident; STRAN=0.0_AbqRK; DSTRAN=0.0_AbqRK; STRESS_t=0.0_AbqRK; STATEV_t=0.0_AbqRK
        DDSDDE=0.0_AbqRK; DDSDDT=0.0_AbqRK; DRPLDE=0.0_AbqRK; PREDEF=0.0_AbqRK; DPRED=0.0_AbqRK; COORDS=0.0_AbqRK; RPL=0.0_AbqRK
        TEMP=0.0_AbqRK; DTEMP=0.0_AbqRK; SSE=0.0_AbqRK; SPD=0.0_AbqRK; SCD=0.0_AbqRK; CELENT=0.0_AbqRK; DRPLDT=0.0_AbqRK
        
        DO WHILE (TIME(2)<=training%Step(o)%time_end)
            
            CALL MOVE_ALLOC(timesteps,timesteps_temp); ALLOCATE(timesteps(SIZE(timesteps_temp)+1)); timesteps(1:SIZE(timesteps_temp))=timesteps_temp
            timesteps(SIZE(timesteps))=TIME(2)
            CALL MOVE_ALLOC(increments,increments_temp); ALLOCATE(increments(SIZE(increments_temp)+1)); increments(1:SIZE(timesteps_temp))=increments_temp
            increments(SIZE(increments))=0
            
            !signalize beginning of a time increment
            CALL UEXTERNALDB(1,-NPT_max*NOEL-NPT,[TIME(2),TIME(2)],DTIME,0,KINC)
            
            !find the next value of strain resp. stress -> assume multilinear strain (stress) curve
            IF (training%Step(o)%simulation_type=='strain') THEN
                IF (training%Step(o)%strain_measure=='DSTRAN') THEN
                    CALL linear_interpol(TIME(2),training%Step(o)%timesteps,training%Step(o)%STRAN,DSTRAN)
                    DSTRAN=DSTRAN-STRAN
                    IF (JSTEP(3)==1) CALL DSTRAN_DFGRD1_map(DSTRAN,DFGRD0,DFGRD1,training%NDI,training%NSHR) !in geom. nonl. setting compute DFGRD1
                ELSE IF (training%Step(o)%strain_measure=='DFGRD1') THEN
                    CALL linear_interpol(TIME(2),training%Step(o)%timesteps,training%Step(o)%DFGRD1,DFGRD1)
                    DFGRD1=DFGRD1+Ident
                END IF
            ELSE
                CALL linear_interpol(TIME(2),training%Step(o)%timesteps,training%Step(o)%STRESS,STRESS_set)
                DSTRAN=0.0_AbqRK; DFGRD1=DFGRD0 !initialization
            END IF
            
            DTIME=TIME(2)-TIME(1) !time increment
            
            DO m=1,kinc_max
                
                increments(SIZE(increments))=increments(SIZE(increments))+1 !count the number of increments
                
                STRESS=STRESS_t
                STATEV=STATEV_t
                PNEWDT=1.0_AbqRK
                
                ! Call the UMAT like Abaqus
                CALL UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,&
                        DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,&
                        training%Material_Name,training%NDI,training%NSHR,training%NTENS,&
                        training%NSTATV,training%PROPS,size(training%PROPS),COORDS,DROT,PNEWDT,&
                        CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
                
                IF (training%Step(o)%simulation_type=='strain') THEN
                    EXIT
                ELSE !iterate for the DSTRAN resp. DFGRD1 in a stress driven context
                
                    IF (PNEWDT<1.0_AbqRK) EXIT !if convergence in UMAT not reached, don't iterate for STRESS_set
                    
                    STRESS_rel=STRESS_set-STRESS !difference between actual and stress defined by user
                    IF (MAXVAL(ABS(STRESS_set))>0.00000001_AbqRK) THEN
                        STRESS_criterion=MAXVAL(ABS(STRESS_rel))/MAXVAL(ABS(STRESS_set))
                    ELSE
                        STRESS_criterion=MAXVAL(ABS(STRESS_rel))/MAXVAL(ABS(training%Step(o)%STRESS))
                    END IF
                    
                    IF (STRESS_criterion>0.0001_AbqRK) THEN
                        
                        !convergence not reached
                        IF (m==kinc_max) THEN
                            PNEWDT=0.5_AbqRK
                            EXIT
                        END IF
                        
                        !compute DSTRAN
                        ALLOCATE(ipiv(training%NTENS))
                        CALL dgesv(training%NTENS,1,DDSDDE,training%NTENS,ipiv,STRESS_rel,training%NTENS,info)
                        DEALLOCATE(ipiv)
                        DSTRAN=DSTRAN+STRESS_rel
                        
                        IF (JSTEP(3)==1) CALL DSTRAN_DFGRD1_map(DSTRAN,DFGRD0,DFGRD1,training%NDI,training%NSHR) !in geom. nonl. setting compute DFGRD1
                
                    ELSE
                    
                        EXIT
                        
                    END IF
                    
                END IF
                
            END DO
            
            !-----------------------------Time step control-----------------------------------
            
            IF (TIME(2)==0.0_AbqRK) THEN
                TIME(2)=training%Step(o)%time_inc_begin
            ELSE IF (PNEWDT<1.0_AbqRK) THEN !no convergence reached
                IF (((TIME(2)-TIME(1))*0.5_AbqRK)<training%Step(o)%time_inc_min) THEN
                    CALL UEXTERNALDB(3,-NPT_max*NOEL-NPT,[TIME(2),TIME(2)],DTIME,0,KINC) ! end of analysis
                    STEP_finished=.TRUE.
                ELSE
                    TIME(2)=TIME(1)+(TIME(2)-TIME(1))*0.5_AbqRK
                END IF
            ELSE IF ((ABS(TIME(2)-training%Step(o)%time_end)/training%Step(o)%time_end)<0.0000000001_AbqRK) THEN !convergence reached, last timestep reached
                CALL UEXTERNALDB(2,-NPT_max*NOEL-NPT,[TIME(2),TIME(2)],DTIME,0,KINC) !convergence of timestep
                CALL UEXTERNALDB(3,-NPT_max*NOEL-NPT,[TIME(2),TIME(2)],DTIME,0,KINC) !end of analysis
                STEP_finished=.TRUE.
            ELSE !convergence reached
                CALL UEXTERNALDB(2,-NPT_max*NOEL-NPT,[TIME(2),TIME(2)],DTIME,0,KINC) !convergence of timestep
                IF (DTIME*PNEWDT>training%Step(o)%time_inc_max) THEN !max inc reached
                    IF ((TIME(2)+training%Step(o)%time_inc_max)<training%Step(o)%time_end) THEN
                        TIME(2)=TIME(2)+training%Step(o)%time_inc_max
                        TIME(1)=TIME(1)+DTIME
                    ELSE !end of time reached
                        TIME(2)=training%Step(o)%time_end
                        TIME(1)=TIME(1)+DTIME
                    END IF
                ELSE !max inc not reached
                    IF ((TIME(2)+DTIME*PNEWDT)<training%Step(o)%time_end) THEN
                        TIME(2)=TIME(2)+DTIME*PNEWDT
                        TIME(1)=TIME(1)+DTIME
                    ELSE !end of time reached
                        TIME(2)=training%Step(o)%time_end
                        TIME(1)=TIME(1)+DTIME
                    END IF
                END IF
                STRESS_t=STRESS
                STATEV_t=STATEV
                STRAN=STRAN+DSTRAN
                DFGRD0=DFGRD1
                KINC=KINC+1
            END IF
            
            IF (STEP_finished) THEN
                !$omp critical
                CALL write_step_info_to_msg_file(timesteps,increments,o,training%Step(o)%StepName)
                !$omp end critical
                EXIT
            END IF
            
        END DO
    END DO
    !$omp end parallel do
    !call UEXTERNALDB at the very end of the Analysis
    CALL UEXTERNALDB(3,0,[0.0_AbqRK,0.0_AbqRK],DTIME,0,0)
    
    !write the current time and whole simulation time to the Messagefile
    CALL ITIME(simulation_end_time)
    WRITE(7,fmt='(A27,I2,A1,I2,A1,I2)') 'Ending simulation at time: ',simulation_end_time(1),':',simulation_end_time(2),':',simulation_end_time(3)
    WRITE(7,fmt='(A25,I)') 'Elapsed time in seconds: ',(simulation_end_time(1)-simulation_start_time(1))*60**2+&
                                                          (simulation_end_time(2)-simulation_start_time(2))*60+&
                                                          (simulation_end_time(3)-simulation_start_time(3))
    WRITE(7,fmt='(A21,I,A1)') 'Number of used CPUs: ',nCPUs,NEW_LINE('a')
    
    !close files that are opened to be filled with message information
    CLOSE(unit=6)
    CLOSE(unit=7)
    
    STOP 'Job COMPLETED successfully.'
    
    CONTAINS
    
        SUBROUTINE write_step_info_to_msg_file(timesteps,increments,Step_nbr,step_name)
        !at the end of a step write the information of a step to the .dat file
        
        IMPLICIT NONE
        
        REAL(KIND=AbqRK),DIMENSION(:),INTENT(IN)::timesteps
        INTEGER,DIMENSION(:),INTENT(IN)::increments
        CHARACTER(LEN=*),INTENT(IN)::step_name
        INTEGER,INTENT(IN):: Step_nbr
        INTEGER:: i
        
        WRITE(7,fmt='(A25)') '-------------------------'
        WRITE(7,*) 'Information for Step: '//TRIM(ADJUSTL(step_name))
        WRITE(7,fmt='(A4,T10,A4,T20,A4,T30,A4)') 'STEP','INCR','ITER','TIME'
        DO i=1,SIZE(timesteps)
            WRITE(7,fmt='(I4,T10,I4,T20,I4,T30,E11.4)') Step_nbr,i,increments(i),timesteps(i)
        END DO
        WRITE(7,fmt='(A25)') '-------------------------'
        
        END SUBROUTINE write_step_info_to_msg_file
            
END PROGRAM UMAT_Driver
