module global_data implicit none save #include "include/finclude/petsc.h" #include "include/finclude/petscvec.h" #include "include/finclude/petscmat.h" #include "include/finclude/petscksp.h" #include "include/finclude/petscpc.h" #include "include/finclude/petscsys.h" integer :: size_x,size_y,Total_time_step,new_start,interval,gridgen,safe_int integer :: steady,quasi_steady,Total_k,time,mom_solver,poisson_solver,start_time real(8) :: CFL, Re, scheme,Pi, time_sta,act_time,inv_Re real(8), allocatable :: p(:,:), u(:,:), v(:,:), u_old(:,:), v_old(:,:),phi(:,:) real(8), allocatable :: x(:),y(:), xu(:),yv(:),big_A(:,:),q_p(:) real(8), allocatable :: u_ast(:,:), v_ast(:,:) integer, allocatable :: int_impl(:,:),int_mat(:,:),int_A(:,:) real(8) :: Uc,del_t,del_t_old real(8) :: impl_mat_A,q_impl real(8), allocatable :: diff_pres_x(:,:),diff_pres_y(:,:) real(8), allocatable :: impl_u_x_e(:,:,:),impl_u_x_w(:,:,:),impl_u_x_n(:,:,:),impl_u_x_s(:,:,:),impl_u_x_c(:,:,:) real(8), allocatable :: impl_u_y_e(:,:,:),impl_u_y_s(:,:,:),impl_u_y_se(:,:,:),impl_u_y_c(:,:,:) real(8), allocatable :: impl_v_x_e(:,:,:),impl_v_x_w(:,:,:),impl_v_x_n(:,:,:),impl_v_x_s(:,:,:),impl_v_x_c(:,:,:) real(8), allocatable :: impl_v_y_n(:,:,:),impl_v_y_w(:,:,:),impl_v_y_nw(:,:,:),impl_v_y_c(:,:,:) real(8), allocatable :: fc1_prev_r(:),fc1_prev_c(:),fc2_prev_r(:),fc2_prev_c(:) real(8) :: fc1_1,fc1_2,fc1_3,fc1_4,fd1_1,fd1_2,fd1_3,fd1_4,fc2_1,fc2_2,fc2_3,fc2_4,fd2_1,fd2_2,fd2_3,fd2_4 real(8), allocatable :: fc1_prev_r_old(:),fc1_prev_c_old(:),fc2_prev_r_old(:),fc2_prev_c_old(:) real(8), allocatable :: fd1_prev_r(:),fd1_prev_c(:),fd2_prev_r(:),fd2_prev_c(:),p_prev_r(:),p_prev_c(:) type cell real(8) :: cx,cy,area real(8) :: pd_E,pd_W,pd_N,pd_S,Ed_E,Ed_W,Ed_N,Ed_S end type cell type(cell), allocatable :: cp(:,:),cu(:,:),cv(:,:) !c=cells Vec xx,b_rhs,xx_uv,b_rhs_uv Mat A_mat,A_mat_uv ! /* sparse matrix */ KSP ksp,ksp_uv !/* linear solver context */ PC pc,pc_uv PCType ptype KSPType ksptype PetscTruth flgg KSPConvergedReason reason PetscOffset i_vec !parallel integer :: myid,num_procs,jsta,jend,ksta_p,kend_p,jsta_ext,jend_ext !,jsta_ext0 ! integer :: jsta_ext_u,jend_ext_u,jsta_ext_v,jend_ext_v integer :: ksta_m,kend_m,jsta_u,jsta_v,jend_u,jend_v,jsta_m,jend_m contains subroutine allo_var !allocate memory for variables #include "include/finclude/petsc.h" #include "include/finclude/petscvec.h" #include "include/finclude/petscmat.h" #include "include/finclude/petscksp.h" #include "include/finclude/petscpc.h" #include "include/finclude/petscsys.h" integer :: status(2),ierr,k,NN total_k=(size_x+2)*(size_y+2); NN=2*total_k jsta=myid*(size_y/num_procs)+1; jend=(myid+1)*(size_y/num_procs) if (num_procs==1) then jsta_ext=0; jend_ext=jend+1 !; jsta_ext0=0; jsta2=jsta+1; jend2=jend-1; jend3=jend-2 !jsta_ext=0; jend_ext=jend+1; jsta_ext0=0 jsta_u=jsta; jsta_v=jsta; jend_u=jend; jend_v=jend-1 !jsta_ext_u=0; jend_ext_u=size_y+1; jsta_ext_v=0; jend_ext_v=size_y jsta_m=0; jend_m=size_y+1 else if (myid==0) then jsta_ext=0; jend_ext=jend+2 !; jsta_ext0=0 jsta_u=jsta; jsta_v=jsta; jend_u=jend; jend_v=jend ! jsta_ext_u=0; jend_ext_u=jend; jsta_ext_v=0; jend_ext_v=jend jsta_m=0; jend_m=jend else if (myid==num_procs-1) then jsta_ext=jsta-2; jend_ext=jend+1 !; jsta_ext0=jsta_ext jsta_u=jsta; jsta_v=jsta; jend_u=jend; jend_v=jend-1 ! jsta_ext_u=jsta; jend_ext_u=size_y+1; jsta_ext_v=jsta; jend_ext_v=size_y jsta_m=jsta; jend_m=size_y+1 else jsta_ext=jsta-2; jend_ext=jend+2 !; jsta_ext0=jsta_ext jsta_u=jsta; jsta_v=jsta; jend_u=jend; jend_v=jend ! jsta_ext_u=jsta; jend_ext_u=jend; jsta_ext_v=jsta; jend_ext_v=jend jsta_m=jsta; jend_m=jend end if end if allocate (x(-1:size_x+1), STAT = status(1)) allocate (y(-1:size_y+1), STAT = status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (xu(-1:size_x+1), STAT = status(1)) ! allocate (yu(0:size_y), STAT = status(2)) cancel, no need if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" ! allocate (xv(0:size_x), STAT = status(1)) allocate (yv(-1:size_y+1), STAT = status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (p(0:size_x+1,jsta_ext:jend_ext), STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (u(0:size_x+1,jsta_ext:jend_ext), STAT=status(1));allocate (v(0:size_x+1,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (u_old(0:size_x+1,jsta_ext:jend_ext), STAT=status(1));allocate (v_old(0:size_x+1,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (diff_pres_x(0:size_x+1,jsta_ext:jend_ext), STAT=status(1)) allocate (diff_pres_y(0:size_x+1,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (u_ast(0:size_x+1,jsta_ext:jend_ext), STAT=status(1));allocate (v_ast(0:size_x+1,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (phi(0:size_x+1,jsta_ext:jend_ext), STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (big_A((size_x+2)*(size_y+2),5), STAT=status(1)) allocate (q_p((size_x+2)*(size_y+2)), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (int_A((size_x+2)*(size_y+2),5), STAT=status(1)) if (status(1)/=0) STOP "Cannot allocate memory" if (scheme>=1) then allocate (int_impl(2*(size_x+2)*(size_y+2),11), STAT=status(2)) if (status(2)/=0) STOP "Cannot allocate memory" allocate (impl_u_x_e(2,size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (impl_u_x_w(2,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (impl_u_x_n(2,size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (impl_u_x_s(2,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (impl_u_x_c(8,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(2)/=0) STOP "Cannot allocate memory" allocate (impl_u_y_e(1,size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (impl_u_y_s(1,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (impl_u_y_se(1,size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (impl_u_y_c(1,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (impl_v_x_e(2,size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (impl_v_x_w(2,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (impl_v_x_n(2,size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (impl_v_x_s(2,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (impl_v_x_c(8,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(2)/=0) STOP "Cannot allocate memory" allocate (impl_v_y_n(1,size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (impl_v_y_w(1,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (impl_v_y_nw(1,size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (impl_v_y_c(1,size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" end if !stop allocate (fc1_prev_r(0:size_x), STAT=status(1));allocate (fc1_prev_c(jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (fc2_prev_r(0:size_x+1), STAT=status(1));allocate (fc2_prev_c(jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (fc1_prev_r_old(0:size_x), STAT=status(1));allocate (fc1_prev_c_old(jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (fc2_prev_r_old(0:size_x+1), STAT=status(1));allocate (fc2_prev_c_old(jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (fd1_prev_r(0:size_x), STAT=status(1));allocate (fd1_prev_c(jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (fd2_prev_r(0:size_x+1), STAT=status(1));allocate (fd2_prev_c(jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (p_prev_r(0:size_x+1), STAT=status(1));allocate (p_prev_c(jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (cu(0:size_x+1,jsta_ext:jend_ext), STAT=status(1));allocate (cv(0:size_x+1,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (cp(0:size_x+1,jsta_ext:jend_ext), STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" call MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,NN,NN,11,PETSC_NULL_INTEGER,11,PETSC_NULL_INTEGER,A_mat_uv,ierr) call MatSetFromOptions(A_mat_uv,ierr) call MatGetOwnershipRange(A_mat_uv,ksta_m,kend_m,ierr) call KSPCreate(MPI_COMM_WORLD,ksp_uv,ierr) call VecCreateMPI(MPI_COMM_WORLD,PETSC_DECIDE,2*(size_x+2)*(size_y+2),b_rhs_uv,ierr) call VecDuplicate(b_rhs_uv,xx_uv,ierr) call MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,total_k,total_k,5,PETSC_NULL_INTEGER,5,PETSC_NULL_INTEGER,A_mat,ierr) call MatSetFromOptions(A_mat,ierr) call MatGetOwnershipRange(A_mat,ksta_p,kend_p,ierr) call KSPCreate(MPI_COMM_WORLD,ksp,ierr) call VecCreateMPI(MPI_COMM_WORLD,PETSC_DECIDE,total_k,b_rhs,ierr) call VecDuplicate(b_rhs,xx,ierr) end subroutine allo_var subroutine ini_var1 !initialize/allocate memory for variables integer :: status(2) u=0.;v=0.;p=0.;u_old=0.;v_old=0.;diff_pres_x=0.;diff_pres_y=0.;phi=0. big_A=0.;q_p=0.;int_A=0 impl_mat_A=0.;q_impl=0.;int_impl=0 impl_u_x_e=0.;impl_u_x_w=0.;impl_u_x_n=0.;impl_u_x_s=0.;impl_u_x_c=0. impl_u_y_e=0.;impl_u_y_s=0.;impl_u_y_se=0.;impl_u_y_c=0. impl_v_x_e=0.;impl_v_x_w=0.;impl_v_x_n=0.;impl_v_x_s=0.;impl_v_x_c=0. impl_v_y_n=0.;impl_v_y_w=0.;impl_v_y_nw=0.;impl_v_y_c=0. x=0.;y=0.;xu=0.;yv=0.;u_ast=0.;v_ast=0. end subroutine ini_var1 subroutine de_ini_var !de-initialize/deallocate memory for variables integer :: status(2),ierr deallocate (p, STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (u, STAT=status(1));deallocate (v, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (u_old, STAT=status(1));deallocate (v_old, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (u_ast, STAT=status(1));deallocate (v_ast, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (x, STAT=status(1));deallocate (y, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (q_p, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (int_impl, STAT=status(2)) if (status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_u_x_e, STAT=status(1));deallocate (impl_u_x_w, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_u_x_n, STAT=status(1));deallocate (impl_u_x_s, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_u_x_c, STAT=status(2)) if (status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_u_y_e, STAT=status(1));deallocate (impl_u_y_s, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_u_y_se, STAT=status(1));deallocate (impl_u_y_c, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_v_x_e, STAT=status(1));deallocate (impl_v_x_w, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_v_x_n, STAT=status(1));deallocate (impl_v_x_s, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_v_x_c, STAT=status(2)) if (status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_v_y_n, STAT=status(1));deallocate (impl_v_y_w, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (impl_v_y_nw, STAT=status(1));deallocate (impl_v_y_c, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" call MatDestroy(A_mat_uv,ierr) call KSPDestroy(ksp_uv,ierr) call VecDestroy(xx_uv,ierr) call VecDestroy(b_rhs_uv,ierr) deallocate (diff_pres_x, STAT=status(1));deallocate (diff_pres_y, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (fc1_prev_r, STAT=status(1));deallocate (fc1_prev_c, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (fc2_prev_r, STAT=status(1));deallocate (fc2_prev_c, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (fc1_prev_r_old, STAT=status(1));deallocate (fc1_prev_c_old, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (fc2_prev_r_old, STAT=status(1));deallocate (fc2_prev_c_old, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (fd1_prev_r, STAT=status(1));deallocate (fd1_prev_c, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (fd2_prev_r, STAT=status(1));deallocate (fd2_prev_c, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (p_prev_r, STAT=status(1));deallocate (p_prev_c, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (cu, STAT=status(1));deallocate (cv, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (cp, STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" call VecDestroy(xx,ierr) call VecDestroy(b_rhs,ierr) call MatDestroy(A_mat,ierr) call KSPDestroy(ksp,ierr) call PetscFinalize(ierr) end subroutine de_ini_var C C .................................................................. C C SUBROUTINE PNPOLY C C PURPOSE C TO DETERMINE WHETHER A POINT IS INSIDE A POLYGON C C USAGE C CALL PNPOLY (PX, PY, XX, YY, N, INOUT ) C C DESCRIPTION OF THE PARAMETERS C PX - X-COORDINATE OF POINT IN QUESTION. C PY - Y-COORDINATE OF POINT IN QUESTION. C XX - N LONG VECTOR CONTAINING X-COORDINATES OF C VERTICES OF POLYGON. C YY - N LONG VECTOR CONTAING Y-COORDINATES OF C VERTICES OF POLYGON. C N - NUMBER OF VERTICES IN THE POLYGON. C INOUT - THE SIGNAL RETURNED: C -1 IF THE POINT IS OUTSIDE OF THE POLYGON, C 0 IF THE POINT IS ON AN EDGE OR AT A VERTEX, C 1 IF THE POINT IS INSIDE OF THE POLYGON. C C REMARKS C THE VERTICES MAY BE LISTED CLOCKWISE OR ANTICLOCKWISE. C THE FIRST MAY OPTIONALLY BE REPEATED, IF SO N MAY C OPTIONALLY BE INCREASED BY 1. C THE INPUT POLYGON MAY BE A COMPOUND POLYGON CONSISTING C OF SEVERAL SEPARATE SUBPOLYGONS. IF SO, THE FIRST VERTEX C OF EACH SUBPOLYGON MUST BE REPEATED, AND WHEN CALCULATING C N, THESE FIRST VERTICES MUST BE COUNTED TWICE. C INOUT IS THE ONLY PARAMETER WHOSE VALUE IS CHANGED. C THE SIZE OF THE ARRAYS MUST BE INCREASED IF N > MAXDIM C WRITTEN BY RANDOLPH FRANKLIN, UNIVERSITY OF OTTAWA, 7/70. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C A VERTICAL LINE IS DRAWN THRU THE POINT IN QUESTION. IF IT C CROSSES THE POLYGON AN ODD NUMBER OF TIMES, THEN THE C POINT IS INSIDE OF THE POLYGON. C C .................................................................. C SUBROUTINE PNPOLY(PX,PY,XX,YY,N,INOUT) INTEGER I,J,N,INOUT REAL X(N),Y(N),XX(N),YY(N),PX,PY LOGICAL MX,MY,NX,NY INTEGER O C OUTPUT UNIT FOR PRINTED MESSAGES DATA O/6/ GO TO 6 WRITE(O,7) 7 FORMAT('0WARNING:',I5,' TOO GREAT FOR THIS VERSION OF PNPOLY. 1RESULTS INVALID') RETURN 6 DO 1 I=1,N X(I)=XX(I)-PX 1 Y(I)=YY(I)-PY INOUT=-1 DO 2 I=1,N J=1+MOD(I,N) MX=X(I).GE.0.0 NX=X(J).GE.0.0 MY=Y(I).GE.0.0 NY=Y(J).GE.0.0 IF(.NOT.((MY.OR.NY).AND.(MX.OR.NX)).OR.(MX.AND.NX)) GO TO 2 IF(.NOT.(MY.AND.NY.AND.(MX.OR.NX).AND..NOT.(MX.AND.NX))) GO TO 3 INOUT=-INOUT GO TO 2 3 IF((Y(I)*X(J)-X(I)*Y(J))/(X(J)-X(I))) 2,4,5 4 INOUT=0 RETURN 5 INOUT=-INOUT 2 CONTINUE RETURN END SUBROUTINE PNPOLY end module global_data