C C USER SUBROUTINE TO CONTROL INCREMENTATION FOR FREE ROLLING subroutine urdfil(lstop,lovrwrt,kstep,kinc,pnewdt,time) c c include 'aba_param.inc' c dimension array(513),jrray(nprecd,513),time(*),lout(2) equivalence (array(1),jrray(1,1)) parameter (nmany = 999999,zero=0.0d0,small=1.0d-6,cmax=2.d-2, * cmax1=8.d-1,one=1.0d0) parameter (omega=10.d0,ftol=10.d0) common /cfroll1/ omega1, omega0, temp1, temp0 common /cfroll2/ kflag common /cfroll3/ time0, time1 C C ftol - the tolerance for zero force C cmax - max increment in rolling speed allowed C omega - upper bound estimate for the free-rolling velocity C kflag - flag set when free-rolling soln. is obtained C c write(7,*) "Entering urdfil", kstep, kinc, pnewdt, kflag if (kflag .eq. 1) return if (kinc .eq. 0) return C C SPECIFY NODE NUMBER AND DEGREE OF FREEDOM TO MONITOR C ldofusr = 5 nodeusr = 199 otol = cmax*omega write(7,1000) C C INITIALLY INCREMENT SOLUTION WITH SMALL STEPS C if (kstep.lt.3 .or. kinc .eq. 1) then pnewdt = cmax otol1 = cmax1*omega time0 = time(1) time1 = time(1)+cmax1 omega0 = otol1 omega1 = omega0+otol return end if C C FIND PREVIOUS INCREMENT do k1=kinc-1,1,-1 call dbfile(2,array,jrcd) call posfil(kstep,k1,array,jrcd) if (jrcd .eq. 0) goto 10 end do C EXIT IF RESULTS FOR PREVIOUS INCREMENT CANNOT BE FOUND call xit return 10 continue lproc = jrray(1,7) if (lproc .ne. 85) then C EXIT IF THE PROCEDURE IS NOT SST call xit return end if time0 = array(4) kinc0 = jrray(1,9) omega0 = temp0 C C FIND THE TORQUE AT PREVIOUS INCREMENT torq0 = zero lfnd0 = 0 do k1=1,nmany call dbfile(0,array,jrcd) if (jrcd .ne. 0) goto 100 key = jrray(1,2) if (key .eq. 104) then node = jrray(1,3) if (node .eq. nodeusr) then torq0 = array(3+ldofusr) lfnd0 = 1 goto 100 end if else if (key .eq. 2001) then goto 100 end if end do C EXIT IF END-OF-FILE MARKER NOT FOUND call xit return 100 continue if (lfnd0 .eq. 0) then lout(1) = nodeusr lout(2) = ldofusr C EXIT IF CANNOT FIND FORCES FOR NODE call xit return end if C C FIND CURRENT INCREMENT call posfil(kstep,kinc,array,jrcd) if (jrcd .ne. 0) then C EXIT IF RESULTS FOR CURRENT INCREMENT CANNOT BE FOUND call xit return end if time1 = array(4) kinc1 = jrray(1,9) omega1 = temp1 C C FIND THE TORQUE AT THIS INCREMENT torq1 = zero lfnd1 = 0 do k1=1,nmany call dbfile(0,array,jrcd) if (jrcd .ne. 0) goto 200 key = jrray(1,2) if (key .eq. 104) then node = jrray(1,3) if (node .eq. nodeusr) then torq1 = array(3+ldofusr) lfnd1 = 1 goto 200 end if end if end do C EXIT IF END-OF-FILE MARKER NOT FOUND call xit return 200 continue if (lfnd1 .eq. 0) then lout(1) = nodeusr lout(2) = ldofusr C EXIT IF CANNOT FIND FORCES FOR NODE call xit return end if C C ESTIMATE ANGULAR VELOCITY TO OBTAIN FREE ROLLING denom = torq1-torq0 write(7,1001) time0,torq0,time1,torq1 if (abs(torq1) .le. ftol) then write (7,1004) kflag = 1 pnewdt= one - time(1) return else if (abs(denom) .lt. small*abs(torq1)) then write (7,1002) kinc1, kinc0 domega = otol else omegafr = (torq1*omega0-torq0*omega1)/denom write (7,1003) omegafr domega = omegafr-omega1 if (abs(domega) .gt. abs(otol)) then write (7,1005) otol domega = otol end if end if omega0 = omega1 omega1 = omega1+domega pnewdt = cmax time0 = time(1) time1 = time(1)+pnewdt write(7,1006) time0, omega1, omega0 C 1000 format(//,15x,'FREE ROLLING SOLUTION CONTROL') 1001 format(/,2x,'TORQUE AT TIME ',1pg12.5,' = ',1pg12.5, * /,2x,'TORQUE AT TIME ',1pg12.5,' = ',1pg12.5) 1002 format(/,2x,'TORQUE AT CURRENT INCREMENT (',i5, * ') AND PREVIOUS INCREMENT (',i5,') ARE THE SAME.', * /,2x,'CONTINUE WITH REGULAR INCREMENTATION.',/) 1003 format(2x,'FREE ROLLING ESTIMATED TO OCCUR AT SPEED ',1pg12.5) 1004 format(/,2x,'FREE ROLLING SOLUTION OBTAINED.', * ' END STEP.',//) 1005 format(2x,'ESTIMATED INCREMENT TO OBTAIN FREE ROLLING', * ' EXCEEDS ',1pg12.5, * /,2x,'CONTINUE WITH REGULAR INCREMENTATION.',/) 1006 format(/2x,'END CALCULATIONS AT TIME',1pg12.5, * 'NEW ANGULAR VELOCITY ',1pg12.5, * 'PREVIOUS ESTIMATE ',1pg12.5,/) C return end C C User subroutine to define motion for free-rolling subroutine umotion(u,kstep,kinc,time,node,jdof) c include 'aba_param.inc' c dimension time(2) parameter (nmany = 999999,zero=0.0d0,small=1.0d-6,cmax=2.d-2, * cmax1=8.d-1,one=1.0d0) parameter (omega=10.d0,ftol=10.d0) common /cfroll1/ omega1, omega0, temp1, temp0 common /cfroll2/ kflag common /cfroll3/ time0, time1 C if (kinc .eq. 0) return if (kinc .eq. 1) then kflag = 0 u = cmax1*omega temp0=omega temp1=u return end if C C CHECK IF FREE ROLLING SOLUTION IS OBTAINED C if (kflag .eq. 1) then u = omega1 return else C C CHECK FOR CUTBACK OR INCREASE IN INCREMENT C if (time(1) .eq. time1) then u = omega1 else u = omega0+(omega1-omega0)/(time1-time0)*(time(1)-time0) end if temp0=omega0 temp1=u end if C return end subroutine uexternaldb(lop,lrestart,time,dtime,kstep,kinc) C INCLUDE 'ABA_PARAM.INC' INCLUDE 'mpif.h' dimension time(2) integer ABA_COMM_WORLD, GETCOMMUNICATOR integer rank, ierr, myrank, numprocs integer dest, source, length, tag integer stat(MPI_STATUS_SIZE) common /cfroll1/ omega1, omega0, temp1, temp0 common /cfroll2/ kflag common /cfroll3/ time0, time1 if ( lop.eq.1 ) then ABA_COMM_WORLD = GETCOMMUNICATOR() if (ABA_COMM_WORLD.ne.0) then tag1 = 1 tag2 = 2 tag3 = 3 length1 = 1 length2 = 2 length4 = 4 call MPI_COMM_RANK(ABA_COMM_WORLD, myrank,ierr) call MPI_COMM_SIZE(ABA_COMM_WORLD, numprocs, ierr) if (myrank .eq. 0) then do dest = 1, numprocs-1 call MPI_SEND(omega1, length4, MPI_DOUBLE_PRECISION, dest, tag1, * ABA_COMM_WORLD,ierr) call MPI_SEND(time0, length2, MPI_DOUBLE_PRECISION, dest, tag2, * ABA_COMM_WORLD,ierr) call MPI_SEND(kflag, length1, MPI_INTEGER, dest, tag3, * ABA_COMM_WORLD,ierr) end do else source = 0 call MPI_RECV(omega1, length4, MPI_DOUBLE_PRECISION, source, tag1, * ABA_COMM_WORLD,stat,err) call MPI_RECV(time0, length2, MPI_DOUBLE_PRECISION, source, tag2, * ABA_COMM_WORLD,stat,err) call MPI_RECV(kflag, length1, MPI_INTEGER, source, tag3, * ABA_COMM_WORLD,stat,err) end if end if end if return end