tran.2D.f
C#==============================================================================
C Diffusion in a 2-dimensional finite difference grid
C all inputs are vectors/matrices
C#==============================================================================
SUBROUTINE diff2D (Nx, Ny, C, Cxup, Cxdown, Cyup, Cydown, &
& Fxup, Fxdown, Fyup, Fydown, &
& axup, axdown, ayup, aydown, &
& BcXup, BcXdown, BcYup, BcYdown, &
& D_x, D_y, VF_x, VF_y, A_x, A_y, &
& dx, dxaux, dy, dyaux, FluxX, FluxY, dC)
IMPLICIT NONE
INTEGER Nx,Ny ! dimension of C
C input
DOUBLE PRECISION C(Nx,Ny)
C when Bc.. = 2; boundary concentrations should be a vector then, else a scalar
DOUBLE PRECISION Cxup(*),Cxdown(*),Cyup(*),Cydown(*)
C when Bc.. = 1; boundary fluxes should be a vector, else a scalar
DOUBLE PRECISION Fxup(*), Fxdown(*), Fyup(*), Fydown(*)
C when Bc.. = 4; should be a vector, else a scalar
DOUBLE PRECISION axup(*), axdown(*), ayup(*), aydown(*)
C diffusion coefficients, volume fractions, surface Areas
DOUBLE PRECISION D_x(Nx+1,Ny), D_y(Nx,Ny+1), VF_x(Nx+1,Ny), &
& VF_y(Nx,Ny+1), A_x(Nx+1,Ny),A_y(Nx,Ny+1)
C Grid sizes and auxillary distances (from mid to mid)
DOUBLE PRECISION dx(Nx,Ny),dxaux(Nx+1,Ny),dy(Nx,Ny),dyaux(Nx,Ny+1)
C Boundaries: 1= flux, 2 = conc, 3 = 0 grad, 4 = convection
INTEGER BcXup, BcXdown, BcYup, BcYdown
C Output: Fluxes in X- and Y-direction, rate of change
DOUBLE PRECISION FluxX(Nx+1, Ny), FluxY(Nx, Ny+1), dC(Nx,Ny)
C locals
DOUBLE PRECISION AVx, AVy
INTEGER I, J
C -------------------------------------------------------------------------------
C ## transport in X-direction ##
C First diffusion internal cells
DO I = 2, Nx
DO J = 1, Ny
FluxX(I,J) = -VF_x(I,J)*D_x(I,J)*(C(I,J)-C(I-1,J)) /dxaux(I,J)
ENDDO
ENDDO
FluxX(1,:) = 0.D0
FluxX(Nx+1,:) = 0.D0
C Then the outer cells - depending on boundary
C upstream X-boundary
IF (BcXup .EQ. 1) THEN
DO J = 1, Ny
FluxX(1,J) = Fxup(J)
ENDDO
ELSE IF (BcXup .EQ. 2) THEN
DO J = 1, Ny
FluxX(1,J) = -VF_x(1,J)*D_x(1,J)*(C(1,J)-Cxup(J))/dxaux(1,J)
ENDDO
ELSE IF (BcXup .EQ. 3) THEN
FluxX(1,:) = 0.D0
ELSE IF (BcXup .EQ. 4) THEN
DO J = 1, Ny
FluxX(1,J) = -axup(J)*(C(1,J)-Cxup(J))
ENDDO
ENDIF
C downstream X-boundary
IF (BcXdown .EQ. 1) THEN
DO J = 1, Ny
FluxX(Nx+1,J) = Fxdown(J)
ENDDO
ELSE IF (BcXdown .EQ. 2) THEN
DO J = 1, Ny
FluxX(Nx+1,J) = -VF_x(Nx+1,J)*D_x(Nx+1,J)*(Cxdown(J)-C(Nx,J)) &
& /dxaux(Nx+1,J)
ENDDO
ELSE IF (BcXdown .EQ. 3) THEN
FluxX(Nx+1,:) = 0.D0
ELSE IF (BcXdown .EQ. 4) THEN
DO J = 1, Ny
FluxX(Nx+1,J) = -axdown(J)*(Cxdown(J)-C(Nx,J))
ENDDO
ENDIF
C Rate of change in X-direction
DO I = 1,Nx
DO J = 1,Ny
AVx=0.25*(A_x(I,J)+A_x(I+1,J))*(VF_x(I,J)+VF_x(I+1,J))*dx(I,J)
dC(I,J) = -(A_x(I+1,J)*FluxX(I+1,J)-A_x(I,J)*FluxX(I,J))/AVx
ENDDO
ENDDO
C ## transport in Y-direction ##
C First diffusion internal cells
DO I = 1, Nx
DO J = 2, Ny
FluxY(I,J) = -VF_y(I,J)*D_y(I,J)*(C(I,J)-C(I,J-1)) /dyaux(I,J)
ENDDO
ENDDO
FluxY(:,1) = 0.D0
FluxY(:,Ny+1) = 0.D0
C upstream Y-boundary
IF (BcYup.EQ. 1) THEN
DO I = 1, Nx
FluxY(I,1) = Fyup(I)
ENDDO
ELSE IF (BcYup .EQ. 2) THEN
DO I = 1, Nx
FluxY(I,1) = -VF_y(I,1)*D_y(I,1)*(C(I,1)-Cyup(I))/dyaux(I,1)
ENDDO
ELSE IF (BcYup .EQ. 3) THEN
FluxY(:,1) = 0.D0
ELSE IF (BcYup .EQ. 4) THEN
DO I = 1, Nx
FluxY(I,1) = -ayup(I)*(C(I,1)-Cyup(I))
ENDDO
ENDIF
C downstream Y-boundary
IF (BcYdown.EQ. 1) THEN
DO I = 1, Nx
FluxY(I,Ny+1) = Fydown(I)
ENDDO
ELSE IF (BcYdown .EQ. 2) THEN
DO I = 1, Nx
FluxY(I,Ny+1) = -VF_y(I,Ny+1)*D_y(I,Ny+1)*(Cydown(I)-C(I,Ny)) &
& /dyaux(I,Ny+1)
ENDDO
ELSE IF (BcYdown .EQ. 3) THEN
FluxY(:,Ny+1) = 0.D0
ELSE IF (BcYdown .EQ. 4) THEN
DO I = 1, Nx
FluxY(I,Ny+1) = -aydown(I)*(Cydown(I)-C(I,Ny))
ENDDO
ENDIF
C rate of change, negative flux gradient
DO I = 1,Nx
DO J = 1,Ny
AVy=0.25*(A_y(I,J)+A_y(I,J+1))*(VF_y(I,J)+VF_y(I,J+1))*dy(I,J)
dC(I,J) = dC(I,J) &
& -(A_y(I,J+1)*FluxY(I,J+1)-A_y(I,J)*FluxY(I,J))/AVy
ENDDO
ENDDO
RETURN
END SUBROUTINE diff2D
C#==============================================================================
C Diffusion in a 2-dimensional finite difference grid - inputs are scalars
C#==============================================================================
SUBROUTINE diff2Dscal (Nx, Ny, C, Cxup, Cxdown, Cyup, Cydown, &
& Fxup, Fxdown, Fyup, Fydown, &
& axup, axdown, ayup, aydown, &
& BcXup, BcXdown, BcYup, BcYdown, &
& D, VF, dx, dy, FluxX, FluxY, dC)
IMPLICIT NONE
INTEGER Nx, Ny ! dimension of C
C input
DOUBLE PRECISION C(Nx,Ny),Cxup, Cxdown, Cyup, Cydown
DOUBLE PRECISION Fxup, Fxdown, Fyup, Fydown
DOUBLE PRECISION axup, axdown, ayup, aydown
DOUBLE PRECISION D, VF, dx, dy
C boundaries: 1= flux, 2 = conc, 3 = 0 grad, 4=convection
INTEGER BcXup, BcXdown, BcYup, BcYdown
C output
DOUBLE PRECISION FluxX(Nx+1, Ny), FluxY(Nx, Ny+1), dC(Nx,Ny)
C locals
INTEGER I, J
C -------------------------------------------------------------------------------
C ## transport in X-direction ##
C First diffusion internal cells
DO I = 2, Nx
DO J = 1, Ny
FluxX(I,J) = -VF*D*( (C(I,J)-C(I-1,J)) /dx)
ENDDO
ENDDO
FluxX(1,:) = 0.D0
FluxX(Nx+1,:) = 0.D0
C Then the outer cells - depending on boundary
C upstream X-boundary
IF (BcXup .EQ. 1) THEN
DO J = 1, Ny
FluxX(1,J) = Fxup
ENDDO
ELSE IF (BcXup .EQ. 2) THEN
DO J = 1, Ny
FluxX(1,J) = -VF*D*(C(1,J)-Cxup)/dx
ENDDO
ELSE IF (BcXup .EQ. 3) THEN
FluxX(1,:) = 0.D0
ELSE IF (BcXup .EQ. 4) THEN
DO J = 1, Ny
FluxX(1,J) = -axup*(C(1,J)-Cxup)
ENDDO
ENDIF
C downstream X-boundary
IF (BcXdown .EQ. 1) THEN
DO J = 1, Ny
FluxX(Nx+1,J) = Fxdown
ENDDO
ELSE IF (BcXdown .EQ. 2) THEN
DO J = 1, Ny
FluxX(Nx+1,J) = -VF*D*(Cxdown-C(Nx,J)) /dx
ENDDO
ELSE IF (BcXdown .EQ. 3) THEN
FluxX(Nx+1,:) = 0.D0
ELSE IF (BcXdown .EQ. 4) THEN
DO J = 1, Ny
FluxX(Nx+1,J) = -axdown*(Cxdown-C(Nx,J))
ENDDO
ENDIF
C Rate of change = negative flux gradient
DO I = 1,Nx
DO J = 1,Ny
dC(I,J) = -(FluxX(I+1,J)-FluxX(I,J))/VF/dx
ENDDO
ENDDO
C ## transport in Y-direction ##
C first internal cells
DO I = 1, Nx
DO J = 2, Ny
FluxY(I,J) = -VF*D*( (C(I,J)-C(I,J-1)) /dy)
ENDDO
ENDDO
FluxY(:,1) = 0.D0
FluxY(:,Ny+1) = 0.D0
C upstream Y-boundary
IF (BcYup.EQ. 1) THEN
DO I = 1, Nx
FluxY(I,1) = Fyup
ENDDO
ELSE IF (BcYup .EQ. 2) THEN
DO I = 1, Nx
FluxY(I,1) = -VF*D*(C(I,1)-Cyup)/dy
ENDDO
ELSE IF (BcYup .EQ. 3) THEN
FluxY(:,1) = 0.D0
ELSE IF (BcYup .EQ. 4) THEN
DO I = 1, Nx
FluxY(I,1) = -ayup*(C(I,1)-Cyup)
ENDDO
ENDIF
C downstream Y-boundary
IF (BcYdown.EQ. 1) THEN
DO I = 1, Nx
FluxY(I,Ny+1) = Fydown
ENDDO
ELSE IF (BcYdown .EQ. 2) THEN
DO I = 1, Nx
FluxY(I,Ny+1) = -VF*D*(Cydown-C(I,Ny))/dy
ENDDO
ELSE IF (BcYdown .EQ. 3) THEN
FluxY(:,Ny+1) = 0.D0
ELSE IF (BcYdown .EQ. 4) THEN
DO I = 1, Nx
FluxY(I,Ny+1) = -aydown*(Cydown-C(I,Ny))
ENDDO
ENDIF
C Rate of change = negative flux gradient
DO I = 1,Nx
DO J = 1,Ny
dC(I,J) = dC(I,J) -(FluxY(I,J+1)-FluxY(I,J))/VF/dy
ENDDO
ENDDO
RETURN
END SUBROUTINE diff2Dscal
C#==============================================================================
C Diffusion in a 2-dimensional polar coordinate system - all inputs are vectors
C C is matrix
C#==============================================================================
SUBROUTINE diffpolar (Nr, Ntet, C, Crup, Crdown, Ctetup, Ctetdown, &
& Frup, Frdown, Ftetup, Ftetdown, &
& Bcrup, Bcrdown, Bctetup, Bctetdown, &
& D_r, D_tet, r, tet, Fluxr, Fluxtet, dC)
IMPLICIT NONE
C dimension of
INTEGER Nr, Ntet
C concentration
DOUBLE PRECISION C(Nr,Ntet)
C used when Bc.. = 2, should be a vector then, else one number
DOUBLE PRECISION Crup(*), Crdown(*), Ctetup(*), Ctetdown(*)
C used when Bc.. = 1, should be a vectorthen , else one number
DOUBLE PRECISION Frup(*), Frdown(*), Ftetup(*), Ftetdown(*)
DOUBLE PRECISION D_r(Nr+1), D_tet(Ntet+1), &
& r(Nr+1), tet(Ntet+1), rc(Nr), tetc(Ntet)
DOUBLE PRECISION dr(Nr),draux(Nr+1),dtet(Ntet),dtetaux(Ntet+1)
c boundaries: 1= flux, 2 = conc, 3 = zero- grad, 5 = cyclic boundary
INTEGER BcRup, BcRdown, BcTetup, BcTetdown
C output
DOUBLE PRECISION Fluxr(Nr+1,Ntet), Fluxtet(Nr,Ntet+1), dC(Nr,Ntet)
C locals
INTEGER I, J
C -------------------------------------------------------------------------------
C The grid
DO I = 1, Nr
rc(I) = 0.5*(r(I)+r(I+1))
dr(I) = r(I+1) - r(I)
ENDDO
DO I = 1, Nr-1
draux(I+1) = (rc(I+1) - rc(I))
ENDDO
draux(1) = rc(1)-r(1)
draux(Nr+1) = r(Nr+1) - rc(Nr)
DO I = 1, Ntet
tetc(I) = 0.5*(tet(I)+tet(I+1))
dtet(I) = tet(I+1) - tet(I)
ENDDO
DO I = 1, Ntet-1
dtetaux(I+1) = (tetc(I+1) - tetc(I))
ENDDO
dtetaux(1) = tetc(1)-tet(1)
dtetaux(Ntet+1) = tet(Ntet+1) - tetc(Ntet)
C ## transport in r-direction ##
C First diffusion internal cells
DO I = 2, Nr
DO J = 1, Ntet
Fluxr(I,J) = -D_r(I)*(C(I,J)-C(I-1,J)) /draux(I)
ENDDO
ENDDO
Fluxr(1,:) = 0.D0
Fluxr(Nr+1,:) = 0.D0
C Then the outer cells - depending on boundary type
C upstream r-boundary
IF (Bcrup .EQ. 1) THEN
DO J = 1, Ntet
Fluxr(1,J) = Frup(J)
ENDDO
ELSE IF (Bcrup .EQ. 2) THEN
DO J = 1, Ntet
Fluxr(1,J) = -D_r(1) * (C(1,J)-Crup(J)) /draux(1)
ENDDO
ELSE IF (Bcrup .EQ. 3) THEN
Fluxr(1,:) = 0.D0
ELSE IF (Bcrup .EQ. 5) THEN
DO J = 1, Ntet
Crup(J) = (C(1,J)*draux(Nr+1) + C(Nr,J)*draux(1)) / &
& (draux(1)+draux(Nr+1))
Crdown(J) = Crup(J)
Fluxr(1,J) = -D_r(1)*(C(1,J)-Crup(J))/draux(1)
ENDDO
ENDIF
C downstream r-boundary
IF (Bcrdown .EQ. 1) THEN
DO J = 1, Ntet
Fluxr(Nr+1,J) = Frdown(J)
ENDDO
ELSE IF (Bcrdown .EQ. 2 .OR. Bcrdown .EQ. 5) THEN ! if 5: Crdown already estimated
DO J = 1, Ntet
Fluxr(Nr+1,J) = -D_r(Nr+1) * (Crdown(J)-C(Nr,J)) /draux(Nr+1)
ENDDO
ELSE IF (Bcrdown .EQ. 3) THEN
Fluxr(Nr+1,:) = 0.D0
ENDIF
C Rate of change in r-direction
DO I = 1,Nr
DO J = 1,Ntet
dC(I,J) = -(Fluxr(I+1,J)*r(I+1)-Fluxr(I,J)*r(I))/dr(I)/rc(I)
ENDDO
ENDDO
C ## transport in teta -direction ##
C First diffusion internal cells
DO I = 1, Nr
DO J = 2, Ntet
Fluxtet(I,J) = -D_tet(J)*(C(I,J)-C(I,J-1)) /dtetaux(J)/rc(I)
ENDDO
ENDDO
Fluxtet(:,1) = 0.D0
Fluxtet(:,Ntet+1) = 0.D0
C upstream teta-boundary
IF (BcTetup.EQ. 1) THEN
DO I = 1, Nr
Fluxtet(I,1) = Ftetup(I)
ENDDO
ELSE IF (BcTetup .EQ. 2) THEN
DO I = 1, Nr
Fluxtet(I,1) = -D_tet(1)*(C(I,1)-Ctetup(I))/dtetaux(1)/rc(I)
ENDDO
ELSE IF (BcTetup .EQ. 3) THEN
Fluxtet(:,1) = 0.D0
ELSE IF (BcTetup .EQ. 5) THEN
DO I = 1, Nr
Ctetup(I) = (C(I,1)*dtetaux(Ntet+1) + C(I,Ntet)*dtetaux(1)) / &
& (dtetaux(1)+dtetaux(Ntet+1))
Ctetdown(I) = Ctetup(I)
Fluxtet(I,1) = -D_tet(1)*(C(I,1)-Ctetup(I))/dtetaux(1)/rc(I)
ENDDO
ENDIF
C downstream teta-boundary
IF (BcTetdown.EQ. 1) THEN
DO I = 1, Nr
Fluxtet(I,Ntet+1) = Ftetdown(I)
ENDDO
ELSE IF (BcTetdown .EQ. 2 .OR. BcTetdown .EQ. 5) THEN
DO I = 1, Nr
Fluxtet(I,Ntet+1) = -D_tet(Ntet+1)*(Ctetdown(I)-C(I,Ntet)) &
& /dtetaux(Ntet+1)/rc(I)
ENDDO
ELSE IF (BcTetdown .EQ. 3) THEN
Fluxtet(:,Ntet+1) = 0.D0
ENDIF
C rate of change, negative flux gradient
DO I = 1,Nr
DO J = 1,Ntet
dC(I,J) = dC(I,J) &
& -(Fluxtet(I,J+1)-Fluxtet(I,J))/dtet(J)/rc(I)
ENDDO
ENDDO
RETURN
END
