module_poisson_solver.f90
Go to the documentation of this file.
1 !*************************************************************************
2 ! This file is part of HERCULES.
3 !
4 ! HERCULES is free software; you can redistribute it and/or modify
5 ! it under the terms of the GNU General Public License as published by
6 ! the Free Software Foundation; either version 3 of the License, or
7 ! (at your option) any later version.
8 !
9 ! HERCULES is distributed in the hope that it will be useful,
10 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
11 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 ! GNU General Public License for more details.
13 !
14 ! You should have received a copy of the GNU General Public License
15 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
16 !
17 ! Copyright 2015, Ping He, MEAS, NCSU
18 !*************************************************************************
19 
20 
21 
23 
25  implicit none
26 
27 contains
28 
29 
30 
31  ! this function calculates a,b,c for p equation
32  subroutine get_p_eqn_coeff(aa,bb,cc)
33 
34  use mpi
35 
36  implicit none
37 
38  real(mytype), dimension(:), intent(out) :: aa,cc
39  complex(mytype),dimension(:,:,:),intent(out) :: bb
40 
41  integer :: i,j,k,l,m
42  complex(mytype) :: kx,ky
43 
44  aa(:)=0.d0
45  bb(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
46  cc(:)=0.d0
47 
48  ! calculate coefficients a,b,c for the pressure equation
49  do k=kstr3,kend3,1
50 
51  aa(k)=1.d0/dz(k)/dz_b(k)
52  cc(k)=1.d0/dz(k)/dz_t(k)
53 
54  do j=1,zsize(2),1
55  do i=1,zsize(1),1
56 
57  ! deal with index
58  ! convert i to l
59  if ( (i+i_offset) .le. nx/2) then
60  l=i+i_offset-1
61  else
62  l=i+i_offset-nx-1
63  endif
64  ! convert j to m
65  if ( (j+j_offset) .le. ny/2) then
66  m=j+j_offset-1
67  else
68  m=j+j_offset-ny-1
69  endif
70 
71  ! horizontal modified wavenumber,
72  ! check to use 2nd or 4th schemes
73  if ( cds .eq. 1) then
74 
75  kx = ( wx1**l + wx1**(-l) -2.d0 ) /dx2
76 
77  ky = ( wy1**m + wy1**(-m) -2.d0 ) /dy2
78 
79  elseif ( cds .eq. 2) then
80 
81  kx = ( wx1**(3.d0*l) - &
82  54.d0*wx1**(2.d0*l) + &
83  783.d0*wx1**l - &
84  1460.d0 + &
85  783.d0*wx1**(-l) - &
86  54.d0*wx1**(-2.d0*l) + &
87  wx1**(-3.d0*l) ) /576.d0/dx2
88 
89  ky = ( wy1**(3.d0*m) - &
90  54.d0*wy1**(2.d0*m) + &
91  783.d0*wy1**m - &
92  1460.d0 + &
93  783.d0*wy1**(-m) - &
94  54.d0*wy1**(-2.d0*m) + &
95  wy1**(-3.d0*m) ) /576.d0/dy2
96 
97  elseif ( cds .eq. 3) then
98 
99  kx = -(2.d0*pi/lx*l)**2.d0
100  ky = -(2.d0*pi/ly*m)**2.d0
101 
102  endif
103 
104  bb(i,j,k)= kx + ky - aa(k) - cc(k)
105 
106  ! Neumann boundary condition
107  if (k .eq. kstr3) then
108  bb(i,j,k)=bb(i,j,k)+aa(k)
109  aa(k)=0.d0
110  endif
111  if (k .eq. kend3) then
112  bb(i,j,k)=bb(i,j,k)+cc(k)
113  cc(k)=0.d0
114  endif
115 
116  enddo
117  enddo
118  enddo
119 
120  return
121 
122  end subroutine
123 
124 
125 
126  ! this function calculates the right hand side for p equation
127  subroutine get_p_eqn_src(u,v,w,qp)
129  use mpi
130  use module_tools
131 
132  implicit none
133 
134  real(mytype),dimension(:,:,:),intent(in) :: u,v,w
135  real(mytype),dimension(:,:,:),intent(out) :: qp
136 
137  integer :: i,j,k
138  real(mytype) :: dudx,dvdy,dwdz
139 
140  ! calculate qp
141  do k=kstr3,kend3,1
142  do j=jstr3,jend3,1
143  do i=istr3,iend3,1
144 
145  ! check to use 2nd or 4th schemes
146  if ( cds .eq. 1) then
147 
148  dudx = ( u(i,j,k)-u(i-1,j,k) )/dx
149  dvdy = ( v(i,j,k)-v(i,j-1,k) )/dy
150 
151  elseif ( cds .eq. 2) then
152 
153  dudx = ( -u(i+1,j,k) + &
154  27.d0*u(i,j,k) - &
155  27.d0*u(i-1,j,k) + &
156  u(i-2,j,k) ) / 24.d0 / dx
157 
158  dvdy = ( -v(i,j+1,k) + &
159  27.d0*v(i,j,k) - &
160  27.d0*v(i,j-1,k) + &
161  v(i,j-2,k) ) / 24.d0 / dy
162 
163  endif
164 
165  ! z direction is 2nd scheme
166  dwdz = ( w(i,j,k)-w(i,j,k-1) )/dz(k)
167 
168  qp(i,j,k)=(dudx+dvdy+dwdz)/dt
169 
170  enddo
171  enddo
172  enddo
173 
174  return
175 
176  end subroutine
177 
178 
179 
180  ! this function calculates a,b,c for u equation
181  subroutine get_u_eqn_coeff(aa,bb,cc)
183  use mpi
184 
185  implicit none
186 
187  real(mytype), dimension(:), intent(out) :: aa,cc
188  complex(mytype),dimension(:,:,:),intent(out) :: bb
189 
190  integer :: i,j,k,l,m
191  complex(mytype) :: kx,ky
192 
193  aa(:)=0.d0
194  bb(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
195  cc(:)=0.d0
196 
197  ! calculate coefficient for the u equation
198  do k=kstr3,kend3,1
199 
200  aa(k)=-0.5d0*dt*nu/dz(k)/dz_b(k)
201  cc(k)=-0.5d0*dt*nu/dz(k)/dz_t(k)
202 
203  do j=1,zsize(2),1
204  do i=1,zsize(1),1
205 
206  ! deal with index
207  ! convert i to l
208  if ( (i+i_offset) .le. nx/2) then
209  l=i+i_offset-1
210  else
211  l=i+i_offset-nx-1
212  endif
213  ! convert j to m
214  if ( (j+j_offset) .le. ny/2) then
215  m=j+j_offset-1
216  else
217  m=j+j_offset-ny-1
218  endif
219 
220  ! horizontal modified wavenumber,
221  ! check to use 2nd or 4th schemes
222  if ( cds .eq. 1) then
223 
224  kx = 0.5d0 * dt * ( wx1**l + wx1**(-l) - 2.d0 ) *nu/dx2
225 
226  ky = 0.5d0 * dt * ( wy1**m + wy1**(-m) - 2.d0 ) *nu/dy2
227 
228  elseif ( cds .eq. 2) then
229 
230  kx = dt * ( -wx1**(2.d0*l) + &
231  16.d0*wx1**l - &
232  30.d0 + &
233  16.d0*wx1**(-l) - &
234  wx1**(-2.d0*l) ) /24.d0*nu/dx2
235 
236  ky = dt * ( -wy1**(2.d0*m) + &
237  16.d0*wy1**m - &
238  30.d0 + &
239  16.d0*wy1**(-m) - &
240  wy1**(-2.d0*m) ) /24.d0*nu/dy2
241 
242  elseif ( cds .eq. 3) then
243 
244  kx = 0.5d0*dt*nu*( -(2.d0*pi/lx*l)**2.d0 )
245  ky = 0.5d0*dt*nu*( -(2.d0*pi/ly*m)**2.d0 )
246 
247  endif
248 
249  ! note that here we should do 1-kx-ky-aa-cc for bb
250  ! however, to be convenient to implement RK3, we
251  ! add the missing 1 in the subroutine "assign_abc"
252  bb(i,j,k)= - kx - ky - aa(k) - cc(k)
253 
254  ! boundary condition
255  ! we need to add the missing terms to qu
256  ! in Poisson equation solver
257  if (k .eq. kstr3) then
258  bb(i,j,k)=bb(i,j,k)-aa(k) ! non-slip
259  aa(k)=0.d0
260  endif
261 
262  ! we need to add the missing terms to qu
263  ! in Poisson equation solver
264  if (k .eq. kend3) then
265 
266  if (ochannel .eq. 0) then
267  bb(i,j,k)=bb(i,j,k)-cc(k) ! non-slip
268  cc(k)=0.d0
269  elseif (ochannel .eq. 1) then
270  bb(i,j,k)=bb(i,j,k)+cc(k) ! slip
271  cc(k)=0.d0
272  endif
273 
274  endif
275 
276  enddo
277  enddo
278  enddo
279 
280  return
281 
282  end subroutine
283 
284 
285 
286  ! this function calculate a,b,c for w equation
287  subroutine get_w_eqn_coeff(aa,bb,cc)
289  use mpi
290 
291  implicit none
292 
293  real(mytype), dimension(:), intent(out) :: aa,cc
294  complex(mytype),dimension(:,:,:),intent(out) :: bb
295 
296  integer :: i,j,k,l,m
297  complex(mytype) :: kx,ky
298 
299  aa(:)=0.d0
300  bb(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
301  cc(:)=0.d0
302 
303  ! calculate coefficient for the u equation
304  do k=kstr3,kend3-1,1 ! note this index is different
305 
306  aa(k)=-0.5d0*dt*nu/dz(k)/dz_t(k)
307  cc(k)=-0.5d0*dt*nu/dz(k+1)/dz_t(k)
308 
309  do j=1,zsize(2),1
310  do i=1,zsize(1),1
311 
312  ! deal with index
313  ! convert i to l
314  if ( (i+i_offset) .le. nx/2) then
315  l=i+i_offset-1
316  else
317  l=i+i_offset-nx-1
318  endif
319  ! convert j to m
320  if ( (j+j_offset) .le. ny/2) then
321  m=j+j_offset-1
322  else
323  m=j+j_offset-ny-1
324  endif
325 
326  ! horizontal modified wavenumber,
327  ! check to use 2nd or 4th schemes
328  if ( cds .eq. 1) then
329 
330  kx = 0.5d0 * dt * ( wx1**l + wx1**(-l) - 2.d0 ) *nu/dx2
331 
332  ky = 0.5d0 * dt * ( wy1**m + wy1**(-m) - 2.d0 ) *nu/dy2
333 
334  elseif ( cds .eq. 2) then
335 
336  kx = dt * ( -wx1**(2.d0*l) + &
337  16.d0*wx1**l - &
338  30.d0 + &
339  16.d0*wx1**(-l) - &
340  wx1**(-2.d0*l) ) /24.d0*nu/dx2
341 
342  ky = dt * ( -wy1**(2.d0*m) + &
343  16.d0*wy1**m - &
344  30.d0 + &
345  16.d0*wy1**(-m) - &
346  wy1**(-2.d0*m) ) /24.d0*nu/dy2
347 
348  elseif ( cds .eq. 3) then
349 
350  kx = 0.5d0*dt*nu*( -(2.d0*pi/lx*l)**2.d0 )
351  ky = 0.5d0*dt*nu*( -(2.d0*pi/ly*m)**2.d0 )
352 
353  endif
354 
355  ! note that here we should do 1-kx-ky-aa-cc for bb
356  ! however, to be convenient to implement RK3, we
357  ! add the missing 1 in the subroutine "assign_abc"
358  bb(i,j,k)= - kx - ky - aa(k) - cc(k)
359 
360  ! Non-slip condition
361  if (k .eq. kstr3) then
362  aa(k)=0.d0
363  endif
364  if (k .eq. kend3-1) then
365  cc(k)=0.d0
366  endif
367 
368  enddo
369  enddo
370  enddo
371 
372  return
373 
374  end subroutine
375 
376 
377 
378  ! this function calculates a,b,c for t equation
379  subroutine get_t_eqn_coeff(aa,bb,cc)
381  use mpi
382 
383  implicit none
384 
385  real(mytype), dimension(:), intent(out) :: aa,cc
386  complex(mytype),dimension(:,:,:) ,intent(out) :: bb
387 
388  integer :: i,j,k,l,m
389  complex(mytype) :: kx,ky
390 
391  aa(:)=0.d0
392  bb(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
393  cc(:)=0.d0
394 
395  ! calculate coefficients a,b,c for the pressure equation
396  do k=kstr3,kend3,1
397 
398  aa(k)=-0.5d0*dt*alpha/dz(k)/dz_b(k)
399  cc(k)=-0.5d0*dt*alpha/dz(k)/dz_t(k)
400 
401  do j=1,zsize(2),1
402  do i=1,zsize(1),1
403 
404  ! deal with index
405  ! convert i to l
406  if ( (i+i_offset) .le. nx/2) then
407  l=i+i_offset-1
408  else
409  l=i+i_offset-nx-1
410  endif
411  ! convert j to m
412  if ( (j+j_offset) .le. ny/2) then
413  m=j+j_offset-1
414  else
415  m=j+j_offset-ny-1
416  endif
417 
418  ! horizontal modified wavenumber,
419  ! check to use 2nd or 4th schemes
420  if ( cds .eq. 1) then
421 
422  kx = 0.5d0*dt*( wx1**l + wx1**(-l) - 2.d0 )*alpha/dx2
423 
424  ky = 0.5d0*dt*( wy1**m + wy1**(-m) - 2.d0 )*alpha/dy2
425 
426  elseif ( cds .eq. 2) then
427 
428  kx = dt * ( -wx1**(2.d0*l) + &
429  16.d0*wx1**l - &
430  30.d0 + &
431  16.d0*wx1**(-l) - &
432  wx1**(-2.d0*l) ) /24.d0*alpha/dx2
433 
434  ky = dt * ( -wy1**(2.d0*m) + &
435  16.d0*wy1**m - &
436  30.d0 + &
437  16.d0*wy1**(-m) - &
438  wy1**(-2.d0*m) ) /24.d0*alpha/dy2
439 
440  elseif ( cds .eq. 3) then
441 
442  kx = 0.5d0*dt*alpha*( -(2.d0*pi/lx*l)**2.d0 )
443  ky = 0.5d0*dt*alpha*( -(2.d0*pi/ly*m)**2.d0 )
444 
445  endif
446 
447  ! note that here we should do 1-kx-ky-aa-cc for bb
448  ! however, to be convenient to implement RK3, we
449  ! add the missing 1 in the subroutine "assign_abc"
450  bb(i,j,k)= - kx - ky - aa(k) - cc(k)
451 
452  ! Dirichlet condition
453  ! we need to add the missing terms to qt
454  ! in Poisson equation solver
455  if (k .eq. kstr3) then
456 
457  bb(i,j,k)=bb(i,j,k)-aa(k)
458  aa(k)=0.d0
459 
460  endif
461 
462  ! we need to add the missing terms to qt
463  ! in Poisson equation solver
464  if (k .eq. kend3) then
465 
466  bb(i,j,k)=bb(i,j,k)-cc(k)
467  cc(k)=0.d0
468 
469  endif
470 
471  enddo
472  enddo
473  enddo
474 
475  return
476 
477  end subroutine
478 
479 
480  ! Poisson equation fft solver,
481  ! we call it Poisson equation solver,
482  ! although we also use it to solve Helmholtz equation for U,V,W,T
483  subroutine poisson_solver_fft(aa,bb,cc,r_src,r_var,zstagger,delnull,tbnd,ubnd)
485  ! C standard, for using fftw
486  use, intrinsic :: iso_c_binding
487  use module_tools
488  ! these are temporary arrays, their dimensions are: nx/p_row,ny/p_col,nz
489  ! they are on z pencil
490  use module_variables, only : c_src,c_var
491 
492  implicit none
493 
494  ! for using fftw
495  include 'include/fftw3.f03'
496 
497  ! these are Thomas algorithm coefficients
498  real(mytype), dimension(:,:,:),intent(inout) :: aa,cc
499  complex(mytype),dimension(:,:,:),intent(inout) :: bb
500 
501  ! 3d real array, input from source term
502  real(mytype), dimension(:,:,:), intent(in) :: r_src
503 
504  ! these are flags
505  ! zstagger: if is staggered grid in z direction.
506  ! set it to 1 for w equation only
507  ! delnull: whether to delete null space, set it to 1 for p equation only
508  ! tbnd: whether to treat the Dirichlet boundary condition,
509  ! set it to 1 for temperature equation only
510  ! ubnd: special treatment for u boundary condition (for u_mrf),
511  ! this is because in a moving coordinate, the bottom wall u is not zero!!
512  ! set it to 1 for u equation only
513  integer, intent(in) :: zstagger,delnull,tbnd,ubnd
514 
515  ! 3d real array, output for the Poisson fft solution
516  real(mytype),dimension(:,:,:), intent(out) :: r_var
517 
518  ! temporary variable used in Thomas algorithm
519  complex(mytype) :: tmp
520  ! temporary variable used in T and u eqn boundary condition
521  real(mytype) :: aa_t_zstr,cc_t_zend
522  real(mytype) :: aa_u_zstr,cc_u_zend
523  integer :: i,j,k,kmax3
524 
525 
526  ! initiate variables
527  r_var(:,:,:)=0.d0
528  c_var(:,:,:)=cmplx(0.d0,0.d0 , kind=mytype )
529  c_src(:,:,:)=cmplx(0.d0,0.d0 , kind=mytype )
530 
531  ! assign the real source term r_src to complex source term c_src
532  ! note that by default, we are on z pencil
533  do k=1,zsize(3),1
534  do j=1,zsize(2),1
535  do i=1,zsize(1),1
536  c_src(i,j,k)=cmplx(r_src(i+ghst,j+ghst,k+ghst),0.d0, &
537  kind=mytype)
538  enddo
539  enddo
540  enddo
541 
542  ! mpi 2d fft for c_src
543  call mpi_2d_fft(c_src,-1)
544 
545  ! check some flags
546  ! check if it is w equation, the z range will be different
547  if (zstagger .eq. 1) then
548  kmax3=kend3-1
549  elseif (zstagger .eq. 0) then
550  kmax3=kend3
551  else
552  print *, 'Invalid zstagger flag!'
553  stop
554  endif
555  ! check if we need to remove the null space for pressure Poisson equation
556  ! we set the mean component of aa,cc,c_src to 0,
557  ! and therefore the mean pressure at the 1st grid level is 0
558  if (delnull .eq. 1) then
559  ! remove null space
560  if (myid .eq. 0) then
561  aa(1,1,kstr3)=0.d0
562  cc(1,1,kstr3)=0.d0
563  c_src(1,1,1)=cmplx(0.d0,0.d0 , kind=mytype )
564  endif
565  elseif (delnull .eq. 0) then
566 
567  else
568  print *, 'delnull flag error!'
569  stop
570  endif
571  ! check if it is t equation.
572  if (tbnd .eq. 1) then
573  ! add the mising terms to qt for the Dirichlet condition
574  if (myid .eq. 0) then
575  ! applied boundary condition for temperature
576  aa_t_zstr=-0.5d0*rk_c*dt*alpha/dz(kstr3)/dz_b(kstr3)
577  c_src(1,1,1)=c_src(1,1,1)-2.d0*aa_t_zstr*tbot
578  cc_t_zend=-0.5d0*rk_c*dt*alpha/dz(kend3)/dz_t(kend3)
579  c_src(1,1,nz)=c_src(1,1,nz)-2.d0*cc_t_zend*ttop
580  endif
581 
582  elseif (tbnd .eq. 0) then
583 
584  else
585  print *, 'tbnd flag error!'
586  stop
587  endif
588  ! check if it is u equation. We need some treatment for
589  ! the moving reference frame
590  if (ubnd .eq. 1) then
591  ! add the missing terms to qu for the Dirichlet condition
592  if (myid .eq. 0) then
593  ! applied boundary condition for u_mrf
594  aa_u_zstr=-0.5d0*rk_c*dt*nu/dz(kstr3)/dz_b(kstr3)
595  c_src(1,1,1)=c_src(1,1,1)+2.d0*aa_u_zstr*u_mrf
596  if (ochannel .eq. 0) then
597  cc_u_zend=-0.5d0*rk_c*dt*nu/dz(kend3)/dz_t(kend3)
598  c_src(1,1,nz)=c_src(1,1,nz)+2.d0*cc_u_zend*u_mrf
599  endif
600  endif
601 
602  elseif (ubnd .eq. 0) then
603 
604  else
605  print *, 'ubnd flag error!'
606  stop
607  endif
608 
609  ! Thomas algorithm to calculate c_var
610  ! forward loop
611  do k=kstr3+1,kmax3,1
612  do j=1,zsize(2),1
613  do i=1,zsize(1),1
614 
615  tmp=aa(i,j,k)/bb(i,j,k-1)
616  bb(i,j,k)=bb(i,j,k)-cc(i,j,k-1)*tmp
617  c_src(i,j,k-ghst)=c_src(i,j,k-ghst)-c_src(i,j,k-ghst-1)*tmp
618 
619  enddo
620  enddo
621  enddo
622  ! backward loop
623  do j=1,zsize(2),1
624  do i=1,zsize(1),1
625 
626  c_var(i,j,kmax3-ghst)=c_src(i,j,kmax3-ghst)/bb(i,j,kmax3)
627 
628  enddo
629  enddo
630 
631  do k=kmax3-1,kstr3,-1
632  do j=1,zsize(2),1
633  do i=1,zsize(1),1
634 
635  c_var(i,j,k-ghst)=( c_src(i,j,k-ghst)- &
636  cc(i,j,k)*c_var(i,j,k-ghst+1) ) / bb(i,j,k)
637 
638  enddo
639  enddo
640  enddo
641 
642  ! mpi 2d ifft, transfer c_var to r_var
643  call mpi_2d_fft(c_var,1)
644 
645  ! assign c_var to r_var
646  do k=kstr3,kmax3,1
647  do j=1,zsize(2),1
648  do i=1,zsize(1),1
649 
650  r_var(i+ghst,j+ghst,k)=real(c_var(i,j,k-ghst), kind=mytype )
651 
652  enddo
653  enddo
654  enddo
655 
656  return
657 
658  end subroutine
659 
660 
661 
662 end module
integer, save, protected myid
integer, save, protected ochannel
integer, save, protected ny
subroutine mpi_2d_fft(c_inout3, direction)
subroutine get_u_eqn_coeff(aa, bb, cc)
subroutine get_p_eqn_src(u, v, w, qp)
complex(mytype), dimension(:,:,:), allocatable, save c_src
real(mytype), parameter pi
real(mytype), save, protected dt
complex(mytype), save, protected wy1
integer, save, protected jend3
real(mytype), save, protected nu
real(mytype), save, protected ttop
integer, save, protected nz
subroutine get_t_eqn_coeff(aa, bb, cc)
real(mytype), save, protected ly
real(mytype), dimension(:), allocatable, save, protected dz
subroutine poisson_solver_fft(aa, bb, cc, r_src, r_var, zstagger, delnull, tbnd, ubnd)
real(mytype), save, protected dy2
real(mytype), save, protected alpha
subroutine get_p_eqn_coeff(aa, bb, cc)
integer, save, protected istr3
real(mytype), save, protected dx
integer, save, protected jstr3
real(mytype), dimension(:), allocatable, save, protected dz_b
integer, save, protected i_offset
real(mytype), save, protected dx2
integer, save, protected j_offset
real(mytype), save, protected u_mrf
integer, parameter ghst
integer, save, protected cds
real(mytype), save, protected lx
real(mytype), dimension(:), allocatable, save, protected dz_t
real(mytype), save, protected dy
integer, save, protected nx
subroutine get_w_eqn_coeff(aa, bb, cc)
complex(mytype), save, protected wx1
real(mytype), save, protected rk_c
integer, save, protected iend3
real(mytype), save, protected tbot
complex(mytype), dimension(:,:,:), allocatable, save c_var
integer, save, protected kstr3
integer, save, protected kend3