| 1 | subroutine initialize_cfd |
|---|
| 2 | #define DEBUG_FLAG 3 |
|---|
| 3 | |
|---|
| 4 | use grid_info |
|---|
| 5 | use mpi_parameters |
|---|
| 6 | use integration_weights |
|---|
| 7 | use dimensional_scales |
|---|
| 8 | use user_parameters |
|---|
| 9 | use pde_parameters |
|---|
| 10 | use dependent_variables |
|---|
| 11 | use intermediate_variables |
|---|
| 12 | use rhs_variables |
|---|
| 13 | use boundary_information |
|---|
| 14 | use counters_flags_etc |
|---|
| 15 | use stress_parameters |
|---|
| 16 | |
|---|
| 17 | integer :: rank,flag |
|---|
| 18 | real :: end_values(2),gamma,drdz |
|---|
| 19 | real :: dy |
|---|
| 20 | |
|---|
| 21 | character(len=3) :: cid |
|---|
| 22 | character(len=80) :: filename |
|---|
| 23 | |
|---|
| 24 | !!! real :: maxu,minu,maxv,minv |
|---|
| 25 | |
|---|
| 26 | include '../input/problem_size.h' |
|---|
| 27 | include 'mpif.h' |
|---|
| 28 | include '../input/params.h' |
|---|
| 29 | |
|---|
| 30 | !module integration_weights |
|---|
| 31 | allocate ( W1(nx), W2(ny), W3(locnz) ) |
|---|
| 32 | |
|---|
| 33 | !module dependent_variables |
|---|
| 34 | allocate ( U(nx,ny,locnz),v(nx,ny,locnz) ) |
|---|
| 35 | allocate ( W(nx,ny,locnz),pd(nx,ny,locnz) ) |
|---|
| 36 | allocate ( pressure(nx,ny,locnz,2) ) |
|---|
| 37 | allocate ( s1(nx,ny,locnz),s2(nx,ny,locnz) ) |
|---|
| 38 | allocate ( s1_bar(nx,locnz) ) ! not a dep. variable, but goes w/ s1 |
|---|
| 39 | allocate ( s1_bar_xi(nx,locnz) ) |
|---|
| 40 | allocate ( s1_bar_zeta(nx,locnz) ) |
|---|
| 41 | allocate ( s2_bar(nx,locnz) ) ! not a dep. variable, but goes w/ s2 |
|---|
| 42 | allocate ( s2_bar_xi(nx,locnz) ) |
|---|
| 43 | allocate ( s2_bar_zeta(nx,locnz) ) |
|---|
| 44 | allocate ( rho_bar(nx,locnz) ) ! not a dep. variable, but goes w/ pd |
|---|
| 45 | |
|---|
| 46 | |
|---|
| 47 | !module intermediate_variables |
|---|
| 48 | allocate ( div_u(nx,ny,locnz,2),b(nx,ny,locnz) ) |
|---|
| 49 | allocate ( L_of_p(nx,ny,locnz) ) |
|---|
| 50 | allocate ( grad_x(nx,ny,locnz) ) |
|---|
| 51 | allocate ( grad_y(nx,ny,locnz) ) |
|---|
| 52 | allocate ( grad_z(nx,ny,locnz) ) |
|---|
| 53 | allocate ( u_star(nx,ny,locnz) ) |
|---|
| 54 | allocate ( v_star(nx,ny,locnz) ) |
|---|
| 55 | allocate ( w_star(nx,ny,locnz) ) |
|---|
| 56 | allocate ( u_star_old(ny,locnz,2) ) |
|---|
| 57 | allocate ( v_star_old(nx,locnz,2) ) |
|---|
| 58 | allocate ( w_star_old(nx,ny,2) ) |
|---|
| 59 | |
|---|
| 60 | !module rhs_variables |
|---|
| 61 | allocate ( rhs1(nx,ny,locnz,AB_order) ) |
|---|
| 62 | allocate ( rhs2(nx,ny,locnz,AB_order) ) |
|---|
| 63 | allocate ( rhs3(nx,ny,locnz,AB_order) ) |
|---|
| 64 | allocate ( rhs4(nx,ny,locnz,AB_order) ) |
|---|
| 65 | allocate ( rhs5(nx,ny,locnz,AB_order) ) |
|---|
| 66 | |
|---|
| 67 | !module pde_parameters |
|---|
| 68 | allocate ( K_m(nx,ny,locnz) ) |
|---|
| 69 | |
|---|
| 70 | ! module boundary_information |
|---|
| 71 | !!! allocate ( boundary_type1(ny,locnz),boundary_type2(ny,locnz) ) |
|---|
| 72 | !!! allocate ( boundary_type3(nx,locnz),boundary_type4(nx,locnz) ) |
|---|
| 73 | !!! allocate ( boundary_type5(nx,ny),boundary_type6(nx,ny) ) |
|---|
| 74 | !!! allocate ( BC_pressure1(ny,locnz),BC_pressure2(ny,locnz) ) |
|---|
| 75 | !!! allocate ( BC_pressure3(nx,locnz),BC_pressure4(nx,locnz) ) |
|---|
| 76 | !!! allocate ( BC_pressure5(nx,ny),BC_pressure6(nx,ny) ) |
|---|
| 77 | !!! allocate ( BC_vel1(ny,locnz),BC_vel2(ny,locnz) ) |
|---|
| 78 | !!! allocate ( BC_vel3(nx,locnz),BC_vel4(nx,locnz) ) |
|---|
| 79 | !!! allocate ( BC_vel5(nx,ny),BC_vel6(nx,ny) ) |
|---|
| 80 | !!! allocate ( BC_scalar1(ny,locnz,num_scalars),BC_scalar2(ny,locnz,num_scalars) ) |
|---|
| 81 | !!! allocate ( BC_scalar3(nx,locnz,num_scalars),BC_scalar4(nx,locnz,num_scalars) ) |
|---|
| 82 | !!! allocate ( BC_scalar5(nx,ny,num_scalars),BC_scalar6(nx,ny,num_scalars) ) |
|---|
| 83 | allocate ( BC_values1(ny,locnz,3+num_scalars),BC_values2(ny,locnz,3+num_scalars) ) |
|---|
| 84 | allocate ( BC_values3(nx,locnz,3+num_scalars),BC_values4(nx,locnz,3+num_scalars) ) |
|---|
| 85 | allocate ( BC_values5(nx,ny,3+num_scalars),BC_values6(nx,ny,3+num_scalars) ) |
|---|
| 86 | allocate ( phiY(nx,locnz,6,3,5) ) |
|---|
| 87 | allocate ( inflow(nx,locnz,3+num_scalars) ) |
|---|
| 88 | |
|---|
| 89 | |
|---|
| 90 | !******************************************************************* |
|---|
| 91 | write(0,*) 'IN INITIALIZE CFD, DEBUG_FLAG, myid = ', DEBUG_FLAG,myid |
|---|
| 92 | !!!#if DEBUG_FLAG >= 1 |
|---|
| 93 | !!! if(myid .eq. 0 ) then |
|---|
| 94 | !!! write(0,*) ' ' |
|---|
| 95 | write(0,*) 'hello world from initialize_cfd, myid =', myid |
|---|
| 96 | !!! endif |
|---|
| 97 | !!!#endif |
|---|
| 98 | |
|---|
| 99 | ! Check sizing parameters for obvious problems. |
|---|
| 100 | if( mod(nz,numprocs) /= 0 ) then |
|---|
| 101 | write(0,*) 'nz must be divisible by numprocs' |
|---|
| 102 | STOP |
|---|
| 103 | endif |
|---|
| 104 | if( mod(nx-1,numprocs) /= 0 ) then |
|---|
| 105 | write(0,*) 'nx-1 must be divisible by numprocs for xyfft logic' |
|---|
| 106 | STOP |
|---|
| 107 | endif |
|---|
| 108 | if( locnz < 4 ) then |
|---|
| 109 | write(0,*) 'currently, locnz >= 4' |
|---|
| 110 | write(0,*) 'WILL TRY TO CARRY ON, BE CAREFUL!!!' |
|---|
| 111 | endif |
|---|
| 112 | if( AB_order /= 2 .and. AB_order /= 3 .and. AB_order /= 4) then |
|---|
| 113 | write(0,*) 'only 2nd, 3rd or 4th order Adams Bashforth time stepping is allowed' |
|---|
| 114 | STOP |
|---|
| 115 | endif |
|---|
| 116 | |
|---|
| 117 | init_derived_type=0 ! MPI derived data types not yet defined |
|---|
| 118 | subslice=-9999 ! MPI derived data types not yet defined, dummy value |
|---|
| 119 | |
|---|
| 120 | ! Each processor imports the main parameter file: |
|---|
| 121 | !!! xlv definition de pi dans params.h |
|---|
| 122 | !!! pi=acos(-1.0) |
|---|
| 123 | |
|---|
| 124 | ! set some sizing parameters that depend on case xy periodic |
|---|
| 125 | do i=0,Grid(0)%n_mg_levels-1 |
|---|
| 126 | |
|---|
| 127 | if( boundary_type2(1,1)=='periodic' .and. boundary_type2(1,1)=='periodic' ) then |
|---|
| 128 | Grid(i)%n1_it=1 |
|---|
| 129 | Grid(i)%n2_it=1 |
|---|
| 130 | Grid(i)%n3_it=Grid(i)%n3 ! sizes for U&V arrays in matvec |
|---|
| 131 | else |
|---|
| 132 | Grid(i)%n1_it=Grid(i)%n1 |
|---|
| 133 | Grid(i)%n2_it=Grid(i)%n2 |
|---|
| 134 | Grid(i)%n3_it=Grid(i)%locn3 ! sizes for U&V arrays in matvec |
|---|
| 135 | endif |
|---|
| 136 | |
|---|
| 137 | ! setup the PIM data for iterative solves at each grid level |
|---|
| 138 | if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then |
|---|
| 139 | lda = Grid(i)%n3 |
|---|
| 140 | loclen = lda ! iterations local to each processor |
|---|
| 141 | pret = 0 ! preconditioning type 1=left, 0 =none |
|---|
| 142 | stopt = 4 ! stop criteria (4 fastest, GMRES has its own) |
|---|
| 143 | tol = 1.0E-19 ! stopping tolerance eps |
|---|
| 144 | basis=maxits ! " |
|---|
| 145 | else ! iterative method for 3d pressur field |
|---|
| 146 | lda = Grid(i)%n1 * Grid(i)%n2 * Grid(i)%n3 |
|---|
| 147 | loclen = LDA/numprocs |
|---|
| 148 | pret = 0 ! preconditioning type 1=left, 0 =none |
|---|
| 149 | stopt = 4 ! stop criteria (4 fastest, GMRES has its own) |
|---|
| 150 | tol = 1.0E-19 ! stopping tolerance eps |
|---|
| 151 | basis=maxits ! " |
|---|
| 152 | endif |
|---|
| 153 | ! REMOVING BECAUSE WE WILL USE SOLVE_FOR_XYPER_PRESSURE INSTEAD |
|---|
| 154 | ! CALL PIMSSETPAR(Grid(i)%ipar,Grid(i)%spar, & |
|---|
| 155 | ! lda,lda,loclen,loclen,basis,numprocs,myid,pret,stopt,maxits,tol) |
|---|
| 156 | |
|---|
| 157 | enddo |
|---|
| 158 | |
|---|
| 159 | |
|---|
| 160 | !!! Domaine |
|---|
| 161 | call MPI_ALLREDUCE(Grid(0)%x(nx,1)*L_scale ,Lx,1,MPI_REAL,MPI_MAX,comm,ierr) |
|---|
| 162 | call MPI_ALLREDUCE(Grid(0)%y(ny)*L_scale ,Ly,1,MPI_REAL,MPI_MAX,comm,ierr) |
|---|
| 163 | call MPI_ALLREDUCE(Grid(0)%z(1,locnz)*L_scale,Lz,1,MPI_REAL,MPI_MAX,comm,ierr) |
|---|
| 164 | print*,'Domain : Lx,Ly,Lz=',Lx,Ly,Lz |
|---|
| 165 | |
|---|
| 166 | |
|---|
| 167 | #if DEBUG_FLAG >= 1 |
|---|
| 168 | if(myid .eq. 0 ) write(0,*) ' initialize_cfd: user parameter file ../input/params.h included ' |
|---|
| 169 | #endif |
|---|
| 170 | |
|---|
| 171 | write(0,*) 'derive inside initialize_cfd, nu,kappa1,kappa2,myid =', nu,kappa1,kappa2,myid |
|---|
| 172 | ! Derive dimensionless parameters given user input: |
|---|
| 173 | Re = U_scale*L_scale/nu ! bulk Reynolds number |
|---|
| 174 | Pr1 = nu/kappa1 ! "Prandtl" number for scalar 1 |
|---|
| 175 | Pr2 = nu/kappa2 ! "Prandtl" number for scalar 2 |
|---|
| 176 | if( closure_flag .eq. 'LES' ) then |
|---|
| 177 | C_s = 0.17 ! see Kaltenbach et el JFM 94, or |
|---|
| 178 | ! Lesieur & Metais Ann Rev of Fluid Mech vol28 |
|---|
| 179 | Pr = 1.0 |
|---|
| 180 | endif |
|---|
| 181 | Ri = (Nsquared * L_scale**2)/U_scale**2 ! bulk Richardson number |
|---|
| 182 | Rot = f*L_scale/U_scale ! bulk inverse Rossby number |
|---|
| 183 | beta_tilde = beta_dimensional * L_scale**2/U_scale |
|---|
| 184 | allocate ( beta_window(ny) ) |
|---|
| 185 | dy = Grid(0)%y(2)*L_scale - Grid(0)%y(1)*L_scale |
|---|
| 186 | beta_window(:) = beta_tilde * (Grid(0)%y(:)-Ly/(2*L_scale)) * tanh ( (Grid(0)%y(:)*L_scale-5*dy)/(Ly/100) ) |
|---|
| 187 | |
|---|
| 188 | #if DEBUG_FLAG >= 1 |
|---|
| 189 | if(myid .eq. 0 ) write(0,*) ' Re, Pr1, Pr2, Ri, Rot, beta_tilde : ', Re, Pr1, Pr2, Ri, Rot, beta_tilde |
|---|
| 190 | #endif |
|---|
| 191 | |
|---|
| 192 | !!! Enregistrement de f_of_y pour chaque processeur (pour vérification) |
|---|
| 193 | write(unit=cid,fmt=501) myid |
|---|
| 194 | 501 format(I3.3) |
|---|
| 195 | filename='../output/beta_effect_'//cid |
|---|
| 196 | open(905,file=filename,status='unknown') |
|---|
| 197 | do j=1,ny |
|---|
| 198 | write(905,*) Grid(0)%y(j)*L_scale, ( Rot + beta_window(j) ) / (L_scale/U_scale) |
|---|
| 199 | enddo |
|---|
| 200 | close(905) |
|---|
| 201 | |
|---|
| 202 | |
|---|
| 203 | if( restart_flag == 'no' ) then |
|---|
| 204 | |
|---|
| 205 | ! Set the ambient stratification profiles s1_bar(x,z), s2_bar(x,z), d'less values returned |
|---|
| 206 | call set_ambient_stratification |
|---|
| 207 | ! set rho_bar(x,z) to -999.9 (will be overwritten in 1st call to equation_of_state) |
|---|
| 208 | rho_bar(:,:)=-999.9 |
|---|
| 209 | ! Each processor initializes its portion of the dep. variable arrays |
|---|
| 210 | if( ics_flag == 'homogeneous' ) then |
|---|
| 211 | U=0.0 |
|---|
| 212 | v=0.0 |
|---|
| 213 | W=0.0 |
|---|
| 214 | s1=0.0 |
|---|
| 215 | if( num_scalars > 1 ) s2=0.0 |
|---|
| 216 | pressure=0.0 |
|---|
| 217 | |
|---|
| 218 | else |
|---|
| 219 | |
|---|
| 220 | ! read in ics data |
|---|
| 221 | |
|---|
| 222 | endif |
|---|
| 223 | |
|---|
| 224 | ! call equation_of_state to compute d'less pd and rho_bar |
|---|
| 225 | call equation_of_state |
|---|
| 226 | |
|---|
| 227 | print*,'[rhobar]=',minval(rho_bar),maxval(rho_bar) |
|---|
| 228 | |
|---|
| 229 | |
|---|
| 230 | ! Compute xi and zeta derivatives of s1_bar and s2_bar |
|---|
| 231 | flag=1;end_values=-1e32 |
|---|
| 232 | |
|---|
| 233 | #if DEBUG_FLAG >= 1 |
|---|
| 234 | if(myid .eq. 0 ) write(0,*) ' initialize_cfd: about to differentiate s1_bar & s2_bar ' |
|---|
| 235 | #endif |
|---|
| 236 | |
|---|
| 237 | |
|---|
| 238 | call compact_ddx(s1_bar,s1_bar_xi,flag,end_values,nx,1,locnz) |
|---|
| 239 | if(num_scalars>1) call compact_ddx(s2_bar,s2_bar_xi,flag,end_values,nx,1,locnz) |
|---|
| 240 | |
|---|
| 241 | #if DEBUG_FLAG >= 1 |
|---|
| 242 | if(myid .eq. 0 ) write(0,*) ' initialize_cfd: return from ddx ' |
|---|
| 243 | #endif |
|---|
| 244 | |
|---|
| 245 | call compact_ddz_mpi(s1_bar,s1_bar_zeta,flag,end_values,nx,1,locnz) |
|---|
| 246 | if(num_scalars>1) call compact_ddz_mpi(s2_bar,s2_bar_zeta,flag,end_values,nx,1,locnz) |
|---|
| 247 | |
|---|
| 248 | print*,'minval(s1_bar_zeta),maxval(s1_bar_zeta)=',minval(s1_bar_zeta),maxval(s1_bar_zeta) |
|---|
| 249 | print*,'minval(s2_bar_zeta),maxval(s2_bar_zeta)=',minval(s2_bar_zeta),maxval(s2_bar_zeta) |
|---|
| 250 | |
|---|
| 251 | |
|---|
| 252 | #if DEBUG_FLAG >= 1 |
|---|
| 253 | if(myid .eq. 0 ) write(0,*) ' initialize_cfd: return from ddz ' |
|---|
| 254 | #endif |
|---|
| 255 | |
|---|
| 256 | |
|---|
| 257 | init_derived_type=0 ! MPI derived data types not yet defined |
|---|
| 258 | subslice=-9999 ! MPI derived data types not yet defined, dummy value |
|---|
| 259 | |
|---|
| 260 | #if DEBUG_FLAG >= 1 |
|---|
| 261 | if(myid .eq. 0 ) write(0,*) ' initialize_cfd: differentiated s1/s2_bar ' |
|---|
| 262 | #endif |
|---|
| 263 | |
|---|
| 264 | |
|---|
| 265 | endif ! end restart_flag == 'no' block |
|---|
| 266 | |
|---|
| 267 | ! Each processor fills weight arrays for spatial integration dxi,deta&dzeta: |
|---|
| 268 | call set_integration_weights(Grid(0)%dxi,Grid(0)%deta,Grid(0)%dzeta) |
|---|
| 269 | |
|---|
| 270 | |
|---|
| 271 | ! Set integer BC flags |
|---|
| 272 | BC_flag(:,:)=1 ! default=no end conditions specified |
|---|
| 273 | if( boundary_type1(1,1) == 'periodic' ) BC_flag(:,1)=5 |
|---|
| 274 | if( boundary_type3(1,1) == 'periodic' ) BC_flag(:,2)=5 |
|---|
| 275 | if( boundary_type5(1,1) == 'periodic' ) BC_flag(:,3)=5 |
|---|
| 276 | |
|---|
| 277 | ! Initial values of counters etc |
|---|
| 278 | if( restart_flag == 'no' ) then |
|---|
| 279 | istart=0 |
|---|
| 280 | iend=istart+nsteps |
|---|
| 281 | new=1;old=2 |
|---|
| 282 | MM0=1;MM1=2;MM2=3;MM3=4 |
|---|
| 283 | step_flag='EULER' |
|---|
| 284 | write3D_counter=1 |
|---|
| 285 | |
|---|
| 286 | ! Intermediate variables initial data values |
|---|
| 287 | b=0.0 |
|---|
| 288 | div_u=0.0 |
|---|
| 289 | u_star=0.0 |
|---|
| 290 | v_star=0.0 |
|---|
| 291 | w_star=0.0 |
|---|
| 292 | u_star_old=0.0 |
|---|
| 293 | v_star_old=0.0 |
|---|
| 294 | w_star_old=0.0 |
|---|
| 295 | grad_x=0.0 |
|---|
| 296 | grad_y=0.0 |
|---|
| 297 | grad_z=0.0 |
|---|
| 298 | phiY=0.0 |
|---|
| 299 | elseif( restart_flag == 'yes' ) then |
|---|
| 300 | call restart |
|---|
| 301 | !!! maxu = MAXVAL(u(:,:,:)) |
|---|
| 302 | !!! minu = MINVAL(u(:,:,:)) |
|---|
| 303 | !!! maxv = MAXVAL(v(:,:,:)) |
|---|
| 304 | !!! minv = MINVAL(v(:,:,:)) |
|---|
| 305 | !!! write(0,*) 'MAXMIN',myid,maxu,minu,maxv,minv |
|---|
| 306 | endif |
|---|
| 307 | |
|---|
| 308 | |
|---|
| 309 | ! Dimensionless eddy viscosity |
|---|
| 310 | K_m(:,:,:)=1./Re |
|---|
| 311 | |
|---|
| 312 | #if DEBUG_FLAG >= 1 |
|---|
| 313 | if(myid .eq. 0 ) write(0,*) ' initialize_cfd: calling initialize immersed bdry' |
|---|
| 314 | #endif |
|---|
| 315 | |
|---|
| 316 | call initialize_immersed_bdry |
|---|
| 317 | |
|---|
| 318 | end subroutine initialize_cfd |
|---|
