84 = 2.50662827463100050241576528481104525300698674
86 REAL(c_double),
DIMENSION(:,:),
POINTER :: original
87 COMPLEX(C_DOUBLE_COMPLEX),
DIMENSION(:,:),
POINTER :: transformed
92 TYPE(c_ptr) :: plan_r2c
93 TYPE(c_ptr) :: plan_c2r
94 TYPE(
field_typ),
DIMENSION(:),
ALLOCATABLE & !< z-layered density/potential field
96 REAL(c_double),
DIMENSION(:,:),
POINTER,
CONTIGUOUS &
98 COMPLEX(C_DOUBLE_COMPLEX),
POINTER,
CONTIGUOUS &
99 :: fmass3d(:,:,:) => null(),&
100 fsum3d(:,:,:) => null(), &
103 REAL,
DIMENSION(:,:,:),
POINTER,
CONTIGUOUS &
104 :: fmass3d_real => null(), &
107 INTEGER(C_INTPTR_T) :: local_joff
108 REAL,
DIMENSION(:),
POINTER &
109 :: kx => null(),ky => null()
112 REAL,
DIMENSION(:),
POINTER &
113 :: joff=>null(),jrem=>null()
116 REAL,
DIMENSION(:,:,:),
POINTER :: mpi_buf => null()
117 INTEGER(C_INTPTR_T) :: c_inum, c_jnum
118 INTEGER(C_INTPTR_T) :: alloc_local, local_jnum
119 TYPE(c_ptr) :: fftw_real_pointer, &
163 TYPE(
dict_typ),
POINTER :: config,IO
165 INTEGER :: err, valwrite, i,j,k
166#if defined(HAVE_FFTW) && defined(PARALLEL)
167 INTEGER(C_INTPTR_T) :: C_INUM,C_JNUM
170 CALL this%InitGravity(mesh,physics,
"shearingbox spectral solver",config,io)
173#if !defined(HAVE_FFTW)
174 CALL this%Error(
"InitGravity_sboxspectral", &
175 "No fftw package could be loaded.")
177 CALL getattr(config,
"order", this%order, 2)
178 SELECT CASE(this%order)
182 IF (mesh%KNUM.GT.1) &
183 CALL this%Error(
"InitGravity_sboxspectral", &
184 "order 4 is currently not supported in 3D")
186 CALL this%Error(
"InitGravity_sboxspectral", &
187 "Order must be one of {2,4}.")
197 IF (mesh%ROTSYM.NE.0) &
198 CALL this%Error(
"InitGravity_sboxspectral", &
199 "Rotational symmetry not supported.")
202 IF (mesh%NDIMS.EQ.1) &
203 CALL this%Error(
"InitGravity_sboxspectral", &
204 "Only 2D (shearingsheet) and 3D (shearingbox) simulations allowed with this module.")
207 SELECT TYPE (phys => physics)
211 CALL this%Error(
"InitGravity_sboxspectral", &
212 "Physics modules with density necessary for this module.")
220 IF (this%GetNumProcs().GT.1 .AND. mesh%shear_dir.EQ.2)
THEN
221 CALL this%Error(
"InitGravity_sboxspectral", &
222 "Execution of the parallel Fourier solver is only possible with pencil "// &
223 "decomposition. The first dimension should be fully accessible by the "//&
224 "first core. The shear thereto needs to vary along the second " // &
225 "dimension. In order to achieve this set the boundaries accordingly. "// &
226 "Check InitBoundary where the shear direction is set.")
228 IF (mesh%dims(1).GT.1 .AND. mesh%shear_dir.EQ.1)
THEN
229 CALL this%Error(
"InitGravity_sboxspectral", &
230 "The first dimension needs to be fully accessible by each process. "// &
231 "Please change domain decomposition to pencil decompositon. " // &
232 "This needs to be done during initialization. In the Mesh dictionary " // &
233 "use the key 'decompositon' and the value (/ 1,-1, 1/) to force decomposition " // &
234 "along second direction")
239 IF (mesh%KNUM.GT.1)
THEN
240 IF (mesh%zmin.NE.-mesh%zmax) &
241 CALL this%Error(
"InitGravity_sboxspectral",
"z_min != -zmax not allowed: " // &
242 achar(13) //
" computational domain must be symmetric " // &
243 "with respect to z=0 plane")
245 CALL this%Warning(
"InitGravity_sboxspectral",
"sboxspectral 3D is experimental")
250 this%phi(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
251 this%qk(mesh%IMIN/2+1:mesh%IMAX/2+1,mesh%JMIN:mesh%JMAX,0:mesh%KNUM+1),&
252 this%field(mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1), &
256 this%Fmass3D(mesh%IMIN/2+1:mesh%IMAX/2+1,mesh%JMIN:mesh%JMAX,mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1), &
259 this%Fsum3D(mesh%IMIN/2+1:mesh%IMAX/2+1,mesh%JMIN:mesh%JMAX,mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1), &
261 this%Kxy2(mesh%IMIN/2+1:mesh%IMAX/2+1,mesh%JMIN:mesh%JMAX), &
262 this%kx(mesh%INUM), &
263 this%ky(mesh%JNUM), &
264 this%den_ip(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1), &
267 CALL this%Error(
"InitGravity_sboxspectral",
"Memory allocation failed.")
276 ALLOCATE(this%mpi_buf(mesh%IMIN:mesh%IMAX,mesh%GJNUM,mesh%KMIN:mesh%KMAX), &
279 CALL this%Error(
"InitGravity_sboxspectral",
"Memory allocation failed for mpi_buf.")
281 this%alloc_local = fftw_mpi_local_size_2d(c_jnum, c_inum, &
282 mpi_comm_world, this%local_JNUM, this%local_joff)
283 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
284 this%fftw_real_pointer = fftw_alloc_real(2*this%alloc_local)
285 this%fftw_complex_pointer = fftw_alloc_complex(this%alloc_local)
286 CALL c_f_pointer(this%fftw_real_pointer,this%field(k)%original, &
287 [2*(c_inum/2+1),this%local_JNUM])
288 CALL c_f_pointer(this%fftw_complex_pointer,this%field(k)%transformed, &
289 [c_inum/2+1,this%local_JNUM])
298 CALL this%Info(
" initializing FFTW: " // &
299#if !defined(PARALLEL)
301 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
302 this%field(k)%original => this%den_ip(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k)
303 this%field(k)%transformed => this%Fmass3D(mesh%IMIN/2+1:mesh%IMAX/2+1,mesh%JMIN:mesh%JMAX,k)
305 this%plan_r2c = fftw_plan_dft_r2c_2d(mesh%JMAX-mesh%JMIN+1, &
306 mesh%IMAX-mesh%IMIN+1,this%field(mesh%KMIN)%original, &
307 this%field(mesh%KMIN)%transformed,fftw_measure)
308 this%plan_c2r = fftw_plan_dft_c2r_2d(mesh%JMAX-mesh%JMIN+1, &
309 mesh%IMAX-mesh%IMIN+1, this%field(mesh%KMIN)%transformed, &
310 this%field(mesh%KMIN)%original, fftw_measure)
311 IF ((.NOT. c_associated(this%plan_r2c)) .OR. (.NOT. c_associated(this%plan_c2r)))
THEN
312 CALL this%Error(
"InitGravity_sboxspectral",
"FFT plan could not be created.")
314#elif defined(PARALLEL)
316 this%plan_r2c = fftw_mpi_plan_dft_r2c_2d(c_jnum,c_inum,this%field(mesh%KMIN)%original, &
317 this%field(mesh%KMIN)%transformed, &
318 mpi_comm_world, fftw_measure)
319 this%plan_c2r = fftw_mpi_plan_dft_c2r_2d(c_jnum,c_inum, this%field(mesh%KMIN)%transformed, &
320 this%field(mesh%KMIN)%original, &
321 mpi_comm_world, fftw_measure)
323 CALL getattr(config,
"print_plans", valwrite, 0)
324 IF (valwrite.GT.0)
THEN
326 IF (this%GetRank().EQ.0)
THEN
327 CALL fftw_print_plan(this%plan_r2c)
329 CALL fftw_print_plan(this%plan_c2r)
333 CALL this%Warning(
"InitGravity_sboxspectral",
"fftw_print_plan currently not supported on SX-Aurora",0)
340 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
341 this%field(k)%transformed(:,:) = cmplx(0.,0.)
343 this%Fsum3D(:,:,k) = cmplx(0.,0)
350 this%kx(:) = cshift((/(i-(mesh%INUM+1)/2,i=0,mesh%INUM-1)/),(mesh%INUM+1)/2) &
351 *2.*pi/(mesh%XMAX-mesh%XMIN)
352 this%ky(:) = cshift((/(j-(mesh%JNUM+1)/2,j=0,mesh%JNUM-1)/),(mesh%JNUM+1)/2) &
353 *2.*pi/(mesh%YMAX-mesh%YMIN)
356 this%maxKxy2 = 0.5*min(pi/mesh%dx,pi/mesh%dy)**2
358 SELECT CASE (mesh%shear_dir)
360 this%shiftconst = mesh%Q*mesh%OMEGA*(mesh%ymax-mesh%ymin)/mesh%dx
362 this%joff(mesh%JMIN:mesh%JMAX), &
363 this%jrem(mesh%JMIN:mesh%JMAX), &
366 this%shiftconst = mesh%Q*mesh%OMEGA*(mesh%xmax-mesh%xmin)/mesh%dy
368 this%joff(mesh%IMIN:mesh%IMAX), &
369 this%jrem(mesh%IMIN:mesh%IMAX), &
372 CALL this%Error(
"InitGravity_sboxspectral", &
373 "Shear direction must be one of WE_shear (x-direction) or SN_shear (y-direction).")
376 CALL this%Error(
"InitGravity_sboxspectral",
"Memory allocation failed for joff & jrem.")
381 CALL this%SetOutput(mesh,physics,config,io)
391 TYPE(
dict_typ),
POINTER :: config,IO
393 INTEGER :: valwrite,err
397 CALL getattr(config,
"output/potential", valwrite, 0)
398 IF (valwrite .EQ. 1)
THEN
399 IF (
ASSOCIATED(this%phi)) &
400 CALL setattr(io,
"potential", &
401 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
404 CALL getattr(config,
"output/Fmass3D", valwrite, 0)
405 IF (valwrite .EQ. 1)
THEN
406 ALLOCATE(this%Fmass3D_real(mesh%IMIN:mesh%IMAX+2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
409 CALL this%Error(
"sboxspectral::SetOutput",
"Memory allocation failed for Fmass3D_real.")
410 this%Fmass3D_real(:,:,:) = 0.0
411 CALL setattr(io,
"Fmass3D_real", &
412 this%Fmass3D_real(mesh%IMIN:mesh%IMAX+2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
432 REAL,
INTENT(IN) :: time,dt
439 CALL this%CalcPotential(mesh,physics,time,pvar)
441 IF (this%order.EQ.2)
THEN
443 DO k = mesh%KMIN,mesh%KMAX
445 DO j = mesh%JMIN,mesh%JMAX
447 DO i = mesh%IMIN,mesh%IMAX
449 this%accel%data4d(i,j,k,1) = -1.0*(this%phi(i+1,j,k)-this%phi(i-1,j,k))/ &
450 (2*mesh%dlx%data3d(i,j,k))
451 this%accel%data4d(i,j,k,2) = -1.0*(this%phi(i,j+1,k)-this%phi(i,j-1,k))/ &
452 (2*mesh%dly%data3d(i,j,k))
453 IF (physics%VDIM.EQ.3) &
454 this%accel%data4d(i,j,k,3) = -1.0*(this%phi(i,j,k+1)-this%phi(i,j,k-1))/ &
455 (2*mesh%dlz%data3d(i,j,k))
463 DO k = mesh%KMIN,mesh%KMAX
465 DO j = mesh%JMIN,mesh%JMAX
467 DO i = mesh%IMIN,mesh%IMAX
469 this%accel%data4d(i,j,k,1) = -1.0*(w1*this%phi(i-2,j,k)-w2*this%phi(i-1,j,k)+ &
470 w2*this%phi(i+1,j,k)-w1*this%phi(i+2,j,k))/ &
471 (mesh%dlx%data3d(i,j,k))
472 this%accel%data4d(i,j,k,2) = -1.0*(w1*this%phi(i,j-2,k)-w2*this%phi(i,j-1,k)+ &
473 w2*this%phi(i,j+1,k)-w1*this%phi(i,j+2,k)) / &
474 (mesh%dly%data3d(i,j,k))
475 IF (physics%VDIM.EQ.3) &
476 this%accel%data4d(i,j,k,3) = -1.0*(w1*this%phi(i,j,k-2)-w2*this%phi(i,j,k-1)+ &
477 w2*this%phi(i,j,k+1)-w1*this%phi(i,j,k+2)) / &
478 (mesh%dlz%data3d(i,j,k))
529 REAL,
INTENT(IN) :: time
539 IF (mesh%OMEGA .EQ. 0)
THEN
541 ELSE IF (mesh%shear_dir.EQ.2)
THEN
542 time0 = nint(time*mesh%Q*mesh%OMEGA*(mesh%xmax-mesh%xmin)/ &
543 (mesh%ymax-mesh%ymin))*(mesh%ymax-mesh%ymin)/(mesh%Q*mesh%OMEGA &
544 *(mesh%xmax-mesh%xmin))
545 ELSE IF (mesh%shear_dir.EQ.1)
THEN
546 time0 = nint(time*mesh%Q*mesh%OMEGA*(mesh%ymax-mesh%ymin)/ &
547 (mesh%xmax-mesh%xmin))*(mesh%xmax-mesh%xmin)/(mesh%Q*mesh%OMEGA &
548 *(mesh%ymax-mesh%ymin))
553 SELECT TYPE(p => pvar)
555 CALL this%FieldShift(mesh,physics,delt, &
556 p%density%data3d(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
557 this%field(mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1))
559 CALL this%Error(
"gravity_sboxspectral::CalcPotential",
"unsupported state vector")
565CALL ftrace_region_begin(
"forward_fft")
571 CALL this%FFT_Forward(mesh,physics,this%field(mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1))
575CALL ftrace_region_end(
"forward_fft")
579 IF (mesh%shear_dir.EQ.2)
THEN
581 DO j = mesh%JMIN,mesh%JMAX
583 DO i = mesh%IMIN/2+1,mesh%IMAX/2+1
584 this%Kxy2(i,j) = (this%kx(i) + mesh%Q*mesh%OMEGA*this%ky(j)*delt)**2 + this%ky(j)**2
587 ELSE IF (mesh%shear_dir.EQ.1)
THEN
589 DO j = mesh%JMIN,mesh%JMAX
591 DO i = mesh%IMIN/2+1,mesh%IMAX/2+1
592 this%Kxy2(i,j) = this%kx(i)**2 + (this%ky(j)-mesh%Q*mesh%OMEGA*this%kx(i)*delt)**2
598 IF (abs(mesh%zmax-mesh%zmin).LT.tiny(mesh%zmin))
THEN
601 DO j = mesh%JMIN,mesh%JMAX
603 DO i = mesh%IMIN/2+1,mesh%IMAX/2+1
604 IF ((this%Kxy2(i,j).LT.this%maxKxy2).AND.(this%Kxy2(i,j).GT. 0))
THEN
605 this%field(k)%transformed(i,j-this%local_joff) = -2.*pi*physics%Constants%GN &
606 *this%field(k)%transformed(i,j-this%local_joff)/sqrt(this%Kxy2(i,j))
608 this%field(k)%transformed(i,j-this%local_joff) = 0.0
621 this%qk(:,:,0) = exp(-0.5*mesh%dz*sqrt(this%Kxy2(:,:)))
623 IF (mesh%KNUM.EQ.1)
THEN
627 DO j = mesh%JMIN,mesh%JMAX
629 DO i = mesh%IMIN/2+1,mesh%IMAX/2+1
630 IF ((this%Kxy2(i,j).LT.this%maxKxy2).AND.(this%Kxy2(i,j).GT. 0))
THEN
631 this%field(k)%transformed(i,j-this%local_joff) = -4.*pi*physics%Constants%GN &
632 *(1.-this%qk(i,j,0))/this%Kxy2(i,j) * this%field(k)%transformed(i,j-this%local_joff)
634 this%field(k)%transformed(i,j-this%local_joff) = 0.0
640 this%qk(:,:,1) = this%qk(:,:,0)*this%qk(:,:,0)
642 this%qk(:,:,k) = this%qk(:,:,k-1) * this%qk(:,:,1)
644 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
645 this%Fsum3D(:,:,k) = 0.0
648 DO kk=mesh%KMIN-mesh%KP1,k-1
650 DO j = mesh%JMIN,mesh%JMAX
652 DO i = mesh%IMIN,mesh%IMAX/2+1
653 this%Fsum3D(i,j,k) = this%Fsum3D(i,j,k) + this%qk(i,j,abs(k-kk)) &
654 * this%field(kk)%transformed(i,j-this%local_joff)
660 DO j = mesh%JMIN,mesh%JMAX
662 DO i = mesh%IMIN/2+1,mesh%IMAX/2+1
663 DO kk=k+1,mesh%KMAX+mesh%KP1
664 this%Fsum3D(i,j,k) = this%Fsum3D(i,j,k) + this%qk(i,j,abs(k-kk)) &
665 * this%field(kk)%transformed(i,j-this%local_joff)
670 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
672 DO j = mesh%JMIN,mesh%JMAX
674 DO i = mesh%IMIN/2+1,mesh%IMAX/2+1
675 IF ((this%Kxy2(i,j).LT.this%maxKxy2).AND.(this%Kxy2(i,j).GT. 0))
THEN
676 this%field(k)%transformed(i,j-this%local_joff) = -4.*pi*physics%Constants%GN &
677 * (1.-this%qk(i,j,0))/this%Kxy2(i,j) &
678 * (this%field(k)%transformed(i,j-this%local_joff) + 0.5*(1.0+1./this%qk(i,j,0)) &
679 * this%Fsum3D(i,j,k))
681 this%field(k)%transformed(i,j-this%local_joff) = 0.0
691CALL ftrace_region_begin(
"backward_fft")
696 CALL this%FFT_Backward(mesh,physics,this%field(mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1))
699CALL ftrace_region_end(
"backward_fft")
704 DO k = mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
706 DO j = mesh%JMIN,mesh%JMAX
708 DO i = mesh%IMIN,mesh%IMAX
709 this%phi(i,j,k) = this%field(k)%original(i,j-this%local_joff)/ &
710 (mesh%JNUM*mesh%INUM)
716 CALL this%FieldShift(mesh,physics,-delt, &
717 this%phi(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
718 this%field(mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1))
721 DO k = mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
723 DO j = mesh%JMIN,mesh%JMAX
725 DO i = mesh%IMIN,mesh%IMAX
726 this%phi(i,j,k) = this%field(k)%original(i,j-this%local_joff)
731 CALL this%SetBoundaryData(mesh,physics,delt)
744 REAL,
INTENT(IN) :: delt
749 INTEGER :: status(MPI_STATUS_SIZE),ierror
752 IF(mesh%shear_dir.EQ.2)
THEN
754 DO k = mesh%KMIN,mesh%KMAX
757 DO i = mesh%IMIN,mesh%IMAX
758 this%phi(i,mesh%JMIN-j,k) = this%phi(i,mesh%JMAX-j+1,k)
759 this%phi(i,mesh%JMAX+j,k) = this%phi(i,mesh%JMIN+j-1,k)
765 joff = -this%shiftconst*delt
766 jrem = joff - floor(joff)
767 DO k = mesh%KMIN,mesh%KMAX
770 DO j = mesh%JMIN,mesh%JMAX
772 this%phi(mesh%IMIN-i,j,k) = this%phi(mesh%IMAX-i+1,j,k)
775 this%phi(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,k) = &
776 cshift(this%phi(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,k),floor(joff))
777 this%phi(mesh%IMIN-i,mesh%JMAX+1,k) = this%phi(mesh%IMIN-i,mesh%JMIN,k)
779 DO j = mesh%JMIN,mesh%JMAX
780 this%phi(mesh%IMIN-i,j,k) = &
781 (1.0 - jrem)*this%phi(mesh%IMIN-i,j,k) + jrem*this%phi(mesh%IMIN-i,j+1,k)
786 joff = this%shiftconst*delt
787 jrem = joff - floor(joff)
788 DO k = mesh%KMIN,mesh%KMAX
791 DO j = mesh%JMIN,mesh%JMAX
793 this%phi(mesh%IMAX+i,j,k) = this%phi(mesh%IMIN+i-1,j,k)
796 this%phi(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,k) = &
797 cshift(this%phi(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,k),floor(joff))
798 this%phi(mesh%IMAX+i,mesh%JMAX+1,:) = this%phi(mesh%IMAX+i,mesh%JMIN,:)
800 DO j = mesh%JMIN,mesh%JMAX
801 this%phi(mesh%IMAX+i,j,k) = &
802 (1.0 - jrem)*this%phi(mesh%IMAX+i,j,k) + jrem*this%phi(mesh%IMAX+i,j+1,k)
806 ELSE IF(mesh%shear_dir.EQ.1)
THEN
808 DO k = mesh%KMIN,mesh%KMAX
809 DO j = mesh%JMIN,mesh%JMAX
812 this%phi(mesh%IMIN-i,j,k) = this%phi(mesh%IMAX-i+1,j,k)
813 this%phi(mesh%IMAX+i,j,k) = this%phi(mesh%IMIN+i-1,j,k)
821 IF(mesh%dims(2).GT.1)
THEN
822 this%mpi_buf(mesh%IMIN:mesh%IMAX,1:mesh%GJNUM,mesh%KMIN:mesh%KMAX) = &
823 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX-mesh%GJNUM+1:mesh%JMAX,mesh%KMIN:mesh%KMAX)
824 CALL mpi_sendrecv_replace(&
826 SIZE(this%mpi_buf), &
829 mesh%neighbor(
south), mpi_any_tag, &
830 mesh%comm_cart, status, ierror)
832 this%phi(mesh%IMIN:mesh%IMAX,mesh%JGMIN:mesh%JMIN-1,mesh%KMIN:mesh%KMAX) = &
833 this%mpi_buf(mesh%IMIN:mesh%IMAX,1:mesh%GJNUM,mesh%KMIN:mesh%KMAX)
835 this%mpi_buf(mesh%IMIN:mesh%IMAX,1:mesh%GJNUM,mesh%KMIN:mesh%KMAX) = &
836 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMIN+mesh%GJNUM-1,mesh%KMIN:mesh%KMAX)
838 CALL mpi_sendrecv_replace(&
840 SIZE(this%mpi_buf), &
843 mesh%neighbor(
north), mpi_any_tag, &
844 mesh%comm_cart, status, ierror)
845 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+1:mesh%JGMAX,mesh%KMIN:mesh%KMAX) = &
846 this%mpi_buf(mesh%IMIN:mesh%IMAX,1:mesh%GJNUM,mesh%KMIN:mesh%KMAX)
851 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX)
852 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX)
858 IF (mesh%mycoords(2).EQ.0)
THEN
864 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX)
866 joff = this%shiftconst*delt
867 jrem = joff - floor(joff)
869 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX) = &
870 cshift(this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX),floor(joff))
872 this%phi(mesh%IMAX+1,mesh%JMIN-j,:) = this%phi(mesh%IMIN,mesh%JMIN-j,:)
873 DO k = mesh%KMIN,mesh%KMAX
874 DO i=mesh%IMIN,mesh%IMAX
875 this%phi(i,mesh%JMIN-j,k) = &
876 (1.0 - jrem)*this%phi(i,mesh%JMIN-j,k) + jrem*this%phi(i+1,mesh%JMIN-j,k)
885 IF (mesh%mycoords(2).EQ.mesh%dims(2)-1)
THEN
890 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX)
893 joff = -this%shiftconst*delt
894 jrem = joff - floor(joff)
896 this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX) = &
897 cshift(this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX),floor(joff))
899 this%phi(mesh%IMAX+1,mesh%JMAX+j,:) = this%phi(mesh%IMIN,mesh%JMAX+j,:)
900 DO k = mesh%KMIN,mesh%KMAX
901 DO i = mesh%IMIN,mesh%IMAX
902 this%phi(i,mesh%JMAX+j,k) = &
903 (1.0 - jrem)*this%phi(i,mesh%JMAX+j,k) + jrem*this%phi(i+1,mesh%JMAX+j,k)
920 TYPE(
field_typ),
DIMENSION(Mesh%KMIN-Mesh%KP1:Mesh%KMAX+Mesh%KP1), &
921 INTENT(INOUT) :: field
927CALL ftrace_region_begin(
"forward FFT")
930 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
932 CALL fftw_mpi_execute_dft_r2c(this%plan_r2c,field(k)%original, &
933 this%field(k)%transformed)
935 CALL fftw_execute_dft_r2c(this%plan_r2c,this%field(k)%original, &
936 this%field(k)%transformed)
940 IF (
ASSOCIATED(this%Fmass3D_real))
THEN
942 DO k = mesh%KMIN,mesh%KMAX
944 DO j = mesh%JMIN,mesh%JMAX
946 DO i = mesh%IMIN/2+1,mesh%IMAX/2+1
947 this%Fmass3D_real(2*i-1,j,k) = real(real(this%Fmass3D(i,j-this%local_joff,k)))
948 this%Fmass3D_real(2*i,j,k) = real(aimag(this%Fmass3D(i,j-this%local_joff,k)))
955CALL ftrace_region_end(
"foward FFT")
967 TYPE(
field_typ),
DIMENSION(Mesh%KMIN-Mesh%KP1:Mesh%KMAX+Mesh%KP1), &
968 INTENT(INOUT) :: field
972 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
974 CALL fftw_mpi_execute_dft_c2r(this%plan_c2r,field(k)%transformed, &
977 CALL fftw_execute_dft_c2r(this%plan_c2r,field(k)%transformed, &
1027 TYPE(
marray_base),
INTENT(INOUT) :: bccsound,h_ext,height
1030 REAL :: cs2,p,q,rho_c
1033 IF (abs(mesh%zmax-mesh%zmin).LT.tiny(mesh%zmin))
THEN
1037 DO j=mesh%JGMIN,mesh%JGMAX
1039 DO i=mesh%IGMIN,mesh%IGMAX
1040 cs2 = bccsound%data3d(i,j,k)*bccsound%data3d(i,j,k)
1041 p = -
sqrttwopi*physics%Constants%GN*pvar%data4d(i,j,k,physics%DENSITY) &
1043 q = cs2/mesh%OMEGA**2.
1045 height%data3d(i,j,k) = p+sqrt(q+p*p)
1051 IF (k.LT.mesh%KMIN.OR.k.GT.mesh%KMAX) &
1052 CALL this%Error(
"gravity_sboxspectral::CalcDiskHeight_single",
"not yet parallelized")
1055 DO j=mesh%JGMIN,mesh%JGMAX
1057 DO i=mesh%IGMIN,mesh%IGMAX
1059 IF (mod(mesh%KNUM,2).EQ.0)
THEN
1061 rho_c = 0.5*(pvar%data4d(i,j,k,physics%DENSITY)+pvar%data4d(i,j,k+1,physics%DENSITY))
1064 rho_c = pvar%data4d(i,j,k+1,physics%DENSITY)
1067 height%data3d(i,j,mesh%KGMIN:mesh%KGMAX) = bccsound%data3d(i,j,k) &
1068 / sqrt(4*pi*physics%Constants%GN * rho_c + mesh%OMEGA**2)
1084 SUBROUTINE fieldshift(this,Mesh,Physics,delt,field,shifted_field)
1090 REAL,
INTENT(IN) :: delt
1091 REAL,
DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
1093 TYPE(
field_typ),
DIMENSION(Mesh%KMIN-Mesh%KP1:Mesh%KMAX+Mesh%KP1), &
1094 INTENT(OUT) :: shifted_field
1096 INTEGER :: i,j,k,local_joff
1098 local_joff = this%local_joff
1100 IF (mesh%shear_dir.EQ.1)
THEN
1102 DO k = mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
1104 DO j = mesh%JMIN,mesh%JMAX
1105 this%joff(j) = mesh%Q*mesh%OMEGA*mesh%bcenter(mesh%IMIN,j,k,2)* &
1107 this%jrem(j) = this%joff(j) - floor(this%joff(j))
1109 DO i = mesh%IMIN,mesh%IMAX
1110 shifted_field(k)%original(i,j-local_joff) = &
1111 (1.0-this%jrem(j))*field(1+modulo(i-1+floor(this%joff(j)), &
1112 mesh%IMAX-mesh%IMIN+1),j,k) + this%jrem(j)*field(1+modulo(i+ &
1113 floor(this%joff(j)),mesh%IMAX-mesh%IMIN+1),j,k)
1119 DO k = mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
1121 DO i = mesh%IMIN,mesh%IMAX
1122 this%joff(i) = -mesh%Q*mesh%OMEGA*mesh%bcenter(i,mesh%JMIN,k,1)* &
1124 this%jrem(i) = this%joff(i) - floor(this%joff(i))
1126 DO j = mesh%JMIN,mesh%JMAX
1128 DO i = mesh%IMIN,mesh%IMAX
1129 shifted_field(k)%original(i,j-local_joff) = &
1130 (1.0-this%jrem(i))*field(i,1+modulo(j-1+floor(this%joff(i)), &
1131 mesh%JMAX-mesh%JMIN+1),k) + this%jrem(i)*field(i,1+modulo(j+ &
1132 floor(this%joff(i)),mesh%JMAX-mesh%JMIN+1),k)
1148 CALL fftw_destroy_plan(this%plan_r2c)
1149 CALL fftw_destroy_plan(this%plan_c2r)
1150#if defined(PARALLEL)
1151 CALL fftw_free(this%fftw_real_pointer)
1152 CALL fftw_free(this%fftw_complex_pointer)
1153 DEALLOCATE(this%mpi_buf)
1155 DEALLOCATE(this%Fmass3D,this%Fsum3D,this%field)
1161 IF (
ASSOCIATED(this%Fmass3D_real))
DEALLOCATE(this%Fmass3D_real)
1164 this%kx, this%ky, this%Kxy2, &
1171 CALL this%Finalize_base()
Dictionary for generic data types.
base module for numerical flux functions
generic gravity terms module providing functionaly common to all gravity terms
Poisson solver using spectral methods within the shearingbox.
subroutine fft_backward(this, Mesh, Physics, field)
Calculates the FFT backward transformation.
real, parameter sqrttwopi
subroutine updategravity_single(this, Mesh, Physics, Fluxes, pvar, time, dt)
Calculates the acceleration from potential.
subroutine initgravity_sboxspectral(this, Mesh, Physics, config, IO)
Constructor of gravity sboxspectral module.
subroutine fft_forward(this, Mesh, Physics, field)
Calculates the FFT forward.
subroutine calcpotential(this, Mesh, Physics, time, pvar)
Computes the potential with FFT method within a shearingsheet.
subroutine fieldshift(this, Mesh, Physics, delt, field, shifted_field)
Shifts the whole field to the next periodic point.
subroutine setboundarydata(this, Mesh, Physics, delt)
Set ghost cell data for the potential.
subroutine calcdiskheight_single(this, Mesh, Physics, pvar, bccsound, h_ext, height)
Compute disk pressure scale height for geometrically thin self-gravitating shearingsheet.
2D poisson solver using spectral methods for direct integration
subroutine finalize(this)
Destructor of logging_base class.
integer, save default_mpi_real
default real type for MPI
base class for mesh arrays
derived class for compound of mesh arrays
integer, parameter south
named constant for southern boundary
subroutine setoutput(this, config, IO)
Setup mesh fields for i/o.
integer, parameter north
named constant for northern boundary
physics module for 1D,2D and 3D isothermal Euler equations