- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi, I got this error when I write my own code. It's look strange because the syntax I think is correct. The error exist in this part
auE_h=(-uCe,(uDx-uCe)/2,0)
Any help will be great thank you and here I posted the whole code
PROGRAM BACKWARD
IMPLICIT NONE ! Declaration parameter ! INTEGER TYPE INTEGER, PARAMETER :: NX=300 INTEGER, PARAMETER :: NY=30 INTEGER dx, dy, i,j, maxIter, maxNMiter, n ! REAL TYPE REAL, PARAMETER :: L=15 REAL, PARAMETER :: H=1 REAL(8) alpha_u, alpha_p, Cs, err, Ho, Hs, nu, Re, uMax REAL(8) uCe, uCw, uCn, uCs !CONVECTIVE TERM PARAMETER FOR u REAL(8) vCe, vCw, vCn, vCs !CONVECTIVE TERM PARAMETER FOR v REAL(8) dudy, dudx, nutu !SGS CLASSIC SMAGORINSKY PARAMETER for u REAL(8) dvdy, dvdx, nutv !SGS CLASSIC SMAGORINSKY PARAMETER for v REAL(8) uDx, uDy !DIFFUISIVE TERM PARAMETER for u REAL(8) vDx, vDy !DIFFUISIVE TERM PARAMETER for v ! REAL TYPE WITH DIMENSION REAL, DIMENSION(NX+1,NY+2) :: u, uold, ustar, uprime REAL, DIMENSION(NX+2,NY+1) :: v, vold, vstar, vprime REAL, DIMENSION(NX+2,NY+2) :: p, pstar, pprime REAL, DIMENSION(1,NY+2) :: uLeft REAL, DIMENSION(1,NY+2) :: vLeft REAL(8) auE_h, auW_h, auN_h, auS_h REAL, DIMENSION(NX,NY) ::auE REAL, DIMENSION(NX,NY) ::auW REAL, DIMENSION(NX,NY) ::auN REAL, DIMENSION(NX,NY) ::auS REAL, DIMENSION(NX,NY) ::auP REAL, DIMENSION(NX,NY) ::dU REAL, DIMENSION(NX,NY) ::dV REAL(8) avE_h, avW_h, avN_h, avS_h REAL, DIMENSION(NX,NY) ::avE REAL, DIMENSION(NX,NY) ::avW REAL, DIMENSION(NX,NY) ::avN REAL, DIMENSION(NX,NY) ::avS REAL, DIMENSION(NX,NY) ::avP REAL(8) apE_h, apW_h, apN_h, apS_h REAL, DIMENSION(NX,NY) ::apE REAL, DIMENSION(NX,NY) ::apW REAL, DIMENSION(NX,NY) ::apN REAL, DIMENSION(NX,NY) ::apS REAL, DIMENSION(NX,NY) ::apP REAL, DIMENSION(NX,NY) ::bpP REAL, DIMENSION(NX,NY) ::ux ! Convergence coeff REAL, DIMENSION(NX,NY) ::uxold REAL, DIMENSION(NX,NY) ::vx REAL, DIMENSION(NX,NY) ::vxold REAL(8) cmax_1, cmax_2, convergence !CONVERGENCE REAL, DIMENSION(1,NY) ::Y ! FILL THE VALUE maxIter=50000 !INTEGER VALUE maxNMiter=1 alpha_u=0.3 !REAL VALUE alpha_p=0.7 Cs=0.07 err=1e-5 Ho=1.0 Re=100 nu=2e-3 dx=L/NX dy=H/NY !Step size Hs=H-Ho uMax=Re*nu/Hs ! ALLOCATE THE GUESS VALUE AND INITIAL CONDITION u(:,:)=0 uold(:,:)=0 !INITIAL CONDITION ustar(:,:)=0 uprime(:,:)=0 dU(:,:)=0 v(:,:)=0 vold(:,:)=0 !INITIAL CONDITION vstar(:,:)=0 vprime(:,:)=0 dV(:,:)=0 p(:,:)=0 pstar(:,:)=0 pprime(:,:)=0 uLeft(:,:)=0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!! SIMPLE SCHEME APPLY FROM HERE & !!!! CLASSIC SMAGORINSKY SGS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO n=1,maxIter ! STEP 1a, solve x-momentum as ustar DO i=2,NX DO j=2,NY+1 ! Convective flux terms use CENTRAL DIFFERENCE uCe=dy*0.5*(uold(i,j)+uold(i+1,j)) uCw=dy*0.5*(uold(i,j)+uold(i-1,j)) uCn=dx*0.5*(vold(i,j)+vold(i,j+1)) uCs=dx*0.5*(vold(i,j-1)+vold(i+1,j-1)) ! Calculation for rate of strain tensor dudy=(0.5*(uold(i,j)+uold(i,j+1)))/dy dudy=(0.5*(uold(i,j)+uold(i+1,j)))/dx nutu=(Cs*FLOAT(dx))**2*((0.5*(dudx+dudy)**2)**0.5)*(dudx+dudy) ! Diffusive flux term uDx=(FLOAT(dy)/FLOAT(dx)*((1/Re)+nutu)) uDy=(FLOAT(dx)/FLOAT(dy))*((1/Re)+nutu) ! ASSEMBLY COEFFICIENT FROM HYBRID SCHEMEauE_h=(-uCe,(uDx-uCe)/2,0) auW_h=(uCw, (uDx+0.5*uCw), 0) auN_h=(-uCn, (uDy-0.5*uCn), 0) auS_h=(uCs, (uDy+0.5*uCs), 0)auE(i,j)=MAX(auE_h) auW(i,j)=MAX(uCw, (uDx+0.5*uCw), 0) auN(i,j)=MAX(-uCn, (uDy-0.5*uCn), 0) auS(i,j)=MAX(uCs, (uDy+0.5*uCs), 0) !Central coeff & under relaxation auP(i,j)=(auE(i,j)+auW(i,j)+auN(i,j)+auS(i,j))/alpha_u !Velocity component for PCE dU(i,j)=FLOAT(dy)/auP(i,j) END END !Use uold for calculate ustar ustar=uold DO iter=maxNMiter DO i=2, NX DO j=2, NY+1 ustar(i,j)=(auE(i,j)*ustar(i+1,j)+auW(i,j)*ustar(i-1,j)+auN(i,j)*ustar(i,j+1)+auS(i,j)*ustar(i,j-1)-FLOAT(dy)*(pstar(i+1,j)-pstar(i,j))+(1-alpha_u)*auP(i,j)*uold(i,j))/auP(i,j) END END END ! ustar FOR X-MOMENTUM HAS BEEN SOLVE !!!!!!!!!! STEP 1B BEGIN ! Solve y-momentum as vstar DO i=2, NX+1 DO j=2, NY ! Convective flux terms use CENTRAL DIFFERENCE vCe=dy*0.5*(uold(i,j)+uold(i,j+1)) vCw=dy*0.5*(uold(i-1,j)+uold(i-1,j+1)) vCn=dx*0.5*(vold(i,j)+vold(i,j+1)) vCs=dx*0.5*(vold(i,j)+vold(i,j-1)) ! Calculation for rate of strain tensor dvdy=(0.5*(vold(i,j)+vold(i,j+1)))/dy dvdy=(0.5*(vold(i,j)+vold(i+1,j)))/dx nutv=(Cs*FLOAT(dx))**2*((0.5*(dvdx+dvdy)**2)**0.5)*(dvdx+dvdy) ! Diffusive flux term vDx=(FLOAT(dy)/FLOAT(dx))*((1/Re)+nutv) vDy=(FLOAT(dx)/FLOAT(dy))*((1/Re)+nutv) ! ASSEMBLY COEFFICIENT FROM HYBRID SCHEME avE_h=(-vCe, (vDx-0.5*vCe), 0) avW_h=(vCw, (vDx+0.5*vCw), 0) avN_h=(-vCn, (vDy-0.5*vCn), 0) avS_h=(vCs, (vDy+0.5*vCs), 0) avE(i,j)=MAX(avE_h) avW(i,j)=MAX(avW_h) avN(i,j)=MAX(avN_h) avS(i,j)=MAX(avS_h) !Central coeff & under relaxation avP(i,j)=(avE(i,j)+avW(i,j)+avN(i,j)+avS(i,j))/alpha_u !Velocity component for PCE dV(i,j)=FLOAT(dx)/avP(i,j) END END ! Calculate vstar vstar=vold DO iter=maxNMiter DO i=2, NX+1 DO j=2, NY vstar(i,j)=(avE(i,j)*vstar(i+1,j)+avW(i,j)*vstar(i-1,j)+avN(i,j)*vstar(i,j+1)+avS(i,j)*vstar(i,j-1)-FLOAT(dy)*(pstar(i+1,j)-pstar(i,j))+(1-alpha_u)*avP(i,j)*vold(i,j))/avP(i,j) END END END ! vstar has been calculated, move to next step !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! STEP 2 CALCULATE PRESSURE CORRECTION DO i=2, NX+1 DO j=2, NY+1 ! Boundary coefficient apE(i,j) = dU(i,j)*dy apW(i,j) = dU(i-1,j)*dy apN(i,j) = dV(i,j)*dx apS(i,j) = dV(i,j-1)*dx ! Central coefficient apP(i,j) = apE(i,j)+apW(i,j)+apN(i,j)+apS(i,j) ! RHS value bpP(i,j)= (ustar(i-1,j)-ustar(i,j))*dy+(vstar(i,j-1)-vstar(i,j))*dx END END !! APPLY BOUNDARY CONDITION FOR PRESSURE IF (n==1 .AND. i=2 .AND. j=2) THEN apP(i,j)=1.0 apE(i,j)=0 apW(i,j)=0 apN(i,j)=0 apS(i,j)=0 bpP(i,j)=0 END !!Calculate pprime pprime(:,:)=0 DO iter=maxNMiter DO i=2, NX+1 DO j=2, NY+1 pprime(i,j)= pprime(i,j)=(apE(i,j)*pprime(i+1,j)+apW(i,j)*pprime(i-1,j)+apN(i,j)*pprime(i,j+1)+apS(i,j)*pprime(i,j-1)+bpP(i,j)/apP(i,j) END END END !!!!!!!!!! pprime has been calculated move to next section !!!!!!!!!! STEP 3 CALCULATE VELOCITY CORRECTION !! u-velocity correction DO i=2, NX DO j=2, NY+1 uprime(i,j)=dU(i,j)*(pprime(i,j)-pprime(i+1,j)) END END !! v-velocity correction DO i=2, NX+1 DO j=2, NY vprime(i,j)=dV(i,j)*(pprime(i,j)-pprime(i,j+1)) END END !!!!!! STEP 3 FINISH MOVE TO STEP 4 !!!!!! STEP 4 CALCULATE REAL u,v, and p !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! CALCULATE REAL u Do i=2, NX Do j=2, NY+1 u(i,j)=ustar(i,j)+uprime(i,j) END END !! CALCULATE REAL v Do i=2, NX+1 Do j=2, NY v(i,j)=vstar(i,j)+vprime(i,j) END END !! CALCULATE REAL p Do i=2, NX+1 Do j=2, NY+1 p(i,j)=pstar(i,j)+pprime(i,j)*alpha_p END END !!!!!!! STEP 4 FINISH !!!!!!!!!! !!!!!!! STEP 5 ITERATION UNTIL CONVERGENCE !!!!!!! ux=u(1:NX+1, 1:NY+1) uxold=uold(1:NX+1, 1:NY+1) vx=v(1:NX+1, 1:NY+1) vxold=vold(1:NX+1, 1:NY+1) cmax1=MAX(MAX(ux-uxold)) cmax2=MAX(MAX(vx-vxold)) convergence=MAX(cmax1, cmax2) IF convergence<err THEN STOP ELSE uold=u vold=v pstar=p !!!!!!!! CONVERGENCE CRITERIA IS SET !!!!!!!! STEP 6 PROBLEM SET UP !!!!!! APPLY INLET VELOCITY FUNCTION DO j=NY+1, 2, -1 Y(1,j)=(1, H:-dy:0) IF Y(1,j)<Ho THEN uLeft(1,j)=4*umax*(Y(i,j)/Ho)*(1-(Y(1,j)/Ho)) vLeft(1,j)=0 ELSE IF Y(1,j)>Ho THEN uLeft(1,j)=0 vLeft(1,j)=0 END END !!!!!! APPLY BOUNDARY CONDITION FOR VELOCITY uold(:,1)=0.0 !bottom bound vold(:,1)=0.0 uOld(:,Ny+2) = 0.0 vOld(:,Ny+1) = 0.0 uOld(1,:) = uLeft(1,:) vOld(:,1) = vLeft(:,1) uOld(Nx+1,:) = uOld(Nx,:) vOld(Nx+2,:) = vOld(Nx+1,:) END END END END END PROGRAM BACKWARD
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
avE_h, for example, is a REAL(8) two-dimensional array. You have:
avE_h=(-vCe, (vDx-0.5*vCe), 0)
The syntax of values in parentheses separated by commas looks like a complex constant, but there are three values, so that doesn't fit. There's no other Fortran syntax that matches this.
What do you want this code to do? Were you trying for an array constructor?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
In looking at the rest of the code, I think there is a missing MAX or MIN to produce a scalar. But his is a speculation on my part.
IOW: What do you want this code to do?
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The code is not yet Fortran. For example, it has a number of places where an apparent DO block is terminated by END instead of END DO. Similarly, IF blocks have no matching ENDIF or ELSE. The lack of indentation makes the code hard to follow. Note that, as far as the compiler is concerned, an END statement terminates the program (on line 134, and many more times again!). Since that is not followed by a SUBROUTINE, FUNCTION or MODULE statement, everything after line 134 is invalid.
Line 122 starts with <s>, which is not meaningful in Fortran. Line 125 ends with </s>. Did you, perhaps, copy and paste the Fortran source code from a Web page? That might also explain why the parenthesis around the conditional clauses of IF statements were lost.
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page