compart.f

对工业生产过程结晶过程的一个仿真程序软件包

     &         +f(i+1,j+1)*dble(i))/4.0D0
     &         *delR1**2*delR2
          mu01(1)=mu01(1)
     &         +(f(i,j)*dble(j-1)
     &         +f(i,j+1)*dble(j)
     &         +f(i+1,j)*dble(j-1)
     &         +f(i+1,j+1)*dble(j))/4.0D0
     &         *delR1*delR2**2
          mu11(1)=mu11(1)
     &         +(f(i,j)*dble(i-1)*dble(j-1)
     &         +f(i,j+1)*dble(i-1)*dble(j)
     &         +f(i+1,j)*dble(i)*dble(j-1)
     &         +f(i+1,j+1)*dble(i)*dble(j))/4.0D0
     &         *delR1**2*delR2**2
          mu20(1)=mu20(1)
     &         +(f(i,j)*dble(i-1)**2
     &         +f(i,j+1)*dble(i-1)**2
     &         +f(i+1,j)*dble(i)**2
     &         +f(i+1,j+1)*dble(i)**2)/4.0D0
     &         *delR1**3*delR2
          mu02(1)=mu02(1)
     &         +(f(i,j)*dble(j-1)**2
     &         +f(i,j+1)*dble(j)**2
     &         +f(i+1,j)*dble(j-1)**2
     &         +f(i+1,j+1)*dble(j)**2)/4.0D0
     &         *delR1*delR2**3
          mu21(1)=mu21(1)
     &         +(f(i,j)*dble(i-1)**2*dble(j-1)
     &         +f(i,j+1)*dble(i-1)**2*dble(j)
     &         +f(i+1,j)*dble(i)**2*dble(j-1)
     &         +f(i+1,j+1)*dble(i)**2*dble(j))/4.0D0
     &         *delR1**3*delR2**2 
          mu12(1)=mu12(1)
     &         +(f(i,j)*dble(i-1)*dble(j-1)**2
     &         +f(i,j+1)*dble(i-1)*dble(j)**2
     &         +f(i+1,j)*dble(i)*dble(j-1)**2
     &         +f(i+1,j+1)*dble(i)*dble(j)**2)/4.0D0
     &         *delR1**2*delR2**3
          mu22(1)=mu22(1)
     &         +(f(i,j)*dble(i-1)**2*dble(j-1)**2
     &         +f(i,j+1)*dble(i-1)**2*dble(j)**2
     &         +f(i+1,j)*dble(i)**2*dble(j-1)**2
     &         +f(i+1,j+1)*dble(i)**2*dble(j)**2)/4.0D0
     &         *delR1**3*delR2**3     
          mu30(1)=mu30(1)
     &         +(f(i,j)*dble(i-1)**3
     &         +f(i,j+1)*dble(i-1)**3
     &         +f(i+1,j)*dble(i)**3
     &         +f(i+1,j+1)*dble(i)**3)/4.0D0
     &         *delR1**4*delR2      
          mu31(1)=mu31(1)
     &         +(f(i,j)*dble(i-1)**3*dble(j-1)
     &         +f(i,j+1)*dble(i-1)**3*dble(j)
     &         +f(i+1,j)*dble(i)**3*dble(j-1)
     &         +f(i+1,j+1)*dble(i)**3*dble(j))/4.0D0
     &         *delR1**4*delR2**2
*          volume(1)=volume(1)+(dble(i)**3/3.0D0+(dble(j)-dble(i))*
*     &         dble(i)**2)*f(i,j)          
       ENDDO
      ENDDO    
      DO i=meshi2,meshi3
         DO j=meshj2,meshj3
            mu00(1)=mu00(1)+(f(i,j)+f(i,j+1)+f(i+1,j)+f(i+1,j+1))/4.0D0
     &           *delR1*delR2     
            mu10(1)=mu10(1)
     &           +(f(i,j)*dble(i-1)
     &           +f(i,j+1)*dble(i-1)
     &           +f(i+1,j)*dble(i)
     &           +f(i+1,j+1)*dble(i))/4.0D0
     &           *delR1**2*delR2
            mu01(1)=mu01(1)
     &           +(f(i,j)*dble(j-1)
     &           +f(i,j+1)*dble(j)
     &           +f(i+1,j)*dble(j-1)
     &           +f(i+1,j+1)*dble(j))/4.0D0
     &           *delR1*delR2**2
            mu11(1)=mu11(1)
     &           +(f(i,j)*dble(i-1)*dble(j-1)
     &           +f(i,j+1)*dble(i-1)*dble(j)
     &           +f(i+1,j)*dble(i)*dble(j-1)
     &           +f(i+1,j+1)*dble(i)*dble(j))/4.0D0
     &           *delR1**2*delR2**2
            mu20(1)=mu20(1)
     &           +(f(i,j)*dble(i-1)**2
     &           +f(i,j+1)*dble(i-1)**2
     &           +f(i+1,j)*dble(i)**2
     &           +f(i+1,j+1)*dble(i)**2)/4.0D0
     &           *delR1**3*delR2
            mu02(1)=mu02(1)
     &           +(f(i,j)*dble(j-1)**2
     &           +f(i,j+1)*dble(j)**2
     &           +f(i+1,j)*dble(j-1)**2
     &           +f(i+1,j+1)*dble(j)**2)/4.0D0
     &           *delR1*delR2**3
            mu21(1)=mu21(1)
     &           +(f(i,j)*dble(i-1)**2*dble(j-1)
     &           +f(i,j+1)*dble(i-1)**2*dble(j)
     &           +f(i+1,j)*dble(i)**2*dble(j-1)
     &           +f(i+1,j+1)*dble(i)**2*dble(j))/4.0D0
     &           *delR1**3*delR2**2
            mu12(1)=mu12(1)
     &           +(f(i,j)*dble(i-1)*dble(j-1)**2
     &           +f(i,j+1)*dble(i-1)*dble(j)**2
     &           +f(i+1,j)*dble(i)*dble(j-1)**2
     &           +f(i+1,j+1)*dble(i)*dble(j)**2)/4.0D0
     &           *delR1**2*delR2**3
            mu22(1)=mu22(1)
     &           +(f(i,j)*dble(i-1)**2*dble(j-1)**2
     &           +f(i,j+1)*dble(i-1)**2*dble(j)**2
     &           +f(i+1,j)*dble(i)**2*dble(j-1)**2
     &           +f(i+1,j+1)*dble(i)**2*dble(j)**2)/4.0D0
     &           *delR1**3*delR2**3     
            mu30(1)=mu30(1)
     &           +(f(i,j)*dble(i-1)**3
     &           +f(i,j+1)*dble(i-1)**3
     &           +f(i+1,j)*dble(i)**3
     &           +f(i+1,j+1)*dble(i)**3)/4.0D0
     &           *delR1**4*delR2      
            mu31(1)=mu31(1)
     &           +(f(i,j)*dble(i-1)**3*dble(j-1)
     &           +f(i,j+1)*dble(i-1)**3*dble(j)
     &           +f(i+1,j)*dble(i)**3*dble(j-1)
     &           +f(i+1,j+1)*dble(i)**3*dble(j))/4.0D0
     &           *delR1**4*delR2**2
*          volume(1)=volume(1)+(dble(i)**3/3.0D0+(dble(j)-dble(i))*
*     &         dble(i)**2)*f(i,j)
         ENDDO
      ENDDO
*     Initial volume
      volume(1)=mu21(1)-2.0D0/3.0D0*mu30(1)
*     Current nucleation rate
      Brate(1)=birth(Conc(1),Concs(1), volume(1))
*     Current growth rate
      G1rate(1)=growth1(Conc(1), Concs(1))
      G2rate(1)=growth2(Conc(1), Concs(1))
*     Mesh boundary
      mratei=mratei+G1rate(1)*delT
      mratej=mratej+G2rate(1)*delT
      meshi1=starti1+DNINT(mratei)
      meshi2=starti2+DNINT(mratei)
      meshi3=starti3+DNINT(mratei)
      meshj1=startj1+DNINT(mratej)
      meshj2=startj2+DNINT(mratej)
      meshj3=startj3+DNINT(mratej)

*     Write f to file
      nfile=10
      nUnit=nfile+myrank
      WRITE(fname(2:3), '(I2)')nfile
      WRITE(fname(5:6), '(I2)')myrank+10
      WRITE(logname(5:6),'(I2)')myrank+10
      OPEN(UNIT=nUnit, FILE=fname, FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=200+myrank, FILE=logname, FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      DO i=1, totR1pts, skipR1
         DO j=1, totR2pts, skipR2
               WRITE(nUnit,30) i, j, time(1), f(i,j)
         ENDDO
      ENDDO
      CLOSE(UNIT=nUnit)      
      m=1
      WRITE(myrank+200,120)'tau is',tau

*     Simulation      
      DO k=2,totTpts
*        Current time
         time(k)=time(k-1)+delT
*        Current temperature
         T(k)=temp(time(k))
*        Current saturated Conc
         Concs(k)=Csat(T(k))
         IF(myrank.eq.0)THEN
            PRINT*, 'k=',k
         END IF         
*        Modify the mesh boundary
         IF(numprocs.GT.1)then
            CALL MPI_ALLREDUCE(meshi3,tmpij,1,MPI_INTEGER,MPI_MAX,
     &           MPI_COMM_WORLD,ierr)
            meshi3=tmpij
            CALL MPI_ALLREDUCE(meshj3,tmpij,1,MPI_INTEGER,MPI_MAX,
     &           MPI_COMM_WORLD,ierr)
            meshj3=tmpij
*           Compute the output from each compartment
            DO j=1, meshj3-1
               DO i=1, meshi3-1
                  Wu=1.0D0-alpha1/DMAX1(dble(meshi3),dble(meshj3))
     &                 *dble(i)-alpha2/DMAX1(dble(meshi3),dble(meshj3))
     &                 *dble(j)
                  Wd=1.0D0+alpha1/DMAX1(dble(meshi3),dble(meshj3))
     &                 *dble(i)+alpha2/DMAX1(dble(meshi3),dble(meshj3))
     &                 *dble(j)     
                  IF(myrank.EQ.numprocs-1) THEN
                     IF(f(i,j).GT.0.0D0)THEN
                        ftemp5(i,j)=Wd*f(i,j)*tau*delT
                        IF(f(i,j).LT.ftemp5(i,j))ftemp5(i,j)=f(i,j)
                     ELSE
                        ftemp5(i,j)=0.0D0
                     END IF
                     f(i,j)=f(i,j)-ftemp5(i,j)
                  ELSE IF(myrank.EQ.0) THEN
                        IF(f(i,j).GT.0.0D0)THEN
                           ftemp4(i,j)=Wu*f(i,j)*tau*delT
                           IF(f(i,j).LT.ftemp4(i,j))ftemp4(i,j)=f(i,j)        
                        ELSE
                           ftemp4(i,j)=0.0D0
                        END IF
                        f(i,j)=f(i,j)-ftemp4(i,j)
                  ELSE
                     IF(f(i,j).GT.0.0D0)THEN
                        ftemp5(i,j)=Wd*f(i,j)*tau*delT
                        ftemp4(i,j)=Wu*f(i,j)*tau*delT
                        IF(f(i,j).LT.(ftemp5(i,j)+ftemp4(i,j)))THEN
                           ftemp5(i,j)=ftemp5(i,j)/(ftemp5(i,j)+
     &                          ftemp4(i,j))*f(i,j)
                           ftemp4(i,j)=ftemp4(i,j)/(ftemp5(i,j)+
     &                          ftemp4(i,j))*f(i,j)
                        END IF
                     ELSE
                        ftemp4(i,j)=0.0D0
                        ftemp5(i,j)=0.0D0
                     END IF
                     f(i,j)=f(i,j)-ftemp5(i,j)-ftemp4(i,j)
                  END IF
               END DO
            END DO 
*           Exchange the concentration between nbrs
            CALL MPI_IRECV(b_conc,1,MPI_DOUBLE_PRECISION,nbrbottom,
     &           0,com1d,req(1),ierr)
            CALL MPI_IRECV(t_conc,1,MPI_DOUBLE_PRECISION,nbrtop,
     &           1,com1d,req(2),ierr)
            CALL MPI_ISEND(conc(k-1),1,MPI_DOUBLE_PRECISION,nbrtop,
     &           0,com1d,req(3),ierr)
            CALL MPI_ISEND(conc(k-1),1,MPI_DOUBLE_PRECISION,nbrbottom,
     &           1,com1d,req(4),ierr)
            CALL MPI_WAITALL(4,req,status_array,ierr)
*           Exchange f between nbrs
            DO i=1,iter
               CALL MPI_IRECV(ftemp2(-1,(i-1)*numcols-1),maxf1,
     &              MPI_DOUBLE_PRECISION, nbrbottom,2,com1d,req(1),ierr)
               CALL MPI_IRECV(ftemp3(-1,(i-1)*numcols-1),maxf1,
     &              MPI_DOUBLE_PRECISION,nbrtop,3,com1d,req(2),ierr)
               CALL MPI_ISEND(ftemp4(-1,(i-1)*numcols-1),maxf1,
     &              MPI_DOUBLE_PRECISION,nbrtop,2,com1d,req(3),ierr)
               CALL MPI_ISEND(ftemp5(-1,(i-1)*numcols-1),maxf1,
     &              MPI_DOUBLE_PRECISION,nbrbottom,3,com1d,req(4),ierr)
               CALL MPI_WAITALL(4,req,status_array,ierr)
            ENDDO  
         END IF 
*        Compute current concentration
         IF(numprocs.GT.1)THEN
            IF(myrank.EQ.numprocs-1) THEN
               conc(k)=conc(k-1)-dble(2)*alpha*delT*G1rate(k-1)
     &              *(mu11(k-1)-mu20(k-1))
     &              -alpha*delT*G2rate(k-1)*mu20(k-1)
     &              +(b_conc-conc(k-1))*tau*delT
            ELSE IF(myrank.EQ.0) THEN
               conc(k)=conc(k-1)-dble(2)*alpha*delT*G1rate(k-1)
     &              *(mu11(k-1)-mu20(k-1))
     &              -alpha*delT*G2rate(k-1)*mu20(k-1)
     &              +(t_conc-conc(k-1))*tau*delT 
            ELSE
               conc(k)=conc(k-1)-dble(2)*alpha*delT*G1rate(k-1)
     &              *(mu11(k-1)-mu20(k-1))
     &              -alpha*delT*G2rate(k-1)*mu20(k-1)  
     &              +(t_conc-2.0D0*conc(k-1)+b_conc)*tau*delT
            END IF
         ELSE
            conc(k)=conc(k-1)-dble(2)*alpha*delT*G1rate(k-1)
     &           *(mu11(k-1)-mu20(k-1))
     &           -alpha*delT*G2rate(k-1)*mu20(k-1)

         END IF
         IF(conc(k).GE.1.or.conc(k).LE.0)then 
            PRINT*, 'myrank=', myrank
            PRINT*, 'Concentration is negative. Abort'
            WRITE(200+myrank,100)'conc is negative abort' 
            call MPI_ABORT(MPI_COMM_WORLD,0,ierr)
            GOTO 1000
         END IF
*       Current super concentration
        superConc(k)=(conc(k)-Concs(k))/Concs(k)

        IF(numprocs.GT.1)THEN
*          Add the inputs from nbrs
           DO j=1, meshj3
              DO i=1, meshi3
                 IF(myrank.EQ.numprocs-1) THEN
                    f(i,j)=f(i,j)+ftemp2(i,j) 
                 ELSE IF(myrank.EQ.0) THEN
                  f(i,j)=f(i,j)+ftemp3(i,j)
                 ELSE
                  f(i,j)=f(i,j)+ftemp2(i,j)+ftemp3(i,j)
                 END IF       
                 IF (f(i,j).LT.-0.5D0)THEN
                    PRINT*, 'myrank=', myrank
                    PRINT*, 'f is negative. See log file'
                    WRITE(200+myrank,100)'f is negative '
                    write(200+myrank,30)i,j,time(k), f(i,j)
                 END IF
              ENDDO
           ENDDO          
        END IF        
        f(1,1)=f(1,1)+Brate(k-1)*delT/delR1/delR2
*       Compute Current population density f   
        vei=delT*G1rate(k-1)/delR1
        vej=delT*G2rate(k-1)/delR2
        IF(numprocs.GT.1)THEN
*         compute f in j-direction          
           DO i=1, meshi3-1
              DO j=1, meshj3-1              
                 if((f(i,j+1)-f(i,j)).EQ.0.0D0)then
                    j2theta=1.0D10*(f(i,j)-f(i,j-1))
                 else
                    j2theta=(f(i,j)-f(i,j-1))/(f(i,j+1)-f(i,j)) 
                 endif
                 if(((f(i,j)-f(i,j-1)).EQ.0.0D0))then
                    j1theta=1.0D10*(f(i,j-1)-f(i,j-2))
                 else
                    j1theta=(f(i,j-1)-f(i,j-2))/(f(i,j)-f(i,j-1))
                 endif
                 j1phi=(DABS(j1theta) +j1theta)/(1.0D0+DABS(j1theta))
                 j2phi=(DABS(j2theta) +j2theta)/(1.0D0+DABS(j2theta))                
*                High resolution             
                 ftemp1(i,j)=f(i,j)-(vej-0.5D0*vej*(1.0D0-vej)*j1phi)
     &                *(f(i,j)-f(i,j-1))-0.5D0*vej*(1.0D0-vej)*j2phi
     &                *(f(i,j+1)-f(i,j))