[petsc-users] Can't compile code after upgrading to VS2017 and Intel One API + PETSc 3.15
TAY wee-beng
zonexo at gmail.com
Tue Apr 13 19:58:26 CDT 2021
Hi,
I have just upgraded to VS2017 and Intel One API. Previously, it was
VS2015 + Intel 2018 compilers + PETSc 3.9.3.
Now, I realise that I can't even compile. I tried both PETSc 3.9.3 and
3.15.0. I didn't change other things. I have attached the code. The
first few lines are:
/*module global_data*//*
*//*
*//*#include "petsc/finclude/petsc.h"*//*
*//*
*//*use petsc*//*
*//*
*//*use kdtree2_module*//*
*//*
*//*implicit none*//*
*//*
*//*save*//*
*//*
*//*PetscInt :: ...*//*
*//*
*//*PetscErrorCode :: ierr*/
The error is:
/*Compiling with Intel® Fortran Compiler Classic 2021.2.0 [Intel(R)
64]...*//*
*//*global.F90*//*
*//*global.F90(3): #error: can't find include file:
petsc/finclude/petsc.h*//*
*//*global.F90(937): #error: can't find include file:
petsc/finclude/petsc.h*/
--
Thank you very much.
Yours sincerely,
================================================
TAY Wee-Beng 郑伟明 (Zheng Weiming)
Personal research webpage: _http://tayweebeng.wixsite.com/website_
Youtube research showcase: _https://goo.gl/PtvdwQ_
linkedin: _https://www.linkedin.com/in/tay-weebeng_
================================================
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210414/523cf383/attachment-0001.html>
-------------- next part --------------
module global_data
#include "petsc/finclude/petsc.h"
use petsc
use kdtree2_module
implicit none
save
!grid variables
PetscInt :: size_x,size_y,size_z,grid_type,status(4),inout_test,terminate !size_x1,size_x2,size_x3,size_y1,size_y2,size_y3
PetscErrorCode :: ierr
PetscInt :: sta_ijk_uniform(3),end_ijk_uniform(3),one,three
real(8), parameter :: pi = 3.141592653589793238462643383279502884197d0
real(8) :: min_delta,min_delta_x,min_delta_x2,min_delta_y,min_delta_z,top_bottom_x_dist,mid_x_dist,mid_x_dist2,top_bottom_y_dist,mid_y_dist,top_x_dist,bottom_x_dist
real(8) :: upstream_z_dist,midstream_z_dist,downstream_z_dist,total_x_dist,uniform_z_sta
real(8), allocatable :: x(:),y(:),z(:),xu(:),yu(:),zu(:),xv(:),yv(:),zv(:),xw(:),yw(:),zw(:),c_cx(:),cu_cx(:),c_cy(:),cv_cy(:),c_cz(:),cw_cz(:)
!solver variables
logical :: tecplot_grid_write
PetscInt :: Total_time_step,new_start,interval,safe_int,safe_int2,OS,file_read,rec_start,rec_quasi,vort_plot,q_criterion_plot,avg_q_criterion_cycle,start_avg_q_criterion_cycle
PetscInt :: steady,quasi_steady,Total_ijk,time,poisson_solver,ibm_choice,start_time,time_scheme,x_symmetry,y_symmetry,boundary_order,reuse,convex_body
PetscInt :: total_snapshots,file_write_no,file_write_no2,file_write_no3,cycle_no,snapshots_no,avg_vorticity,spatial_scheme,spatial_scheme_initial,fixed_time,fixed_interval
PetscInt :: count_rate_int,start_simulation,end_simulation,no_test,mpi_no_test,cycles_initial,new_start_initial,vort_plot_initial,body3_unsteady_initial
PetscInt :: total_txyz_uvwp_snapshots
real(8) :: count_rate_real,simulation_runtime,start_simulation_sec,end_simulation_sec,del_t_initial,coef_tol
!boundary_order = 1/2/3/4/5/6 = 1st order, normal/ 1st order, simple/ 1st order, simple (no RHS update)/
!2nd order, normal/ 2nd order, simple/ 2nd order, simple (no RHS update)
!x_symmetry = 0/1/2/3/4/5/6 = no symmetry/sym at west/sym at east,west/sym at west,wall at east/
!sym at west,wall at east but larger x dist/sym at west and south/same as 4, except sym at east
PetscInt :: interpolation_choice,IB_force_choice,debug_lvl,CFL_check,cycles,write_xyz_check,write_misc_check
real(8) :: f_coef_ratio,vorticity_interval_time,quasi_steady_time,beta_m,inv_beta_m,save_txyz_uvwp_interval_time
!ibm_choice = 1 (Kim), 2 (Yang)
real(8) :: CFL, Re,inv_Re,u_inf_initial,u_inf,time_sta,act_time,Uc,Uc_old,del_t,del_t_old,last_period
real(8) :: t1_Eldredge,t2_Eldredge,t3_Eldredge,t3_Eldredge_initial,t4_Eldredge,a_const_Eldredge,a_max_Eldredge,max_body_vel(3)
!last_period = number of period to record, can be half etc
!max_b_vel = max body velocity
real(8) :: Uc_e,Uc_w,Uc_n,Uc_s !convective velocity for east, west, north, south
!FFT variables
PetscInt :: fft_Modes
real(8), allocatable :: fft_coeff_a(:,:,:),fft_coeff_b(:,:,:),cosF(:),sinF(:),discrete_time(:),fft_coeff_a2(:,:,:),fft_coeff_b2(:,:,:),fft_coeff_a3(:,:,:),fft_coeff_b3(:,:,:)
real(8), allocatable :: n_fft_coeff_a(:,:,:,:),n_fft_coeff_b(:,:,:,:)
!Kriging variables
PetscInt :: krig_span_no,krig_chord_no,krig_raw_data_no,krig_raw_pts_no,krig_skip_no
type krig_model
PetscInt :: T
real(8) :: mu(3),sig(3),theta(3),dt
real(8), allocatable :: scti(:,:),sct(:,:),y0(:,:),xyzi(:,:)
end type krig_model
type(krig_model) :: krigmodel
!airfoil variables
PetscInt :: airfoil_config,hover,motion_type,wing_config,body_input_type,no_body,deformation,body_multi_no,cross_sect_sym
real(8) :: freq_max,deform_s_max,deform_c_max
!deform_s/c_max = spanwise/chordwise deformation max
!>motion_type = 1,2 = linear,rotation
!1st body
PetscInt :: surface_normal_correct1,no_line_pts1
real(8) :: h1y0,k0,freq1,phase_ang1,theta10_x,theta10_y,theta10_z,loc_rot1,leading_edge_pos1 !center of rotation wrt body
!>theta0_x,theta0_z correspond to max flapping,pitching angles
!>airfoil_config=1 for 1 airfoil, 2 for tandem, 3 for biplane, 4 for biplane + tail
real(8) :: plate_h_tk1,plate_len1a,plate_len1b,plate_wid1
!Plate half thickness, length, width for body_input_type = 3
real(8) :: h1x0,h1z0,ini_h1x0,ini_h1y0,ini_theta10_x,ini_theta10_y,ini_theta10_z,body_cg1_ini(3),body_cg1(3),body_cg1_old(3),body_cg1_old2(3)
real(8) :: rot_body1(3),rot_body1_old(3),fixed_rot_body1(3) !center of rotation wrt body, fixed center of rotation wrt body at initial starting
real(8) :: pvm_xyz1(9),cd_cl_cs_mom_implicit1(6),cd_cl_cs_mom_explicit1(6),cd_cl_cs_mom_power1(8)
real(8) :: cd_tmp,cl_tmp,cs_tmp,body_vel1(3),body_vel1_old(3),body_displacement1(3)
PetscInt :: IIB_equal_cell_no_u1_max,I_equal_cell_no_u1_max,IIB_cell_no_u1_max,I_cell_no_u1_max,IIB_global_cell_no_u1_max,I_global_cell_no_u1_max
PetscInt :: IIB_equal_cell_no_u2_max,I_equal_cell_no_u2_max,IIB_cell_no_u2_max,I_cell_no_u2_max,IIB_global_cell_no_u2_max,I_global_cell_no_u2_max
PetscInt :: IIB_equal_cell_no_u3_max,I_equal_cell_no_u3_max,IIB_cell_no_u3_max,I_cell_no_u3_max,IIB_global_cell_no_u3_max,I_global_cell_no_u3_max
PetscInt :: IIB_equal_cell_no_u1_max_max,I_equal_cell_no_u1_max_max,IIB_cell_no_u1_max_max,I_cell_no_u1_max_max,IIB_global_cell_no_u1_max_max,I_global_cell_no_u1_max_max
PetscInt :: final_IIB_equal_no1,final_I_equal_no1,final_IIB_no1,final_I_no1,final_IIB_global_no1,final_I_global_no1
real(8) :: shift_x
!for JFM2007 paper's kinematics
real(8) :: d2r,azi_amp,azi_ph,azi_ampoff,kshape,lat_amp,lat_ph,lat_ampoff,nshape,pit_amp,pit_ph,pit_ampoff,cshape
!>fx_p,fy_p,fz_p = pvm_xyz(1:3), pvm = pressure,viscous,moment
!>fx_v,fy_v,fz_v = pvm_xyz(4:6)
!>mom_x,mom_y,mom_z = pvm_xyz(7:9)
!2nd body
PetscInt :: surface_normal_correct2,no_line_pts2,body2_unsteady
real(8) :: h2y0,freq2,phase_ang2,theta20_x,theta20_y,theta20_z,loc_rot2,leading_edge_pos2 !center of rotation wrt body
real(8) :: plate_h_tk2,plate_len2a,plate_len2b,plate_wid2
real(8) :: shift_x2,z_1_to_2,phase_ang_1_to_2 !center of rotation wrt body
!>theta0_x,theta0_z correspond to max flapping,pitching angles
real(8) :: h2x0,h2z0,ini_h2x0,ini_h2y0,ini_theta20_x,ini_theta20_y,ini_theta20_z,body_cg2_ini(3),body_cg2(3),body_cg2_old(3),body_cg2_old2(3)
real(8) :: rot_body2(3),rot_body2_old(3),fixed_rot_body2(3) !center of rotation wrt body, fixed center of rotation wrt body at initial starting
real(8) :: pvm_xyz2(9),cd_cl_cs_mom_implicit2(6),cd_cl_cs_mom_explicit2(6),cd_cl_cs_mom_power2(8)
real(8) :: cd_tmp2,cl_tmp2,cs_tmp2,body_vel2(3),body_vel2_old(3),body_displacement2(3)
PetscInt :: IIB_equal_cell_no_u2_max_max,I_equal_cell_no_u2_max_max,IIB_cell_no_u2_max_max,I_cell_no_u2_max_max,IIB_global_cell_no_u2_max_max,I_global_cell_no_u2_max_max
PetscInt :: final_IIB_equal_no2,final_I_equal_no2,final_IIB_no2,final_I_no2,final_IIB_global_no2,final_I_global_no2
!3rd body
PetscInt :: surface_normal_correct3,no_line_pts3,body3_unsteady
real(8) :: h3y0,freq3,phase_ang3,theta30_x,theta30_y,theta30_z,loc_rot3,leading_edge_pos3 !center of rotation wrt body
real(8) :: plate_h_tk3,plate_len3a,plate_len3b,plate_wid3
real(8) :: shift_x3,z_1_to_3,phase_ang_1_to_3 !center of rotation wrt body
!>theta0_x,theta0_z correspond to max flapping,pitching angles
real(8) :: h3x0,h3z0,ini_h3x0,ini_h3y0,ini_theta30_x,ini_theta30_y,ini_theta30_z,body_cg3_ini(3),body_cg3(3),body_cg3_old(3),body_cg3_old2(3)
real(8) :: rot_body3(3),rot_body3_old(3),fixed_rot_body3(3) !center of rotation wrt body, fixed center of rotation wrt body at initial starting
real(8) :: pvm_xyz3(9),cd_cl_cs_mom_implicit3(6),cd_cl_cs_mom_explicit3(6),cd_cl_cs_mom_power3(8)
real(8) :: cd_tmp3,cl_tmp3,cs_tmp3,body_vel3(3),body_vel3_old(3),body_displacement3(3)
PetscInt :: IIB_equal_cell_no_u3_max_max,I_equal_cell_no_u3_max_max,IIB_cell_no_u3_max_max,I_cell_no_u3_max_max,IIB_global_cell_no_u3_max_max,I_global_cell_no_u3_max_max
PetscInt :: final_IIB_equal_no3,final_I_equal_no3,final_IIB_no3,final_I_no3,final_IIB_global_no3,final_I_global_no3
!n bodies
PetscInt, allocatable :: n_surface_normal_correct(:),n_no_line_pts(:)
real(8), allocatable :: n_hy0(:),n_freq(:),n_phase_ang(:),n_theta0_x(:),n_theta0_y(:),n_theta0_z(:),n_loc_rot(:),n_leading_edge_pos(:) !center of rotation wrt body
real(8), allocatable :: n_plate_h_tk(:),n_plate_lena(:),n_plate_lenb(:),n_plate_wid(:)
real(8), allocatable :: n_shift_x(:),z_1_to_n(:),phase_ang_1_to_n(:) !center of rotation wrt body
real(8), allocatable :: n_hx0(:),n_hz0(:),n_ini_hx0(:),n_ini_hy0(:),n_ini_theta0_x(:),n_ini_theta0_y(:),n_ini_theta0_z(:),n_body_cg_ini(:,:),n_body_cg(:,:),n_body_cg_old(:,:),n_body_cg_old2(:,:)
real(8), allocatable :: n_rot_body(:,:),n_rot_body_old(:,:),n_fixed_rot_body(:,:) !center of rotation wrt body, fixed center of rotation wrt body at initial starting
real(8), allocatable :: n_pvm_xyz(:,:),n_cd_cl_cs_mom_implicit(:,:),n_cd_cl_cs_mom_explicit(:,:),n_cd_cl_cs_mom_power(:,:)
real(8), allocatable :: n_cd_tmp(:),n_cl_tmp(:),n_cs_tmp(:),n_body_vel(:,:),n_body_vel_old(:,:),n_body_displacement(:,:)
!real(8), allocatable :: n_hx0(no_body),n_hz0(no_body),n_ini_hx0(no_body),n_ini_hy0(no_body),n_ini_theta0_x(no_body),n_ini_theta0_y(no_body),n_ini_theta0_z(no_body),n_body_cg_ini(3,no_body),n_body_cg(3,no_body),n_body_cg_old(3,no_body),n_body_cg_old2(3,no_body)
!real(8), allocatable :: n_rot_body(3,no_body),n_rot_body_old(3,no_body),n_fixed_rot_body(3,no_body) !center of rotation wrt body, fixed center of rotation wrt body at initial starting
!real(8), allocatable :: n_pvm_xyz(9,no_body),n_cd_cl_cs_mom_implicit(6,no_body),n_cd_cl_cs_mom_explicit(6,no_body),n_cd_cl_cs_mom_power(8,no_body)
!real(8), allocatable :: n_cd_tmp(no_body),n_cl_tmp(no_body),n_cs_tmp(no_body),n_body_vel(3,no_body),n_body_vel_old(3,no_body),n_body_displacement(3,no_body)
!PetscInt, allocatable :: n_IIB_equal_cell_no_u_max_max(:),n_I_equal_cell_no_u_max_max(:),n_IIB_cell_no_u_max_max(:),n_I_cell_no_u_max_max(:),n_IIB_global_cell_no_u_max_max(:),n_I_global_cell_no_u_max_max(:)
PetscInt :: n_IIB_equal_cell_no_u_max_max,n_I_equal_cell_no_u_max_max,n_IIB_cell_no_u_max_max,n_I_cell_no_u_max_max,n_IIB_global_cell_no_u_max_max,n_I_global_cell_no_u_max_max
PetscInt :: n_final_IIB_equal_no,n_final_I_equal_no,n_final_IIB_no,n_final_I_no,n_final_IIB_global_no,n_final_I_global_no
PetscInt :: n_final_IIB_equal_no_cur,n_final_I_equal_no_cur,n_final_IIB_no_cur,n_final_I_no_cur,n_final_IIB_global_no_cur,n_final_I_global_no_cur
!ibm variables
PetscInt :: near_2_body_check,near_2_body12_global(4),near_2_body12_global_tmp(4),near_2_body23_global(4),near_2_body23_global_tmp(4),IIB_IB_sf_global,IIB_IB_sf,IIB_IB_sf_equal,initial_check,type_option
!near_2_body12_global = between body 1 & 2, near_2_body23_global = between body 2 & 3
PetscInt :: body_line,delfly_time_steps,check_IIB_I_cell_no_size_int !body=1,line=2
PetscInt :: IIB_I_sta_domain(3),IIB_I_end_domain(3) !body=1,line=2
!IIB_I_sta_domain etc = min/max global domain index
real(8) :: min_element_length,max_element_length
!flow variables
real(8), allocatable :: p(:,:,:),u(:,:,:),v(:,:,:),w(:,:,:),u_old(:,:,:),v_old(:,:,:),w_old(:,:,:),phi(:,:,:)
!real(8), allocatable :: u_old(:,:,:),v_old(:,:,:),w_old(:,:,:),phi(:,:,:)
real(8), allocatable :: u_ast(:,:,:), v_ast(:,:,:), w_ast(:,:,:),u_ast2(:,:,:), v_ast2(:,:,:), w_ast2(:,:,:)
real(4), allocatable :: Q_Soria(:,:,:),Q_criterion(:,:,:),Q_Jeong(:,:,:),vorticity_avg(:,:,:,:),vorticity_tmp(:,:,:) !,q_crit(:,:,:)
real(4), allocatable :: Q_criterion_avg(:,:,:,:),Q_criterion_tmp(:,:,:)
!iterative solver variables
real(8), allocatable :: big_A(:,:),qk(:) !,q_p(:)
PetscInt, allocatable :: int_A(:,:),int_semi_xyz(:,:),int_index_xyz(:)
real(8), allocatable :: diff_pres_x(:,:,:),diff_pres_y(:,:,:),diff_pres_z(:,:,:)
!real(8), allocatable :: semi_diag_x(:),semi_diag_y(:),semi_diag_z(:)
real(8), allocatable :: semi_mat_xyz(:,:),semi_mat_xyz_da(:) !,q_semi_vect_xyz(:) !,q_semi_x_force(:,:,:),q_semi_y_force(:,:,:),q_semi_z_force(:,:,:)
!real(8), allocatable :: fc1_prev_r(:,:),fc1_prev_c(:,:),fc1_prev_rc(:,:),fc2_prev_r(:,:),fc2_prev_c(:,:),fc2_prev_rc(:,:)
real(8) :: fc1(6),fc2(6),fc3(6),fc1_old(6),fc2_old(6),fc3_old(6),fd1(6),fd2(6),fd3(6),fc_prev_e(3),fc_prev_e_old(3),fd_prev_e(3)
real(8), allocatable :: fc_prev_n(:,:),fc_prev_n_old(:,:),fd_prev_n(:,:),fc_prev_f(:,:,:),fc_prev_f_old(:,:,:),fd_prev_f(:,:,:)
!real(8), allocatable :: fd1_prev_r(:,:),fd1_prev_c(:,:),fd1_prev_rc(:,:),fd2_prev_r(:,:),fd2_prev_c(:,:),fd2_prev_rc(:,:)
!real(8), allocatable :: fd1_prev_r_old(:,:),fd1_prev_c_old(:,:),fd1_prev_rc_old(:,:),fd2_prev_r_old(:,:),fd2_prev_c_old(:,:),fd2_prev_rc_old(:,:)
real(8), allocatable :: s_cf_u_e(:,:,:),s_cf_v_e(:,:,:),s_cf_u_n(:,:,:),s_cf_v_n(:,:,:)
real(8), allocatable :: s_cf_w_e(:,:,:),s_cf_w_n(:,:,:),s_cf_u_f(:,:,:),s_cf_v_f(:,:,:),s_cf_w_f(:,:,:)
!s_cf = spatial upwind scheme coefficients
!FSI variables
PetscInt :: fsi_run,fsi_iter_no,fsi_iteration_completed
real(8) :: fsi_tol,fsi_max_tol,spring_const,theta_rot_fsi,theta_rot_fsi_old,theta_rot_fsi_old2,theta_rot_fsi_old3,initial_origin_pt(3),origin_pt(3),span_pt(3),thick_pt(3),trail_pt(3)
real(8) :: omega_rot_fsi,omega_rot_fsi_old,alpha_x_rot_fsi,alpha_x_rot_fsi_old,theta_rot_fsi_inner_iter,theta_rot_fsi_inner_iter_old,theta_rot_fsi_inner_iter_old2,theta_rot_fsi_iter1
real(8) :: akin_parameter,inertia_const,fsi_qt(2),fsi_qt1(2),fsi_qt2(2),fsi_qt3(2),fsi_qt4(2),fsi_qk1(2),fsi_qk(2)
real(8) :: fsi_zt(2),fsi_zt1(2),fsi_zt2(2),fsi_zt3(2),fsi_zt4(2),fsi_pz(2),fsi_zk(2),fsi_zk1(2),fsi_et(2),fsi_et1(2)
real(8) :: fsi_z_tol(2)
!k1 = k - 1, t1 = t - 1
!Optimization flexibility study
real(8) :: flex_center,nom_lead_flex,nom_trail_flex,lead_flex,trail_flex,phase_flex_ang
!Lua's FF SYR motion
PetscInt :: lua_data_no,lua_data_no_time,lua_prev_time_location !remember old time position for faster access
!lua_data_no = no. of data, including repeated time, lua_data_no_time = actual no. of data w/o repeated time
real(8), allocatable :: ff_syr_time_period(:),ff_syr_theta_x_rad(:),ff_syr_theta_z_rad(:)
type IIB_I_cell
logical :: check
!PetscInt(1) :: types,types_old
PetscInt :: types,types_old
PetscInt :: s_no,e_no,v_no,s_no2,v_no2 !type01 = 0 outside, 1 inside, s/e/v_no = surface/edge/vertex no
end type IIB_I_cell
type(IIB_I_cell), allocatable :: IIB_I_cp(:,:,:),IIB_I_cu(:,:,:),IIB_I_cv(:,:,:),IIB_I_cw(:,:,:)
type cell
!PetscInt :: types,types_old,check,s_no,v_no !type01 = 0 outside, 1 inside, s/v_no = surface/vertex no
PetscInt :: types
real(8) :: vol,dr_E,dr_N,dr_F,dr_W,dr_S,dr_B !,vel,area_cv_percent,ne_x,nw_x,se_x,sw_x,ne_y,nw_y,se_y,sw_y
end type cell
type(cell), allocatable :: cp(:,:,:),cu(:,:,:),cv(:,:,:),cw(:,:,:)
type cellx
real(8) :: cx
real(8) :: pd_E,pd_W,ratio_E,ratio_W
end type cellx
type(cellx), allocatable :: cp_x(:),cu_x(:),cv_x(:),cw_x(:)
type celly
real(8) :: cy
real(8) :: pd_N,pd_S,ratio_N,ratio_S
end type celly
type(celly), allocatable :: cp_y(:),cu_y(:),cv_y(:),cw_y(:)
type cellz
real(8) :: cz
real(8) :: pd_F,pd_B,ratio_F,ratio_B
end type cellz
type(cellz), allocatable :: cp_z(:),cu_z(:),cv_z(:),cw_z(:)
type cellxy
real(8) :: fc_F !for z cells use
end type cellxy
type(cellxy), allocatable :: cp_xy(:,:),cu_xy(:,:),cv_xy(:,:),cw_xy(:,:)
type cellyz
real(8) :: fc_E !for x cells use
end type cellyz
type(cellyz), allocatable :: cp_yz(:,:),cu_yz(:,:),cv_yz(:,:),cw_yz(:,:)
type cellzx
real(8) :: fc_N !for y cells use
end type cellzx
type(cellzx), allocatable :: cp_zx(:,:),cu_zx(:,:),cv_zx(:,:),cw_zx(:,:)
!c=cells
!For IBM
!real(8), allocatable :: ibm_force(:,:,:,:)
PetscInt, allocatable :: IIB_index_u(:,:,:),IIB_index_v(:,:,:),IIB_index_w(:,:,:),IIB_index_p(:,:,:)
!IIB/I_equal_cell_no_global_u/v/w/p1 represent no. of equal cells across diff procs, known to all procs
!1st body
!cells variables
PetscInt :: IIB_cell_no_u1,IIB_cell_no_v1,IIB_cell_no_w1,IIB_cell_no_p1,I_cell_no_u1,I_cell_no_v1,I_cell_no_w1,I_cell_no_p1,pIB_cell_no_u1,pIB_cell_no_v1,pIB_cell_no_w1,pIB_cell_no_p1
PetscInt :: IIB_global_cell_no_u1,IIB_global_cell_no_v1,IIB_global_cell_no_w1,IIB_global_cell_no_p1,I_global_cell_no_u1,I_global_cell_no_v1,I_global_cell_no_w1,I_global_cell_no_p1
PetscInt :: IIB_equal_cell_no_u1,IIB_equal_cell_no_v1,IIB_equal_cell_no_w1,IIB_equal_cell_no_p1,I_equal_cell_no_u1,I_equal_cell_no_v1,I_equal_cell_no_w1,I_equal_cell_no_p1
PetscInt :: IIB_I_cell_no_uvw_total1(6),IIB_I_cell_no_p_total1(2)
PetscInt, allocatable :: IIB_I_cell_no_uvw_global1(:),IIB_I_cell_no_p_global1(:)
PetscInt, allocatable :: IIB_equal_cell_no_global_u1(:),IIB_equal_cell_no_global_v1(:),IIB_equal_cell_no_global_w1(:),IIB_equal_cell_no_global_p1(:)
PetscInt, allocatable :: I_equal_cell_no_global_u1(:),I_equal_cell_no_global_v1(:),I_equal_cell_no_global_w1(:),I_equal_cell_no_global_p1(:)
!body variables
PetscInt :: sta_x1,end_x1,sta_y1,end_y1,sta_z1,end_z1,sta_xu1,end_xu1,sta_yv1,end_yv1,sta_zw1,end_zw1,no_vertices1,no_surfaces1,normal_inout
PetscInt :: start_ijk_IIB1(3),end_ijk_IIB1(3),start_ijk_int_IIB1(3),end_ijk_int_IIB1(3),start_ijk_ghost_IIB1(3),end_ijk_ghost_IIB1(3)
!sta_x1 = starting index for 1st body (of entire domain)
!start/end_ijk_IIB1 -starting / ending indices for 1st body scan, based on IIB domain
!start_ijk_int_IIB refers to internal
PetscInt, allocatable :: sta_yy1(:),end_yy1(:)
PetscInt :: no_body_pts1
real(8), allocatable :: body_pts1(:,:),body_pts1_ini(:,:),body_pts1_old(:,:),body_length_seg1(:)
real(8) :: sta_cept_xyz1(3),end_cept_xyz1(3)
real(8), allocatable :: IIB_IB_vel1(:,:),IIB_IB_vel1_old(:,:)
!for 2nd body
!cells variables
PetscInt :: IIB_cell_no_u2,IIB_cell_no_v2,IIB_cell_no_w2,IIB_cell_no_p2,I_cell_no_u2,I_cell_no_v2,I_cell_no_w2,I_cell_no_p2
PetscInt :: IIB_global_cell_no_u2,IIB_global_cell_no_v2,IIB_global_cell_no_w2,IIB_global_cell_no_p2,I_global_cell_no_u2,I_global_cell_no_v2,I_global_cell_no_w2,I_global_cell_no_p2
PetscInt :: IIB_equal_cell_no_u2,IIB_equal_cell_no_v2,IIB_equal_cell_no_w2,IIB_equal_cell_no_p2,I_equal_cell_no_u2,I_equal_cell_no_v2,I_equal_cell_no_w2,I_equal_cell_no_p2
PetscInt :: IIB_I_cell_no_uvw_total2(6),IIB_I_cell_no_p_total2(2)
PetscInt, allocatable :: IIB_I_cell_no_uvw_global2(:),IIB_I_cell_no_p_global2(:)
PetscInt, allocatable :: IIB_equal_cell_no_global_u2(:),IIB_equal_cell_no_global_v2(:),IIB_equal_cell_no_global_w2(:),IIB_equal_cell_no_global_p2(:)
PetscInt, allocatable :: I_equal_cell_no_global_u2(:),I_equal_cell_no_global_v2(:),I_equal_cell_no_global_w2(:),I_equal_cell_no_global_p2(:)
!body variables
PetscInt :: sta_x2,end_x2,sta_y2,end_y2,sta_z2,end_z2,sta_xu2,end_xu2,sta_yv2,end_yv2,sta_zw2,end_zw2,no_vertices2,no_surfaces2
PetscInt :: start_ijk_IIB2(3),end_ijk_IIB2(3),start_ijk_int_IIB2(3),end_ijk_int_IIB2(3),start_ijk_ghost_IIB2(3),end_ijk_ghost_IIB2(3)
PetscInt, allocatable :: sta_yy2(:),end_yy2(:)
PetscInt :: no_body_pts2
real(8), allocatable :: body_pts2(:,:),body_pts2_ini(:,:),body_pts2_old(:,:),body_length_seg2(:)
real(8) :: sta_cept_xyz2(3),end_cept_xyz2(3)
real(8), allocatable :: IIB_IB_vel2(:,:),IIB_IB_vel2_old(:,:)
!for 3rd body
!cells variables
PetscInt :: IIB_cell_no_u3,IIB_cell_no_v3,IIB_cell_no_w3,IIB_cell_no_p3,I_cell_no_u3,I_cell_no_v3,I_cell_no_w3,I_cell_no_p3
PetscInt :: IIB_global_cell_no_u3,IIB_global_cell_no_v3,IIB_global_cell_no_w3,IIB_global_cell_no_p3,I_global_cell_no_u3,I_global_cell_no_v3,I_global_cell_no_w3,I_global_cell_no_p3
PetscInt :: IIB_equal_cell_no_u3,IIB_equal_cell_no_v3,IIB_equal_cell_no_w3,IIB_equal_cell_no_p3,I_equal_cell_no_u3,I_equal_cell_no_v3,I_equal_cell_no_w3,I_equal_cell_no_p3
PetscInt :: IIB_I_cell_no_uvw_total3(6),IIB_I_cell_no_p_total3(2)
PetscInt, allocatable :: IIB_I_cell_no_uvw_global3(:),IIB_I_cell_no_p_global3(:)
PetscInt, allocatable :: IIB_equal_cell_no_global_u3(:),IIB_equal_cell_no_global_v3(:),IIB_equal_cell_no_global_w3(:),IIB_equal_cell_no_global_p3(:)
PetscInt, allocatable :: I_equal_cell_no_global_u3(:),I_equal_cell_no_global_v3(:),I_equal_cell_no_global_w3(:),I_equal_cell_no_global_p3(:)
!body variables
PetscInt :: sta_x3,end_x3,sta_y3,end_y3,sta_z3,end_z3,sta_xu3,end_xu3,sta_yv3,end_yv3,sta_zw3,end_zw3,no_vertices3,no_surfaces3
PetscInt :: start_ijk_IIB3(3),end_ijk_IIB3(3),start_ijk_int_IIB3(3),end_ijk_int_IIB3(3),start_ijk_ghost_IIB3(3),end_ijk_ghost_IIB3(3)
PetscInt, allocatable :: sta_yy3(:),end_yy3(:)
PetscInt :: no_body_pts3
real(8), allocatable :: body_pts3(:,:),body_pts3_ini(:,:),body_pts3_old(:,:),body_length_seg3(:)
real(8) :: sta_cept_xyz3(3),end_cept_xyz3(3)
real(8), allocatable :: IIB_IB_vel3(:,:),IIB_IB_vel3_old(:,:)
!n bodies
!cells variables
PetscInt, allocatable :: n_IIB_cell_no_u(:),n_IIB_cell_no_v(:),n_IIB_cell_no_w(:),n_IIB_cell_no_p(:),n_I_cell_no_u(:),n_I_cell_no_v(:),n_I_cell_no_w(:),n_I_cell_no_p(:)
PetscInt, allocatable :: n_IIB_global_cell_no_u(:),n_IIB_global_cell_no_v(:),n_IIB_global_cell_no_w(:),n_IIB_global_cell_no_p(:),n_I_global_cell_no_u(:),n_I_global_cell_no_v(:),n_I_global_cell_no_w(:),n_I_global_cell_no_p(:)
PetscInt, allocatable :: n_IIB_equal_cell_no_u(:),n_IIB_equal_cell_no_v(:),n_IIB_equal_cell_no_w(:),n_IIB_equal_cell_no_p(:),n_I_equal_cell_no_u(:),n_I_equal_cell_no_v(:),n_I_equal_cell_no_w(:),n_I_equal_cell_no_p(:)
PetscInt, allocatable :: n_IIB_equal_cell_no_u_max(:),n_I_equal_cell_no_u_max(:),n_IIB_cell_no_u_max(:),n_I_cell_no_u_max(:),n_IIB_global_cell_no_u_max(:),n_I_global_cell_no_u_max(:)
PetscInt, allocatable:: n_IIB_I_cell_no_uvw_total(:,:),n_IIB_I_cell_no_p_total(:,:)
!PetscInt :: IIB_I_cell_no_uvw_total3(6),IIB_I_cell_no_p_total3(2)
PetscInt, allocatable :: n_IIB_I_cell_no_uvw_global(:,:),n_IIB_I_cell_no_p_global(:,:)
PetscInt, allocatable :: n_IIB_equal_cell_no_global_u(:,:),n_IIB_equal_cell_no_global_v(:,:),n_IIB_equal_cell_no_global_w(:,:),n_IIB_equal_cell_no_global_p(:,:)
PetscInt, allocatable :: n_I_equal_cell_no_global_u(:,:),n_I_equal_cell_no_global_v(:,:),n_I_equal_cell_no_global_w(:,:),n_I_equal_cell_no_global_p(:,:)
!body variables
PetscInt, allocatable :: n_sta_x(:),n_end_x(:),n_sta_y(:),n_end_y(:),n_sta_z(:),n_end_z(:)
PetscInt, allocatable :: n_no_vertices(:),n_no_surfaces(:)
!n_sta_x etc = starting index of bodies along x direction in global domain
PetscInt, allocatable :: n_body_sta_ijk(:,:),n_body_end_ijk(:,:),n_body_sta_ijk_ghost(:,:),n_body_end_ijk_ghost(:,:)
!n_body_sta_ijk etc = starting index along x direction in global MPI domain
PetscInt, allocatable :: n_body_sta_ijk_IIB(:,:),n_body_end_ijk_IIB(:,:),n_body_sta_ijk_IIB_ghost(:,:),n_body_end_ijk_IIB_ghost(:,:)
!n_body_sta_ijk_IIB etc = starting index along x direction in IIB MPI domain
PetscInt, allocatable :: n_body_sta_ijk_IIB_array(:,:),n_body_end_ijk_IIB_array(:,:)
!n_body_sta_ijk_IIB_array etc = starting index along x direction in IIB MPI domain, used only for creating and specifying limit/range of array
!used becos array must be created or checked at range bigger than n_body_sta_ijk_IIB. only difference for edge cells
!PetscInt, allocatable :: n_body_sta_ijk_IIB_ext(:,:),n_body_end_ijk_IIB_ext(:,:)
!n_body_sta_ijk_IIB_ext etc = starting index along x direction in IIB MPI domain, used only for checking in/out cells
!used becos array must be created or checked at range bigger than n_body_sta_ijk_IIB. only difference for edge cells
PetscInt, allocatable :: n_sta_yy(:,:),n_end_yy(:,:)
PetscInt, allocatable :: n_no_body_pts(:)
real(8), allocatable :: n_body_pts(:,:,:),n_body_pts_ini(:,:,:),n_body_pts_old(:,:,:),n_body_length_seg(:,:)
real(8), allocatable :: n_sta_cept_xyz(:,:),n_end_cept_xyz(:,:)
!real(8) :: sta_cept_xyz3(3),end_cept_xyz3(3)
real(8), allocatable :: n_IIB_IB_vel(:,:,:),n_IIB_IB_vel_old(:,:,:)
!>normal_inout = normal inout, if normal pt in/outside body, normal_inout = -1/1
!>for some bodies which have errors
!real(8), allocatable :: vol_percent(:,:,:,:),pIB_IB_vel(:,:),pIB_IB_vel_old(:,:),IIB_IB_vel1(:,:),IIB_IB_vel1_old(:,:)
real(8), allocatable :: pIB_IB_vel(:,:),pIB_IB_vel_old(:,:)
!pIB_IB_vel = IB u,v,w velocity for pIB p cells
type surface_polygon
!surface triangle polygon defn
PetscInt :: e_connect_no,connect(3),position !element neighbour no, connection to 3 vertices
!position = 0/1/-1 = center/top/bottom
PetscInt, allocatable :: ngh_surfaces(:) !element ngh_surfaces
real(8) :: pts(9),nxyz(3),center(3),area,le_dist
!le_dist = distance from leading edige
end type surface_polygon
!3 bodies
!type(surface_polygon), allocatable :: surface1(:),surface2(:),surface3(:),surface1_ini(:),surface2_ini(:),surface3_ini(:) !,surface1_multi(:,:),surface2_multi(:,:)
type(surface_polygon), allocatable :: surface1(:) !,surface1_ini(:)
type(surface_polygon), allocatable :: surface2(:) !,surface2_ini(:)
type(surface_polygon), allocatable :: surface3(:) !,surface3_ini(:)
type(surface_polygon), allocatable :: n_surface(:,:) !,surface3_ini(:)
real(8), allocatable :: surface1_pts(:,:),surface2_pts(:,:),surface3_pts(:,:),n_surface_pts(:,:,:)
type vertex_pts
!vertex pts defn
PetscInt :: s_connect_no,v_connect_no,v_local_conn_no,conn_vertex(15),conn_surface(15) !s/v_connect_no = surfaces/ngh vertices connect no, connection to max 15 other vertices/surfaces
!v_local_conn_no = only local connected vertices (excluding own vertex no)
PetscInt, allocatable :: ngh_vertics(:) !neighbouring vertices
real(8) :: xyz(3),pres,vel(3),center_y !,thickness !,ang_vel_flex
end type vertex_pts
!3 bodies
type(vertex_pts), allocatable :: vertex1(:),vertex1_ini(:),vertex1_multi(:,:),vertex1_old(:),vertex1_old2(:)
type(vertex_pts), allocatable :: vertex2(:),vertex2_ini(:),vertex2_multi(:,:),vertex2_old(:),vertex2_old2(:)
type(vertex_pts), allocatable :: vertex3(:),vertex3_ini(:),vertex3_multi(:,:),vertex3_old(:),vertex3_old2(:)
type(vertex_pts), allocatable :: n_vertex(:,:),n_vertex_ini(:,:),n_vertex_multi(:,:,:),n_vertex_old(:,:),n_vertex_old2(:,:)
type IIB_cell
!PetscInt :: ijk(3),int_template_ijk(9),stencil_no,v_no,e_no,s_no,v_no2,s_no2,near_2_body !ij index, lower left ij bilinear interpolation template for external pt,
!req v_no,e_no,s_no at both IIB_cells and
!stencil_no= no. of values to be used in template, normal, trilinear / 2 boundary pts / 1 = 3 / 2 / only bilinear + 2 boundary pts
PetscInt :: ijk(3),int_template_ijk(75),stencil_no,v_no,e_no,s_no,v_no2,s_no2,near_2_body
real(8) :: xyz(3)
!real(8) :: mirror_xyz(3) !mirror coordinates from IIB along boundary
real(8) :: IB_xyz(3),IB_xyz_old(3) !IB coordinates - perpendicular intersection to boundary from IIB
real(8) :: IB_xyz2(3),IB_xyz3(3),IB_xyz4(3)
!real(8) :: int_template_coef(4) !bilinear interpolation coefficient template for external pt, arranged as 3 4
!use int_template_coef as vol
!real(8) :: int_template_coef(4) !using IVD, max is 26 pts
real(8) :: int_template_coef(26) !using IVD, max is 26 pts
real(8) :: vel_b,vel_b2,vel_b3,vel_b4 !can be u or v ! 1 2
!vel_b2/3/4 indicates velocity at 2nd,3rd and 4th boundary pts. Boundary pts can be only own or other body.
!real(8) :: int_template_coef_matrix(4,4),ibm_force
real(8) :: int_template_coef_matrix(8,8),ibm_force
end type IIB_cell
type tmp_IIB_cell
PetscInt :: ijk(3)
end type tmp_IIB_cell
!3 bodies
type(IIB_cell), allocatable :: IIB_cell_u1(:),IIB_cell_v1(:),IIB_cell_w1(:),IIB_cell_p1(:),pIB_cell_u(:),pIB_cell_v(:),pIB_cell_w(:),pIB_cell_p(:)
type(IIB_cell), allocatable :: IIB_cell_u2(:),IIB_cell_v2(:),IIB_cell_w2(:),IIB_cell_p2(:)
type(IIB_cell), allocatable :: IIB_cell_u3(:),IIB_cell_v3(:),IIB_cell_w3(:),IIB_cell_p3(:)
type(IIB_cell), allocatable :: n_IIB_cell_u(:,:),n_IIB_cell_v(:,:),n_IIB_cell_w(:,:),n_IIB_cell_p(:,:)
!PetscInt, allocatable :: IIB_cell_p1_s_no(:),IIB_cell_p2_s_no(:),IIB_cell_p3_s_no(:),IIB_cell_p1_v_no(:),IIB_cell_p2_v_no(:),IIB_cell_p3_v_no(:)
!specifically for pressure IIB cells
!type IIB_global_cell
! PetscInt :: ijk(3),int_template_ijk(9) !ij index, lower left ij bilinear interpolation template for external pt,
!stencil_no= no. of values to be used in template, normal, trilinear / 2 boundary pts / 1 = 3 / 2 / only bilinear + 2 boundary pts
!with new even global distribution
! real(8) :: xyz(3)
!real(8) :: mirror_xyz(3) !mirror coordinates from IIB along boundary
! real(8) :: IB_xyz(3)
! real(8) :: int_template_coef(4) !bilinear interpolation coefficient template for external pt, arranged as 3 4
!use int_template_coef as vol
! real(8) :: vel_b
!vel_b2/3/4 indicates velocity at 2nd,3rd and 4th boundary pts. Boundary pts can be only own or other body.
! real(8) :: int_template_coef_matrix(4,4)
!end type IIB_global_cell
type(IIB_cell), allocatable :: IIB_global_cell_u1(:),IIB_global_cell_v1(:),IIB_global_cell_w1(:),IIB_global_cell_p1(:)
type(IIB_cell), allocatable :: IIB_global_cell_u2(:),IIB_global_cell_v2(:),IIB_global_cell_w2(:),IIB_global_cell_p2(:)
type(IIB_cell), allocatable :: IIB_global_cell_u3(:),IIB_global_cell_v3(:),IIB_global_cell_w3(:),IIB_global_cell_p3(:)
type(IIB_cell), allocatable :: n_IIB_global_cell_u(:,:),n_IIB_global_cell_v(:,:),n_IIB_global_cell_w(:,:),n_IIB_global_cell_p(:,:)
!type IIB_equal_cell
! PetscInt :: ijk(3),int_template_ijk(9) !ij index, lower left ij bilinear interpolation template for external pt,
!stencil_no= no. of values to be used in template, normal, trilinear / 2 boundary pts / 1 = 3 / 2 / only bilinear + 2 boundary pts
!with new even equal distribution
! real(8) :: xyz(3)
!real(8) :: mirror_xyz(3) !mirror coordinates from IIB along boundary
! real(8) :: IB_xyz(3)
! real(8) :: int_template_coef(4) !bilinear interpolation coefficient template for external pt, arranged as 3 4
!use int_template_coef as vol
! real(8) :: vel_b
!vel_b2/3/4 indicates velocity at 2nd,3rd and 4th boundary pts. Boundary pts can be only own or other body.
! real(8) :: int_template_coef_matrix(4,4)
!end type IIB_equal_cell
type(IIB_cell), allocatable :: IIB_equal_cell_u1(:),IIB_equal_cell_v1(:),IIB_equal_cell_w1(:),IIB_equal_cell_p1(:)
type(IIB_cell), allocatable :: IIB_equal_cell_u2(:),IIB_equal_cell_v2(:),IIB_equal_cell_w2(:),IIB_equal_cell_p2(:)
type(IIB_cell), allocatable :: IIB_equal_cell_u3(:),IIB_equal_cell_v3(:),IIB_equal_cell_w3(:),IIB_equal_cell_p3(:)
type(IIB_cell), allocatable :: n_IIB_equal_cell_u(:,:),n_IIB_equal_cell_v(:,:),n_IIB_equal_cell_w(:,:),n_IIB_equal_cell_p(:,:)
type I_cell !internal cells
PetscInt :: ijk(3),s_no,e_no,v_no !ij index
real(8) :: vel,ibm_force
end type I_cell
type tmp_I_cell
PetscInt :: ijk(3)
end type tmp_I_cell
type(I_cell), allocatable :: I_cell_u1(:),I_cell_v1(:),I_cell_w1(:),I_cell_p1(:)
type(I_cell), allocatable :: I_cell_u2(:),I_cell_v2(:),I_cell_w2(:),I_cell_p2(:)
type(I_cell), allocatable :: I_cell_u3(:),I_cell_v3(:),I_cell_w3(:),I_cell_p3(:)
type(I_cell), allocatable :: n_I_cell_u(:,:),n_I_cell_v(:,:),n_I_cell_w(:,:),n_I_cell_p(:,:)
!type I_global_cell !internal cells
! PetscInt :: ijk(3) !ij index
! real(8) :: vel
!end type I_global_cell
type(I_cell), allocatable :: I_global_cell_u1(:),I_global_cell_v1(:),I_global_cell_w1(:),I_global_cell_p1(:)
type(I_cell), allocatable :: I_global_cell_u2(:),I_global_cell_v2(:),I_global_cell_w2(:),I_global_cell_p2(:)
type(I_cell), allocatable :: I_global_cell_u3(:),I_global_cell_v3(:),I_global_cell_w3(:),I_global_cell_p3(:)
type(I_cell), allocatable :: n_I_global_cell_u(:,:),n_I_global_cell_v(:,:),n_I_global_cell_w(:,:),n_I_global_cell_p(:,:)
!type I_equal_cell !internal cells
! PetscInt :: ijk(3) !ij index
! real(8) :: vel
!end type I_equal_cell
type(I_cell), allocatable :: I_equal_cell_u1(:),I_equal_cell_v1(:),I_equal_cell_w1(:),I_equal_cell_p1(:)
type(I_cell), allocatable :: I_equal_cell_u2(:),I_equal_cell_v2(:),I_equal_cell_w2(:),I_equal_cell_p2(:)
type(I_cell), allocatable :: I_equal_cell_u3(:),I_equal_cell_v3(:),I_equal_cell_w3(:),I_equal_cell_p3(:)
type(I_cell), allocatable :: n_I_equal_cell_u(:,:),n_I_equal_cell_v(:,:),n_I_equal_cell_w(:,:),n_I_equal_cell_p(:,:)
!kdtree
type(kdtree2), pointer :: tree_surface1,tree_surface2,tree_surface3
type(kdtree2_result) :: nearest_pt_kdtree(1)
!PETSc variables
Vec xx,b_rhs,b_rhs_semi_xyz,xx_semi_xyz
Mat A_mat,A_semi_xyz
MatNullSpace nullsp
KSP ksp,ksp_semi_xyz
DM da_u,da_v,da_w,da_p,da_uvw_ast,da_q_semi_vect_xyz,da_q_p !,da_uvw_old !DM stuffs
DM da_cu_types !,da_cv_types,da_cw_types,da_cp_types
Vec u_local,u_global,v_local,v_global,w_local,w_global,p_local,p_global,uvw_ast_local,uvw_ast_global,q_semi_vect_xyz_local,q_semi_vect_xyz_global,q_p_local,q_p_global !,uvw_old_local,uvw_old_global
Vec cu_types_local,cu_types_global !,cv_types_local,cv_types_global,cw_types_local,cw_types_global,cp_types_local,cp_types_global
PetscInt :: ione,iseven
PetscScalar,pointer :: u_array(:,:,:),v_array(:,:,:),w_array(:,:,:),p_array(:,:,:),uvw_ast_array(:,:,:,:),q_semi_vect_xyz_array(:,:,:,:),q_p_array(:,:,:) !,uvw_old_array(:,:,:,:)
PetscScalar,pointer :: cu_types_array(:,:,:),cv_types_array(:,:,:),cw_types_array(:,:,:),cp_types_array(:,:,:)
PetscScalar :: value_insert(7)
MatStencil :: row(4,1),col(4,7)
PC pc,pc_semi_xyz
PCType ptype
KSPType ksptype
PetscBool flgg
KSPConvergedReason reason
PetscOffset i_vec
PetscViewer viewer
PetscLogEvent Poisson_log
!PetscLogStage stage1
!PetscScalar,pointer :: b_rhs_semi_local_array(:,:,:,:)
!PetscScalar,pointer :: b_rhs_semi_local_array(:,:,:) !dof = 1
!real(8), allocatable :: u_array(:,:,:),v_array(:,:,:),w_array(:,:,:),p_array(:,:,:) !,uvw_old_array(:,:,:,:)
!parallel
PetscMPIInt :: myid,num_procs,num_procs_old
PetscInt :: ksta,kend,ksta_ext,kend_ext,ksta2,kend2,procs_x,procs_y,procs_z
PetscInt :: jsta,jend,jsta_ext,jend_ext,jsta2,jend2 !,ista,iend,ista_ext,iend_ext,ista2,iend2
!3 bodies
!for overall domain encompassing all bodies
!not applicable for > 3 bodies since it will not be using overall domain. each body will only use its own domain, since no overlap
PetscInt :: ista_IIB,iend_IIB,jsta_IIB,jend_IIB,ksta_IIB,kend_IIB !i/j/k_sta/end for IIB only
PetscInt :: IIB_ksta1,IIB_kend1,IIB_ksta2,IIB_kend2,IIB_ksta3,IIB_kend3
PetscInt :: IIB_jsta1,IIB_jend1,IIB_jsta2,IIB_jend2,IIB_jsta3,IIB_jend3
PetscInt :: IIB_ista1,IIB_iend1,IIB_ista2,IIB_iend2,IIB_ista3,IIB_iend3
PetscInt, allocatable :: n_IIB_ista(:),n_IIB_iend(:),n_IIB_jsta(:),n_IIB_jend(:),n_IIB_ksta(:),n_IIB_kend(:)
PetscInt :: ijk_sta_p,ijk_end_p,ijk_sta_m,ijk_end_m,ijk_sta_mx,ijk_end_mx,ijk_sta_my,ijk_end_my,ijk_sta_mz,ijk_end_mz
PetscInt :: ijk_sta,ijk_end,ijk_dm_sta,ijk_dm_end,ijk_xyz_sta,ijk_xyz_end !>includes extension
PetscInt :: num_procs_xyz(3),start_ijk(3),width_ijk(3),width_ijk_ghost(3),end_ijk(3),start_ijk_ghost(3),end_ijk_ghost(3),stencil_width
PetscInt :: num_procs_xyz_IIB(3),width_ijk_IIB(3),width_ijk_ghost_IIB(3),stencil_width_IIB
PetscInt :: start_ijk_IIB(3),end_ijk_IIB(3),start_ijk_int_IIB(3),end_ijk_int_IIB(3),start_ijk_ghost_IIB(3),end_ijk_ghost_IIB(3)
!IIB_i/j/ksta,IIB_i/j/kend = standard IIB boundary
!start_ijk_IIB,end_ijk_IIB = start/end index for IIB domain of all bodies
!start_ijk_int_IIB,end_ijk_int_IIB = intermediate IIB boundary +/- 1
!start_ijk_ghost_IIB,end_ijk_ghost_IIB = standard IIB boundary +/- 2
PetscInt, allocatable :: num_procs_ksta_kend(:,:),num_procs_jsta_jend(:,:),num_procs_size(:),temp_procs(:)
!PetscInt(8), allocatable :: num_procs_size_start_index(:)
PetscInt, allocatable :: num_procs_size_start_index(:)
!>num_procs_size_start_index represents (starting index - 1) of each procs
character(2) :: disk_loc
character(5) :: procs,file_write_no_char,body_no_char
contains
subroutine allo_var
!allocate memory for variables
#include "petsc/finclude/petsc.h"
PetscInt :: ijk,status(3)
PetscInt :: lx(1),ly(1) !,lz(0:num_procs-1)
!>Grid partitioning, using PETSc DM
!grid variables
!ksta=myid*(size_z/num_procs)+1; kend=(myid+1)*(size_z/num_procs)
one = 1; three = 3
sta_ijk_uniform(1) = 1; end_ijk_uniform(1) = size_x
allocate (num_procs_jsta_jend(2,0:num_procs-1), STAT = status(1))
if (status(1)/=0) stop "1 Cannot allocate memory"
allocate (num_procs_ksta_kend(2,0:num_procs-1), STAT = status(1))
if (status(1)/=0) stop "2 Cannot allocate memory"
allocate (num_procs_size(0:num_procs-1), STAT = status(1))
if (status(1)/=0) stop "3 Cannot allocate memory"
allocate (num_procs_size_start_index(0:num_procs-1), STAT = status(1))
if (status(1)/=0) stop "4 Cannot allocate memory"
allocate (temp_procs(3*num_procs), STAT = status(1))
if (status(1)/=0) stop "5 Cannot allocate memory"
!DM stuffs
!procs_x = 1; procs_y = 2; procs_z = 1
stencil_width = 2; stencil_width_IIB = 1
lx(1) = size_x ; ly(1) = size_y
!do ijk = 0,num_procs-1
! call mpi_para_range(1,size_z,num_procs,ijk,num_procs_ksta_kend(1,ijk),num_procs_ksta_kend(2,ijk))
! lz(ijk) = num_procs_ksta_kend(2,ijk) - num_procs_ksta_kend(1,ijk) + 1 !DM stuffs
!end do
call DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&
size_z,one,PETSC_DECIDE,PETSC_DECIDE,one,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_u,ierr)
call DMSetFromOptions(da_u,ierr)
call DMSetUp(da_u,ierr)
if (inout_test <= 1) then
call DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&
size_z,one,PETSC_DECIDE,PETSC_DECIDE,one,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_v,ierr)
call DMSetFromOptions(da_v,ierr)
call DMSetUp(da_v,ierr)
call DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&
size_z,one,PETSC_DECIDE,PETSC_DECIDE,one,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_w,ierr)
call DMSetFromOptions(da_w,ierr)
call DMSetUp(da_w,ierr)
call DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&
size_z,one,PETSC_DECIDE,PETSC_DECIDE,one,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_p,ierr)
call DMSetFromOptions(da_p,ierr)
call DMSetUp(da_p,ierr)
call DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&
size_z,one,PETSC_DECIDE,PETSC_DECIDE,three,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_uvw_ast,ierr)
call DMSetFromOptions(da_uvw_ast,ierr)
call DMSetUp(da_uvw_ast,ierr)
call DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&
size_z,one,PETSC_DECIDE,PETSC_DECIDE,three,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_q_semi_vect_xyz,ierr)
call DMSetFromOptions(da_q_semi_vect_xyz,ierr)
call DMSetUp(da_q_semi_vect_xyz,ierr)
call DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&
size_z,one,PETSC_DECIDE,PETSC_DECIDE,one,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_q_p,ierr)
call DMSetFromOptions(da_q_p,ierr)
call DMSetUp(da_q_p,ierr)
end if
!u_global = tVec(1)
call DMCreateGlobalVector(da_u,u_global,ierr)
call DMCreateLocalVector(da_u,u_local,ierr)
if (inout_test <= 1) then
call DMCreateGlobalVector(da_v,v_global,ierr)
call DMCreateLocalVector(da_v,v_local,ierr)
call DMCreateGlobalVector(da_w,w_global,ierr)
call DMCreateLocalVector(da_w,w_local,ierr)
call DMCreateGlobalVector(da_p,p_global,ierr)
call DMCreateLocalVector(da_p,p_local,ierr)
call DMCreateGlobalVector(da_uvw_ast,uvw_ast_global,ierr)
call DMCreateLocalVector(da_uvw_ast,uvw_ast_local,ierr)
call DMCreateGlobalVector(da_q_semi_vect_xyz,q_semi_vect_xyz_global,ierr)
call DMCreateLocalVector(da_q_semi_vect_xyz,q_semi_vect_xyz_local,ierr)
call DMCreateGlobalVector(da_q_p,q_p_global,ierr)
call DMCreateLocalVector(da_q_p,q_p_local,ierr)
!call DMCreateGlobalVector(da_uvw_old,uvw_old_global,ierr)
!call DMCreateLocalVector(da_uvw_old,uvw_old_local,ierr)
end if
call DMDAGetInfo(da_u,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,num_procs_xyz(1),num_procs_xyz(2),&
num_procs_xyz(3),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
!call DMDAGetInfo(da_v,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,num_procs_xyz(1),num_procs_xyz(2),& !same for all
!num_procs_xyz(3),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
!call DMDAGetInfo(da_w,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,num_procs_xyz(1),num_procs_xyz(2),&
!num_procs_xyz(3),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
!call DMDAGetInfo(da_p,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,num_procs_xyz(1),num_procs_xyz(2),&
!num_procs_xyz(3),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
!call DMDAGetInfo(da_uvw_ast,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,num_procs_xyz(1),num_procs_xyz(2),&
!num_procs_xyz(3),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
call DMDAGetCorners(da_u,start_ijk(1),start_ijk(2),start_ijk(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)
!call DMDAGetCorners(da_v,start_ijk(1),start_ijk(2),start_ijk(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)
!call DMDAGetCorners(da_w,start_ijk(1),start_ijk(2),start_ijk(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)
!call DMDAGetCorners(da_p,start_ijk(1),start_ijk(2),start_ijk(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)
if (inout_test <= 1) call DMDAGetCorners(da_uvw_ast,start_ijk(1),start_ijk(2),start_ijk(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)
call DMDAGetGhostCorners(da_u,start_ijk_ghost(1),start_ijk_ghost(2),start_ijk_ghost(3),width_ijk_ghost(1),width_ijk_ghost(2),width_ijk_ghost(3),ierr)
!dec$ if defined(__linux)
!call DMView(da_u,PETSC_VIEWER_DRAW_WORLD,ierr)
!dec$ endif
!call DMDAGetGhostCorners(da_uvw_ast,start_ijk_ghost(1),start_ijk_ghost(2),start_ijk_ghost(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)
!if (myid == 0) print *, "start_ijk_ghost,width_ijk_ghost x", start_ijk_ghost,width_ijk_ghost
!call DMDAGetGhostCorners(da_v,start_ijk_ghost(1),start_ijk_ghost(2),start_ijk_ghost(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)
!if (myid == 0) print *, "start_ijk_ghost,width_ijk y", start_ijk_ghost,width_ijk
!call DMDAGetGhostCorners(da_w,start_ijk_ghost(1),start_ijk_ghost(2),start_ijk_ghost(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)
!if (myid == 0) print *, "start_ijk_ghost,width_ijk z", start_ijk_ghost,width_ijk
!call DMDAGetGhostCorners(da_p,start_ijk_ghost(1),start_ijk_ghost(2),start_ijk_ghost(3),width_ijk(1),width_ijk(2),width_ijk(3),ierr)
!if (myid == 0) print *, "start_ijk_ghost,width_ijk p", start_ijk_ghost,width_ijk
!Here we shift the starting indices up by one so that we can easily
!use the Fortran convention of 1-based indices (rather 0-based indices)
!start/end_ijk_ghost to shift during intent(in)
start_ijk = start_ijk + 1
end_ijk = start_ijk + width_ijk - 1
start_ijk_ghost = start_ijk_ghost + 1
end_ijk_ghost = start_ijk_ghost + width_ijk_ghost - 1
!do ijk = 0,num_procs - 1 !to verify correct
! call MPI_Barrier(MPI_COMM_WORLD,ierr)
!if (ijk == myid) print *, "myid,num_procs_xyz,start_ijk,width_ijk,end_ijk", myid,num_procs_xyz,start_ijk,width_ijk,end_ijk
!if (ijk == myid) print *, "myid,num_procs_xyz,start_ijk_ghost,width_ijk,end_ijk_ghost", myid,num_procs_xyz,start_ijk_ghost,width_ijk,end_ijk_ghost
!if (ijk == myid) print *, "myid,ksta,kend,ksta_ext,kend_ext",myid,ksta,kend,ksta_ext,kend_ext
!end do
call MPI_ALLGATHER(start_ijk,3,MPIU_INTEGER,temp_procs,3,MPIU_INTEGER,MPI_COMM_WORLD,ierr)
do ijk = 0,num_procs-1
num_procs_jsta_jend(1,ijk) = temp_procs((ijk + 1)*3 - 1)
num_procs_ksta_kend(1,ijk) = temp_procs((ijk + 1)*3)
end do
call MPI_ALLGATHER(end_ijk,3,MPIU_INTEGER,temp_procs,3,MPIU_INTEGER,MPI_COMM_WORLD,ierr)
do ijk = 0,num_procs-1
num_procs_jsta_jend(2,ijk) = temp_procs((ijk + 1)*3 - 1)
num_procs_ksta_kend(2,ijk) = temp_procs((ijk + 1)*3)
end do
do ijk = 0,num_procs-1
!call mpi_para_range(1,size_y,num_procs,ijk,num_procs_jsta_jend(1,ijk),num_procs_jsta_jend(2,ijk))
!call mpi_para_range(1,size_z,num_procs,ijk,num_procs_ksta_kend(1,ijk),num_procs_ksta_kend(2,ijk))
!num_procs_size(ijk) = (num_procs_ksta_kend(2,ijk) - num_procs_ksta_kend(1,ijk) + 1) * size_x * size_y
num_procs_size(ijk) = (num_procs_ksta_kend(2,ijk) - num_procs_ksta_kend(1,ijk) + 1) * size_x * (num_procs_jsta_jend(2,ijk) - num_procs_jsta_jend(1,ijk) + 1)
!num_procs_size_start_index(ijk) = (num_procs_ksta_kend(1,ijk)-1) * size_x * size_y
!lz(ijk) = num_procs_ksta_kend(2,ijk) - num_procs_ksta_kend(1,ijk) + 1 !DM stuffs
end do
num_procs_size_start_index = 1
do ijk = 1,num_procs-1
num_procs_size_start_index(ijk) = num_procs_size_start_index(ijk - 1) + num_procs_size(ijk - 1)
!print *, "ijk,num_procs_size_start_index(ijk)",ijk,num_procs_size_start_index(ijk)
end do
!if (myid == 0) then
! print *, "num_procs_ksta_kend,num_procs_size,num_procs_size_start_index" !reuse debug
! print *, num_procs_ksta_kend,num_procs_size,num_procs_size_start_index
!end if
jsta = start_ijk(2); jend = end_ijk(2)
ksta = start_ijk(3); kend = end_ijk(3)
jsta_ext = start_ijk_ghost(2); jend_ext = end_ijk_ghost(2)
ksta_ext = start_ijk_ghost(3); kend_ext = end_ijk_ghost(3)
!call mpi_para_range(1,size_z,num_procs,myid,ksta,kend)
!ijk_sta = size_x*size_y*(ksta-1); ijk_end = size_x*size_y*kend ! 0 based indices for ijk_sta only
!ijk_xyz_sta = 3*size_x*size_y*(ksta-1); ijk_xyz_end = 3*size_x*size_y*kend
!ijk_sta = size_x*(jsta-1)*(ksta-1); ijk_end = size_x*jend*kend ! 0 based indices for ijk_sta only
!ijk_xyz_sta = 3*size_x*(jsta-1)*(ksta-1); ijk_xyz_end = 3*size_x*jend*kend
ijk_sta = num_procs_size_start_index(myid); ijk_end = num_procs_size_start_index(myid) + num_procs_size(myid) - 1
ijk_xyz_sta = 3*ijk_sta - 2; ijk_xyz_end = 3*ijk_end
!do ijk = 0,num_procs - 1 !to verify correct
!call MPI_Barrier(MPI_COMM_WORLD,ierr)
!if (ijk == myid) print *, "myid,jsta,jend,ksta,kend,ijk_sta,ijk_end",myid,jsta,jend,ksta,kend,ijk_sta,ijk_end
!if (ijk == myid) print *, "myid,num_procs_xyz,start_ijk_ghost,width_ijk,end_ijk_ghost", myid,num_procs_xyz,start_ijk_ghost,width_ijk,end_ijk_ghost
!if (ijk == myid) print *, "myid,ksta,kend,ksta_ext,kend_ext",myid,ksta,kend,ksta_ext,kend_ext
!end do
if (kend - ksta + 1 <= 2) then
print *, "z grid divid too small!"
print *, "myid,each procs z size", myid, kend - ksta + 1
end if
if (jend - jsta + 1 <= 2) then
print *, "y grid divid too small!"
print *, "myid,each procs y size", myid, jend - jsta + 1
end if
if (ksta == 1) then
ksta2 = ksta+1
else
ksta2 = ksta
end if
if (kend == size_z) then
kend2 = kend - 1
else
kend2 = kend
end if
if (jsta == 1) then
jsta2 = jsta+1
else
jsta2 = jsta
end if
if (jend == size_y) then
jend2 = jend - 1
else
jend2 = jend
end if
if (num_procs < 50) print *, "myid,jsta,jend,ksta,kend",myid,jsta,jend,ksta,kend !quick debug
!print *, "myid,jsta2,jend2,start_ijk_ghost,end_ijk_ghost",jsta2,jend2,start_ijk_ghost,end_ijk_ghost
!print *, "myid,ksta,kend,ksta_ext,kend_ext,ksta2,kend2",myid,ksta,kend,ksta_ext,kend_ext,ksta2,kend2 !reuse debug
!allocate (x(0:size_x), STAT = status(1))
!allocate (y(0:size_y), STAT = status(2))
!allocate (z(0:size_z), STAT = status(3))
!if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "Cannot allocate memory"
allocate (xu(0:size_x), STAT = status(1))
allocate (yu(0:size_y), STAT = status(2))
allocate (zu(0:size_z), STAT = status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "6 Cannot allocate memory"
allocate (xv(0:size_x), STAT = status(1))
allocate (yv(0:size_y), STAT = status(2))
allocate (zv(0:size_z), STAT = status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "7 Cannot allocate memory"
allocate (xw(0:size_x), STAT = status(1))
allocate (yw(0:size_y), STAT = status(2))
allocate (zw(0:size_z), STAT = status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "8 Cannot allocate memory"
allocate (c_cx(size_x), STAT = status(1))
allocate (c_cy(size_y), STAT = status(2))
allocate (c_cz(size_z), STAT = status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "9 Cannot allocate memory"
allocate (cu_cx(size_x), STAT = status(1))
allocate (cv_cy(size_y), STAT = status(2))
allocate (cw_cz(size_z), STAT = status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "Cannot allocate memory"
!solver variables
Total_ijk=size_x*size_y*size_z
!body variables
!ibm variables
if (inout_test <= 1) then
!flow variables
!allocate (p(size_x,size_y,ksta_ext:kend_ext), STAT=status(1))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
!allocate (u(size_x,size_y,ksta_ext:kend_ext), STAT=status(1))
!allocate (v(size_x,size_y,ksta_ext:kend_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,ksta_ext:kend_ext), STAT=status(1));allocate (v_old(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "10 Cannot allocate memory"
!allocate (w(size_x,size_y,ksta_ext:kend_ext), STAT=status(1));
allocate (w_old(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "11 Cannot allocate memory"
if (q_criterion_plot >= 1) then
!if (.not. allocated(Q_Soria)) allocate (Q_Soria(size_x,size_y,ksta_ext:kend_ext), STAT=status(1))
!if (.not. allocated(q_crit)) allocate (q_crit(size_x,size_y,ksta_ext:kend_ext), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
if (.not. allocated(Q_criterion)) allocate (Q_criterion(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(1))
if (status(1)/=0) stop "12 Cannot allocate memory"
!if (.not. allocated(Q_Jeong)) allocate (Q_Jeong(size_x,size_y,ksta_ext:kend_ext), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
if (q_criterion_plot == 2) then
allocate (Q_criterion_avg(size_x,jsta:jend,ksta:kend,total_snapshots), STAT=status(1));allocate (Q_criterion_tmp(size_x,jsta:jend,ksta:kend), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "13 Cannot allocate memory"
end if
end if
if (avg_vorticity == 1) then
allocate (vorticity_avg(size_x,jsta:jend,ksta:kend,total_snapshots), STAT=status(1));allocate (vorticity_tmp(size_x,jsta:jend,ksta:kend), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "14 Cannot allocate memory"
end if
allocate (diff_pres_x(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(1))
allocate (diff_pres_y(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
allocate (diff_pres_z(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "15 Cannot allocate memory"
!allocate (u_ast(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(1));allocate (v_ast(size_x,jsta_ext:jend_ext,ksta_ext:kend_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,ksta_ext:kend_ext), STAT=status(1));allocate (v_ast2(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "16 Cannot allocate memory"
!allocate (w_ast(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(1))
!if (status(1)/=0) stop "Cannot allocate memory"
allocate (phi(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(1));allocate (w_ast2(size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "17 Cannot allocate memory"
!ijk_sta_p = size_x*size_y*(ksta-1); ijk_end_p = size_x*size_y*kend
ijk_sta_p = ijk_sta; ijk_end_p = ijk_end
!allocate (IIB_index_u(size_x,size_y,ksta_ext:kend_ext), STAT=status(1));allocate (IIB_index_v(size_x,size_y,ksta_ext:kend_ext), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
!allocate (IIB_index_w(size_x,size_y,ksta_ext:kend_ext), STAT=status(1));allocate (IIB_index_p(size_x,size_y,ksta_ext:kend_ext), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
end if
if (no_body <= 3) then
!body1
allocate (IIB_I_cell_no_uvw_global1(0:6*num_procs-1), STAT=status(1));allocate (IIB_I_cell_no_p_global1(0:2*num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "18 Cannot allocate memory"
allocate (IIB_equal_cell_no_global_u1(0:num_procs-1), STAT=status(1));allocate (IIB_equal_cell_no_global_v1(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "19 Cannot allocate memory"
allocate (IIB_equal_cell_no_global_w1(0:num_procs-1), STAT=status(1));allocate (I_equal_cell_no_global_u1(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "20 Cannot allocate memory"
allocate (I_equal_cell_no_global_v1(0:num_procs-1), STAT=status(1));allocate (I_equal_cell_no_global_w1(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "21 Cannot allocate memory"
allocate (IIB_equal_cell_no_global_p1(0:num_procs-1), STAT=status(1));allocate (I_equal_cell_no_global_p1(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "22 Cannot allocate memory"
!body2
if (no_body >= 2) then
allocate (IIB_I_cell_no_uvw_global2(0:6*num_procs-1), STAT=status(1));allocate (IIB_I_cell_no_p_global2(0:2*num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "23 Cannot allocate memory"
allocate (IIB_equal_cell_no_global_u2(0:num_procs-1), STAT=status(1));allocate (IIB_equal_cell_no_global_v2(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "24 Cannot allocate memory"
allocate (IIB_equal_cell_no_global_w2(0:num_procs-1), STAT=status(1));allocate (I_equal_cell_no_global_u2(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "25 Cannot allocate memory"
allocate (I_equal_cell_no_global_v2(0:num_procs-1), STAT=status(1));allocate (I_equal_cell_no_global_w2(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "26 Cannot allocate memory"
allocate (IIB_equal_cell_no_global_p2(0:num_procs-1), STAT=status(1));allocate (I_equal_cell_no_global_p2(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "27 Cannot allocate memory"
if (no_body == 3) then
allocate (IIB_I_cell_no_uvw_global3(0:6*num_procs-1), STAT=status(1));allocate (IIB_I_cell_no_p_global3(0:2*num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "28 Cannot allocate memory"
allocate (IIB_equal_cell_no_global_u3(0:num_procs-1), STAT=status(1));allocate (IIB_equal_cell_no_global_v3(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "29 Cannot allocate memory"
allocate (IIB_equal_cell_no_global_w3(0:num_procs-1), STAT=status(1));allocate (I_equal_cell_no_global_u3(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "30 Cannot allocate memory"
allocate (I_equal_cell_no_global_v3(0:num_procs-1), STAT=status(1));allocate (I_equal_cell_no_global_w3(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "31 Cannot allocate memory"
allocate (IIB_equal_cell_no_global_p3(0:num_procs-1), STAT=status(1));allocate (I_equal_cell_no_global_p3(0:num_procs-1), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "32 Cannot allocate memory"
end if
end if
end if
if (inout_test <= 1) then
!PETSc variables
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 MatCreateAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,size_x*size_y*size_z,size_x*size_y*size_z,7,PETSC_NULL_INTEGER,7,PETSC_NULL_INTEGER,A_mat,ierr)
!call MatCreateAIJ(MPI_COMM_WORLD,ijk_end_p - ijk_sta_p + 1,ijk_end_p - ijk_sta_p + 1,PETSC_DECIDE,PETSC_DECIDE,7,PETSC_NULL_INTEGER,7,PETSC_NULL_INTEGER,A_mat,ierr)
call DMSetMatType(da_p,MATAIJ,ierr)
call DMCreateMatrix(da_p,A_mat,ierr)
call MatSetFromOptions(A_mat,ierr)
call MatGetOwnershipRange(A_mat,ijk_sta_p,ijk_end_p,ierr)
ijk_sta_p = ijk_sta_p + 1 !start from index 1
!print *, "myid,ijk_sta_p,ijk_end_p",myid,ijk_sta_p,ijk_end_p
call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
!call VecCreateMPI(MPI_COMM_WORLD,ijk_end_p - ijk_sta + 1,PETSC_DECIDE,b_rhs,ierr)
!call VecCreateMPI(MPI_COMM_WORLD,PETSC_DECIDE,size_x*size_y*size_z,b_rhs,ierr)
call DMCreateGlobalVector(da_p,b_rhs,ierr)
call VecDuplicate(b_rhs,xx,ierr)
end if
!iterative solver variables
allocate (big_A(ijk_sta_p:ijk_end_p,7), STAT=status(1))
!allocate (q_p(ijk_sta_p:ijk_end_p), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "33 Cannot allocate memory"
allocate (int_A(ijk_sta_p:ijk_end_p,7), STAT=status(1))
allocate (qk(ijk_sta_p:ijk_end_p), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "34 Cannot allocate memory"
!allocate (fc1_prev_r(size_x), STAT=status(1));allocate (fc1_prev_c(ksta_ext:kend_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(ksta_ext:kend_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(ksta_ext:kend_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(ksta_ext:kend_ext), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
allocate (fc_prev_n(3,size_x), STAT=status(1));allocate (fc_prev_n_old(3,size_x), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "35 Cannot allocate memory"
allocate (fd_prev_n(3,size_x), STAT=status(1)) ;allocate (fc_prev_f(3,size_x,size_y), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "36 Cannot allocate memory"
allocate (fc_prev_f_old(3,size_x,size_y), STAT=status(1));allocate (fd_prev_f(3,size_x,size_y), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
!allocate (fd2_prev_r_old(size_x), STAT=status(1));allocate (fd2_prev_c_old(ksta_ext:kend_ext), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
if (spatial_scheme == 3 .or. spatial_scheme == 6) then
allocate (s_cf_u_e(3,2,size_x), STAT=status(1)); allocate (s_cf_v_e(3,2,size_x), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "37 Cannot allocate memory"
allocate (s_cf_u_n(3,2,jsta2 - 1:jend), STAT=status(1)); allocate (s_cf_v_n(3,2,jsta2 - 1:jend), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "38 Cannot allocate memory"
allocate (s_cf_w_e(3,2,size_x), STAT=status(1)); allocate (s_cf_w_n(3,2,jsta2 - 1:jend), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "39 Cannot allocate memory"
allocate (s_cf_u_f(3,2,ksta2 - 1:kend), STAT=status(1)); allocate (s_cf_v_f(3,2,ksta2 - 1:kend), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "40 Cannot allocate memory"
allocate (s_cf_w_f(3,2,ksta2 - 1:kend), STAT=status(1))
if (status(1)/=0) STOP "41 Cannot allocate memory"
end if
end if
!cell variables
allocate (cu(1:size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(1)) !inout_test_cell_uv
if (status(1)/=0) stop "42 Cannot allocate memory"
!allocate (IIB_I_cu(size_x,size_y,size_z), STAT=status(1)) !inout_test_cell_uv
!if (status(1)/=0) stop "Cannot allocate memory"
if (inout_test <= 1) then
allocate (cv(1:size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2)) !inout_test
if (status(2)/=0) stop "43 Cannot allocate memory"
allocate (cp(1:size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(1)) !inout_test
allocate (cw(1:size_x,jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2)) !inout_test
if (status(1)/=0 .or. status(2)/=0) stop "44 Cannot allocate memory"
!allocate (IIB_I_cv(size_x,size_y,size_z), STAT=status(2)) !inout_test
!if (status(2)/=0) stop "Cannot allocate memory"
!allocate (IIB_I_cp(size_x,size_y,size_z), STAT=status(1)) !inout_test
!allocate (IIB_I_cw(size_x,size_y,size_z), STAT=status(2)) !inout_test
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
end if
allocate (cu_x(size_x), STAT=status(1));allocate (cv_x(size_x), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "45 Cannot allocate memory"
allocate (cp_x(size_x), STAT=status(1));allocate (cw_x(size_x), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "46 Cannot allocate memory"
allocate (cu_y(size_y), STAT=status(1));allocate (cv_y(size_y), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "47 Cannot allocate memory"
allocate (cp_y(size_y), STAT=status(1));allocate (cw_y(size_y), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "48 Cannot allocate memory"
allocate (cu_z(size_z), STAT=status(1));allocate (cv_z(size_z), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "49 Cannot allocate memory"
allocate (cp_z(size_z), STAT=status(1));allocate (cw_z(size_z), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "50 Cannot allocate memory"
if (inout_test <= 1) then
allocate (cp_xy(size_x,jsta_ext:jend_ext), STAT=status(1));allocate (cp_yz(jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "51 Cannot allocate memory"
allocate (cp_zx(size_x,ksta_ext:kend_ext), STAT=status(1))
if (status(1)/=0) stop "52 Cannot allocate memory"
allocate (cu_xy(size_x,jsta_ext:jend_ext), STAT=status(1));allocate (cu_yz(jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "53 Cannot allocate memory"
allocate (cu_zx(size_x,ksta_ext:kend_ext), STAT=status(1))
if (status(1)/=0) stop "54 Cannot allocate memory"
allocate (cv_xy(size_x,jsta_ext:jend_ext), STAT=status(1));allocate (cv_yz(jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "55 Cannot allocate memory"
allocate (cv_zx(size_x,ksta_ext:kend_ext), STAT=status(1))
if (status(1)/=0) stop "56 Cannot allocate memory"
allocate (cw_xy(size_x,jsta_ext:jend_ext), STAT=status(1));allocate (cw_yz(jsta_ext:jend_ext,ksta_ext:kend_ext), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "57 Cannot allocate memory"
allocate (cw_zx(size_x,ksta_ext:kend_ext), STAT=status(1))
if (status(1)/=0) stop "58 Cannot allocate memory"
end if
if (body_input_type == 2 .and. no_body <= 3) then
allocate (sta_yy1(sta_ijk_uniform(1):end_ijk_uniform(1)), STAT=status(1));allocate (end_yy1(sta_ijk_uniform(1):end_ijk_uniform(1)), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "59 Cannot allocate memory"
if (no_body >= 2) then
allocate (sta_yy2(sta_ijk_uniform(1):end_ijk_uniform(1)), STAT=status(1));allocate (end_yy2(sta_ijk_uniform(1):end_ijk_uniform(1)), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "60 Cannot allocate memory"
if (no_body == 3) then
allocate (sta_yy3(sta_ijk_uniform(1):end_ijk_uniform(1)), STAT=status(1));allocate (end_yy3(sta_ijk_uniform(1):end_ijk_uniform(1)), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "61 Cannot allocate memory"
end if
end if
end if
call PetscLogEventRegister('Poisson log',0,Poisson_log,ierr)
!ibm variables
!allocate (ibm_force(size_x,size_y,ksta_ext:kend_ext,3), STAT=status(2))
!if (status(2)/=0) stop "Cannot allocate memory"
!if (ibm_choice==3) then !Ravoux - use int_template_coef as vol
! allocate (vol_percent(size_x,size_y,ksta_ext:kend_ext,3), STAT=status(2))
! if (status(2)/=0) stop "Cannot allocate memory"
!end if
!allocate (IIB_cell_u(size_x*(kend-ksta)), STAT=status(1))
!allocate (IIB_cell_v(size_x*(kend-ksta)), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
!allocate (IIB_cell_p(size_x*(kend-ksta)), STAT=status(1))
!if (status(1)/=0) stop "Cannot allocate memory"
!allocate (pIB_cell_u(size_x*(kend-ksta)), STAT=status(1))
!allocate (pIB_cell_v(size_x*(kend-ksta)), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
!allocate (pIB_cell_p(size_x*(kend-ksta)), STAT=status(1))
!if (status(1)/=0) stop "Cannot allocate memory"
!allocate (I_cell_u(size_x*(kend-ksta)), STAT=status(1))
!allocate (I_cell_v(size_x*(kend-ksta)), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
!allocate (I_cell_p(size_x*(kend-ksta)), STAT=status(1))
!if (status(1)/=0) stop "Cannot allocate memory"
!print *,myid,num_procs,ksta,kend,ksta_ext,kend_ext,ijk_sta_m,ijk_end_m,ksta2,kend2,kend3
!if (new_start==1 .or. (new_start==0 .and. time_scheme==3)) then
end subroutine allo_var
subroutine semi_allo_var
!allocate memory for variables, from time_step>2 for semi-implicit
PetscInt :: i,NN
NN=size_x*size_y*size_z
if (inout_test <= 1) then
!call MatCreateAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,NN,NN,7,PETSC_NULL_INTEGER,7,PETSC_NULL_INTEGER,A_semi_x,ierr)
!call MatCreateAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,NN,NN,7,PETSC_NULL_INTEGER,7,PETSC_NULL_INTEGER,A_semi_y,ierr)
!call MatCreateAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,NN,NN,7,PETSC_NULL_INTEGER,7,PETSC_NULL_INTEGER,A_semi_z,ierr)
!call MatCreateAIJ(MPI_COMM_WORLD,ijk_xyz_end-ijk_xyz_sta,ijk_xyz_end-ijk_xyz_sta,PETSC_DECIDE,PETSC_DECIDE,7,PETSC_NULL_INTEGER,7,PETSC_NULL_INTEGER,A_semi_xyz,ierr)
!call MatCreateAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,3*NN,3*NN,7,PETSC_NULL_INTEGER,7,PETSC_NULL_INTEGER,A_semi_xyz,ierr)
call DMSetMatType(da_uvw_ast,MATAIJ,ierr)
call DMCreateMatrix(da_uvw_ast,A_semi_xyz,ierr)
call MatSetFromOptions(A_semi_xyz,ierr)
!call MatGetOwnershipRange(A_semi_xyz,ijk_xyz_sta_mz,ijk_xyz_end_mz,ierr)
!print *, "myid,ijk_xyz_sta_mz,ijk_xyz_end_mz,ijk_xyz_sta,ijk_xyz_end",myid,ijk_xyz_sta_mz,ijk_xyz_end_mz,ijk_xyz_sta,ijk_xyz_end
!call VecGetOwnershipRange(b_rhs_semi_global,ijk_dm_sta,ijk_dm_end,ierr)
!ijk_xyz_sta_mz = ijk_xyz_sta_mz + 1 !start from 1 index
!call VecCreateMPI(MPI_COMM_WORLD,PETSC_DECIDE,3*NN,b_rhs_semi_xyz,ierr)
call DMCreateGlobalVector(da_uvw_ast,b_rhs_semi_xyz,ierr)
call VecDuplicate(b_rhs_semi_xyz,xx_semi_xyz,ierr)
call KSPCreate(MPI_COMM_WORLD,ksp_semi_xyz,ierr)
!allocate (semi_mat_xyz(ijk_xyz_sta:ijk_xyz_end,7), STAT=status(1))
allocate (semi_mat_xyz_da(ijk_xyz_sta:ijk_xyz_end), STAT=status(1))
!allocate (int_semi_xyz(ijk_xyz_sta:ijk_xyz_end,7), STAT=status(2))
!allocate (q_semi_vect_xyz(ijk_xyz_sta_ext:ijk_xyz_end_ext), STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "62 Cannot allocate memory"
allocate (int_index_xyz(ijk_xyz_end - ijk_xyz_sta + 1), STAT=status(1))
if (status(1)/=0) stop "63 Cannot allocate memory"
do i = 1,ijk_xyz_end - ijk_xyz_sta + 1
int_index_xyz(i) = ijk_xyz_sta + i - 1 - 1
end do
end if
end subroutine semi_allo_var
subroutine n_bodies_allo_var1
!allocate memory for n bodies, part 1, only essential
if (.not. allocated(n_hy0)) then
allocate (n_hy0(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "qCannot allocate memory"
allocate (n_freq(no_body), STAT=status(1)); allocate (n_phase_ang(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "wCannot allocate memory"
allocate (n_theta0_x(no_body), STAT=status(1)); allocate (n_theta0_y(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "eCannot allocate memory"
allocate (n_theta0_z(no_body), STAT=status(1)); allocate (n_loc_rot(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "rCannot allocate memory"
allocate (n_leading_edge_pos(no_body), STAT=status(1)); allocate (n_plate_h_tk(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "tCannot allocate memory"
allocate (n_plate_lena(no_body), STAT=status(1)); allocate (n_plate_lenb(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "yCannot allocate memory"
allocate (n_plate_wid(no_body), STAT=status(1)); allocate (n_hx0(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "uCannot allocate memory"
allocate (n_shift_x(no_body), STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "iCannot allocate memory"
allocate (z_1_to_n(no_body), STAT=status(1)); allocate (phase_ang_1_to_n(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "oCannot allocate memory"
allocate (n_hz0(no_body), STAT=status(1)); allocate (n_ini_hx0(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "oCannot allocate memory"
allocate (n_ini_hy0(no_body), STAT=status(1)); allocate (n_ini_theta0_x(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "pCannot allocate memory"
allocate (n_ini_theta0_y(no_body), STAT=status(1)); allocate (n_ini_theta0_z(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "aCannot allocate memory"
allocate (n_body_cg_ini(3,no_body), STAT=status(1)); allocate (n_body_cg(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "sCannot allocate memory"
allocate (n_body_cg_old(3,no_body), STAT=status(1)); allocate (n_body_cg_old2(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "dCannot allocate memory"
allocate (n_rot_body(3,no_body), STAT=status(1)); allocate (n_rot_body_old(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "fCannot allocate memory"
allocate (n_fixed_rot_body(3,no_body), STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "fCannot allocate memory"
end if
end subroutine n_bodies_allo_var1
subroutine n_bodies_allo_var2
!allocate memory for n bodies part 2
if (.not. allocated(n_IIB_I_cell_no_uvw_global)) then
allocate (n_IIB_I_cell_no_uvw_global(0:6*num_procs-1,no_body), STAT=status(1)); allocate (n_IIB_I_cell_no_p_global(0:2*num_procs-1,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1834 Cannot allocate memory"
allocate (n_IIB_equal_cell_no_global_u(0:num_procs-1,no_body), STAT=status(1));allocate (n_IIB_equal_cell_no_global_v(0:num_procs-1,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "193 Cannot allocate memory"
allocate (n_IIB_equal_cell_no_global_w(0:num_procs-1,no_body), STAT=status(1));allocate (n_I_equal_cell_no_global_u(0:num_procs-1,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "203 Cannot allocate memory"
allocate (n_I_equal_cell_no_global_v(0:num_procs-1,no_body), STAT=status(1));allocate (n_I_equal_cell_no_global_w(0:num_procs-1,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "213 Cannot allocate memory"
allocate (n_IIB_equal_cell_no_global_p(0:num_procs-1,no_body), STAT=status(1));allocate (n_I_equal_cell_no_global_p(0:num_procs-1,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "223 Cannot allocate memory"
allocate (n_sta_yy(sta_ijk_uniform(1):end_ijk_uniform(1),no_body), STAT=status(1));allocate (n_end_yy(sta_ijk_uniform(1):end_ijk_uniform(1),no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "593 Cannot allocate memory"
allocate (n_surface_normal_correct(no_body), STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "qCannot allocate memory"
allocate (n_pvm_xyz(9,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "gCannot allocate memory"
allocate (n_cd_cl_cs_mom_implicit(6,no_body), STAT=status(1)); allocate (n_cd_cl_cs_mom_explicit(6,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "hCannot allocate memory"
allocate (n_cd_cl_cs_mom_power(8,no_body), STAT=status(1)); allocate (n_cd_tmp(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "jCannot allocate memory"
allocate (n_cl_tmp(no_body), STAT=status(1)); allocate (n_cs_tmp(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "kCannot allocate memory"
allocate (n_body_vel(3,no_body), STAT=status(1)); allocate (n_body_vel_old(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "lCannot allocate memory"
allocate (n_body_displacement(3,no_body), STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "zCannot allocate memory"
allocate (n_IIB_cell_no_u(no_body), STAT=status(1)); allocate (n_IIB_cell_no_v(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "xCannot allocate memory"
allocate (n_IIB_cell_no_w(no_body), STAT=status(1)); allocate (n_IIB_cell_no_p(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "cCannot allocate memory"
allocate (n_I_cell_no_u(no_body), STAT=status(1)); allocate (n_I_cell_no_v(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "vCannot allocate memory"
allocate (n_I_cell_no_w(no_body), STAT=status(1)); allocate (n_I_cell_no_p(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "bCannot allocate memory"
allocate (n_IIB_global_cell_no_u(no_body), STAT=status(1)); allocate (n_IIB_global_cell_no_v(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "nCannot allocate memory"
allocate (n_IIB_global_cell_no_w(no_body), STAT=status(1)); allocate (n_IIB_global_cell_no_p(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "mCannot allocate memory"
allocate (n_I_global_cell_no_u(no_body), STAT=status(1)); allocate (n_I_global_cell_no_v(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1Cannot allocate memory"
allocate (n_I_global_cell_no_w(no_body), STAT=status(1)); allocate (n_I_global_cell_no_p(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "2Cannot allocate memory"
allocate (n_IIB_equal_cell_no_u(no_body), STAT=status(1)); allocate (n_IIB_equal_cell_no_v(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "3Cannot allocate memory"
allocate (n_IIB_equal_cell_no_w(no_body), STAT=status(1)); allocate (n_IIB_equal_cell_no_p(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "4Cannot allocate memory"
allocate (n_I_equal_cell_no_u(no_body), STAT=status(1)); allocate (n_I_equal_cell_no_v(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "5Cannot allocate memory"
allocate (n_I_equal_cell_no_w(no_body), STAT=status(1)); allocate (n_I_equal_cell_no_p(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "6Cannot allocate memory"
allocate (n_IIB_I_cell_no_uvw_total(6,no_body), STAT=status(1)); allocate (n_IIB_I_cell_no_p_total(2,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "7Cannot allocate memory"
allocate (n_IIB_equal_cell_no_u_max(no_body), STAT=status(1)); allocate (n_I_equal_cell_no_u_max(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "71Cannot allocate memory"
allocate (n_IIB_cell_no_u_max(no_body), STAT=status(1)); allocate (n_I_cell_no_u_max(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "72Cannot allocate memory"
allocate (n_IIB_global_cell_no_u_max(no_body), STAT=status(1)); allocate (n_I_global_cell_no_u_max(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "73Cannot allocate memory"
allocate (n_sta_x(no_body), STAT=status(1)); allocate (n_end_x(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "8Cannot allocate memory"
allocate (n_sta_y(no_body), STAT=status(1)); allocate (n_end_y(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "9Cannot allocate memory"
allocate (n_sta_z(no_body), STAT=status(1)); allocate (n_end_z(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "0Cannot allocate memory"
allocate (n_no_vertices(no_body), STAT=status(1)); allocate (n_no_surfaces(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "aaCannot allocate memory"
allocate (n_body_sta_ijk_IIB(3,no_body), STAT=status(1)); allocate (n_body_end_ijk_IIB(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "ssCannot allocate memory"
allocate (n_body_sta_ijk_IIB_ghost(3,no_body), STAT=status(1)); allocate (n_body_end_ijk_IIB_ghost(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "ddCannot allocate memory"
!allocate (n_body_sta_ijk_IIB_ext(3,no_body), STAT=status(1)); allocate (n_body_end_ijk_IIB_ext(3,no_body), STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot allocate memory"
allocate (n_body_sta_ijk_IIB_array(3,no_body), STAT=status(1)); allocate (n_body_end_ijk_IIB_array(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "ffCannot allocate memory"
allocate (n_body_sta_ijk(3,no_body), STAT=status(1)); allocate (n_body_end_ijk(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "ggCannot allocate memory"
allocate (n_body_sta_ijk_ghost(3,no_body), STAT=status(1)); allocate (n_body_end_ijk_ghost(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "hhCannot allocate memory"
allocate (n_no_body_pts(no_body), STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "jjCannot allocate memory"
allocate (n_sta_cept_xyz(3,no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "kkCannot allocate memory"
allocate (n_end_cept_xyz(3,no_body), STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "llCannot allocate memory"
allocate (n_IIB_ista(no_body), STAT=status(1)); allocate (n_IIB_iend(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "zzCannot allocate memory"
allocate (n_IIB_jsta(no_body), STAT=status(1)); allocate (n_IIB_jend(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "xxCannot allocate memory"
allocate (n_IIB_ksta(no_body), STAT=status(1)); allocate (n_IIB_kend(no_body), STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "ccCannot allocate memory"
end if
end subroutine n_bodies_allo_var2
subroutine ini_var
!initialize/allocate memory for variables
!grid variables
!#include "petsc/finclude/petsc.h"
!x=0.;y=0.;z=0.
if (new_start == 1) then
act_time = 0.d0; time_sta = 0.d0
end if
xu = 0.; yu = 0.; zu = 0.; xv = 0.; yv = 0.; zv = 0.; xw = 0.; yw = 0.; zw = 0.
near_2_body12_global = 0; near_2_body23_global = 0
if (inout_test <= 1) then
!flow variables
!DM stuffs
!if (myid==0) print *,"0"
call DMDAVecGetArrayF90(da_u,u_local,u_array,ierr) !cannot directly access u_local/global. must use a 'normal' array
call DMDAVecGetArrayF90(da_v,v_local,v_array,ierr)
call DMDAVecGetArrayF90(da_w,w_local,w_array,ierr)
call DMDAVecGetArrayF90(da_p,p_local,p_array,ierr)
call DMDAVecGetArrayF90(da_uvw_ast,uvw_ast_local,uvw_ast_array,ierr)
call DMDAVecGetArrayF90(da_q_semi_vect_xyz,q_semi_vect_xyz_local,q_semi_vect_xyz_array,ierr)
call DMDAVecGetArrayF90(da_q_p,q_p_local,q_p_array,ierr)
u_array = 0.d0
v_array = 0.d0
w_array = 0.d0
p_array = 0.d0
uvw_ast_array = 0.d0
q_semi_vect_xyz_array = 0.d0
q_p_array = 0.d0
call DMDAVecRestoreArrayF90(da_p,p_local,p_array,ierr)
call DMDAVecRestoreArrayF90(da_w,w_local,w_array,ierr)
call DMDAVecRestoreArrayF90(da_v,v_local,v_array,ierr)
call DMDAVecRestoreArrayF90(da_u,u_local,u_array,ierr)
call DMDAVecRestoreArrayF90(da_uvw_ast,uvw_ast_local,uvw_ast_array,ierr)
call DMDAVecRestoreArrayF90(da_q_semi_vect_xyz,q_semi_vect_xyz_local,q_semi_vect_xyz_array,ierr)
call DMDAVecRestoreArrayF90(da_q_p,q_p_local,q_p_array,ierr)
!call DMDAVecRestoreArrayF90(da_u,u_local,u_array,ierr) !cannot directly access u_local/global. must use a 'normal' array
!call DMDAVecRestoreArrayF90(da_v,v_local,v_array,ierr); call DMDAVecRestoreArrayF90(da_w,w_local,w_array,ierr)
!call DMDAVecRestoreArrayF90(da_p,p_local,p_array,ierr)
!if (myid==0) print *,"02"
u_old=0.;v_old=0.;w_old=0.;diff_pres_x=0.;diff_pres_y=0.;diff_pres_z=0.;u_ast2=0.;v_ast2=0.;w_ast2=0.
!phi=0.
!iterative solver variables
big_A=0.d0; int_A = 0; qk=0.d0 !;q_p=0.d0
!semi_mat_xyz = 0.d0; int_semi_xyz = 0;
semi_mat_xyz_da = 0.d0
if (avg_vorticity == 1) then
vorticity_avg = 0.;vorticity_tmp = 0.
end if
if (q_criterion_plot == 2) then
Q_criterion_avg = 0.;Q_criterion_tmp = 0.
end if
!ibm variables
!ibm_force=0.
!q_semi_x_force=0.;q_semi_y_force=0.;q_semi_z_force=0.
if (spatial_scheme == 3 .or. spatial_scheme == 6) then
s_cf_u_e = 0.d0; s_cf_v_e = 0.d0; s_cf_u_n = 0.d0; s_cf_w_e = 0.d0; s_cf_v_f = 0.d0
s_cf_v_n = 0.d0; s_cf_w_n = 0.d0; s_cf_u_f = 0.d0; s_cf_w_f = 0.d0
end if
end if
if (no_body <= 3) then
IIB_cell_no_u1=0; IIB_cell_no_v1=0; IIB_cell_no_w1=0; IIB_cell_no_p1=0
IIB_I_cell_no_uvw_total1 = 0; IIB_I_cell_no_p_total1 = 0
if (no_body >= 2) then
IIB_cell_no_u2=0; IIB_cell_no_v2=0; IIB_cell_no_w2=0; IIB_cell_no_p2=0
IIB_I_cell_no_uvw_total2 = 0; IIB_I_cell_no_p_total2 = 0
if (no_body == 3) then
IIB_cell_no_u3=0; IIB_cell_no_v3=0; IIB_cell_no_w3=0; IIB_cell_no_p3=0
IIB_I_cell_no_uvw_total3 = 0; IIB_I_cell_no_p_total3 = 0
end if
end if
else
n_IIB_cell_no_u(:) = 0; n_IIB_cell_no_v(:) = 0; n_IIB_cell_no_w(:) = 0; n_IIB_cell_no_p(:) = 0
n_I_cell_no_u(:) = 0; n_I_cell_no_v(:) = 0; n_I_cell_no_w(:) = 0; n_I_cell_no_p(:) = 0
n_IIB_global_cell_no_u = 0; n_IIB_global_cell_no_v = 0; n_IIB_global_cell_no_w = 0; n_IIB_global_cell_no_p = 0
n_I_global_cell_no_u = 0; n_I_global_cell_no_v = 0; n_I_global_cell_no_w = 0; n_I_global_cell_no_p = 0
n_IIB_equal_cell_no_u = 0; n_IIB_equal_cell_no_v = 0; n_IIB_equal_cell_no_w = 0; n_IIB_equal_cell_no_p = 0
n_I_equal_cell_no_u = 0; n_I_equal_cell_no_v = 0; n_I_equal_cell_no_w = 0; n_I_equal_cell_no_p = 0
n_IIB_equal_cell_no_u_max = 0; n_I_equal_cell_no_u_max = 0; n_IIB_cell_no_u_max = 0; n_I_cell_no_u_max = 0; n_IIB_global_cell_no_u_max = 0; n_I_global_cell_no_u_max = 0
n_IIB_I_cell_no_uvw_total(:,:) = 0; n_IIB_I_cell_no_p_total(:,:) = 0
end if
!IIB_I_cu(1:size_x,jsta:jend,ksta:kend)%types=0 !inout_test_cell_uv
!IIB_I_cu(1:size_x,jsta:jend,ksta:kend)%types_old = 0
!cu(1:size_x,jsta:jend,ksta:kend)%types_temp = 0
!IIB_index_u = 0; IIB_index_v = 0; IIB_index_w = 0; IIB_index_p = 0
!pIB_IB_vel = 0.; pIB_IB_vel_old = 0.
!if (inout_test <= 1) then
! IIB_I_cp(1:size_x,1:size_y,1:size_z)%types=0
! IIB_I_cv(1:size_x,1:size_y,1:size_z)%types=0
! IIB_I_cw(1:size_x,1:size_y,1:size_z)%types=0
! IIB_I_cp(1:size_x,1:size_y,1:size_z)%types_old = 0
! IIB_I_cv(1:size_x,1:size_y,1:size_z)%types_old = 0
! IIB_I_cw(1:size_x,1:size_y,1:size_z)%types_old = 0
!cp(1:size_x,jsta:jend,ksta:kend)%types_temp = 0; cv(1:size_x,jsta:jend,ksta:kend)%types_temp = 0; cw(1:size_x,jsta:jend,ksta:kend)%types_temp = 0
!end if
!cp(1:size_x,jsta:jend,ksta:kend)%types01=0; cu(1:size_x,jsta:jend,ksta:kend)%types01=0; cv(1:size_x,jsta:jend,ksta:kend)%types01=0; cw(1:size_x,jsta:jend,ksta:kend)%types01=0
!cp(1:size_x,jsta:jend,ksta:kend)%vel=0.; cu(1:size_x,jsta:jend,ksta:kend)%vel=0.; cv(1:size_x,jsta:jend,ksta:kend)%vel=0.; cw(1:size_x,jsta:jend,ksta:kend)%vel=0.
!airfoil variables
if (no_body <= 3) then
sta_x1 = size_x/2; end_x1 = size_x/2; sta_xu1 = size_x/2; end_xu1 = size_x/2 !guess
sta_y1 = size_y/2; end_y1 = size_y/2; sta_yv1 = size_y/2; end_yv1 = size_y/2 !guess
sta_z1 = (kend - ksta + 1)/2; end_z1 = (kend - ksta + 1)/2; sta_zw1 = (kend - ksta + 1)/2; end_zw1 = (kend - ksta + 1)/2 !guess
body_vel1 = 0.d0; body_vel1_old = 0.d0; body_displacement1 = 0.d0
if (no_body >= 2) then
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
sta_z2 = (kend - ksta + 1)/2; end_z2 = (kend - ksta + 1)/2; sta_zw2 = (kend - ksta + 1)/2; end_zw2 = (kend - ksta + 1)/2 !guess
body_vel2 = 0.d0; body_vel2_old = 0.d0; body_displacement2 = 0.d0
if (no_body == 3) then
sta_x3 = size_x/2; end_x3 = size_x/2; sta_xu3 = size_x/2; end_xu3 = size_x/2 !guess
sta_y3 = size_y/2; end_y3 = size_y/2; sta_yv3 = size_y/2; end_yv3 = size_y/2 !guess
sta_z3 = (kend - ksta + 1)/2; end_z3 = (kend - ksta + 1)/2; sta_zw3 = (kend - ksta + 1)/2; end_zw3 = (kend - ksta + 1)/2 !guess
body_vel3 = 0.d0; body_vel3_old = 0.d0; body_displacement3 = 0.d0
end if
end if
else
n_sta_x(:) = size_x/2; n_end_x(:) = size_x/2
n_sta_y(:) = size_y/2; n_end_y(:) = size_y/2
n_sta_z(:) = (kend - ksta + 1)/2; n_end_z(:) = (kend - ksta + 1)/2
n_body_vel(:,:) = 0.d0; n_body_vel_old(:,:) = 0.d0; n_body_displacement(:,:) = 0.d0
end if
!already done at creation
!if (inout_test <= 2) IIB_cell_u1(:)%stencil_no = 4
!if (inout_test <= 1) then
! IIB_cell_v1(:)%stencil_no = 4; IIB_cell_w1(:)%stencil_no = 4; IIB_cell_p1(:)%stencil_no = 4
!end if
!if (no_body >= 2) then
! if (inout_test <= 2) IIB_cell_u2(:)%stencil_no = 4
! if (inout_test <= 1) then
! IIB_cell_v2(:)%stencil_no = 4; IIB_cell_w2(:)%stencil_no = 4; IIB_cell_p2(:)%stencil_no = 4
! end if
! if (no_body == 3) then
! if (inout_test <= 2) IIB_cell_u3(:)%stencil_no = 4
! if (inout_test <= 1) then
! IIB_cell_v3(:)%stencil_no = 4; IIB_cell_w3(:)%stencil_no = 4; IIB_cell_p3(:)%stencil_no = 4
! end if
! end if
!end if
end subroutine ini_var
subroutine reset_update_var
!reset variables
!do ij=1,I_cell_no_u
! IIB_cell_u(ij)%int_template_ij(:)=0
!end do
!do ij=1,I_cell_no_v
! IIB_cell_v(ij)%int_template_ij(:)=0
!end do
!do ij=1,I_cell_no_p
! IIB_cell_p(ij)%int_template_ij(:)=0
!end do
!I_cell_no_u=0; I_cell_no_v=0; I_cell_no_w=0; IIB_cell_no_u=0; IIB_cell_no_v=0; IIB_cell_no_w=0
!I_cell_no_p=0; IIB_cell_no_p=0
!IIB_I_cell_no_uvw_total = 0; IIB_I_cell_no_p_total=0
!if (body2_unsteady == 1) then
! IIB_I_cell_no_uvw_total2 = 0; IIB_I_cell_no_p_total2 = 0
!end if
!pIB_cell_no_u=0; pIB_cell_no_v=0; pIB_cell_no_w=0; pIB_cell_no_p=0
!pIB_cell_no_u1=0; pIB_cell_no_v1=0; pIB_cell_no_p1=0
!cu(1:size_x,jsta:jend,ksta:kend)%types_old = cu(1:size_x,jsta:jend,ksta:kend)%types
!IIB_cell_u(:)%vel_b=0.
!I_cell_u(:)%vel=0.
!if (fsi_run == 0 .or. (fsi_run >= 1 .and. fsi_iteration_completed == 1)) then
! body_cg1_old2 = body_cg1_old !req to update body_cg
! body_cg1_old = body_cg1
! rot_body1_old = rot_body1
!end if
near_2_body12_global = 0; near_2_body23_global = 0; near_2_body12_global_tmp = 0; near_2_body23_global_tmp = 0
if (no_body <= 3 .and. near_2_body_check == 1 .and. inout_test == 0) then !only for max 3 bodies at the moment
IIB_cell_u1(:)%near_2_body = 0; IIB_cell_v1(:)%near_2_body = 0; IIB_cell_w1(:)%near_2_body = 0; IIB_cell_p1(:)%near_2_body = 0
if (no_body >= 2) then
IIB_cell_u2(:)%near_2_body = 0; IIB_cell_v2(:)%near_2_body = 0; IIB_cell_w2(:)%near_2_body = 0; IIB_cell_p2(:)%near_2_body = 0
if (no_body == 3) then
IIB_cell_u3(:)%near_2_body = 0; IIB_cell_v3(:)%near_2_body = 0; IIB_cell_w3(:)%near_2_body = 0; IIB_cell_p3(:)%near_2_body = 0
end if
end if
end if
if (no_body <= 3) then
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. motion_type == 8 .or. deformation /= 0 .or. (fsi_run >= 1 .and. fsi_iter_no == 1)) then
!if (deformation /= 0 .or. fsi_run >= 1) then !will use old vertex position to get velocity / acc
!if (time == start_time) then
! if (myid == 0) print *, "body_cg1_old,body_cg1,rot_body1",body_cg1_old,body_cg1,rot_body1
!end if
body_cg1_old2 = body_cg1_old !req to update body_cg
body_cg1_old = body_cg1
rot_body1_old = rot_body1
!call kdtree2_destroy(tree_surface1) !destroy old one 1st !kdtree_2
end if
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0 .or. (fsi_run >= 1 .and. fsi_iter_no == 1)) then
vertex1_old2 = vertex1_old
vertex1_old = vertex1
end if
if (body_input_type == 1) body_pts1_old = body_pts1
!if (inout_test <= 1) then
!cv(1:size_x,jsta:jend,ksta:kend)%types_old=cv(1:size_x,jsta:jend,ksta:kend)%types; cw(1:size_x,jsta:jend,ksta:kend)%types_old=cw(1:size_x,jsta:jend,ksta:kend)%types; cp(1:size_x,jsta:jend,ksta:kend)%types_old=cp(1:size_x,jsta:jend,ksta:kend)%types
!cu(1:size_x,jsta:jend,ksta:kend)%types=0; cv(1:size_x,jsta:jend,ksta:kend)%types=0; cw(1:size_x,jsta:jend,ksta:kend)%types=0; cp(1:size_x,jsta:jend,ksta:kend)%types=0
!cp(1:size_x,jsta:jend,ksta:kend)%types01=0; cu(1:size_x,jsta:jend,ksta:kend)%types01=0; cv(1:size_x,jsta:jend,ksta:kend)%types01=0; cw(1:size_x,jsta:jend,ksta:kend)%types01=0
!IIB_cell_v1(:)%stencil_no = 4; IIB_cell_w1(:)%stencil_no = 4; IIB_cell_p1(:)%stencil_no = 4
!cp(1:size_x,jsta:jend,ksta:kend)%vel=0.; cu(1:size_x,jsta:jend,ksta:kend)%vel=0.; cv(1:size_x,jsta:jend,ksta:kend)%vel=0.; cw(1:size_x,jsta:jend,ksta:kend)%vel=0.
!ibm_force=0.
!q_semi_x_force=0.; q_semi_y_force=0.; q_semi_z_force=0.
!IIB_cell_v(:)%vel_b=0.; IIB_cell_w(:)%vel_b=0.; IIB_cell_p(:)%vel_b=0.
!I_cell_v(:)%vel=0.; I_cell_w(:)%vel=0.; I_cell_p(:)%vel=0.
!end if
if (no_body >= 2) then
!if (inout_test <= 2) IIB_cell_u2(:)%stencil_no = 4
!if (inout_test <= 1) then
! IIB_cell_v2(:)%stencil_no = 4; IIB_cell_w2(:)%stencil_no = 4; IIB_cell_p2(:)%stencil_no = 4
!end if
!if (deformation /= 0 .or. fsi_run >= 1) then !will use old vertex position to get velocity / acc
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0) then
body_cg2_old2 = body_cg2_old !req to update body_cg
body_cg2_old = body_cg2
rot_body2_old = rot_body2
end if
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0) then
vertex2_old2 = vertex2_old
vertex2_old = vertex2
end if
if (body_input_type == 1) body_pts2_old = body_pts2
if (no_body == 3) then
!if (inout_test <= 2) IIB_cell_u3(:)%stencil_no = 4
!if (inout_test <= 1) then
! IIB_cell_v3(:)%stencil_no = 4; IIB_cell_w3(:)%stencil_no = 4; IIB_cell_p3(:)%stencil_no = 4
!end if
!if (deformation /= 0 .or. fsi_run >= 1) then !will use old vertex position to get velocity / acc
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0) then
body_cg3_old2 = body_cg3_old !req to update body_cg
body_cg3_old = body_cg3
rot_body3_old = rot_body3
end if
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0) then
vertex3_old2 = vertex3_old
vertex3_old = vertex3
end if
if (body_input_type == 1) body_pts3_old = body_pts3
end if
end if
else
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. motion_type == 8 .or. deformation /= 0 .or. (fsi_run >= 1 .and. fsi_iter_no == 1)) then
n_body_cg_old2 = n_body_cg_old !req to update body_cg
n_body_cg_old = n_body_cg
n_rot_body_old = n_rot_body
end if
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0 .or. (fsi_run >= 1 .and. fsi_iter_no == 1)) then
n_vertex_old2 = n_vertex_old
n_vertex_old = n_vertex
end if
if (body_input_type == 1) n_body_pts_old = n_body_pts
end if
end subroutine reset_update_var
subroutine de_ini_var
!de-initialize/deallocate memory for variables
!#include "petsc/finclude/petsc.h"
PetscInt :: status(3)
deallocate (num_procs_jsta_jend, STAT=status(1));deallocate (temp_procs, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1 Cannot deallocate memory"
deallocate (num_procs_ksta_kend, STAT=status(1));deallocate (num_procs_size, STAT=status(2));deallocate (num_procs_size_start_index, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "2 Cannot deallocate memory"
deallocate (xu, STAT=status(1));deallocate (yu, STAT=status(2));deallocate (zu, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "3 Cannot deallocate memory"
deallocate (xv, STAT=status(1));deallocate (yv, STAT=status(2));deallocate (zv, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "4 Cannot deallocate memory"
deallocate (xw, STAT=status(1));deallocate (yw, STAT=status(2));deallocate (zw, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "5 Cannot deallocate memory"
deallocate (c_cx, STAT=status(1));deallocate (c_cy, STAT=status(2));deallocate (c_cz, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "6 Cannot deallocate memory"
deallocate (cu_cx, STAT=status(1));deallocate (cv_cy, STAT=status(2));deallocate (cw_cz, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "7 Cannot deallocate memory"
deallocate (x, STAT=status(1));deallocate (y, STAT=status(2));deallocate (z, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "8 Cannot deallocate memory"
deallocate (cu_x, STAT=status(1));deallocate (cv_x, STAT=status(2));deallocate (cw_x, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "9 Cannot deallocate memory"
deallocate (cu_y, STAT=status(1));deallocate (cv_y, STAT=status(2));deallocate (cw_y, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "10 Cannot deallocate memory"
deallocate (cu_z, STAT=status(1));deallocate (cv_z, STAT=status(2));deallocate (cw_z, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "11 Cannot deallocate memory"
deallocate (cp_x, STAT=status(1));deallocate (cp_y, STAT=status(1));deallocate (cp_z, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "12 Cannot deallocate memory"
deallocate (cu, STAT=status(1))
if (status(1)/=0) stop "13 Cannot deallocate memory"
deallocate (IIB_I_cu, STAT=status(1))
if (status(1)/=0) stop "14 Cannot deallocate memory"
if (no_body <= 3) then
if (inout_test <= 2) then
deallocate (IIB_cell_u1, STAT=status(1))
if (status(1)/=0) stop "15 Cannot deallocate memory"
deallocate (I_cell_u1, STAT=status(1))
if (status(1)/=0) stop "16 Cannot deallocate memory"
deallocate (IIB_global_cell_u1, STAT=status(1))
if (status(1)/=0) stop "17 Cannot deallocate memory"
deallocate (I_global_cell_u1, STAT=status(1))
if (status(1)/=0) stop "18 Cannot deallocate memory"
deallocate (IIB_equal_cell_u1, STAT=status(1))
if (status(1)/=0) stop "19 Cannot deallocate memory"
deallocate (I_equal_cell_u1, STAT=status(1))
if (status(1)/=0) stop "20 Cannot deallocate memory"
if (no_body >= 2) then
deallocate (IIB_cell_u2, STAT=status(1));deallocate (I_cell_u2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "21 Cannot deallocate memory"
deallocate (IIB_global_cell_u2, STAT=status(1));deallocate (I_global_cell_u2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "22 Cannot deallocate memory"
deallocate (IIB_equal_cell_u2, STAT=status(1));deallocate (I_equal_cell_u2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "23 Cannot deallocate memory"
if (no_body == 3) then
deallocate (IIB_cell_u3, STAT=status(1));deallocate (I_cell_u3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "24 Cannot deallocate memory"
deallocate (IIB_global_cell_u3, STAT=status(1));deallocate (I_global_cell_u3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "25 Cannot deallocate memory"
deallocate (IIB_equal_cell_u3, STAT=status(1));deallocate (I_equal_cell_u3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "26 Cannot deallocate memory"
end if
end if
end if
deallocate (IIB_I_cell_no_uvw_global1, STAT=status(1)); deallocate (IIB_I_cell_no_p_global1, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "27 Cannot deallocate memory"
deallocate (IIB_equal_cell_no_global_u1, STAT=status(1));deallocate (IIB_equal_cell_no_global_v1, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "28 Cannot deallocate memory"
deallocate (IIB_equal_cell_no_global_w1, STAT=status(1));deallocate (I_equal_cell_no_global_u1, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "29 Cannot deallocate memory"
deallocate (I_equal_cell_no_global_v1, STAT=status(1));deallocate (I_equal_cell_no_global_w1, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "30 Cannot deallocate memory"
deallocate (IIB_equal_cell_no_global_p1, STAT=status(1));deallocate (I_equal_cell_no_global_p1, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "31 Cannot deallocate memory"
if (no_body >= 2) then
deallocate (IIB_I_cell_no_uvw_global2, STAT=status(1)); deallocate (IIB_I_cell_no_p_global2, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "32 Cannot deallocate memory"
deallocate (IIB_equal_cell_no_global_u2, STAT=status(1));deallocate (IIB_equal_cell_no_global_v2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "33 Cannot deallocate memory"
deallocate (IIB_equal_cell_no_global_w2, STAT=status(1));deallocate (I_equal_cell_no_global_u2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "34 Cannot deallocate memory"
deallocate (I_equal_cell_no_global_v2, STAT=status(1));deallocate (I_equal_cell_no_global_w2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "35 Cannot deallocate memory"
deallocate (IIB_equal_cell_no_global_p2, STAT=status(1));deallocate (I_equal_cell_no_global_p2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "36 Cannot deallocate memory"
if (no_body == 3) then
deallocate (IIB_I_cell_no_uvw_global3, STAT=status(1)); deallocate (IIB_I_cell_no_p_global3, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "37 Cannot deallocate memory"
deallocate (IIB_equal_cell_no_global_u3, STAT=status(1));deallocate (IIB_equal_cell_no_global_v3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "38 Cannot deallocate memory"
deallocate (IIB_equal_cell_no_global_w3, STAT=status(1));deallocate (I_equal_cell_no_global_u3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "39 Cannot deallocate memory"
deallocate (I_equal_cell_no_global_v3, STAT=status(1));deallocate (I_equal_cell_no_global_w3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "40 Cannot deallocate memory"
deallocate (IIB_equal_cell_no_global_p3, STAT=status(1));deallocate (I_equal_cell_no_global_p3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "41 Cannot deallocate memory"
end if
end if
deallocate (surface1, STAT=status(1)) ; deallocate (vertex1, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "42 Cannot deallocate memory"
if (allocated(vertex1_ini)) then
deallocate (vertex1_ini, STAT=status(1))
if (status(1)/=0) stop "43 Cannot deallocate memory"
end if
if (body_input_type /= 3) then
!deallocate (vertex1_ini, STAT=status(1))
deallocate (surface1_pts, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "43 Cannot deallocate memory"
end if
if (body_input_type == 3) then
deallocate (vertex1_multi, STAT=status(1))
if (status(1)/=0) stop "45 Cannot deallocate memory"
if (deformation == 0) then
deallocate (fft_coeff_a, STAT=status(1)) ; deallocate (fft_coeff_b, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "46 Cannot deallocate memory"
deallocate (cosF, STAT=status(1)) ; deallocate (sinF, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "47 Cannot deallocate memory"
deallocate (discrete_time, STAT=status(1))
if (status(1)/=0) stop "48 Cannot deallocate memory"
end if
if (no_body >= 2) then
if (deformation == 0) then
deallocate (fft_coeff_a2, STAT=status(1)) ; deallocate (fft_coeff_b2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "49 Cannot deallocate memory"
if (no_body == 3) then
deallocate (fft_coeff_a3, STAT=status(1)) ; deallocate (fft_coeff_b3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "50 Cannot deallocate memory"
end if
end if
end if
end if
!if (deformation /= 0 .or. fsi_run >= 1) then
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0 .or. fsi_run >= 1) then
deallocate (vertex1_old, STAT=status(1)) !;
deallocate (vertex1_old2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "51 Cannot deallocate memory"
end if
if (no_body >= 2) then
deallocate (surface2, STAT=status(1)) ; deallocate (vertex2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "52 Cannot deallocate memory"
deallocate (vertex2_ini, STAT=status(1))
if (status(1)/=0) stop "53 Cannot deallocate memory"
if (body_input_type /= 3) then
deallocate (surface2_pts, STAT=status(2))
if (status(2)/=0) stop "54 Cannot deallocate memory"
end if
if (body_input_type == 3) then
deallocate (vertex2_multi, STAT=status(1))
if (status(1)/=0) stop "55 Cannot deallocate memory"
end if
!if (deformation /= 0) then
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0) then
deallocate (vertex2_old, STAT=status(1)) !; deallocate (vertex2_old2, STAT=status(2))
deallocate (vertex2_old2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "56 Cannot deallocate memory"
end if
if (no_body == 3) then
deallocate (surface3, STAT=status(1)) ; deallocate (vertex3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "57 Cannot deallocate memory"
deallocate (vertex3_ini, STAT=status(1))
if (status(1)/=0) stop "58 Cannot deallocate memory"
if (body_input_type /= 3) then
deallocate (surface3_pts, STAT=status(2))
if (status(2)/=0) stop "59 Cannot deallocate memory"
end if
if (body_input_type == 3) then
deallocate (vertex3_multi, STAT=status(1))
if (status(1)/=0) stop "60 Cannot deallocate memory"
end if
!if (deformation /= 0) then
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0) then
deallocate (vertex3_old, STAT=status(1)) !; deallocate (vertex3_old2, STAT=status(2))
deallocate (vertex3_old2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "61 Cannot deallocate memory"
end if
end if
end if
if (body_input_type == 1) then
deallocate (body_pts1, STAT=status(1));deallocate (body_length_seg1, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "62 Cannot deallocate memory"
deallocate (body_pts1_old, STAT=status(1)); deallocate (body_pts1_ini, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "63 Cannot deallocate memory"
if (no_body >= 2) then
deallocate (body_pts2, STAT=status(1));deallocate (body_length_seg2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "64 Cannot deallocate memory"
deallocate (body_pts2_old, STAT=status(1)); deallocate (body_pts2_ini, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "65 Cannot deallocate memory"
if (no_body == 3) then
deallocate (body_pts3, STAT=status(1));deallocate (body_length_seg3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "66 Cannot deallocate memory"
deallocate (body_pts3_old, STAT=status(1)); deallocate (body_pts3_ini, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "67 Cannot deallocate memory"
end if
end if
end if
if (body_input_type == 2) then
deallocate (sta_yy1, STAT=status(1));deallocate (end_yy1, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "77 Cannot deallocate memory"
if (no_body >= 2) then
deallocate (sta_yy2, STAT=status(1));deallocate (end_yy2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "78 Cannot deallocate memory"
if (no_body == 3) then
deallocate (sta_yy3, STAT=status(1));deallocate (end_yy3, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "79 Cannot deallocate memory"
end if
end if
end if
end if
call DMDestroy(da_u,ierr)
call VecDestroy(u_local,ierr)
call VecDestroy(u_global,ierr)
call DMDestroy(da_cu_types,ierr)
if (inout_test <= 1) then
!deallocate (p, STAT=status(1))
!if (status(1)/=0) stop "Cannot deallocate memory"
!deallocate (u, STAT=status(1));deallocate (v, STAT=status(2));deallocate (w, STAT=status(3))
!if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "Cannot deallocate memory"
deallocate (u_old, STAT=status(1));deallocate (v_old, STAT=status(2));deallocate (w_old, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "83 Cannot deallocate memory"
!deallocate (u_ast, STAT=status(1));deallocate (v_ast, STAT=status(2));deallocate (w_ast, STAT=status(3))
!if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "Cannot deallocate memory"
deallocate (u_ast2, STAT=status(1));deallocate (v_ast2, STAT=status(2));deallocate (w_ast2, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "84 Cannot deallocate memory"
deallocate (phi, STAT=status(1))
if (status(1)/=0) stop "85 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 (qk, STAT=status(2)) !;deallocate (big_A, STAT=status(1))
if (status(2)/=0) stop "86 Cannot deallocate memory"
!deallocate (int_A, STAT=status(1))
!if (status(1)/=0) stop "Cannot deallocate memory"
!deallocate (semi_mat_xyz, STAT=status(1));deallocate (int_semi_xyz, STAT=status(2))
deallocate (semi_mat_xyz_da, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "87 Cannot deallocate memory"
deallocate (int_index_xyz, STAT=status(1))
if (status(1)/=0) stop "88 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_force, STAT=status(1));deallocate (q_semi_y_force, STAT=status(2));deallocate (q_semi_z_force, STAT=status(3))
!if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "Cannot deallocate memory"
deallocate (diff_pres_x, STAT=status(1));deallocate (diff_pres_y, STAT=status(2));deallocate (diff_pres_z, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "89 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 (fc_prev_n, STAT=status(1));deallocate (fc_prev_n_old, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "90 Cannot deallocate memory"
deallocate (fd_prev_n, STAT=status(1)) ;deallocate (fc_prev_f, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "91 Cannot deallocate memory"
deallocate (fc_prev_f_old, STAT=status(1));deallocate (fd_prev_f, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "92Cannot deallocate memory"
!deallocate (fd2_prev_r_old, STAT=status(1));deallocate (fd2_prev_c_old, STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot deallocate memory"
deallocate (cp, STAT=status(1));deallocate (cv, STAT=status(2));deallocate (cw, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "93 Cannot deallocate memory"
deallocate (IIB_I_cp, STAT=status(1));deallocate (IIB_I_cv, STAT=status(2));deallocate (IIB_I_cw, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "94 Cannot deallocate memory"
deallocate (cp_xy, STAT=status(1));deallocate (cp_yz, STAT=status(2));deallocate (cp_zx, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "95 Cannot deallocate memory"
deallocate (cu_xy, STAT=status(1));deallocate (cu_yz, STAT=status(2));deallocate (cu_zx, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "96 Cannot deallocate memory"
deallocate (cv_xy, STAT=status(1));deallocate (cv_yz, STAT=status(2));deallocate (cv_zx, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "97 Cannot deallocate memory"
deallocate (cw_xy, STAT=status(1));deallocate (cw_yz, STAT=status(2));deallocate (cw_zx, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "98 Cannot deallocate memory"
if (no_body <= 3) then
deallocate (IIB_cell_p1, STAT=status(1));deallocate (IIB_cell_v1, STAT=status(2));deallocate (IIB_cell_w1, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "99 Cannot deallocate memory"
deallocate (I_cell_p1, STAT=status(1));deallocate (I_cell_v1, STAT=status(2));deallocate (I_cell_w1, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "100 Cannot deallocate memory"
deallocate (IIB_global_cell_p1, STAT=status(1));deallocate (IIB_global_cell_v1, STAT=status(2));deallocate (IIB_global_cell_w1, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "101Cannot deallocate memory"
deallocate (I_global_cell_p1, STAT=status(1));deallocate (I_global_cell_v1, STAT=status(2));deallocate (I_global_cell_w1, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "102 Cannot deallocate memory"
deallocate (IIB_equal_cell_p1, STAT=status(1));deallocate (IIB_equal_cell_v1, STAT=status(2));deallocate (IIB_equal_cell_w1, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "103 Cannot deallocate memory"
deallocate (I_equal_cell_p1, STAT=status(1));deallocate (I_equal_cell_v1, STAT=status(2));deallocate (I_equal_cell_w1, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "104 Cannot deallocate memory"
if (no_body >= 2) then
deallocate (IIB_cell_p2, STAT=status(1));deallocate (IIB_cell_v2, STAT=status(2));deallocate (IIB_cell_w2, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "105 Cannot deallocate memory"
deallocate (I_cell_p2, STAT=status(1));deallocate (I_cell_v2, STAT=status(2));deallocate (I_cell_w2, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "106 Cannot deallocate memory"
deallocate (IIB_global_cell_p2, STAT=status(1));deallocate (IIB_global_cell_v2, STAT=status(2));deallocate (IIB_global_cell_w2, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "107 Cannot deallocate memory"
deallocate (I_global_cell_p2, STAT=status(1));deallocate (I_global_cell_v2, STAT=status(2));deallocate (I_global_cell_w2, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "108 Cannot deallocate memory"
deallocate (IIB_equal_cell_p2, STAT=status(1));deallocate (IIB_equal_cell_v2, STAT=status(2));deallocate (IIB_equal_cell_w2, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "109 Cannot deallocate memory"
deallocate (I_equal_cell_p2, STAT=status(1));deallocate (I_equal_cell_v2, STAT=status(2));deallocate (I_equal_cell_w2, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "110 Cannot deallocate memory"
if (no_body == 3) then
deallocate (IIB_cell_p3, STAT=status(1));deallocate (IIB_cell_v3, STAT=status(2));deallocate (IIB_cell_w3, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "111 Cannot deallocate memory"
deallocate (I_cell_p3, STAT=status(1));deallocate (I_cell_v3, STAT=status(2));deallocate (I_cell_w3, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "112 Cannot deallocate memory"
deallocate (IIB_global_cell_p3, STAT=status(1));deallocate (IIB_global_cell_v3, STAT=status(2));deallocate (IIB_global_cell_w3, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "113 Cannot deallocate memory"
deallocate (I_global_cell_p3, STAT=status(1));deallocate (I_global_cell_v3, STAT=status(2));deallocate (I_global_cell_w3, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "114 Cannot deallocate memory"
deallocate (IIB_equal_cell_p3, STAT=status(1));deallocate (IIB_equal_cell_v3, STAT=status(2));deallocate (IIB_equal_cell_w3, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "115 Cannot deallocate memory"
deallocate (I_equal_cell_p3, STAT=status(1));deallocate (I_equal_cell_v3, STAT=status(2));deallocate (I_equal_cell_w3, STAT=status(3))
if (status(1)/=0 .or. status(2)/=0 .or. status(3)/=0) stop "116 Cannot deallocate memory"
end if
end if
!if (time>steady .and. field_ext==2) then
! deallocate (pIB_cell_u, STAT=status(1));deallocate (pIB_cell_v, STAT=status(2))
! if (status(1)/=0 .or. status(2)/=0) stop "Cannot deallocate memory"
! deallocate (pIB_cell_w, STAT=status(1));deallocate (pIB_cell_p, STAT=status(2))
! if (status(1)/=0 .or. status(2)/=0) stop "Cannot deallocate memory"
! deallocate (pIB_IB_vel, STAT=status(1));deallocate (pIB_IB_vel_old, STAT=status(2))
! if (status(1)/=0 .or. status(2)/=0) stop "Cannot deallocate memory"
!end if
if (IB_force_choice == 2 .or. vort_plot == 1) then
deallocate (IIB_IB_vel1, STAT=status(1));deallocate (IIB_IB_vel1_old, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "117 Cannot deallocate memory"
!deallocate (IIB_cell_p1_s_no, STAT=status(1));deallocate (IIB_cell_p1_v_no, STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot deallocate memory"
if (no_body >= 2) then
deallocate (IIB_IB_vel2, STAT=status(1));deallocate (IIB_IB_vel2_old, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "118 Cannot deallocate memory"
!deallocate (IIB_cell_p2_s_no, STAT=status(1));deallocate (IIB_cell_p2_v_no, STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot deallocate memory"
if (no_body == 3) then
deallocate (IIB_IB_vel3, STAT=status(1));deallocate (IIB_IB_vel3_old, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "119 Cannot deallocate memory"
!deallocate (IIB_cell_p3_s_no, STAT=status(1));deallocate (IIB_cell_p3_v_no, STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot deallocate memory"
end if
end if
end if
end if
!deallocate (ibm_force, STAT=status(2))
!if (status(2)/=0) stop "Cannot deallocate memory"
!deallocate (IIB_index_u, STAT=status(1));deallocate (IIB_index_v, STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "Cannot deallocate memory"
!deallocate (IIB_index_w, STAT=status(1));;deallocate (IIB_index_p, STAT=status(2))
!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)
!call PCDestroy(pc,ierr)
end if
if (allocated(Q_Soria)) deallocate (Q_Soria, STAT=status(1))
!if (.not. allocated(q_crit)) allocate (q_crit(size_x,size_y,size_z), STAT=status(2))
if (status(1)/=0) stop "120 Cannot deallocate memory"
if (allocated(Q_criterion)) deallocate (Q_criterion, STAT=status(1))
if (allocated(Q_Jeong)) deallocate (Q_Jeong, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "121 Cannot deallocate memory"
if (allocated(Q_criterion_avg)) deallocate (Q_criterion_avg, STAT=status(1))
if (allocated(Q_criterion_tmp)) deallocate (Q_criterion_tmp, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "122 Cannot deallocate memory"
if (allocated(vorticity_tmp)) deallocate (vorticity_tmp, STAT=status(1))
if (allocated(vorticity_avg)) deallocate (vorticity_avg, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "123 Cannot deallocate memory"
if (spatial_scheme == 3 .or. spatial_scheme == 6) then
deallocate (s_cf_u_e, STAT=status(1));deallocate (s_cf_v_e, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "124 Cannot deallocate memory"
deallocate (s_cf_u_n, STAT=status(1));deallocate (s_cf_v_n, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "125 Cannot deallocate memory"
deallocate (s_cf_w_e, STAT=status(1));deallocate (s_cf_w_n, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "126 Cannot deallocate memory"
deallocate (s_cf_u_f, STAT=status(1))
if (status(1)/=0) STOP "127 Cannot deallocate memory"
deallocate (s_cf_v_f, STAT=status(1));deallocate (s_cf_w_f, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "128 Cannot deallocate memory"
end if
call VecDestroy(xx_semi_xyz,ierr)
call VecDestroy(b_rhs_semi_xyz,ierr)
!call VecDestroy(velocity_global,ierr); call VecDestroy(velocity_local,ierr)
!call VecDestroy(b_rhs_semi_global,ierr); call VecDestroy(b_rhs_semi_local,ierr)
call MatDestroy(A_semi_xyz,ierr)
call KSPDestroy(ksp_semi_xyz,ierr)
!call PCDestroy(pc_semi_x,ierr)
!call PCDestroy(pc_semi_y,ierr)
!call PCDestroy(pc_semi_z,ierr)
!call PCDestroy(pc_semi_xyz,ierr)
!DM stuffs
call DMDestroy(da_v,ierr)
call VecDestroy(v_local,ierr)
call VecDestroy(v_global,ierr)
call DMDestroy(da_w,ierr)
call VecDestroy(w_local,ierr)
call VecDestroy(w_global,ierr)
call DMDestroy(da_p,ierr)
call VecDestroy(p_local,ierr)
call VecDestroy(p_global,ierr)
call DMDestroy(da_uvw_ast,ierr)
call DMDestroy(da_q_semi_vect_xyz,ierr)
call DMDestroy(da_q_p,ierr)
call VecDestroy(uvw_ast_local,ierr)
call VecDestroy(uvw_ast_global,ierr)
call VecDestroy(q_semi_vect_xyz_local,ierr)
call VecDestroy(q_semi_vect_xyz_global,ierr)
call VecDestroy(q_p_local,ierr)
call VecDestroy(q_p_global,ierr)
!call DMDestroy(da_uvw_old,ierr)
!call VecDestroy(uvw_old_local,ierr)
!call VecDestroy(uvw_old_global,ierr)
end if
if (motion_type == 13) then
deallocate (ff_syr_time_period, STAT=status(1));deallocate (ff_syr_theta_x_rad, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) STOP "129 Cannot deallocate memory"
deallocate (ff_syr_theta_z_rad, STAT=status(1))
if (status(1)/=0) STOP "130 Cannot deallocate memory"
end if
!dec$ define flag = 3
!dec$ if (flag < 2)
! print *, "abc"
!dec$ else if (flag < 1)
! print *, "abc"
!dec$ else
! print *, "abc"
!dec$ endif
end subroutine de_ini_var
subroutine de_ini_var_n_bodies1
!de-initialize/deallocate memory for variables of n bodies
!part 1
!#include "petsc/finclude/petsc.h"
PetscInt :: status(3)
deallocate (n_hy0, STAT=status(1));deallocate (n_freq, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "13 Cannot deallocate memory"
deallocate (n_phase_ang, STAT=status(1));deallocate (n_theta0_x, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "14 Cannot deallocate memory"
deallocate (n_theta0_y, STAT=status(1));deallocate (n_theta0_z, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "15 Cannot deallocate memory"
deallocate (n_loc_rot, STAT=status(1));deallocate (n_leading_edge_pos, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "16 Cannot deallocate memory"
deallocate (n_plate_h_tk, STAT=status(1));deallocate (n_plate_lena, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "17 Cannot deallocate memory"
deallocate (n_plate_lenb, STAT=status(1));deallocate (n_plate_wid, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "18 Cannot deallocate memory"
deallocate (n_shift_x, STAT=status(1));deallocate (z_1_to_n, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "19 Cannot deallocate memory"
deallocate (phase_ang_1_to_n, STAT=status(1));deallocate (n_hx0, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "11 Cannot deallocate memory"
deallocate (n_hz0, STAT=status(1));deallocate (n_ini_hx0, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "21 Cannot deallocate memory"
deallocate (n_ini_hy0, STAT=status(1));deallocate (n_ini_theta0_x, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "31 Cannot deallocate memory"
deallocate (n_ini_theta0_y, STAT=status(1));deallocate (n_ini_theta0_z, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "41 Cannot deallocate memory"
deallocate (n_body_cg_ini, STAT=status(1));deallocate (n_body_cg, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "51 Cannot deallocate memory"
deallocate (n_body_cg_old, STAT=status(1));deallocate (n_body_cg_old2, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "61 Cannot deallocate memory"
deallocate (n_rot_body, STAT=status(1));deallocate (n_rot_body_old, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "71 Cannot deallocate memory"
deallocate (n_fixed_rot_body, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "81 Cannot deallocate memory"
end subroutine de_ini_var_n_bodies1
subroutine de_ini_var_n_bodies2
!de-initialize/deallocate memory for variables of n bodies
!part 2
!#include "petsc/finclude/petsc.h"
PetscInt :: status(3)
if (inout_test <= 2) then
deallocate (n_IIB_cell_u, STAT=status(1))
if (status(1)/=0) stop "15 Cannot deallocate memory"
deallocate (n_I_cell_u, STAT=status(1))
if (status(1)/=0) stop "16 Cannot deallocate memory"
deallocate (n_IIB_global_cell_u, STAT=status(1))
if (status(1)/=0) stop "17 Cannot deallocate memory"
deallocate (n_I_global_cell_u, STAT=status(1))
if (status(1)/=0) stop "18 Cannot deallocate memory"
deallocate (n_IIB_equal_cell_u, STAT=status(1))
if (status(1)/=0) stop "19 Cannot deallocate memory"
deallocate (n_I_equal_cell_u, STAT=status(1))
if (status(1)/=0) stop "20 Cannot deallocate memory"
if (inout_test <= 1) then
deallocate (n_IIB_cell_v, STAT=status(1)); deallocate (n_IIB_cell_w, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "15 Cannot deallocate memory"
deallocate (n_I_cell_v, STAT=status(1)); deallocate (n_I_cell_w, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "16 Cannot deallocate memory"
deallocate (n_IIB_global_cell_v, STAT=status(1)); deallocate (n_IIB_global_cell_w, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "17 Cannot deallocate memory"
deallocate (n_I_global_cell_v, STAT=status(1)); deallocate (n_I_global_cell_w, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "18 Cannot deallocate memory"
deallocate (n_IIB_equal_cell_v, STAT=status(1)); deallocate (n_IIB_equal_cell_w, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "19 Cannot deallocate memory"
deallocate (n_I_equal_cell_v, STAT=status(1)); deallocate (n_I_equal_cell_w, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "20 Cannot deallocate memory"
end if
end if
deallocate (n_IIB_I_cell_no_uvw_global, STAT=status(1)); deallocate (n_IIB_I_cell_no_p_global, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "37 Cannot deallocate memory"
deallocate (n_IIB_equal_cell_no_global_u, STAT=status(1));deallocate (n_IIB_equal_cell_no_global_v, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "38 Cannot deallocate memory"
deallocate (n_IIB_equal_cell_no_global_w, STAT=status(1));deallocate (n_I_equal_cell_no_global_u, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "39 Cannot deallocate memory"
deallocate (n_I_equal_cell_no_global_v, STAT=status(1));deallocate (n_I_equal_cell_no_global_w, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "40 Cannot deallocate memory"
deallocate (n_IIB_equal_cell_no_global_p, STAT=status(1));deallocate (n_I_equal_cell_no_global_p, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "41 Cannot deallocate memory"
deallocate (n_body_sta_ijk, STAT=status(1));deallocate (n_body_end_ijk, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "41 Cannot deallocate memory"
deallocate (n_body_sta_ijk_ghost, STAT=status(1));deallocate (n_body_end_ijk_ghost, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "41 Cannot deallocate memory"
deallocate (n_body_sta_ijk_IIB, STAT=status(1));deallocate (n_body_end_ijk_IIB, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "41 Cannot deallocate memory"
deallocate (n_body_sta_ijk_IIB_ghost, STAT=status(1));deallocate (n_body_end_ijk_IIB_ghost, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "41 Cannot deallocate memory"
!deallocate (n_body_sta_ijk_IIB_ext, STAT=status(1));deallocate (n_body_end_ijk_IIB_ext, STAT=status(2))
!if (status(1)/=0 .or. status(2)/=0) stop "41 Cannot deallocate memory"
deallocate (n_body_sta_ijk_IIB_array, STAT=status(1));deallocate (n_body_end_ijk_IIB_array, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "41 Cannot deallocate memory"
if (allocated(n_surface)) then
deallocate (n_surface, STAT=status(1))
if (status(1) /= 0) then
! print *, "myid,n_surface deallocate success",myid
!else
print *, "myid,n_surface deallocate not success",myid
stop
end if
end if
if (allocated(n_vertex)) then
deallocate (n_vertex, STAT=status(1))
if (status(1) /= 0) then
! print *, "myid,n_vertex deallocate success",myid
!else
print *, "myid,n_vertex deallocate not success",myid
stop
end if
end if
if (allocated(n_vertex_ini)) then
deallocate (n_vertex_ini, STAT=status(1))
if (status(1)/=0) stop "43 Cannot deallocate memory"
end if
if (body_input_type /= 3) then
!deallocate (vertex1_ini, STAT=status(1))
deallocate (n_surface_pts, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "43 Cannot deallocate memory"
end if
if (body_input_type == 3) then
deallocate (n_vertex_multi, STAT=status(1))
if (status(1)/=0) stop "45 Cannot deallocate memory"
deallocate (n_fft_coeff_a, STAT=status(1)) ; deallocate (n_fft_coeff_b, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "49 Cannot deallocate memory"
end if
if (body_input_type == 0 .or. body_input_type == 3 .or. motion_type == 2 .or. deformation /= 0 .or. fsi_run >= 1) then
if (allocated(n_vertex_old)) then
deallocate (n_vertex_old, STAT=status(1))
if (status(1) /= 0) then
! print *, "myid,n_vertex_old deallocate success",myid
!else
print *, "myid,n_vertex_old deallocate not success",myid
stop
end if
end if
if (allocated(n_vertex_old2)) then
deallocate (n_vertex_old2, STAT=status(1))
if (status(1) /= 0) then
! print *, "myid,n_vertex_old2 deallocate success",myid
!else
print *, "myid,n_vertex_old2 deallocate not success",myid
stop
end if
end if
end if
deallocate (n_sta_yy, STAT=status(1)) ;deallocate (n_end_yy, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1 Cannot deallocate memory"
deallocate (n_surface_normal_correct, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "12 Cannot deallocate memory"
deallocate (n_pvm_xyz, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "81 Cannot deallocate memory"
deallocate (n_cd_cl_cs_mom_implicit, STAT=status(1));deallocate (n_cd_cl_cs_mom_explicit, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "91 Cannot deallocate memory"
deallocate (n_cd_cl_cs_mom_power, STAT=status(1));deallocate (n_cd_tmp, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1a Cannot deallocate memory"
deallocate (n_cl_tmp, STAT=status(1));deallocate (n_cs_tmp, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1s Cannot deallocate memory"
deallocate (n_body_vel, STAT=status(1));deallocate (n_body_vel_old, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1d Cannot deallocate memory"
deallocate (n_body_displacement, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "1f Cannot deallocate memory"
deallocate (n_IIB_cell_no_u, STAT=status(1)); deallocate (n_IIB_cell_no_v, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1g Cannot deallocate memory"
deallocate (n_IIB_cell_no_w, STAT=status(1));deallocate (n_IIB_cell_no_p, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1h Cannot deallocate memory"
deallocate (n_I_cell_no_u, STAT=status(1)); deallocate (n_I_cell_no_v, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1j Cannot deallocate memory"
deallocate (n_I_cell_no_w, STAT=status(1));deallocate (n_I_cell_no_p, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1k Cannot deallocate memory"
deallocate (n_IIB_global_cell_no_u, STAT=status(1));deallocate (n_IIB_global_cell_no_v, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1l Cannot deallocate memory"
deallocate (n_IIB_global_cell_no_w, STAT=status(1));deallocate (n_IIB_global_cell_no_p, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1; Cannot deallocate memory"
deallocate (n_I_global_cell_no_u, STAT=status(1));deallocate (n_I_global_cell_no_v, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1q Cannot deallocate memory"
deallocate (n_I_global_cell_no_w, STAT=status(1));deallocate (n_I_global_cell_no_p, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1w Cannot deallocate memory"
deallocate (n_IIB_equal_cell_no_u, STAT=status(1));deallocate (n_IIB_equal_cell_no_v, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1e Cannot deallocate memory"
deallocate (n_IIB_equal_cell_no_w, STAT=status(1));deallocate (n_IIB_equal_cell_no_p, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1r Cannot deallocate memory"
deallocate (n_I_equal_cell_no_u, STAT=status(1));deallocate (n_I_equal_cell_no_v, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1t Cannot deallocate memory"
deallocate (n_I_equal_cell_no_w, STAT=status(1));deallocate (n_I_equal_cell_no_p, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1t Cannot deallocate memory"
deallocate (n_IIB_I_cell_no_uvw_total, STAT=status(1));deallocate (n_IIB_I_cell_no_p_total, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1t Cannot deallocate memory"
deallocate (n_IIB_equal_cell_no_u_max, STAT=status(1));deallocate (n_I_equal_cell_no_u_max, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1t Cannot deallocate memory"
deallocate (n_IIB_cell_no_u_max, STAT=status(1));deallocate (n_I_cell_no_u_max, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1t Cannot deallocate memory"
deallocate (n_IIB_global_cell_no_u_max, STAT=status(1));deallocate (n_I_global_cell_no_u_max, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1t Cannot deallocate memory"
deallocate (n_sta_x, STAT=status(1));deallocate (n_end_x, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1y Cannot deallocate memory"
deallocate (n_sta_y, STAT=status(1));deallocate (n_end_y, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1u Cannot deallocate memory"
deallocate (n_sta_z, STAT=status(1));deallocate (n_end_z, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1i Cannot deallocate memory"
deallocate (n_no_vertices, STAT=status(1));deallocate (n_no_surfaces, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1o Cannot deallocate memory"
deallocate (n_no_body_pts, STAT=status(1))
if (status(1)/=0 .or. status(2)/=0) stop "1b Cannot deallocate memory"
deallocate (n_sta_cept_xyz, STAT=status(1));deallocate (n_end_cept_xyz, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1n Cannot deallocate memory"
deallocate (n_IIB_ista, STAT=status(1));deallocate (n_IIB_iend, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1m Cannot deallocate memory"
deallocate (n_IIB_jsta, STAT=status(1));deallocate (n_IIB_jend, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1, Cannot deallocate memory"
deallocate (n_IIB_ksta, STAT=status(1));deallocate (n_IIB_kend, STAT=status(2))
if (status(1)/=0 .or. status(2)/=0) stop "1> Cannot deallocate memory"
end subroutine de_ini_var_n_bodies2
subroutine mpi_para_range(n1, n2, nprocs, irank, ksta, kend)
!>block distribution
PetscInt n1 !The lowest value of the iteration variable (IN)
PetscInt n2 !The highest value of the iteration variable (IN)
PetscInt nprocs !The number of processes (IN)
PetscInt irank !The rank for which you want to know the range of iterations(IN)
PetscInt ksta !The lowest value of the iteration variable that process irank executes (OUT)
PetscInt kend !The highest value of the iteration variable that process irank executes (OUT)
PetscInt iwork1,iwork2
iwork1 = (n2 - n1 + 1) / nprocs
iwork2 = mod(n2 - n1 + 1, nprocs)
ksta = irank * iwork1 + n1 + min(irank, iwork2)
kend = ksta + iwork1 - 1
if (iwork2 > irank) kend = kend + 1
end subroutine mpi_para_range
end module global_data
!
!(c) Matthew Kennel, Institute for Nonlinear Science (2004)
!
! Licensed under the Academic Free License version 2.0
!
! Distributed as part of AxiSEM
! It is distributed from the webpage <http://www.axisem.info>
More information about the petsc-users
mailing list