[4] | 1 | !!==================================================================================== |
---|
| 2 | module particles |
---|
| 3 | !!==================================================================================== |
---|
| 4 | logical :: user_defines_particles=.TRUE. |
---|
| 5 | character(len=80) :: particle_step_flag |
---|
| 6 | |
---|
| 7 | real(kind=8),allocatable :: positions(:,:) |
---|
| 8 | real(kind=8),allocatable :: uvels (:,:) |
---|
| 9 | real(kind=8),allocatable :: vvels (:,:) |
---|
| 10 | real(kind=8),allocatable :: wvels (:,:) |
---|
| 11 | |
---|
| 12 | integer :: nparticles |
---|
| 13 | integer :: j_particle_time |
---|
| 14 | integer :: istart_trajectories=0 |
---|
| 15 | |
---|
| 16 | !!==================================================================================== |
---|
| 17 | end module particles |
---|
| 18 | !!==================================================================================== |
---|
| 19 | |
---|
| 20 | |
---|
| 21 | subroutine InitializeParticles |
---|
| 22 | use user_routines |
---|
| 23 | use particles |
---|
| 24 | use user_data |
---|
| 25 | use mpi_params |
---|
| 26 | use etc, only: docfile |
---|
| 27 | |
---|
| 28 | implicit none |
---|
| 29 | integer :: i, j, npts |
---|
| 30 | real(kind=8) :: urv |
---|
| 31 | real(kind=8) :: sides(3) |
---|
| 32 | |
---|
| 33 | call random_seed !! I'm using the f90 intrinsic functions |
---|
| 34 | !! random_seed, random_number here |
---|
| 35 | j_particle_time = 1 |
---|
| 36 | npts = nparticles |
---|
| 37 | |
---|
| 38 | allocate( positions(npts,3) ) |
---|
| 39 | allocate( uvels(npts,4) ) !! 4 time levels |
---|
| 40 | allocate( vvels(npts,4) ) |
---|
| 41 | allocate( wvels(npts,4) ) |
---|
| 42 | |
---|
| 43 | sides(:)=(/Lx,Ly,Lz/) |
---|
| 44 | user_defines_particles = .TRUE. |
---|
| 45 | if (user_defines_particles) then |
---|
| 46 | |
---|
| 47 | call user_particle_positions(positions,npts,Lx,Ly,Lz) |
---|
| 48 | |
---|
| 49 | else !! I'll distribute them uniformly |
---|
| 50 | |
---|
| 51 | do i=1,npts |
---|
| 52 | do j=1,3 |
---|
| 53 | call RANDOM_NUMBER( urv ) |
---|
| 54 | positions(i,j) = sides(j)*urv |
---|
| 55 | enddo |
---|
| 56 | enddo |
---|
| 57 | |
---|
| 58 | endif |
---|
| 59 | |
---|
| 60 | positions=positions/Lz |
---|
| 61 | uvels(:,:) = 0.d0 |
---|
| 62 | vvels(:,:) = 0.d0 |
---|
| 63 | wvels(:,:) = 0.d0 |
---|
| 64 | |
---|
| 65 | if(myid==0) then |
---|
| 66 | open(1,file=docfile,position='append') |
---|
| 67 | write(1,*) ' -----> InitializeParticles routine exiting normally <---------- ' |
---|
| 68 | write(0,*) ' -----> InitializeParticles routine exiting normally <---------- ' |
---|
| 69 | close(1) |
---|
| 70 | endif |
---|
| 71 | |
---|
| 72 | end subroutine InitializeParticles |
---|
| 73 | !============================================================================ |
---|
| 74 | !============================================================================ |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | |
---|
| 78 | |
---|
| 79 | !=============================================================== |
---|
| 80 | !=============================================================== |
---|
| 81 | subroutine advance_particles |
---|
| 82 | use etc, only: PM0,PM1,PM2,PM3, & |
---|
| 83 | my_z_limits,t_secs, & |
---|
| 84 | numzplanes_below,M_oldest |
---|
| 85 | use dependent_variables, only: velocity |
---|
| 86 | use intermediate_variables, only: explicit_rhs |
---|
| 87 | use user_data, only: dt,U0,Lx,Ly,Lz, & |
---|
| 88 | xdim_periodic, & |
---|
| 89 | ydim_periodic, & |
---|
| 90 | zdim_periodic,pi,ny |
---|
| 91 | use independent_variables |
---|
| 92 | use particles |
---|
| 93 | use mpi_params |
---|
| 94 | implicit none |
---|
| 95 | include 'mpif.h' |
---|
| 96 | |
---|
| 97 | integer :: i,j,k,kk,ngrdpts,ierr |
---|
| 98 | real(kind=8) :: xx,yy,zz,ans1,ans2,ans3 |
---|
| 99 | integer, SAVE :: m(3) |
---|
| 100 | integer, allocatable,SAVE :: my_labels(:) |
---|
| 101 | real(kind=8),allocatable,SAVE :: tmp(:) |
---|
| 102 | real(kind=8), save :: dt2,dt12,dt24 |
---|
| 103 | logical, save :: first_entry = .TRUE. |
---|
| 104 | |
---|
| 105 | if(nparticles<=0) return |
---|
| 106 | |
---|
| 107 | !!================================================================= |
---|
| 108 | !!================================================================= |
---|
| 109 | if( first_entry ) then |
---|
| 110 | dt2=dt/2.d0 |
---|
| 111 | dt12=dt/12.d0 |
---|
| 112 | dt24=dt/24.d0 |
---|
| 113 | m(:)=0 !! interpolate function vals, not derivs |
---|
| 114 | allocate( my_labels(nparticles),tmp(nparticles) ) |
---|
| 115 | first_entry = .FALSE. |
---|
| 116 | endif |
---|
| 117 | !!================================================================= |
---|
| 118 | !!================================================================= |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | !!================================================================= |
---|
| 123 | !!================================================================= |
---|
| 124 | !! periodic wrapping etc |
---|
| 125 | !!================================================================= |
---|
| 126 | !!================================================================= |
---|
| 127 | if(xdim_periodic) then |
---|
| 128 | do i=1,nparticles |
---|
| 129 | if( positions(i,1) < 0 ) then |
---|
| 130 | positions(i,1) = positions(i,1) + Lx/Lz |
---|
| 131 | elseif( positions(i,1) > Lx/Lz ) then |
---|
| 132 | positions(i,1) = positions(i,1) - Lx/Lz |
---|
| 133 | endif |
---|
| 134 | enddo |
---|
| 135 | endif |
---|
| 136 | |
---|
| 137 | if(ydim_periodic .and. ny>1) then |
---|
| 138 | do i=1,nparticles |
---|
| 139 | if( positions(i,2) < 0 ) then |
---|
| 140 | positions(i,2) = positions(i,2) + Ly/Lz |
---|
| 141 | elseif( positions(i,2) > Ly/Lz ) then |
---|
| 142 | positions(i,2) = positions(i,2) - Ly/Lz |
---|
| 143 | endif |
---|
| 144 | enddo |
---|
| 145 | endif |
---|
| 146 | |
---|
| 147 | if(zdim_periodic) then |
---|
| 148 | do i=1,nparticles |
---|
| 149 | if( positions(i,3) < 0 ) then |
---|
| 150 | positions(i,3) = positions(i,3) + Lz/Lz |
---|
| 151 | elseif( positions(i,3) > Lz/Lz ) then |
---|
| 152 | positions(i,3) = positions(i,3) - Lz/Lz |
---|
| 153 | endif |
---|
| 154 | enddo |
---|
| 155 | endif |
---|
| 156 | !!================================================================= |
---|
| 157 | !!================================================================= |
---|
| 158 | |
---|
| 159 | |
---|
| 160 | !!================================================================= |
---|
| 161 | !! mark the indices of particles residing on processor myid |
---|
| 162 | !!================================================================= |
---|
| 163 | |
---|
| 164 | do i=1,nparticles |
---|
| 165 | if( positions(i,3) >= my_z_limits(1) .and. & |
---|
| 166 | positions(i,3) < my_z_limits(2) ) then |
---|
| 167 | |
---|
| 168 | my_labels(i) = 1 |
---|
| 169 | |
---|
| 170 | else |
---|
| 171 | |
---|
| 172 | my_labels(i) = 0 |
---|
| 173 | positions(i,1:3) = 0.d0 |
---|
| 174 | |
---|
| 175 | endif |
---|
| 176 | enddo |
---|
| 177 | !!================================================================= |
---|
| 178 | !!================================================================= |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | !!================================================================= |
---|
| 182 | !!================================================================= |
---|
| 183 | !! get the u,v,w values at the positions of the particles on myid |
---|
| 184 | !! arrays are dimensioned like |
---|
| 185 | !! positions(npts,3) !! 3 space dims |
---|
| 186 | !! uvels(npts,4) !! 4 time levels |
---|
| 187 | !!================================================================= |
---|
| 188 | !!================================================================= |
---|
| 189 | |
---|
| 190 | !! get interpolated u values ( explicit_rhs are scratch arrays ) |
---|
| 191 | j=M_oldest !! oldest of the MM0,MM1,MM2,MM3 set, can be trashed |
---|
| 192 | call spline_interp_value(positions(1,1), & |
---|
| 193 | positions(1,2), & |
---|
| 194 | positions(1,3), & |
---|
| 195 | nparticles, & |
---|
| 196 | velocity(1), & |
---|
| 197 | m(1), & |
---|
| 198 | uvels(1,PM0), & |
---|
| 199 | my_labels, & |
---|
| 200 | explicit_rhs(1,j)%SlabVals, & |
---|
| 201 | explicit_rhs(2,j)%SlabVals, & |
---|
| 202 | explicit_rhs(3,j)%SlabVals ) |
---|
| 203 | |
---|
| 204 | !! get interpolated v values |
---|
| 205 | call spline_interp_value(positions(1,1), & |
---|
| 206 | positions(1,2), & |
---|
| 207 | positions(1,3), & |
---|
| 208 | nparticles, & |
---|
| 209 | velocity(2), & |
---|
| 210 | m(1), & |
---|
| 211 | vvels(1,PM0), & |
---|
| 212 | my_labels, & |
---|
| 213 | explicit_rhs(1,j)%SlabVals, & |
---|
| 214 | explicit_rhs(2,j)%SlabVals, & |
---|
| 215 | explicit_rhs(3,j)%SlabVals ) |
---|
| 216 | |
---|
| 217 | !! get interpolated w values |
---|
| 218 | call spline_interp_value(positions(1,1), & |
---|
| 219 | positions(1,2), & |
---|
| 220 | positions(1,3), & |
---|
| 221 | nparticles, & |
---|
| 222 | velocity(3), & |
---|
| 223 | m(1), & |
---|
| 224 | wvels(1,PM0), & |
---|
| 225 | my_labels, & |
---|
| 226 | explicit_rhs(1,j)%SlabVals, & |
---|
| 227 | explicit_rhs(2,j)%SlabVals, & |
---|
| 228 | explicit_rhs(3,j)%SlabVals ) |
---|
| 229 | |
---|
| 230 | |
---|
| 231 | !!================================================================= |
---|
| 232 | !!================================================================= |
---|
| 233 | !! At this point positions(i,:) are nonzero only if |
---|
| 234 | !! particle i resides on processor myid. Similarly |
---|
| 235 | !! for uvels(i,PM0),vvels(i,PM0) etc. Values for |
---|
| 236 | !! PM1, PM2 etc should be identical for all processors. |
---|
| 237 | !!================================================================= |
---|
| 238 | !!================================================================= |
---|
| 239 | |
---|
| 240 | !! Integrate forward one time step. |
---|
| 241 | do i=1,nparticles |
---|
| 242 | |
---|
| 243 | if( my_labels(i)==1 ) then !! only deal with those on proc. myid |
---|
| 244 | |
---|
| 245 | if( particle_step_flag .eq. 'euler' ) then |
---|
| 246 | positions(i,1)=positions(i,1) + dt*uvels(i,PM0) |
---|
| 247 | positions(i,2)=positions(i,2) + dt*vvels(i,PM0) |
---|
| 248 | positions(i,3)=positions(i,3) + dt*wvels(i,PM0) |
---|
| 249 | particle_step_flag='AB2' |
---|
| 250 | |
---|
| 251 | elseif( particle_step_flag .eq. 'AB2' ) then |
---|
| 252 | positions(i,1)=positions(i,1) + dt2*( 3.*uvels(i,PM0) & |
---|
| 253 | - uvels(i,PM1) ) |
---|
| 254 | |
---|
| 255 | positions(i,2)=positions(i,2) + dt2*( 3.*vvels(i,PM0) & |
---|
| 256 | - vvels(i,PM1) ) |
---|
| 257 | |
---|
| 258 | positions(i,3)=positions(i,3) + dt2*( 3.*wvels(i,PM0) & |
---|
| 259 | - wvels(i,PM1) ) |
---|
| 260 | particle_step_flag='AB3' |
---|
| 261 | |
---|
| 262 | elseif( particle_step_flag .eq. 'AB3' ) then |
---|
| 263 | positions(i,1)=positions(i,1) + dt12*(23.*uvels(i,PM0) & |
---|
| 264 | - 16.*uvels(i,PM1) & |
---|
| 265 | + 5.*uvels(i,PM2) ) |
---|
| 266 | |
---|
| 267 | positions(i,2)=positions(i,2) + dt12*(23.*vvels(i,PM0) & |
---|
| 268 | - 16.*vvels(i,PM1) & |
---|
| 269 | + 5.*vvels(i,PM2) ) |
---|
| 270 | |
---|
| 271 | positions(i,3)=positions(i,3) + dt12*(23.*wvels(i,PM0) & |
---|
| 272 | - 16.*wvels(i,PM1) & |
---|
| 273 | + 5.*wvels(i,PM2) ) |
---|
| 274 | particle_step_flag='AB4' |
---|
| 275 | |
---|
| 276 | elseif( particle_step_flag .eq. 'AB4' ) then |
---|
| 277 | positions(i,1)=positions(i,1) + dt24*(55.*uvels(i,PM0) & |
---|
| 278 | - 59.*uvels(i,PM1) & |
---|
| 279 | + 37.*uvels(i,PM2) & |
---|
| 280 | - 9.*uvels(i,PM3) ) |
---|
| 281 | |
---|
| 282 | positions(i,2)=positions(i,2) + dt24*(55.*vvels(i,PM0) & |
---|
| 283 | - 59.*vvels(i,PM1) & |
---|
| 284 | + 37.*vvels(i,PM2) & |
---|
| 285 | - 9.*vvels(i,PM3) ) |
---|
| 286 | |
---|
| 287 | positions(i,3)=positions(i,3) + dt24*(55.*wvels(i,PM0) & |
---|
| 288 | - 59.*wvels(i,PM1) & |
---|
| 289 | + 37.*wvels(i,PM2) & |
---|
| 290 | - 9.*wvels(i,PM3) ) |
---|
| 291 | endif |
---|
| 292 | |
---|
| 293 | endif |
---|
| 294 | |
---|
| 295 | enddo |
---|
| 296 | |
---|
| 297 | !!================================================================= |
---|
| 298 | !!================================================================= |
---|
| 299 | !! Now globally_update_particle_data |
---|
| 300 | !!================================================================= |
---|
| 301 | !!================================================================= |
---|
| 302 | tmp(:) = uvels(:,PM0) |
---|
| 303 | call mpi_allreduce(tmp,uvels(1,PM0),nparticles, & |
---|
| 304 | mpi_double_precision,mpi_sum,comm,ierr ) |
---|
| 305 | |
---|
| 306 | tmp(:) = vvels(:,PM0) |
---|
| 307 | call mpi_allreduce(tmp,vvels(1,PM0),nparticles, & |
---|
| 308 | mpi_double_precision,mpi_sum,comm,ierr ) |
---|
| 309 | |
---|
| 310 | tmp(:) = wvels(:,PM0) |
---|
| 311 | call mpi_allreduce(tmp,wvels(1,PM0),nparticles, & |
---|
| 312 | mpi_double_precision,mpi_sum,comm,ierr ) |
---|
| 313 | |
---|
| 314 | do i=1,3 |
---|
| 315 | tmp(:) = positions(:,i) |
---|
| 316 | call mpi_allreduce(tmp,positions(1,i),nparticles, & |
---|
| 317 | mpi_double_precision,mpi_sum,comm,ierr ) |
---|
| 318 | enddo |
---|
| 319 | !!================================================================= |
---|
| 320 | !!================================================================= |
---|
| 321 | !! All processors now have updated positions and velocities |
---|
| 322 | !! for all particles, whether they reside locally or not |
---|
| 323 | !!================================================================= |
---|
| 324 | !!================================================================= |
---|
| 325 | end subroutine advance_particles |
---|
| 326 | |
---|
| 327 | |
---|
| 328 | |
---|
| 329 | |
---|
| 330 | |
---|
| 331 | !=============================================================== |
---|
| 332 | !=============================================================== |
---|
| 333 | subroutine spline_interp_value(xvec,yvec,zvec,npts,FF, & |
---|
| 334 | m,ans,my_labels, & |
---|
| 335 | s_x,s_y,s_z, & |
---|
| 336 | nhat_1,nhat_2,nhat_3) |
---|
| 337 | !! |
---|
| 338 | !! (1) Given a vector of (dimensionless) x-y-z locations |
---|
| 339 | !! calculate the corresponding parameter values s. |
---|
| 340 | !! e.g. x=f(s) given x* find s* such that f(s*)=x* |
---|
| 341 | !! Do this in each coordinate direction x-y-z. |
---|
| 342 | !! |
---|
| 343 | !! (2) Call "spline" to get the interpolated field value |
---|
| 344 | !! at each point. This part of the problem is done in |
---|
| 345 | !! the space {s1, s2, s3}, the equally spaced parameters |
---|
| 346 | !! for coordinates x,y,z. The s* values found in (1) define |
---|
| 347 | !! the current particle positions in this space. |
---|
| 348 | !! if |
---|
| 349 | !! m(i) = 0,1,2 return F, F' or F'' wrt the "ith" dimension |
---|
| 350 | !!============================================================= |
---|
| 351 | |
---|
| 352 | use class_ScalarVariable |
---|
| 353 | use user_routines |
---|
| 354 | use user_data, only: Lx,Ly,Lz, & |
---|
| 355 | nx,ny,nz, & |
---|
| 356 | xdim_periodic, & |
---|
| 357 | ydim_periodic, & |
---|
| 358 | zdim_periodic |
---|
| 359 | implicit none |
---|
| 360 | type( ScalarVariable ) :: FF |
---|
| 361 | integer :: m(3) |
---|
| 362 | integer :: npts,i |
---|
| 363 | integer :: my_labels(npts) |
---|
| 364 | real(kind=8) :: xvec(npts) |
---|
| 365 | real(kind=8) :: yvec(npts) |
---|
| 366 | real(kind=8) :: zvec(npts) |
---|
| 367 | real(kind=8) :: s_x (npts) |
---|
| 368 | real(kind=8) :: s_y (npts) |
---|
| 369 | real(kind=8) :: s_z (npts) |
---|
| 370 | real(kind=8),optional :: nhat_1(npts) |
---|
| 371 | real(kind=8),optional :: nhat_2(npts) |
---|
| 372 | real(kind=8),optional :: nhat_3(npts) |
---|
| 373 | real(kind=8) :: locs(3) |
---|
| 374 | real(kind=8) :: svals(3) |
---|
| 375 | real(kind=8) :: ans(npts) |
---|
| 376 | |
---|
| 377 | |
---|
| 378 | !! Get the parametric (s) values corresponding to xyz positions. |
---|
| 379 | do i=1,npts |
---|
| 380 | |
---|
| 381 | if( my_labels(i)==1 ) then |
---|
| 382 | !! get dimensional x,y,z for trajectory i |
---|
| 383 | locs(:)=Lz*(/xvec(i),yvec(i),zvec(i)/) |
---|
| 384 | |
---|
| 385 | !! get the corresponding parameter values s_x,s_y,s_z |
---|
| 386 | call point_map_x2s(locs(1),1,xdim_periodic,Lx,nx,s_x(i)) |
---|
| 387 | call point_map_x2s(locs(2),2,ydim_periodic,Ly,ny,s_y(i)) |
---|
| 388 | call point_map_x2s(locs(3),3,zdim_periodic,Lz,nz,s_z(i)) |
---|
| 389 | endif |
---|
| 390 | |
---|
| 391 | enddo |
---|
| 392 | |
---|
| 393 | if( present(nhat_1) .and. present(nhat_2) .and. present(nhat_3) ) then |
---|
| 394 | call spline(FF,ans,s_x,s_y,s_z,m,npts,my_labels,nhat_1,nhat_2,nhat_3) |
---|
| 395 | else |
---|
| 396 | call spline(FF,ans,s_x,s_y,s_z,m,npts,my_labels) |
---|
| 397 | endif |
---|
| 398 | |
---|
| 399 | |
---|
| 400 | end subroutine spline_interp_value |
---|
| 401 | !=============================================================== |
---|
| 402 | !=============================================================== |
---|
| 403 | |
---|
| 404 | !!=============================================================== |
---|
| 405 | !!=============================================================== |
---|
| 406 | subroutine spline(F,g,x,y,z,m,npts,my_labels,nhat_1,nhat_2,nhat_3) |
---|
| 407 | !! This function interpolates three-dimensional discrete |
---|
| 408 | !! data to a set of arbitary, i.e. non-gridpoint, |
---|
| 409 | !! locations x(1:npts),y(1:npts),z(1:npts). |
---|
| 410 | !! given the discrete gridpoint values F(:,:,:). |
---|
| 411 | !! Results are returned in g(1:npts). |
---|
| 412 | !! |
---|
| 413 | !! m(idim) return m-th derivative, in the idim direction, of F. |
---|
| 414 | !! m can be 0, 1 or 2. |
---|
| 415 | !! |
---|
| 416 | !! Required subroutines |
---|
| 417 | !! db3ink: computes coefficients of the 3D tensor basis functions |
---|
| 418 | !! db3val: evaluates the tensor product spline |
---|
| 419 | !! db2ink: computes coefficients of the 2D tensor basis functions |
---|
| 420 | !! db2val: evaluates the correspondingtensor product spline |
---|
| 421 | !! |
---|
| 422 | !! routine requires no interprocessor communication |
---|
| 423 | !! but, rather, uses the 'not-a-knot' end conditions |
---|
| 424 | !! and computes near edge values locally. Similarly, |
---|
| 425 | !! the spline basis does not explicitly impose periodicity, |
---|
| 426 | !! again relying on the 'not-a-knot' end condition. |
---|
| 427 | !!=============================================================== |
---|
| 428 | !!=============================================================== |
---|
| 429 | use Class_ScalarVariable |
---|
| 430 | use mpi_params |
---|
| 431 | use independent_variables |
---|
| 432 | use user_routines, only: point_map_s2x |
---|
| 433 | use user_data, only: Lx,Ly,Lz, & |
---|
| 434 | xdim_periodic, & |
---|
| 435 | ydim_periodic, & |
---|
| 436 | zdim_periodic |
---|
| 437 | use etc, only: numzplanes, & |
---|
| 438 | numzplanes_below |
---|
| 439 | implicit none |
---|
| 440 | type( ScalarVariable ) :: F |
---|
| 441 | integer :: i,j,k |
---|
| 442 | integer,SAVE :: n1,n2,n3 |
---|
| 443 | integer,SAVE :: nx,ny,nz,globalnz |
---|
| 444 | integer :: m(3),npts,k0,k1 |
---|
| 445 | integer,SAVE :: p,order(3),nwrk,iflag,idim |
---|
| 446 | integer :: my_labels(npts) |
---|
| 447 | real(kind=8) :: x(npts),xval |
---|
| 448 | real(kind=8) :: y(npts),yval |
---|
| 449 | real(kind=8) :: z(npts),zval |
---|
| 450 | real(kind=8) :: g(npts) |
---|
| 451 | real(kind=8),allocatable,SAVE :: Fext(:,:,:) |
---|
| 452 | real(kind=8),allocatable,SAVE :: xgrid(:) |
---|
| 453 | real(kind=8),allocatable,SAVE :: ygrid(:) |
---|
| 454 | real(kind=8),allocatable,SAVE :: zgrid(:) |
---|
| 455 | real(kind=8),allocatable,SAVE :: tx(:) |
---|
| 456 | real(kind=8),allocatable,SAVE :: ty(:) |
---|
| 457 | real(kind=8),allocatable,SAVE :: tz(:) |
---|
| 458 | real(kind=8),allocatable,SAVE :: e3(:,:,:) |
---|
| 459 | real(kind=8),allocatable,SAVE :: wrk(:) |
---|
| 460 | real(kind=8),optional :: nhat_1(npts) |
---|
| 461 | real(kind=8),optional :: nhat_2(npts) |
---|
| 462 | real(kind=8),optional :: nhat_3(npts) |
---|
| 463 | real(kind=8),external :: db3val,db2val |
---|
| 464 | real(kind=8) :: dx,dy,dz,h |
---|
| 465 | real(kind=8) :: f_s,f_x,f_y,f_z,x_s,y_s,z_s |
---|
| 466 | real(kind=8) :: Fmin,Fmax |
---|
| 467 | logical,SAVE :: periodic(3) |
---|
| 468 | logical,SAVE :: limits_test |
---|
| 469 | logical,SAVE :: first_entry=.TRUE. |
---|
| 470 | |
---|
| 471 | if( first_entry ) then |
---|
| 472 | limits_test = .FALSE. |
---|
| 473 | p=3 ! cubic |
---|
| 474 | order(:)=p+1 |
---|
| 475 | periodic=(/xdim_periodic, & |
---|
| 476 | ydim_periodic, & |
---|
| 477 | zdim_periodic/) |
---|
| 478 | |
---|
| 479 | call GetLocalArraySize(F,n1,n2,n3) !! i.e. nx,ny,locnz for F |
---|
| 480 | |
---|
| 481 | !! determine the size of the data array we'll actually work with |
---|
| 482 | if(xdim_periodic) then |
---|
| 483 | nx=n1+1 |
---|
| 484 | else |
---|
| 485 | nx=n1 |
---|
| 486 | endif |
---|
| 487 | if(ydim_periodic .and. n2>1) then |
---|
| 488 | ny=n2+1 |
---|
| 489 | else |
---|
| 490 | ny=n2 |
---|
| 491 | endif |
---|
| 492 | if(zdim_periodic .and. myid==numprocs-1) then |
---|
| 493 | nz=n3+1 |
---|
| 494 | elseif( .not. zdim_periodic .and. myid==numprocs-1 ) then |
---|
| 495 | nz=n3 |
---|
| 496 | endif |
---|
| 497 | if( myid<numprocs-1 ) then |
---|
| 498 | nz=n3+1 |
---|
| 499 | endif |
---|
| 500 | |
---|
| 501 | |
---|
| 502 | nwrk = nx*ny*nz + 2*max(order(1)*(nx+1),order(2)*(ny+1),order(3)*(nz+1)) |
---|
| 503 | !! nwrk=2*nwrk ! I have occaisionally seen errors that go away |
---|
| 504 | !! when the work space is increased, above is adv. size |
---|
| 505 | allocate( wrk(nwrk), e3(nx,ny,nz) ) |
---|
| 506 | allocate( xgrid(nx), ygrid(ny), zgrid(nz) ) |
---|
| 507 | allocate( tx(nx+order(1)), ty(ny+order(2)), tz(nz+order(3)) ) |
---|
| 508 | allocate( Fext(nx,ny,nz) ) |
---|
| 509 | |
---|
| 510 | !! Equally spaced s points where gridded data values reside |
---|
| 511 | !! assumes s in 0 to 1, and endpts now have values |
---|
| 512 | h = 1./(nx-1.) |
---|
| 513 | xgrid(1:nx)=(/( (i-1.)*h,i=1,nx )/) |
---|
| 514 | |
---|
| 515 | if(ny>1) then |
---|
| 516 | h = 1./(ny-1.) |
---|
| 517 | ygrid(1:ny)=(/( (i-1.)*h,i=1,ny )/) |
---|
| 518 | elseif(ny==1) then |
---|
| 519 | ygrid(1)=1.d0 |
---|
| 520 | endif |
---|
| 521 | |
---|
| 522 | !! these are the s positions on myid |
---|
| 523 | k0 = numzplanes_below(myid) + 1 |
---|
| 524 | k1 = k0 + numzplanes(myid) - 1 |
---|
| 525 | zgrid(1:n3) = Grid(0)%SpaceDims(3)%s(k0:k1) |
---|
| 526 | globalnz = ubound(Grid(0)%SpaceDims(3)%s,1) |
---|
| 527 | |
---|
| 528 | !! if we are extending myid's data upward, get the next value as well |
---|
| 529 | if(nz>n3) zgrid(nz) = Grid(0)%SpaceDims(3)%s(k1+1) |
---|
| 530 | |
---|
| 531 | iflag=0 !! use default knot selection scheme |
---|
| 532 | |
---|
| 533 | if( order(1) .ge. nx ) order(1)=nx-1 |
---|
| 534 | if( order(2) .ge. ny ) order(2)=ny-1 |
---|
| 535 | if( order(3) .ge. nz ) order(3)=nz-1 |
---|
| 536 | |
---|
| 537 | first_entry = .FALSE. |
---|
| 538 | endif |
---|
| 539 | |
---|
| 540 | |
---|
| 541 | !! fill extended array |
---|
| 542 | call AssignFextValues(F,n1,n2,n3,Fext,nx,ny,nz,periodic) |
---|
| 543 | |
---|
| 544 | |
---|
| 545 | !!======================================================== |
---|
| 546 | !! Compute coefficients of the 3D tensor basis functions |
---|
| 547 | !! using deBoor routine from nist/gams/cmlib |
---|
| 548 | !!======================================================== |
---|
| 549 | if(ny>1) then |
---|
| 550 | call db3ink(xgrid,nx,ygrid,ny,zgrid,nz, & |
---|
| 551 | Fext(1,1,1),nx,ny, & |
---|
| 552 | order(1),order(2),order(3), & |
---|
| 553 | tx(1),ty(1),tz(1), & |
---|
| 554 | e3,wrk,iflag) |
---|
| 555 | if(iflag .ne. 1 ) stop 'error condition detected by db3ink' |
---|
| 556 | elseif(ny==1) then |
---|
| 557 | call db2ink(xgrid,nx,zgrid,nz,Fext,nx,order(1),order(3), & |
---|
| 558 | tx(1),tz(1),e3,wrk,iflag) |
---|
| 559 | if(iflag .ne. 1 ) stop 'error condition detected by db2ink' |
---|
| 560 | endif |
---|
| 561 | |
---|
| 562 | !!=========================================================== |
---|
| 563 | !! Now loop through each position & interpolate the function |
---|
| 564 | !!=========================================================== |
---|
| 565 | |
---|
| 566 | !!=========================================================== |
---|
| 567 | !! CASE 1: m(:)=0 evaluate function only, |
---|
| 568 | !! and only when my_labels = 1 |
---|
| 569 | !! do 2d ny=1 case separately |
---|
| 570 | !!=========================================================== |
---|
| 571 | if( m(1)==0 .and. m(2)==0 .and. m(3)==0 ) then |
---|
| 572 | if(ny>1) then |
---|
| 573 | do i=1,npts |
---|
| 574 | !! note, positions detected as outside the appropriate |
---|
| 575 | !! range will trigger a fatal error |
---|
| 576 | if( my_labels(i)==1 ) then |
---|
| 577 | g(i) = db3val(x(i),y(i),z(i), & |
---|
| 578 | m(1),m(2),m(3), & |
---|
| 579 | tx,ty,tz, & |
---|
| 580 | nx,ny,nz, & |
---|
| 581 | order(1),order(2),order(3), & |
---|
| 582 | e3,wrk) |
---|
| 583 | else |
---|
| 584 | g(i)=0.d0 |
---|
| 585 | endif |
---|
| 586 | enddo |
---|
| 587 | elseif(ny==1) then |
---|
| 588 | do i=1,npts |
---|
| 589 | !! note, positions detected as outside the appropriate |
---|
| 590 | !! range will trigger a fatal error |
---|
| 591 | if( my_labels(i)==1 ) then |
---|
| 592 | g(i) = db2val(x(i),z(i),m(1),m(3),tx,tz,nx,nz, & |
---|
| 593 | order(1),order(3),e3,wrk) |
---|
| 594 | else |
---|
| 595 | g(i)=0.d0 |
---|
| 596 | endif |
---|
| 597 | enddo |
---|
| 598 | endif |
---|
| 599 | endif |
---|
| 600 | |
---|
| 601 | !!=========================================================== |
---|
| 602 | !! CASE 2: m(1)=0 evaluate Fx,Fy,Fz and dot result with n_hat |
---|
| 603 | !! NB don't bother with my_labels here |
---|
| 604 | !! handle 2d case as above |
---|
| 605 | !! Also, note that values returned from point_map_s2x are |
---|
| 606 | !! DIMENSIONAL and so need to be scaled by Lz |
---|
| 607 | !! Also, though the inputs to point_map_s2x are stored in |
---|
| 608 | !! arrays called x,y,z they are really the s values in |
---|
| 609 | !! each of the 3 dimensions. |
---|
| 610 | !! |
---|
| 611 | !! {x_s}=x_s/Lz df/dx = {df/ds} * {ds/dx} = {df/ds}/{x_s} |
---|
| 612 | !! etc. |
---|
| 613 | !!=========================================================== |
---|
| 614 | if( m(1)==1 .and. m(2)==1 .and. m(3)==1 .and. present(nhat_1) & |
---|
| 615 | .and. present(nhat_2) .and. present(nhat_3) ) then |
---|
| 616 | if(ny>1) then |
---|
| 617 | do i=1,npts |
---|
| 618 | |
---|
| 619 | idim = 1 |
---|
| 620 | call point_map_s2x(x(i),idim,periodic(idim),Lx,n1,xval,x_s) |
---|
| 621 | f_s = db3val(x(i),y(i),z(i), & |
---|
| 622 | 1,0,0, & |
---|
| 623 | tx,ty,tz, & |
---|
| 624 | nx,ny,nz, & |
---|
| 625 | order(1),order(2),order(3), & |
---|
| 626 | e3,wrk) |
---|
| 627 | f_x = f_s / (x_s/Lz) |
---|
| 628 | |
---|
| 629 | idim = 2 |
---|
| 630 | call point_map_s2x(y(i),idim,periodic(idim),Ly,n2,yval,y_s) |
---|
| 631 | f_s = db3val(x(i),y(i),z(i), & |
---|
| 632 | 0,1,0, & |
---|
| 633 | tx,ty,tz, & |
---|
| 634 | nx,ny,nz, & |
---|
| 635 | order(1),order(2),order(3), & |
---|
| 636 | e3,wrk) |
---|
| 637 | f_y = f_s/(y_s/Lz) |
---|
| 638 | |
---|
| 639 | idim = 3 |
---|
| 640 | call point_map_s2x(z(i),idim,periodic(idim),Lz,globalnz,zval,z_s) |
---|
| 641 | f_s = db3val(x(i),y(i),z(i), & |
---|
| 642 | 0,0,1, & |
---|
| 643 | tx,ty,tz, & |
---|
| 644 | nx,ny,nz, & |
---|
| 645 | order(1),order(2),order(3), & |
---|
| 646 | e3,wrk) |
---|
| 647 | f_z = f_s/(z_s/Lz) |
---|
| 648 | |
---|
| 649 | !! dot product |
---|
| 650 | g(i) = nhat_1(i)*f_x + nhat_2(i)*f_y + nhat_3(i)*f_z |
---|
| 651 | |
---|
| 652 | enddo |
---|
| 653 | elseif(ny==1) then |
---|
| 654 | do i=1,npts |
---|
| 655 | |
---|
| 656 | idim = 1 |
---|
| 657 | call point_map_s2x(x(i),idim,periodic(idim),Lx,n1,xval,x_s) |
---|
| 658 | f_s = db2val(x(i),z(i),1,0,tx,tz,nx,nz, & |
---|
| 659 | order(1),order(3),e3,wrk) |
---|
| 660 | f_x = f_s/(x_s/Lz) |
---|
| 661 | |
---|
| 662 | idim = 3 |
---|
| 663 | call point_map_s2x(z(i),idim,periodic(idim),Lz,globalnz,zval,z_s) |
---|
| 664 | f_s = db2val(x(i),z(i),0,1,tx,tz,nx,nz, & |
---|
| 665 | order(1),order(3),e3,wrk) |
---|
| 666 | f_z = f_s/(z_s/Lz) |
---|
| 667 | |
---|
| 668 | !! dot product |
---|
| 669 | g(i) = nhat_1(i)*f_x + nhat_3(i)*f_z |
---|
| 670 | |
---|
| 671 | enddo |
---|
| 672 | endif |
---|
| 673 | endif |
---|
| 674 | |
---|
| 675 | |
---|
| 676 | end subroutine spline |
---|
| 677 | !=============================================================== |
---|
| 678 | !=============================================================== |
---|
| 679 | |
---|
| 680 | |
---|
| 681 | !=============================================================== |
---|
| 682 | !=============================================================== |
---|
| 683 | subroutine AssignFextValues(F,n1,n2,n3,Fext,nx,ny,nz,periodic) |
---|
| 684 | |
---|
| 685 | use class_ScalarVariable |
---|
| 686 | use mpi_params |
---|
| 687 | implicit none |
---|
| 688 | include 'mpif.h' |
---|
| 689 | |
---|
| 690 | type( ScalarVariable ) :: F |
---|
| 691 | logical :: periodic(3) |
---|
| 692 | integer :: n1,n2,n3 !! actual local array size of F%SlabVals |
---|
| 693 | integer :: nx,ny,nz !! size of extended array Fext |
---|
| 694 | integer :: i,j,k |
---|
| 695 | integer :: ierr,reqs(2) |
---|
| 696 | integer :: status_array(mpi_status_size,2) |
---|
| 697 | real(kind=8) :: Fext(nx,ny,nz) |
---|
| 698 | real(kind=8),allocatable,SAVE :: tmp(:,:) |
---|
| 699 | integer, SAVE :: src,count |
---|
| 700 | logical, SAVE :: first_entry=.TRUE. |
---|
| 701 | |
---|
| 702 | |
---|
| 703 | if( first_entry ) then |
---|
| 704 | allocate( tmp(n1,n2) ) |
---|
| 705 | first_entry = .FALSE. |
---|
| 706 | endif |
---|
| 707 | |
---|
| 708 | nreqs=0 |
---|
| 709 | status_array=0 |
---|
| 710 | src=myid+1 |
---|
| 711 | dest=myid-1 |
---|
| 712 | tag=99 |
---|
| 713 | count=n1*n2 |
---|
| 714 | |
---|
| 715 | if(myid /= numprocs-1) then !! all but top processor get ghost data from above |
---|
| 716 | ! recv 1 plane from above |
---|
| 717 | nreqs=nreqs+1 |
---|
| 718 | call MPI_IRECV(tmp,count,mpi_double_precision,src,tag,comm,reqs(nreqs),ierr) |
---|
| 719 | if( ierr .ne. 0 ) then |
---|
| 720 | stop 'MPI IRECV ERROR IN AssignFextValues/particle_routines.f90' |
---|
| 721 | endif |
---|
| 722 | |
---|
| 723 | endif |
---|
| 724 | |
---|
| 725 | !call mpi_barrier(comm,ierr) !! all recvs posted before sends begin |
---|
| 726 | |
---|
| 727 | ! Make these ISENDs, don't worry about buffering (2d) & eliminate the barrier |
---|
| 728 | if( myid /= 0 ) then !! all but bottom processor send ghost data below |
---|
| 729 | ! send bottom plane of F%SlabVals to myid-1 |
---|
| 730 | nreqs=nreqs+1 |
---|
| 731 | call MPI_ISEND(F%SlabVals,count,mpi_double_precision,dest,tag,comm,reqs(nreqs),ierr) |
---|
| 732 | if( ierr .ne. 0 ) then |
---|
| 733 | stop 'MPI IRSEND ERROR IN AssignFextValues/particle_routines.f90' |
---|
| 734 | endif |
---|
| 735 | |
---|
| 736 | endif |
---|
| 737 | |
---|
| 738 | !!================================================== |
---|
| 739 | !! do these copies while waiting for |
---|
| 740 | !! the messages to complete |
---|
| 741 | !!================================================== |
---|
| 742 | Fext(1:n1,1:n2,1:n3) = F%SlabVals(1:n1,1:n2,1:n3) |
---|
| 743 | |
---|
| 744 | if(periodic(1)) then |
---|
| 745 | Fext(nx,1:n2,1:n3)=F%SlabVals(1,1:n2,1:n3) |
---|
| 746 | endif |
---|
| 747 | |
---|
| 748 | if(periodic(2) .and. n2 > 1) then |
---|
| 749 | Fext(1:n1,ny,1:n3)=F%SlabVals(1:n1,1,1:n3) |
---|
| 750 | endif |
---|
| 751 | |
---|
| 752 | if(periodic(1) .and. periodic(2)) then |
---|
| 753 | Fext(nx,ny,1:n3)=F%SlabVals(1,1,1:n3) |
---|
| 754 | endif |
---|
| 755 | !!================================================== |
---|
| 756 | |
---|
| 757 | !!make sure all myid's communications are done before proceeding |
---|
| 758 | if(nreqs>0) then |
---|
| 759 | call mpi_waitall(nreqs,reqs(1),status_array(1,1),ierr) |
---|
| 760 | endif |
---|
| 761 | |
---|
| 762 | !!================================================== |
---|
| 763 | |
---|
| 764 | if(nz>n3) then |
---|
| 765 | Fext(1:n1,1:n2,nz)=tmp(1:n1,1:n2) |
---|
| 766 | Fext(nx,:,nz)=Fext(1,:,nz) |
---|
| 767 | Fext(:,ny,nz)=Fext(:,1,nz) |
---|
| 768 | Fext(nx,ny,nz)=Fext(1,1,nz) |
---|
| 769 | endif |
---|
| 770 | |
---|
| 771 | end subroutine AssignFextValues |
---|