gravity_sboxspectral.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: gravity_sboxspectral.f90 #
5!# #
6!# Copyright (C) 2015-2024 #
7!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
8!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
9!# #
10!# This program is free software; you can redistribute it and/or modify #
11!# it under the terms of the GNU General Public License as published by #
12!# the Free Software Foundation; either version 2 of the License, or (at #
13!# your option) any later version. #
14!# #
15!# This program is distributed in the hope that it will be useful, but #
16!# WITHOUT ANY WARRANTY; without even the implied warranty of #
17!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
18!# NON INFRINGEMENT. See the GNU General Public License for more #
19!# details. #
20!# #
21!# You should have received a copy of the GNU General Public License #
22!# along with this program; if not, write to the Free Software #
23!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
24!# #
25!#############################################################################
35!----------------------------------------------------------------------------!
57!----------------------------------------------------------------------------!
58MODULE gravity_sboxspectral_mod
59 USE gravity_base_mod
60 USE gravity_spectral_mod, ONLY : gravity_spectral
61 USE boundary_base_mod
62 USE fluxes_base_mod
63 USE physics_base_mod
64 USE mesh_base_mod
65 USE logging_base_mod
66 USE marray_base_mod
67 USE marray_compound_mod
68 USE common_dict
69 USE fftw
70#ifdef PARALLEL
71#ifdef HAVE_MPI_MOD
72 USE mpi
73#endif
74#endif
75 IMPLICIT NONE
76#ifdef PARALLEL
77#ifdef HAVE_MPIF_H
78 include 'mpif.h'
79#endif
80#endif
81 !--------------------------------------------------------------------------!
82 PRIVATE
83 REAL, PARAMETER :: sqrttwopi &
84 = 2.50662827463100050241576528481104525300698674
85 TYPE field_typ
86 REAL(c_double), DIMENSION(:,:), POINTER :: original
87 COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:,:), POINTER :: transformed
88 END TYPE field_typ
89
90 TYPE, EXTENDS(gravity_base) :: gravity_sboxspectral
91#ifdef HAVE_FFTW
92 TYPE(c_ptr) :: plan_r2c
93 TYPE(c_ptr) :: plan_c2r
94 TYPE(field_typ), DIMENSION(:), ALLOCATABLE & !< z-layered density/potential field
95 :: field
96 REAL(c_double), DIMENSION(:,:), POINTER, CONTIGUOUS &
97 :: kxy2 => null()
98 COMPLEX(C_DOUBLE_COMPLEX), POINTER, CONTIGUOUS &
99 :: fmass3d(:,:,:) => null(),&
100 fsum3d(:,:,:) => null(), &
102 qk(:,:,:) => null()
103 REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS &
104 :: fmass3d_real => null(), &
105 den_ip => null(), &
106 phi => null()
107 INTEGER(C_INTPTR_T) :: local_joff
108 REAL,DIMENSION(:), POINTER &
109 :: kx => null(),ky => null()
110 REAL :: shiftconst
111 REAL :: maxkxy2
112 REAL, DIMENSION(:), POINTER &
113 :: joff=>null(),jrem=>null()
114 INTEGER :: order
115#ifdef PARALLEL
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, &
120 fftw_complex_pointer
121#endif
122#endif
123
124 CONTAINS
125
126 PROCEDURE :: initgravity_sboxspectral
127 PROCEDURE :: updategravity_single
128! PROCEDURE :: InfoGravity
129 PROCEDURE :: setoutput
130 PROCEDURE :: calcdiskheight_single
131 final :: finalize
132#ifdef HAVE_FFTW
133 PROCEDURE :: calcpotential
134 PROCEDURE :: fft_forward
135 PROCEDURE :: fft_backward
136 PROCEDURE :: setboundarydata
137 PROCEDURE :: fieldshift
138#endif
139 END TYPE
140
141 !--------------------------------------------------------------------------!
142 PUBLIC :: &
143 ! types
144 gravity_sboxspectral
145 !--------------------------------------------------------------------------!
146
147 CONTAINS
148
156 SUBROUTINE initgravity_sboxspectral(this,Mesh,Physics,config,IO)
157 USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
158 IMPLICIT NONE
159 !------------------------------------------------------------------------!
160 CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
161 CLASS(mesh_base), INTENT(IN) :: Mesh
162 CLASS(physics_base), INTENT(IN) :: Physics
163 TYPE(dict_typ), POINTER :: config,IO
164 !------------------------------------------------------------------------!
165 INTEGER :: err, valwrite, i,j,k
166#if defined(HAVE_FFTW) && defined(PARALLEL)
167 INTEGER(C_INTPTR_T) :: C_INUM,C_JNUM
168#endif
169 !------------------------------------------------------------------------!
170 CALL this%InitGravity(mesh,physics,"shearingbox spectral solver",config,io)
171
172 !-------------- checks & warnings for initial conditions ----------------!
173#if !defined(HAVE_FFTW)
174 CALL this%Error("InitGravity_sboxspectral", &
175 "No fftw package could be loaded.")
176#else
177 CALL getattr(config, "order", this%order, 2)
178 SELECT CASE(this%order)
179 CASE(2)
180 ! ok, do nothing
181 CASE(4)
182 IF (mesh%KNUM.GT.1) &
183 CALL this%Error("InitGravity_sboxspectral", &
184 "order 4 is currently not supported in 3D")
185 CASE DEFAULT
186 CALL this%Error("InitGravity_sboxspectral", &
187 "Order must be one of {2,4}.")
188 END SELECT
189
190#ifdef PARALLEL
191 CALL fftw_mpi_init()
192 c_inum = mesh%INUM
193 c_jnum = mesh%JNUM
194#endif
195
196 ! rotational symmetry is not allowed
197 IF (mesh%ROTSYM.NE.0) &
198 CALL this%Error("InitGravity_sboxspectral", &
199 "Rotational symmetry not supported.")
200
201 ! 1D simulations are not allowed
202 IF (mesh%NDIMS.EQ.1) &
203 CALL this%Error("InitGravity_sboxspectral", &
204 "Only 2D (shearingsheet) and 3D (shearingbox) simulations allowed with this module.")
205
206 ! check physics
207 SELECT TYPE (phys => physics)
208 CLASS IS(physics_eulerisotherm)
209 ! do nothing
210 CLASS DEFAULT
211 CALL this%Error("InitGravity_sboxspectral", &
212 "Physics modules with density necessary for this module.")
213 END SELECT
214
215 !------------------------------------------------------------------------!
216
217#ifdef PARALLEL
218 ! check the dimensions if fftw should be used in parallel in order to
219 ! know how to allocate the local arrays
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.")
227 END IF
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")
235 END IF
236#endif
237
238 ! check for vertical extent (3D)
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")
244#ifdef PARALLEL
245 CALL this%Warning("InitGravity_sboxspectral","sboxspectral 3D is experimental")
246#endif
247 END IF
248
249 ALLOCATE( &
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), &
253#ifndef PARALLEL
254 ! allocate this field and set pointer references into z-slices (see below)
255 ! in parallel mode initialization handled by FFTW (see below)
256 this%Fmass3D(mesh%IMIN/2+1:mesh%IMAX/2+1,mesh%JMIN:mesh%JMAX,mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1), &
257#endif
258!!!! TODO: only allocate these for KNUM > 1
259 this%Fsum3D(mesh%IMIN/2+1:mesh%IMAX/2+1,mesh%JMIN:mesh%JMAX,mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1), &
260!!!! END TODO
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), &
265 stat=err)
266 IF (err.NE.0) &
267 CALL this%Error("InitGravity_sboxspectral","Memory allocation failed.")
268
269#ifdef HAVE_FFTW
270 this%local_joff = 0
271#endif
272
273 ! use special allocation pattern from fftw when using MPI in order to
274 ! assure good alignment
275#if defined(PARALLEL)
276 ALLOCATE(this%mpi_buf(mesh%IMIN:mesh%IMAX,mesh%GJNUM,mesh%KMIN:mesh%KMAX), &
277 stat=err)
278 IF (err.NE.0) &
279 CALL this%Error("InitGravity_sboxspectral","Memory allocation failed for mpi_buf.")
280
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])
290 END DO
291#endif
292
293 !------------------- create plans for fftw ------------------------------!
294 ! Pay attention to the argument order of the dimension (JNUM and INUM !
295 ! are switched because of C -> row-major, Fortran -> column-major), !
296 ! BUT ONLY in modern Fortran UNLIKE the legacy version !
297 ! ------------ plans are allocated in dictionary ------------------------!
298 CALL this%Info(" initializing FFTW: " // &
299#if !defined(PARALLEL)
300 "serial mode")
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)
304 END DO
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.")
313 END IF
314#elif defined(PARALLEL)
315 "parallel mode")
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)
322#endif
323 CALL getattr(config, "print_plans", valwrite, 0)
324 IF (valwrite.GT.0) THEN
325#ifndef NECSXAURORA
326 IF (this%GetRank().EQ.0) THEN
327 CALL fftw_print_plan(this%plan_r2c)
328 print *,achar(13)
329 CALL fftw_print_plan(this%plan_c2r)
330 print *,achar(13)
331 END IF
332#else
333 CALL this%Warning("InitGravity_sboxspectral", "fftw_print_plan currently not supported on SX-Aurora",0)
334#endif
335 END IF
336
337 !------------------------------------------------------------------------!
338 ! initialize arrays
339 this%phi(:,:,:) = 0.
340 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
341 this%field(k)%transformed(:,:) = cmplx(0.,0.)
342 ! in serial mode this also zeros this%Fmass3D
343 this%Fsum3D(:,:,k) = cmplx(0.,0)
344 this%qk(:,:,k) = 0.
345 END DO
346 this%kx(:) = 0.
347 this%ky(:) = 0.
348
349 ! precompute wave numbers
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)
354
355 ! cut-off value for wave numbers
356 this%maxKxy2 = 0.5*min(pi/mesh%dx,pi/mesh%dy)**2
357
358 SELECT CASE (mesh%shear_dir)
359 CASE(1)
360 this%shiftconst = mesh%Q*mesh%OMEGA*(mesh%ymax-mesh%ymin)/mesh%dx
361 ALLOCATE(&
362 this%joff(mesh%JMIN:mesh%JMAX), &
363 this%jrem(mesh%JMIN:mesh%JMAX), &
364 stat=err)
365 CASE(2)
366 this%shiftconst = mesh%Q*mesh%OMEGA*(mesh%xmax-mesh%xmin)/mesh%dy
367 ALLOCATE(&
368 this%joff(mesh%IMIN:mesh%IMAX), &
369 this%jrem(mesh%IMIN:mesh%IMAX), &
370 stat=err)
371 CASE DEFAULT
372 CALL this%Error("InitGravity_sboxspectral", &
373 "Shear direction must be one of WE_shear (x-direction) or SN_shear (y-direction).")
374 END SELECT
375 IF (err.NE.0) &
376 CALL this%Error("InitGravity_sboxspectral","Memory allocation failed for joff & jrem.")
377 this%joff(:) = 0.
378 this%jrem(:) = 0.
379
380 ! register data fields for output
381 CALL this%SetOutput(mesh,physics,config,io)
382#endif
383 END SUBROUTINE initgravity_sboxspectral
384
385 SUBROUTINE setoutput(this,Mesh,Physics,config,IO)
386 IMPLICIT NONE
387 !------------------------------------------------------------------------!
388 CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
389 CLASS(mesh_base), INTENT(IN) :: Mesh
390 CLASS(physics_base), INTENT(IN) :: Physics
391 TYPE(dict_typ), POINTER :: config,IO
392 !------------------------------------------------------------------------!
393 INTEGER :: valwrite,err
394 !------------------------------------------------------------------------!
395#ifdef HAVE_FFTW
396 valwrite = 0
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))
402 END IF
403 valwrite = 0
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), &
407 stat=err)
408 IF (err.NE.0) &
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))
413 END IF
414#endif
415 END SUBROUTINE setoutput
416
424 SUBROUTINE updategravity_single(this,Mesh,Physics,Fluxes,pvar,time,dt)
425 IMPLICIT NONE
426 !------------------------------------------------------------------------!
427 CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
428 CLASS(mesh_base), INTENT(IN) :: Mesh
429 CLASS(physics_base), INTENT(IN) :: Physics
430 CLASS(fluxes_base), INTENT(IN) :: Fluxes
431 CLASS(marray_compound), INTENT(INOUT) :: pvar
432 REAL, INTENT(IN) :: time,dt
433 !------------------------------------------------------------------------!
434 INTEGER :: i,j,k
435 REAL :: w1,w2
436 !------------------------------------------------------------------------!
437#ifdef HAVE_FFTW
438 ! calc potential first
439 CALL this%CalcPotential(mesh,physics,time,pvar)
440
441 IF (this%order.EQ.2) THEN
442!NEC$ ivdep
443 DO k = mesh%KMIN,mesh%KMAX
444!NEC$ ivdep
445 DO j = mesh%JMIN,mesh%JMAX
446!NEC$ ivdep
447 DO i = mesh%IMIN,mesh%IMAX
448 ! second order approximation
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))
456 END DO
457 END DO
458 END DO
459 ELSE ! this%order.EQ.4
460 w1 = 3./48.
461 w2 = 30./48.
462!NEC$ ivdep
463 DO k = mesh%KMIN,mesh%KMAX
464!NEC$ ivdep
465 DO j = mesh%JMIN,mesh%JMAX
466!NEC$ ivdep
467 DO i = mesh%IMIN,mesh%IMAX
468 ! fourth order
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))
479 END DO
480 END DO
481 END DO
482 END IF
483#endif
484 END SUBROUTINE updategravity_single
485
521#ifdef HAVE_FFTW
522 SUBROUTINE calcpotential(this,Mesh,Physics,time,pvar)
523 USE physics_eulerisotherm_mod, ONLY : statevector_eulerisotherm
524 IMPLICIT NONE
525 !------------------------------------------------------------------------!
526 CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
527 CLASS(mesh_base), INTENT(IN) :: Mesh
528 CLASS(physics_base), INTENT(IN) :: Physics
529 REAL, INTENT(IN) :: time
530 CLASS(marray_compound), INTENT(INOUT) :: pvar
531 !------------------------------------------------------------------------!
532 INTEGER :: i,j,k,kk
533 REAL :: delt,time0
534 !------------------------------------------------------------------------!
535 !---------------- fourier transformation of density ---------------------!
536 ! calculate the shift of the indice at time t !
537 ! shift density to pretend periodic behavior with interpolation !
538
539 IF (mesh%OMEGA .EQ. 0) THEN
540 time0 = 0.0
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))
549 END IF
550 delt = time - time0
551
552 !----------------- shift field to periodic point ------------------------!
553 SELECT TYPE(p => pvar)
554 CLASS IS(statevector_eulerisotherm)
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))
558 CLASS DEFAULT
559 CALL this%Error("gravity_sboxspectral::CalcPotential","unsupported state vector")
560 END SELECT
561
562 !-------------- fourier transform (forward) of shifted density ----------!
563! for performance checks
564#ifdef _FTRACE
565CALL ftrace_region_begin("forward_fft")
566#endif
567
571 CALL this%FFT_Forward(mesh,physics,this%field(mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1))
572
573! performance checks
574#ifdef _FTRACE
575CALL ftrace_region_end("forward_fft")
576#endif
577
578 !------------- calculate wave number vector squared --------------------!
579 IF (mesh%shear_dir.EQ.2) THEN
580!NEC$ IVDEP
581 DO j = mesh%JMIN,mesh%JMAX
582!NEC$ IVDEP
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
585 END DO
586 END DO
587 ELSE IF (mesh%shear_dir.EQ.1) THEN
588!NEC$ IVDEP
589 DO j = mesh%JMIN,mesh%JMAX
590!NEC$ IVDEP
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
593 END DO
594 END DO
595 END IF
596
597 !------------- calculate potential in fourier domain --------------------!
598 IF (abs(mesh%zmax-mesh%zmin).LT.tiny(mesh%zmin)) THEN
599 ! 2D razor thin disk
600 k=mesh%KMIN
601 DO j = mesh%JMIN,mesh%JMAX
602!NEC$ IVDEP
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))
607 ELSE
608 this%field(k)%transformed(i,j-this%local_joff) = 0.0
609 END IF
610 END DO
611 END DO
612! k=Mesh%KMIN
613! WHERE ((this%Kxy2(:,:).LT.this%maxKxy2).AND.(this%Kxy2(:,:).GT. 0))
614! this%Fmass3D(:,:,k) = -2.*PI*Physics%Constants%GN &
615! *this%Fmass3D(:,:,k)/SQRT(this%Kxy2(:,:))
616! ELSEWHERE
617! this%Fmass3D(:,:,k) = 0.0
618! END WHERE
619 ELSE
620 ! disk with vertical extent
621 this%qk(:,:,0) = exp(-0.5*mesh%dz*sqrt(this%Kxy2(:,:)))
622 ! compute weight factors for disks with more than one vertical cell
623 IF (mesh%KNUM.EQ.1) THEN
624 ! 2D disk with one zone with finite vertical size
625 k=mesh%KMIN
626!NEC$ IVDEP
627 DO j = mesh%JMIN,mesh%JMAX
628!NEC$ IVDEP
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)
633 ELSE
634 this%field(k)%transformed(i,j-this%local_joff) = 0.0
635 END IF
636 END DO
637 END DO
638 ELSE
639 ! 3D disk with more than one vertical zone
640 this%qk(:,:,1) = this%qk(:,:,0)*this%qk(:,:,0)
641 DO k=2,mesh%KNUM+1
642 this%qk(:,:,k) = this%qk(:,:,k-1) * this%qk(:,:,1)
643 END DO
644 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
645 this%Fsum3D(:,:,k) = 0.0
646 ! the following operations are equivalent to a matrix vector multiplication
647 ! where the matrix is a symmetric band matrix with zeros on the main diagonal
648 DO kk=mesh%KMIN-mesh%KP1,k-1
649!NEC$ IVDEP
650 DO j = mesh%JMIN,mesh%JMAX
651!NEC$ IVDEP
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)
655 END DO
656 END DO
657 END DO
658 ! skip k=kk
659!NEC$ IVDEP
660 DO j = mesh%JMIN,mesh%JMAX
661!NEC$ IVDEP
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)
666 END DO
667 END DO
668 END DO
669 END DO
670 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
671!NEC$ IVDEP
672 DO j = mesh%JMIN,mesh%JMAX
673!NEC$ IVDEP
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))
680 ELSE
681 this%field(k)%transformed(i,j-this%local_joff) = 0.0
682 END IF
683 END DO
684 END DO
685 END DO
686 END IF
687 END IF
688
689 !----------- fourier transform (backward) of shifted potential ---------!
690#ifdef _FTRACE
691CALL ftrace_region_begin("backward_fft")
692#endif
693
694 ! take the fourier transformed potential this%field(:)%transformed,
695 ! apply the inverse transforms and store the result in this%field(:)%original
696 CALL this%FFT_Backward(mesh,physics,this%field(mesh%KMIN-mesh%KP1:mesh%KMAX+mesh%KP1))
697
698#ifdef _FTRACE
699CALL ftrace_region_end("backward_fft")
700#endif
701
702 !------ calculate final potential with backshift and normalization ------!
703!NEC$ IVDEP
704 DO k = mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
705!NEC$ IVDEP
706 DO j = mesh%JMIN,mesh%JMAX
707!NEC$ IVDEP
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) ! no norm. by FFTW !
711 END DO
712 END DO
713 END DO
714 !-------------------------- shift field back ----------------------------!
715!!NEC$ IEXPAND
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))
719
720!NEC$ IVDEP
721 DO k = mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
722!NEC$ IVDEP
723 DO j = mesh%JMIN,mesh%JMAX
724!NEC$ IVDEP
725 DO i = mesh%IMIN,mesh%IMAX
726 this%phi(i,j,k) = this%field(k)%original(i,j-this%local_joff)
727 END DO
728 END DO
729 END DO
730
731 CALL this%SetBoundaryData(mesh,physics,delt)
732 END SUBROUTINE calcpotential
733
738 SUBROUTINE setboundarydata(this,Mesh,Physics,delt)
739 IMPLICIT NONE
740 !------------------------------------------------------------------------!
741 CLASS(gravity_sboxspectral), INTENT(INOUT):: this
742 CLASS(mesh_base), INTENT(IN) :: Mesh
743 CLASS(physics_base), INTENT(IN) :: Physics
744 REAL, INTENT(IN) :: delt
745 !------------------------------------------------------------------------!
746 INTEGER :: i,j,k
747 REAL :: joff,jrem
748#ifdef PARALLEL
749 INTEGER :: status(MPI_STATUS_SIZE),ierror
750#endif
751 !------------------------------------------------------------------------!
752 IF(mesh%shear_dir.EQ.2) THEN
753 ! southern / northern boundary: periodic
754 DO k = mesh%KMIN,mesh%KMAX
755!NEC$ SHORTLOOP
756 DO j = 1,mesh%GJNUM
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)
760 END DO
761 END DO
762 END DO
763
764 ! western / eastern boundary: shifted periodic
765 joff = -this%shiftconst*delt
766 jrem = joff - floor(joff)
767 DO k = mesh%KMIN,mesh%KMAX
768!NEC$ SHORTLOOP
769 DO i = 1,mesh%GINUM
770 DO j = mesh%JMIN,mesh%JMAX
771 ! copy eastern data
772 this%phi(mesh%IMIN-i,j,k) = this%phi(mesh%IMAX-i+1,j,k)
773 END DO
774 ! apply integer shift
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)
778 ! apply residual shift (interpolation)
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)
782 END DO
783 END DO
784 END DO
785
786 joff = this%shiftconst*delt
787 jrem = joff - floor(joff)
788 DO k = mesh%KMIN,mesh%KMAX
789!NEC$ SHORTLOOP
790 DO i = 1,mesh%GINUM
791 DO j = mesh%JMIN,mesh%JMAX
792 ! copy western data
793 this%phi(mesh%IMAX+i,j,k) = this%phi(mesh%IMIN+i-1,j,k)
794 END DO
795 ! apply integer shift
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,:)
799 ! apply residual shift (interpolation)
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)
803 END DO
804 END DO
805 END DO
806 ELSE IF(mesh%shear_dir.EQ.1) THEN
807 ! western / eastern boundary: periodic (no MPI communication)
808 DO k = mesh%KMIN,mesh%KMAX
809 DO j = mesh%JMIN,mesh%JMAX
810!NEC$ SHORTLOOP
811 DO i = 1,mesh%GJNUM
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)
814 END DO
815 END DO
816 END DO
817#ifdef PARALLEL
818 ! ATTENTION: send order is important for MPI:
819 ! 1. send non-shifted direction
820 ! 2. send shifted direction
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(&
825 this%mpi_buf,&
826 SIZE(this%mpi_buf), &
827 default_mpi_real, &
828 mesh%neighbor(north), 53+north, &
829 mesh%neighbor(south), mpi_any_tag, &
830 mesh%comm_cart, status, ierror)
831
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)
834
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)
837
838 CALL mpi_sendrecv_replace(&
839 this%mpi_buf,&
840 SIZE(this%mpi_buf), &
841 default_mpi_real, &
842 mesh%neighbor(south), 53+south, &
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)
847 ELSE
848!NEC$ SHORTLOOP
849 DO j = 1,mesh%GJNUM
850 ! southern northern (periodic in first step - further shift-treatment below)
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)
853 END DO
854 END IF
855#endif
856
857#ifdef PARALLEL
858 IF (mesh%mycoords(2).EQ.0) THEN
859#endif
860!NEC$ SHORTLOOP
861 DO j = 1,mesh%GJNUM
862 ! southern (shorn periodic) - residual and integer shift
863#ifndef PARALLEL
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)
865#endif
866 joff = this%shiftconst*delt
867 jrem = joff - floor(joff)
868 !------- integral shift ---------------------------------------------!
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))
871 !------- residual shift ---------------------------------------------!
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)
877 END DO
878 END DO
879 END DO
880#ifdef PARALLEL
881 END IF
882#endif
883
884#ifdef PARALLEL
885 IF (mesh%mycoords(2).EQ.mesh%dims(2)-1) THEN
886#endif
887!NEC$ SHORTLOOP
888 DO j = 1,mesh%GJNUM
889#ifndef PARALLEL
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)
891#endif
892 ! northern (shorn periodic) - residual and integer shift
893 joff = -this%shiftconst*delt
894 jrem = joff - floor(joff)
895 !------- integral shift ---------------------------------------------!
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))
898 !------- residual shift ---------------------------------------------!
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)
904 END DO
905 END DO
906 END DO
907#ifdef PARALLEL
908 END IF
909#endif
910 END IF
911 END SUBROUTINE setboundarydata
912
914 SUBROUTINE fft_forward(this,Mesh,Physics,field)
915 IMPLICIT NONE
916 !------------------------------------------------------------------------!
917 CLASS(gravity_sboxspectral), INTENT(INOUT):: this
918 CLASS(mesh_base), INTENT(IN) :: Mesh
919 CLASS(physics_base), INTENT(IN) :: Physics
920 TYPE(field_typ), DIMENSION(Mesh%KMIN-Mesh%KP1:Mesh%KMAX+Mesh%KP1), &
921 INTENT(INOUT) :: field
922 !------------------------------------------------------------------------!
923 INTEGER :: i,j,k
924 !------------------------------------------------------------------------!
925
926#ifdef _FTRACE
927CALL ftrace_region_begin("forward FFT")
928#endif
929
930 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
931#if defined(PARALLEL)
932 CALL fftw_mpi_execute_dft_r2c(this%plan_r2c,field(k)%original, &
933 this%field(k)%transformed)
934#else
935 CALL fftw_execute_dft_r2c(this%plan_r2c,this%field(k)%original, &
936 this%field(k)%transformed)
937#endif
938 END DO
939
940 IF (ASSOCIATED(this%Fmass3D_real)) THEN
941!NEC$ IVDEP
942 DO k = mesh%KMIN,mesh%KMAX
943!NEC$ IVDEP
944 DO j = mesh%JMIN,mesh%JMAX
945!NEC$ IVDEP
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)))
949 END DO
950 END DO
951 END DO
952 END IF
953
954#ifdef _FTRACE
955CALL ftrace_region_end("foward FFT")
956#endif
957
958 END SUBROUTINE
959
961 SUBROUTINE fft_backward(this,Mesh,Physics,field)
962 IMPLICIT NONE
963 !------------------------------------------------------------------------!
964 CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
965 CLASS(mesh_base), INTENT(IN) :: Mesh
966 CLASS(physics_base), INTENT(IN) :: Physics
967 TYPE(field_typ), DIMENSION(Mesh%KMIN-Mesh%KP1:Mesh%KMAX+Mesh%KP1), &
968 INTENT(INOUT) :: field
969 !------------------------------------------------------------------------!
970 INTEGER :: k
971 !------------------------------------------------------------------------!
972 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
973#if defined(PARALLEL)
974 CALL fftw_mpi_execute_dft_c2r(this%plan_c2r,field(k)%transformed, &
975 field(k)%original)
976#else
977 CALL fftw_execute_dft_c2r(this%plan_c2r,field(k)%transformed, &
978 field(k)%original)
979#endif
980 END DO
981 END SUBROUTINE
982#endif
983
984
1020 SUBROUTINE calcdiskheight_single(this,Mesh,Physics,pvar,bccsound,h_ext,height)
1021 IMPLICIT NONE
1022 !------------------------------------------------------------------------!
1023 CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
1024 CLASS(mesh_base), INTENT(IN) :: Mesh
1025 CLASS(physics_base), INTENT(IN) :: Physics
1026 CLASS(marray_compound), INTENT(INOUT) :: pvar
1027 TYPE(marray_base), INTENT(INOUT) :: bccsound,h_ext,height
1028 !------------------------------------------------------------------------!
1029 INTEGER :: i,j,k
1030 REAL :: cs2,p,q,rho_c
1031 !------------------------------------------------------------------------!
1032 ! pure self-gravitating shearing sheet with external point mass potential
1033 IF (abs(mesh%zmax-mesh%zmin).LT.tiny(mesh%zmin)) THEN
1034 ! 2D razor thin disk
1035 k = mesh%KMIN
1036!NEC$ IVDEP
1037 DO j=mesh%JGMIN,mesh%JGMAX
1038!NEC$ IVDEP
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) &
1042 /mesh%OMEGA**2.
1043 q = cs2/mesh%OMEGA**2.
1044 ! return the new disk height
1045 height%data3d(i,j,k) = p+sqrt(q+p*p)
1046 END DO
1047 END DO
1048 ELSE
1049 k = mesh%KNUM / 2
1050#ifdef PARALLEL
1051 IF (k.LT.mesh%KMIN.OR.k.GT.mesh%KMAX) &
1052 CALL this%Error("gravity_sboxspectral::CalcDiskHeight_single","not yet parallelized")
1053#endif
1054!NEC$ IVDEP
1055 DO j=mesh%JGMIN,mesh%JGMAX
1056!NEC$ IVDEP
1057 DO i=mesh%IGMIN,mesh%IGMAX
1058 ! check for even number grid cells in z-direction
1059 IF (mod(mesh%KNUM,2).EQ.0) THEN
1060 ! average density around midplane
1061 rho_c = 0.5*(pvar%data4d(i,j,k,physics%DENSITY)+pvar%data4d(i,j,k+1,physics%DENSITY))
1062 ELSE
1063 ! take midplane density
1064 rho_c = pvar%data4d(i,j,k+1,physics%DENSITY)
1065 END IF
1066 ! return the new disk height
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)
1069 END DO
1070 END DO
1071 END IF
1072 END SUBROUTINE calcdiskheight_single
1073
1083#ifdef HAVE_FFTW
1084 SUBROUTINE fieldshift(this,Mesh,Physics,delt,field,shifted_field)
1085 IMPLICIT NONE
1086 !------------------------------------------------------------------------!
1087 CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
1088 CLASS(mesh_base), INTENT(IN) :: Mesh
1089 CLASS(physics_base), INTENT(IN) :: Physics
1090 REAL, INTENT(IN) :: delt
1091 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
1092 INTENT(IN) :: field
1093 TYPE(field_typ), DIMENSION(Mesh%KMIN-Mesh%KP1:Mesh%KMAX+Mesh%KP1), &
1094 INTENT(OUT) :: shifted_field
1095 !------------------------------------------------------------------------!
1096 INTEGER :: i,j,k,local_joff
1097 !------------------------------------------------------------------------!
1098 local_joff = this%local_joff
1099
1100 IF (mesh%shear_dir.EQ.1) THEN
1101!NEC$ IVDEP
1102 DO k = mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
1103!NEC$ IVDEP
1104 DO j = mesh%JMIN,mesh%JMAX
1105 this%joff(j) = mesh%Q*mesh%OMEGA*mesh%bcenter(mesh%IMIN,j,k,2)* &
1106 delt/mesh%dx
1107 this%jrem(j) = this%joff(j) - floor(this%joff(j))
1108!NEC$ IVDEP
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)
1114 END DO
1115 END DO
1116 END DO
1117 ELSE ! must be Mesh%shear_dir.EQ.2, because otherwise initialization would raise an error
1118!NEC$ IVDEP
1119 DO k = mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
1120!NEC$ IVDEP
1121 DO i = mesh%IMIN,mesh%IMAX
1122 this%joff(i) = -mesh%Q*mesh%OMEGA*mesh%bcenter(i,mesh%JMIN,k,1)* &
1123 delt/mesh%dy
1124 this%jrem(i) = this%joff(i) - floor(this%joff(i))
1125 END DO
1126 DO j = mesh%JMIN,mesh%JMAX
1127!NEC$ IVDEP
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)
1133 END DO
1134 END DO
1135 END DO
1136 END IF
1137 END SUBROUTINE fieldshift
1138#endif
1139
1141 SUBROUTINE finalize(this)
1142 IMPLICIT NONE
1143 !------------------------------------------------------------------------!
1144 TYPE(gravity_sboxspectral), INTENT(INOUT) :: this
1145 !------------------------------------------------------------------------!
1146#ifdef HAVE_FFTW
1147 ! Destroy plans
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)
1154#else
1155 DEALLOCATE(this%Fmass3D,this%Fsum3D,this%field)
1156#endif
1157 ! Free all temporary memory of FFTW
1158 CALL fftw_cleanup()
1159
1160 ! Free memory
1161 IF (ASSOCIATED(this%Fmass3D_real)) DEALLOCATE(this%Fmass3D_real)
1162 DEALLOCATE(&
1163 this%phi, &
1164 this%kx, this%ky, this%Kxy2, &
1165 this%qk, &
1166 this%joff, &
1167 this%jrem, &
1168 this%den_ip &
1169 )
1170#endif
1171 CALL this%Finalize_base()
1172 END SUBROUTINE finalize
1173
1174END MODULE gravity_sboxspectral_mod
boundary_base_mod
Definition: boundary_base.f90:40
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
fftw
fftw module
Definition: fftw.f90:29
fluxes_base_mod
base module for numerical flux functions
Definition: fluxes_base.f90:39
gravity_base_mod
generic gravity terms module providing functionaly common to all gravity terms
Definition: gravity_base.f90:46
gravity_sboxspectral_mod
Poisson solver using spectral methods within the shearingbox.
Definition: gravity_sboxspectral.f90:58
gravity_sboxspectral_mod::fft_backward
subroutine fft_backward(this, Mesh, Physics, field)
Calculates the FFT backward transformation.
Definition: gravity_sboxspectral.f90:962
gravity_sboxspectral_mod::sqrttwopi
real, parameter sqrttwopi
Definition: gravity_sboxspectral.f90:83
gravity_sboxspectral_mod::updategravity_single
subroutine updategravity_single(this, Mesh, Physics, Fluxes, pvar, time, dt)
Calculates the acceleration from potential.
Definition: gravity_sboxspectral.f90:425
gravity_sboxspectral_mod::initgravity_sboxspectral
subroutine initgravity_sboxspectral(this, Mesh, Physics, config, IO)
Constructor of gravity sboxspectral module.
Definition: gravity_sboxspectral.f90:157
gravity_sboxspectral_mod::fft_forward
subroutine fft_forward(this, Mesh, Physics, field)
Calculates the FFT forward.
Definition: gravity_sboxspectral.f90:915
gravity_sboxspectral_mod::calcpotential
subroutine calcpotential(this, Mesh, Physics, time, pvar)
Computes the potential with FFT method within a shearingsheet.
Definition: gravity_sboxspectral.f90:523
gravity_sboxspectral_mod::fieldshift
subroutine fieldshift(this, Mesh, Physics, delt, field, shifted_field)
Shifts the whole field to the next periodic point.
Definition: gravity_sboxspectral.f90:1085
gravity_sboxspectral_mod::setboundarydata
subroutine setboundarydata(this, Mesh, Physics, delt)
Set ghost cell data for the potential.
Definition: gravity_sboxspectral.f90:739
gravity_sboxspectral_mod::calcdiskheight_single
subroutine calcdiskheight_single(this, Mesh, Physics, pvar, bccsound, h_ext, height)
Compute disk pressure scale height for geometrically thin self-gravitating shearingsheet.
Definition: gravity_sboxspectral.f90:1021
gravity_spectral_mod
2D poisson solver using spectral methods for direct integration
Definition: gravity_spectral.f90:45
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
logging_base_mod::finalize
subroutine finalize(this)
Destructor of logging_base class.
Definition: logging_base.f90:160
logging_base_mod::default_mpi_real
integer, save default_mpi_real
default real type for MPI
Definition: logging_base.f90:80
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
mesh_base_mod::south
integer, parameter south
named constant for southern boundary
Definition: mesh_base.f90:101
mesh_base_mod::setoutput
subroutine setoutput(this, config, IO)
Setup mesh fields for i/o.
Definition: mesh_base.f90:828
mesh_base_mod::north
integer, parameter north
named constant for northern boundary
Definition: mesh_base.f90:101
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
physics_eulerisotherm_mod
physics module for 1D,2D and 3D isothermal Euler equations
Definition: physics_eulerisotherm.f90:44
common_dict::dict_typ
Definition: common_dict.f90:123
fluxes_base_mod::fluxes_base
Definition: fluxes_base.f90:61
gravity_base_mod::gravity_base
Definition: gravity_base.f90:58
gravity_sboxspectral_mod::field_typ
Definition: gravity_sboxspectral.f90:85
gravity_sboxspectral_mod::gravity_sboxspectral
Definition: gravity_sboxspectral.f90:90
gravity_spectral_mod::gravity_spectral
Definition: gravity_spectral.f90:65
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:69
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
physics_base_mod::physics_base
Definition: physics_base.f90:69
physics_eulerisotherm_mod::physics_eulerisotherm
Definition: physics_eulerisotherm.f90:56
physics_eulerisotherm_mod::statevector_eulerisotherm
Definition: physics_eulerisotherm.f90:122
Generated by   doxygen 1.9.4