module global_data implicit none save #include "finclude/petsc.h90" integer :: size_x,size_y,Total_time_step,new_start,interval,gridgen,safe_int,OS,airfoil_no integer :: steady,quasi_steady,Total_k,time,mom_solver,poisson_solver,start_time,motion !size_x must be in multiples of 37/32/36/40/55/41, !size_y must be in multiples of 26/16/36 !gridgen1 - 32x20, gridgen4 - 30x22, gridgen5 - 30x24 real(8) :: CFL, Re, scheme, B, AA, BB,AB, ld,air_centy,Pi,hy0,k0,freq,phase_ang,theta0,loc_rot real(8) :: time_sta,act_time,vel_h,vel_hn,inv_Re real(8) :: hx0, ini_ang,air_centx,air_centx_old,air_centx_old2 !constant "fact" & "B" must give adjacent grid (e.g. bet. front outside n inside cylinder) !of similar spacing to avoid numerical errors !4=2.7, 5=2.5, 6=1.6 !AA=2.7, BB=2.8,AB=2.8 ! AA=2.5, BB=2.17, AB=2.25 ! AA=2.6, BB=2.6,AB=2.7 ! AA=3.8, BB=4.,AB=3.9 !AA=4.2, BB=5.7,AB=5. ! AA=5.2, BB=5.7,AB=5.4 !scheme=0 for explicit, scheme=1 for semi implicit, choice=2 for full implicit !semi implicit cannot start from zero !domain(x/y)=length in x/y !if new_start=1 start from beginning ie time=0 !if new_start=0 start from reading u,v,p values file !pres,vel(new/old) at cell centre ! real(8) :: big_A_1,big_A_2,big_A_3,big_A_4,big_A_5,q_p ! integer :: int_A_1,int_A_2,int_A_3,int_A_4,int_A_5 real(8), allocatable :: p(:,:), u(:,:), v(:,:), u_old(:,:), v_old(:,:),phi(:,:) real(8), allocatable :: x(:),y(:), xu(:),yu(:),xv(:),yv(:),big_A(:,:),q_p(:) real(8), allocatable :: u_ast(:,:), v_ast(:,:),u_ast2(:,:), v_ast2(:,:) integer, allocatable :: int_impl(:,:),int_semi_x(:,:),int_semi_y(:,:),int_mat(:,:),int_A(:,:) ! integer, allocatable :: in_body_u(:,:),in_body_v(:,:) !inside body=1, outside = 0, based on x,y coordinates' position, not cell center real(8) :: Uc,del_t,del_t_old,air_centy_old,air_centy_old2 !full implicit mtd matrix (u & v together) real(8) :: impl_mat_A,q_impl,semi_mat_A,q_semi_xx,q_semi_yy real(8), allocatable :: diff_pres_x(:,:),diff_pres_y(:,:) !,q_impl(:),impl_mat(:,:), 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 :: semi_diag_x(:),semi_diag_y(:) !diff_pres = pres diff in x & y, s_A = 1/ sum (1/A) !semi implicit mtd matrix (u & v separate) real(8), allocatable :: semi_mat_x(:,:),semi_mat_y(:,:) !,q_semi_x(:),q_semi_y(:) 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(:) real(8), allocatable :: area_cv_percent_u(:,:),area_cv_percent_v(:,:) !previous stored values of convective(c), diffusive(d), 1(east), 2(north) type cell !cell no in x and y, types=-1=internal cells, =0=fluid cells, =1=ghost cells, =2=cell center approx equal to boundary value !next= che(end_x-sta_x)*(end_y-sta_y) if cell has ghost cell any of 4 sides - 1 left, 2 right, 3 bottom, 4 top real(8) :: cx,cy,area !,area_cv_percent ,ne_x,nw_x,se_x,sw_x,ne_y,nw_y,se_y,sw_y !center in x and y, area of cell,area of cell covered by body,percentage of area of cell covered by body,corners coord real(8) :: pd_E,pd_W,pd_N,pd_S,Ed_E,Ed_W,Ed_N,Ed_S !perpendicular dist to E,W,N,S from center,length of E,W,N,S edges !note that length of W,S face will be same as E,N face of adjacent cell !real(8) :: d_PE,d_PW,d_PN,d_PS,uf_E,uf_W,uf_N,uf_S,vf_E,vf_W,vf_N,vf_S !real(8) :: uf_E_old,uf_W_old,uf_N_old,uf_S_old,vf_E_old,vf_W_old,vf_N_old,vf_S_old !length from P to E,W,N,S !real(8) :: vel_f_E,vel_f_W,vel_f_N,vel_f_S,vel_f_E_m,vel_f_W_m,vel_f_N_m,vel_f_S_m !real(8) :: vel_f_E_old,vel_f_W_old,vel_f_N_old,vel_f_S_old,vel_f_E_mn,vel_f_W_mn,vel_f_N_mn,vel_f_S_mn !vel normal to face end type cell type(cell), allocatable :: cp(:,:),cu(:,:),cv(:,:) !c=cells !For Ghost cells integer :: sta_x,end_x,sta_y,end_y,gk_u,gk_v,gk_u2,gk_v2,gk2 integer :: sta_xu,end_xu,sta_yv,end_yv,no_body_pts ! integer, allocatable :: g_loc_u(:,:),g_loc_v(:,:) real(8) :: sta_cept_x,end_cept_x,sta_cept_y,end_cept_y real(8), allocatable :: body_pts(:,:),body_pts_ini(:,:),body_pts_old(:,:) type ghost_cell integer :: types,types_c !types=1 if cell partially cut by body, -1 if fully cut by body, else 0 if fully outside !types_c=1 if center inside body, else -1 !(end_x-sta_x)*(end_y-sta_y)=ghost cells count, !ik=internal cells count !g_loc=location of ghost cell (i,j index) !g_near=location of 2 nearest fluid cells (1-4) !sta/end_x/y=start/end x/y coord no. !body_c_line=body center line real(8) :: ub,vb,ub_old,vb_old !g_bc=pt on boundary at which GO is normal to boundary segment, ptxy=location of cut on cell, odd=x, even=y coord !inv_B,Bg=inverse B coeff for dirichlet,neumann BC end type ghost_cell type(ghost_cell), allocatable :: g_c_u(:,:),g_c_v(:,:) Vec xx,b_rhs,xx_uv,b_rhs_uv,xx_semi_x,xx_semi_y,b_rhs_semi_x,b_rhs_semi_y ! /* solution vector, right hand side vector and work vector */ Mat A_mat,A_mat_uv,A_semi_x,A_semi_y ! /* sparse matrix */ KSP ksp,ksp_uv,ksp_semi_x,ksp_semi_y !/* linear solver context */ PC pc,pc_uv,pc_semi_x,pc_semi_y PCType ptype KSPType ksptype PetscTruth flgg KSPConvergedReason reason PetscOffset i_vec !for hypre integer(8) grid_hypre, stencil_hypre, A_hypre, b_hypre , x_hypre integer(8) solver_hypre, precond_hypre integer precond_choice,solver_choice !for 2nd airfoil real(8) :: h2x0,h2y0,freq2,phase_ang2,theta20,loc_rot2,ini_ang2,x_1to2,phase_ang_1to2,ld2 real(8) :: air_centx2,air_centx2_old,air_centx2_old2,air_centy2,air_centy2_old,air_centy2_old2 real(8), allocatable :: body_pts2(:,:),body_pts2_ini(:,:),body_pts2_old(:,:) integer :: sta_x2,end_x2,sta_y2,end_y2 integer :: sta_xu2,end_xu2,sta_yv2,end_yv2,no_body_pts2 real(8) :: sta_cept_x2,end_cept_x2,sta_cept_y2,end_cept_y2 !parallel integer :: myid,num_procs,jsta,jend,ksta_p,kend_p,jsta_ext,jend_ext,jsta_ext0 integer :: sta_yp,end_yp,sta_y2p,end_y2p,sta_yvp,end_yvp,sta_yv2p,end_yv2p integer :: ksta_m,kend_m,jsta2,jend2,jend3,ksta_mx,kend_mx,ksta_my,kend_my integer,allocatable :: jjsta_body_u(:),jjlen_body_u(:),jjsta_body_v(:),jjlen_body_v(:) integer,allocatable :: jjsta_body_u2(:),jjlen_body_u2(:),jjsta_body_v2(:),jjlen_body_v2(:) contains subroutine allo_var !allocate memory for variables integer :: status(2),ierr,k,NN total_k=size_x*size_y; 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=1; jend_ext=jend; jsta_ext0=0; jsta2=jsta+1; jend2=jend-1; jend3=jend-2 else if (myid==0) then jsta_ext=jsta; jend_ext=jend+2; jsta_ext0=0; jsta2=jsta+1; jend2=jend; jend3=jend else if (myid==num_procs-1) then jsta_ext=jsta-2; jend_ext=jend; jsta_ext0=jsta_ext; jsta2=jsta; jend2=jend-1; jend3=jend-2 else jsta_ext=jsta-2; jend_ext=jend+2; jsta_ext0=jsta_ext; jsta2=jsta; jend2=jend; jend3=jend end if end if allocate (jjsta_body_u(0:num_procs-1), STAT = status(1)); allocate (jjlen_body_u(0:num_procs-1), STAT = status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (jjsta_body_v(0:num_procs-1), STAT = status(1)); allocate (jjlen_body_v(0:num_procs-1), STAT = status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" if (airfoil_no==2) then allocate (jjsta_body_u2(0:num_procs-1), STAT = status(1)); allocate (jjlen_body_u2(0:num_procs-1), STAT = status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (jjsta_body_v2(0:num_procs-1), STAT = status(1)); allocate (jjlen_body_v2(0:num_procs-1), STAT = status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" end if allocate (x(0:size_x), STAT = status(1)) allocate (y(0:size_y), STAT = status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (xu(0:size_x), STAT = status(1)) allocate (yu(0:size_y), STAT = status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (xv(0:size_x), STAT = status(1)) allocate (yv(0:size_y), STAT = status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (p(size_x,jsta_ext:jend_ext), STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (u(size_x,jsta_ext:jend_ext), STAT=status(1));allocate (v(size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (u_old(size_x,jsta_ext:jend_ext), STAT=status(1));allocate (v_old(size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (diff_pres_x(size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (diff_pres_y(size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (u_ast(size_x,jsta_ext:jend_ext), STAT=status(1));allocate (v_ast(size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (u_ast2(size_x,jsta_ext:jend_ext), STAT=status(1));allocate (v_ast2(size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" ! allocate (in_body_u(0:size_x,0:size_y), STAT=status(2));allocate (in_body_v(0:size_x,0:size_y), STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (phi(size_x,jsta_ext:jend_ext), STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (area_cv_percent_u(size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (area_cv_percent_v(size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (big_A(size_x*(jsta-1)+1:size_x*jend,5), STAT=status(1)) allocate (q_p(size_x*(jsta-1)+1:size_x*jend), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (int_A(size_x*(jsta-1)+1:size_x*jend,5), STAT=status(1)) if (status(1)/=0) STOP "Cannot allocate memory" if (scheme>=1) then ! allocate (impl_mat(2*total_k,11), STAT=status(1)) ;allocate (q_impl(2*total_k), STAT=status(2)) ! if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (int_impl(2*size_x*size_y,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 allocate (fc1_prev_r(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(size_x), 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(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(size_x), 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(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(size_x), 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(size_x), 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(size_x,jsta_ext:jend_ext), STAT=status(1));allocate (cv(size_x,jsta_ext:jend_ext), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (cp(size_x,jsta_ext:jend_ext), STAT=status(1)) !;allocate (c(size_x,size_y), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" ! allocate (g_c(1:size_x,1:size_y), STAT=status(1)) ! if (status(1)/=0) STOP "Cannot allocate memory" allocate (g_c_u(size_x,jsta_ext:jend_ext), STAT=status(1)) allocate (g_c_v(size_x,jsta_ext:jend_ext), STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" if (mom_solver==2) then if (scheme>=1) then !call MatCreateSeqAIJ(PETSC_COMM_SELF,2*size_x*size_y,2*size_x*size_y,11,PETSC_NULL_INTEGER,A_mat_uv,ierr) 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*size_y,b_rhs_uv,ierr) call VecDuplicate(b_rhs_uv,xx_uv,ierr) !call VecAssemblyBegin(b_rhs_uv,ierr) !mpidis !call VecAssemblyEnd(b_rhs_uv,ierr) !call VecAssemblyBegin(xx_uv,ierr) !mpidis !call VecAssemblyEnd(xx_uv,ierr) end if end if if (poisson_solver==2) then !call MatCreateSeqAIJ(PETSC_COMM_SELF,size_x*size_y,size_x*size_y,5,PETSC_NULL_INTEGER,A_mat,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) !call KSPCreate(PETSC_COMM_SELF,ksp,ierr) !call VecAssemblyBegin(b_rhs,ierr) !call VecAssemblyEnd(b_rhs,ierr) !call VecAssemblyBegin(xx,ierr) !call VecAssemblyEnd(xx,ierr) end if !print *,myid,num_procs,jsta,jend,jsta_ext,jend_ext,ksta_m,kend_m,jsta2,jend2,jend3 end subroutine allo_var subroutine semi_allo_var !allocate memory for variables, from time_step>2 for semi-implicit integer :: status(2),ierr,k,NN NN=size_x*size_y ! deallocate (impl_mat, STAT=status(1)) ;deallocate (q_impl, 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" if (mom_solver==2) then call VecDestroy(xx_uv,ierr) call VecDestroy(b_rhs_uv,ierr) call MatDestroy(A_mat_uv,ierr) call KSPDestroy(ksp_uv,ierr) end if allocate (semi_mat_x(size_x*(jsta-1)+1:size_x*jend,5), STAT=status(1)) allocate (semi_mat_y(size_x*(jsta-1)+1:size_x*jend,5), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" ! allocate (q_semi_x(size_x*size_y), STAT=status(1));allocate (q_semi_y(size_x*size_y), STAT=status(2)) ! if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" allocate (int_semi_x(size_x*(jsta-1)+1:size_x*jend,5), STAT=status(1)) allocate (int_semi_y(size_x*(jsta-1)+1:size_x*jend,5), STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" if (mom_solver==2) then !call MatCreateSeqAIJ(PETSC_COMM_SELF,size_x*size_y,size_x*size_y,5,PETSC_NULL_INTEGER,A_semi_x,ierr) !call MatCreateSeqAIJ(PETSC_COMM_SELF,size_x*size_y,size_x*size_y,5,PETSC_NULL_INTEGER,A_semi_y,ierr) call MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,NN,NN,5,PETSC_NULL_INTEGER,5,PETSC_NULL_INTEGER,A_semi_x,ierr) call MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,NN,NN,5,PETSC_NULL_INTEGER,5,PETSC_NULL_INTEGER,A_semi_y,ierr) call MatSetFromOptions(A_semi_x,ierr) call MatSetFromOptions(A_semi_y,ierr) call MatGetOwnershipRange(A_semi_x,ksta_mx,kend_mx,ierr) call MatGetOwnershipRange(A_semi_y,ksta_my,kend_my,ierr) !call VecCreateSeq(PETSC_COMM_SELF,size_x*size_y,b_rhs_semi_x,ierr) call VecCreateMPI(MPI_COMM_WORLD,PETSC_DECIDE,NN,b_rhs_semi_x,ierr) call VecDuplicate(b_rhs_semi_x,xx_semi_x,ierr) call VecDuplicate(b_rhs_semi_x,b_rhs_semi_y,ierr) call VecDuplicate(b_rhs_semi_x,xx_semi_y,ierr) !call VecAssemblyBegin(b_rhs_semi_x,ierr) !call VecAssemblyEnd(b_rhs_semi_x,ierr) !call VecAssemblyBegin(xx_semi_x,ierr) !call VecAssemblyEnd(xx_semi_x,ierr) !call VecAssemblyBegin(b_rhs_semi_y,ierr) !call VecAssemblyEnd(b_rhs_semi_y,ierr) !call VecAssemblyBegin(xx_semi_y,ierr) !call VecAssemblyEnd(xx_semi_y,ierr) !call KSPCreate(PETSC_COMM_SELF,ksp_semi_x,ierr) !call KSPCreate(PETSC_COMM_SELF,ksp_semi_y,ierr) call KSPCreate(MPI_COMM_WORLD,ksp_semi_x,ierr) call KSPCreate(MPI_COMM_WORLD,ksp_semi_y,ierr) end if semi_mat_x=0.;semi_mat_y=0.;q_semi_xx=0.;q_semi_yy=0.;int_semi_x=0;int_semi_y=0 end subroutine semi_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 if (scheme==1) then semi_mat_x=0.;semi_mat_y=0.;q_semi_xx=0.;q_semi_yy=0.;int_semi_x=0;int_semi_y=0 else if (scheme==2) then 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. end if x=0.;y=0.;xu=0.;yv=0.;u_ast=0.;v_ast=0.;u_ast2=0.;v_ast2=0. !;in_body_u=0;in_body_v=0 g_c_u%types=0;g_c_v%types=0;g_c_u%types_c=-1;g_c_v%types_c=-1 g_c_u%ub=0.;g_c_u%ub_old=0.;g_c_v%vb=0.;g_c_v%vb_old=0. sta_x=size_x/2; end_x=size_x/2; sta_xu=size_x/2; end_xu=size_x/2 !guess sta_y=size_y/2; end_y=size_y/2; sta_yv=size_y/2; end_yv=size_y/2 !guess sta_x2=size_x/2; end_x2=size_x/2; sta_xu2=size_x/2; end_xu2=size_x/2 !guess sta_y2=size_y/2; end_y2=size_y/2; sta_yv2=size_y/2; end_yv2=size_y/2 !guess end subroutine ini_var1 subroutine ini_var2 !initialize/allocate memory for variables integer :: status(2) ! allocate (g_loc_u(int((end_x-sta_x)*(end_y-sta_y)*1.5),2), STAT=status(1)) ! allocate (g_loc_v(int((end_x-sta_x)*(end_y-sta_y)*1.5),2), STAT=status(2)) ! if (status(1)/=0 .or. status(2)/=0) STOP "Cannot allocate memory" ! g_loc_u=0;g_loc_v=0 end subroutine ini_var2 subroutine reset_var !reset variables gk_u=0;gk_v=0 !;g_loc_u=0;g_loc_v=0 g_c_u(:,:)%types=0 ! g_c_u(:,:)%ub=0.;g_c_u(:,:)%ub_old=0. g_c_v(:,:)%types=0 ! g_c_v(:,:)%vb=0.;g_c_v(:,:)%vb_old=0. gk_u2=0;gk_v2=0;gk2=0 g_c_u(:,:)%types_c=-1;g_c_v(:,:)%types_c=-1 area_cv_percent_u(:,:)=0.; area_cv_percent_v(:,:)=0. if (airfoil_no==2) then g_c_u(:,:)%types=0 ! g_c_u(:,:)%ub=0.;g_c_u(:,:)%ub_old=0. g_c_v(:,:)%types=0 ! g_c_v(:,:)%vb=0.;g_c_v(:,:)%vb_old=0. g_c_u(:,:)%types_c=-1;g_c_v(:,:)%types_c=-1 area_cv_percent_u(sta_xu2-2:end_xu2+2,sta_y2-2:end_y2+2)=0. area_cv_percent_v(sta_x2-2:end_x2+2,sta_yv2-2:end_yv2+2)=0. end if end subroutine reset_var 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 (u_ast2, STAT=status(1));deallocate (v_ast2, 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)) !;deallocate (big_A, STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" ! deallocate (int_A, STAT=status(1)) ! if (status(1)/=0) STOP "Cannot deallocate memory" if (scheme==1) then ! deallocate (int_semi_x, STAT=status(1));deallocate (int_semi_y, STAT=status(2)) ! if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (semi_diag_x, STAT=status(1));deallocate (semi_diag_y, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" ! deallocate (q_semi_x, STAT=status(1));deallocate (q_semi_y, STAT=status(2)) ! if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" if (mom_solver==2) then call MatDestroy(A_semi_x,ierr) call MatDestroy(A_semi_y,ierr) call VecDestroy(xx_semi_x,ierr) call VecDestroy(b_rhs_semi_x,ierr) call VecDestroy(xx_semi_y,ierr) call VecDestroy(b_rhs_semi_y,ierr) call KSPDestroy(ksp_semi_x,ierr) call KSPDestroy(ksp_semi_y,ierr) end if else if (scheme==2) then ! deallocate (impl_mat, STAT=status(1)) ;deallocate (q_impl, 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" if (mom_solver==2) then call MatDestroy(A_mat_uv,ierr) call KSPDestroy(ksp_uv,ierr) call VecDestroy(xx_uv,ierr) call VecDestroy(b_rhs_uv,ierr) end if end if 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" deallocate (area_cv_percent_u, STAT=status(1));deallocate (area_cv_percent_v, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (g_c_u, STAT=status(1)) deallocate (g_c_v, STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" ! deallocate (g_loc_u, STAT=status(1)) ! deallocate (g_loc_v, STAT=status(2)) ! if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" ! deallocate (in_body_v, STAT=status(1));deallocate (in_body_u, STAT=status(2)) ! if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (body_pts, STAT=status(1));deallocate (body_pts_ini, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (body_pts_old, STAT=status(1)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" if (poisson_solver==2) then call VecDestroy(xx,ierr) call VecDestroy(b_rhs,ierr) call MatDestroy(A_mat,ierr) call KSPDestroy(ksp,ierr) else if (poisson_solver==3) then call HYPRE_StructGridDestroy(grid_hypre,ierr) call HYPRE_StructStencilDestroy(stencil_hypre,ierr) call HYPRE_StructMatrixDestroy(A_hypre,ierr) call HYPRE_StructVectorDestroy(b_hypre,ierr) call HYPRE_StructVectorDestroy(x_hypre,ierr) if (precond_choice==1) then call HYPRE_StructSMGDestroy(precond_hypre,ierr) else if (precond_choice==2) then call HYPRE_StructPFMGDestroy(precond_hypre,ierr) end if if (solver_choice==1) then call HYPRE_StructGMRESDestroy(solver_hypre,ierr) else if (solver_choice==2) then call HYPRE_StructHybridDestroy(solver_hypre,ierr) else if (solver_choice==3) then call HYPRE_StructBiCGSTABDestroy(solver_hypre,ierr) end if end if ! parallel deallocate (jjsta_body_u, STAT=status(1));deallocate (jjlen_body_u, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (jjsta_body_v, STAT=status(1));deallocate (jjlen_body_v, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" if (airfoil_no==2) then deallocate (jjsta_body_u2, STAT=status(1));deallocate (jjlen_body_u2, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" deallocate (jjsta_body_v2, STAT=status(1));deallocate (jjlen_body_v2, STAT=status(2)) if (status(1)/=0 .or. status(2)/=0) STOP "Cannot deallocate memory" end if call PetscFinalize(ierr) end subroutine de_ini_var ! ! .................................................................. ! ! SUBROUTINE PNPOLY ! ! PURPOSE ! TO DETERMINE WHETHER A POINT IS INSIDE A POLYGON ! ! USAGE ! CALL PNPOLY (PX, PY, XX, YY, N, INOUT ) ! ! DESCRIPTION OF THE PARAMETERS ! PX - X-COORDINATE OF POINT IN QUESTION. ! PY - Y-COORDINATE OF POINT IN QUESTION. ! XX - N LONG VECTOR CONTAINING X-COORDINATES OF ! VERTICES OF POLYGON. ! YY - N LONG VECTOR CONTAING Y-COORDINATES OF ! VERTICES OF POLYGON. ! N - NUMBER OF VERTICES IN THE POLYGON. ! INOUT - THE SIGNAL RETURNED: ! -1 IF THE POINT IS OUTSIDE OF THE POLYGON, ! 0 IF THE POINT IS ON AN EDGE OR AT A VERTEX, ! 1 IF THE POINT IS INSIDE OF THE POLYGON. ! ! REMARKS ! THE VERTICES MAY BE LISTED CLOCKWISE OR ANTICLOCKWISE. ! THE FIRST MAY OPTIONALLY BE REPEATED, IF SO N MAY ! OPTIONALLY BE INCREASED BY 1. ! THE INPUT POLYGON MAY BE A COMPOUND POLYGON CONSISTING ! OF SEVERAL SEPARATE SUBPOLYGONS. IF SO, THE FIRST VERTEX ! OF EACH SUBPOLYGON MUST BE REPEATED, AND WHEN CALCULATING ! N, THESE FIRST VERTICES MUST BE COUNTED TWICE. ! INOUT IS THE ONLY PARAMETER WHOSE VALUE IS CHANGED. ! THE SIZE OF THE ARRAYS MUST BE INCREASED IF N > MAXDIM ! WRITTEN BY RANDOLPH FRANKLIN, UNIVERSITY OF OTTAWA, 7/70. ! ! SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED ! NONE ! ! METHOD ! A VERTICAL LINE IS DRAWN THRU THE POINT IN QUESTION. IF IT ! CROSSES THE POLYGON AN ODD NUMBER OF TIMES, THEN THE ! POINT IS INSIDE OF THE POLYGON. ! ! .................................................................. ! 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 ! 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.RESULTS 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