Context Navigation

source: trunk/00FlowSolve_KW/src/particle_routines.f90 @ 4

Last change on this file since 4 was 4, checked in by xlvlod, 17 years ago
import initial des vraies codes.
File size: 28.4 KB

Line
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 = nxnynz + 2max(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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: