boundary_generic.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: boundary_generic.f90 #
5!# #
6!# Copyright (C) 2006-2016 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Manuel Jung <mjung@astrophysik.uni-kiel.de> #
9!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
10!# #
11!# This program is free software; you can redistribute it and/or modify #
12!# it under the terms of the GNU General Public License as published by #
13!# the Free Software Foundation; either version 2 of the License, or (at #
14!# your option) any later version. #
15!# #
16!# This program is distributed in the hope that it will be useful, but #
17!# WITHOUT ANY WARRANTY; without even the implied warranty of #
18!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19!# NON INFRINGEMENT. See the GNU General Public License for more #
20!# details. #
21!# #
22!# You should have received a copy of the GNU General Public License #
23!# along with this program; if not, write to the Free Software #
24!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25!# #
26!#############################################################################
34!----------------------------------------------------------------------------!
43!----------------------------------------------------------------------------!
61 USE common_dict
62 !--------------------------------------------------------------------------!
63 TYPE, PRIVATE :: boundary_p
64 CLASS(boundary_base), ALLOCATABLE :: p
65 END TYPE
66
69 TYPE(boundary_p) :: boundary(6)
70 CONTAINS
71 PROCEDURE :: initboundary
72 PROCEDURE :: centerboundary
73 PROCEDURE :: setcorneredges
74#ifdef PARALLEL
76 PROCEDURE :: initboundary_mpi
77 PROCEDURE :: mpibuffer2pvar
78 PROCEDURE :: mpipvar2buffer
79#endif
80 PROCEDURE :: finalize
81 END TYPE boundary_generic
82 !--------------------------------------------------------------------------!
83 PUBLIC :: &
84 ! types
86 ! constants
88 !--------------------------------------------------------------------------!
89
90CONTAINS
91
92 SUBROUTINE new_boundary(Boundary,Mesh,Physics,config,IO)
93 IMPLICIT NONE
94 !------------------------------------------------------------------------!
95 CLASS(boundary_generic), ALLOCATABLE :: boundary
96 CLASS(mesh_base), INTENT(INOUT) :: mesh
97 CLASS(physics_base), INTENT(IN) :: physics
98 TYPE(dict_typ), POINTER :: config,io
99 !------------------------------------------------------------------------!
100 ALLOCATE(boundary)
101 CALL boundary%InitBoundary(mesh,physics,config,io)
102 END SUBROUTINE
103
104 SUBROUTINE initboundary(this,Mesh,Physics,config,IO)
105 IMPLICIT NONE
106 !------------------------------------------------------------------------!
107 CLASS(boundary_generic), INTENT(INOUT) :: this
108 CLASS(mesh_base), INTENT(INOUT) :: mesh
109 CLASS(physics_base), INTENT(IN) :: physics
110 TYPE(dict_typ), POINTER :: config,io
111 INTEGER :: western, eastern, southern, northern, bottomer, topper
112 !------------------------------------------------------------------------!
113 INTEGER :: new(6)
114 INTEGER :: western_std, eastern_std, northern_std, southern_std, &
115 topper_std, bottomer_std
116 LOGICAL, DIMENSION(3) :: periods = .false.
117 INTEGER :: dir
118 !------------------------------------------------------------------------!
119 IF (.NOT.physics%Initialized().OR..NOT.mesh%Initialized()) &
120 CALL this%Error("InitBoundary","physics and/or mesh module uninitialized")
121
122 ! check for shearingsheet standard boundaries
123 western_std = 0
124 eastern_std = 0
125 northern_std = 0
126 southern_std = 0
127 bottomer_std = 0
128 topper_std = 0
129 IF (mesh%shear_dir.EQ.1) THEN
130 western_std = periodic
131 eastern_std = periodic
132 southern_std = shearing
133 northern_std = shearing
134 bottomer_std = reflecting
135 topper_std = reflecting
136 ELSE IF (mesh%shear_dir.EQ.2) THEN
137 western_std = shearing
138 eastern_std = shearing
139 southern_std = periodic
140 northern_std = periodic
141 bottomer_std = reflecting
142 topper_std = reflecting
143 END IF
144
145 CALL getattr(config, "western", western, western_std)
146 CALL getattr(config, "eastern", eastern, eastern_std)
147 CALL getattr(config, "southern", southern, northern_std)
148 CALL getattr(config, "northern", northern, southern_std)
149 CALL getattr(config, "bottomer", bottomer, bottomer_std)
150 CALL getattr(config, "topper", topper, topper_std)
151
152 new(west) = western
153 new(east) = eastern
154 new(south) = southern
155 new(north) = northern
156 new(bottom) = bottomer
157 new(top) = topper
158
159#ifdef PARALLEL
160 ! define inner connections, where boundaries are no true ones
161 IF (mesh%mycoords(1).NE.0) new(west) = none
162 IF (mesh%mycoords(1).NE.mesh%dims(1)-1) new(east) = none
163 IF (mesh%mycoords(2).NE.0) new(south) = none
164 IF (mesh%mycoords(2).NE.mesh%dims(2)-1) new(north) = none
165 IF (mesh%mycoords(3).NE.0) new(bottom) = none
166 IF (mesh%mycoords(3).NE.mesh%dims(3)-1) new(top) = none
167#endif
168
169 ! Check for correct shifting and boundaries
170 IF (western.EQ.shearing.AND.eastern.EQ.shearing) THEN
171 IF (.NOT.mesh%shear_dir.EQ.2) &
172 CALL this%Error("InitBoundary", &
173 "Please apply shifting in second dimension, when applying shearing boundaries at western/eastern.")
174#ifdef PARALLEL
175 CALL this%Error("InitBoundary", &
176 "Parallel mode is not allowed with shearing in West-East direction.")
177#endif
178 ELSE IF (southern.EQ.shearing.AND.northern.EQ.shearing) THEN
179 IF (.NOT.mesh%shear_dir.EQ.1) &
180 CALL this%Error("InitBoundary", &
181 "Please apply shifting in first dimension, when applying shearing boundaries at northern/southern.")
182 ELSE IF (bottomer.EQ.shearing.AND.topper.EQ.shearing) THEN
183 CALL this%Error("InitBoundary", &
184 "shifting in topper/bottomer direction not allowed.")
185 END IF
186
187 ! initialize every boundary
188 ! IMPORTANT: do this before anything else
189 DO dir=west,top
190 SELECT CASE(new(dir))
191 CASE(absorbing)
192 ALLOCATE(boundary_absorbing::this%Boundary(dir)%p)
193 CASE(axis)
194 ALLOCATE(boundary_axis::this%Boundary(dir)%p)
195 CASE(custom)
196 ALLOCATE(boundary_custom::this%Boundary(dir)%p)
197 CASE(farfield)
198 ALLOCATE(boundary_farfield::this%Boundary(dir)%p)
199 CASE(fixed)
200 ALLOCATE(boundary_fixed::this%Boundary(dir)%p)
201 CASE(no_gradients)
202 ALLOCATE(boundary_nogradients::this%Boundary(dir)%p)
203 CASE(noslip)
204 ALLOCATE(boundary_noslip::this%boundary(dir)%p)
205 CASE(periodic)
206 ALLOCATE(boundary_periodic::this%Boundary(dir)%p)
207 CASE(reflecting)
208 ALLOCATE(boundary_reflecting::this%Boundary(dir)%p)
209 CASE(shearing)
210 ALLOCATE(boundary_shearing::this%Boundary(dir)%p)
211#ifdef PARALLEL
212 CASE(none)
213 ALLOCATE(boundary_inner::this%Boundary(dir)%p)
214#endif
215 CASE DEFAULT
216 CALL this%Error("boundary_generic::InitBoundary","unknown boundary condition")
217 END SELECT
218
219 SELECT TYPE(obj => this%Boundary(dir)%p)
220 TYPE IS (boundary_absorbing)
221 CALL obj%InitBoundary_absorbing(mesh,physics,dir,config)
222 TYPE IS (boundary_axis)
223 CALL obj%InitBoundary_axis(mesh,physics,dir,config)
224 TYPE IS (boundary_custom)
225 CALL obj%InitBoundary_custom(mesh,physics,dir,config)
226 TYPE IS (boundary_farfield)
227 CALL obj%InitBoundary_farfield(mesh,physics,dir,config)
228 TYPE IS (boundary_fixed)
229 CALL obj%InitBoundary_fixed(mesh,physics,dir,config)
230 TYPE IS (boundary_nogradients)
231 CALL obj%InitBoundary_nogradients(mesh,physics,dir,config)
232 TYPE IS (boundary_noslip)
233 CALL obj%InitBoundary_noslip(mesh,physics,dir,config)
234 TYPE IS (boundary_periodic)
235 CALL obj%InitBoundary_periodic(mesh,physics,dir,config)
236 TYPE IS (boundary_reflecting)
237 CALL obj%InitBoundary_reflecting(mesh,physics,dir,config)
238 TYPE IS (boundary_shearing)
239 CALL obj%InitBoundary_shearing(mesh,physics,dir,config)
240#ifdef PARALLEL
241 TYPE IS (boundary_inner)
242 CALL obj%InitBoundary_inner(mesh,physics,dir,config)
243#endif
244 END SELECT
245 END DO
246
247 ! check periodicity
248 IF ((western.EQ.periodic.AND.eastern.EQ.periodic) .OR. &
249 (western.EQ.shearing.AND.eastern.EQ.shearing)) THEN
250 periods(1) = .true.
251 ELSE IF (western.EQ.periodic.NEQV.eastern.EQ.periodic) THEN
252 CALL this%boundary(west)%p%Error("InitBoundary", &
253 "Opposite boundary should be periodic.")
254 ELSE IF (western.EQ.shearing.NEQV.eastern.EQ.shearing) THEN
255 CALL this%boundary(west)%p%Error("InitBoundary", &
256 "Opposite boundary should be shearing.")
257 END IF
258 IF ((southern.EQ.periodic.AND.northern.EQ.periodic) .OR. &
259 (southern.EQ.shearing.AND.northern.EQ.shearing)) THEN
260 periods(2) = .true.
261 ELSE IF (southern.EQ.periodic.NEQV.northern.EQ.periodic) THEN
262 CALL this%boundary(south)%p%Error("InitBoundary", &
263 "Opposite boundary should be periodic.")
264 ELSE IF (southern.EQ.shearing.NEQV.northern.EQ.shearing) THEN
265 CALL this%boundary(south)%p%Error("InitBoundary", &
266 "Opposite boundary should be shearing.")
267 END IF
268 IF ((bottomer.EQ.periodic.AND.topper.EQ.periodic) .OR. &
269 (bottomer.EQ.shearing.AND.topper.EQ.shearing)) THEN
270 periods(3) = .true.
271 ELSE IF (bottomer.EQ.periodic.NEQV.topper.EQ.periodic) THEN
272 CALL this%boundary(bottom)%p%Error("InitBoundary", &
273 "Opposite boundary should be periodic.")
274 ELSE IF (bottomer.EQ.shearing.NEQV.topper.EQ.shearing) THEN
275 CALL this%boundary(bottom)%p%Error("InitBoundary", &
276 "Opposite boundary should be shearing.")
277 END IF
278
279#ifdef PARALLEL
280 CALL this%InitBoundary_MPI(mesh,physics,periods)
281#endif
282
283 END SUBROUTINE initboundary
284
285
287 SUBROUTINE centerboundary(this,Mesh,Physics,time,pvar,cvar)
290 IMPLICIT NONE
291 !------------------------------------------------------------------------!
292 CLASS(boundary_generic),INTENT(INOUT) :: this
293 CLASS(mesh_base), INTENT(IN) :: mesh
294 CLASS(physics_base), INTENT(IN) :: physics
295 REAL, INTENT(IN) :: time
296 CLASS(marray_compound), INTENT(INOUT) :: pvar, cvar
297 !------------------------------------------------------------------------!
298 CALL physics%Convert2Primitive(mesh%IMIN,mesh%IMAX,mesh%JMIN, &
299 mesh%JMAX,mesh%KMIN,mesh%KMAX,cvar,pvar)
300 SELECT TYPE(c => cvar)
302 SELECT TYPE(p => pvar)
304 p%fargo_transformation_applied = c%fargo_transformation_applied
305 END SELECT
306 END SELECT
307
308 this%err = 0
309 ! set physical boundary conditions at western and eastern boundaries
310 IF (mesh%INUM.GT.1) THEN
311 CALL this%Boundary(west)%p%SetBoundaryData(mesh,physics,time,pvar)
312 IF (this%err.GT.0) CALL this%Error("boundary_generic::CenterBoundary", &
313 "western boundary condition failed")
314 CALL this%Boundary(east)%p%SetBoundaryData(mesh,physics,time,pvar)
315 IF (this%err.GT.0) CALL this%Error("boundary_generic::CenterBoundary", &
316 "eastern boundary condition failed")
317 END IF
318
319 ! set physical boundary conditions at southern and northern boundaries
320 !
321 ! Here an extra case for parallel execution is handled, when shearing
322 ! boundaries are applied in northern/southern direction. The application
323 ! of the boundary conditions is done after the MPI communication, because
324 ! the shearing boundaries need to first apply periodic boundary conditions.
325 ! This is done MPI below across the computational domains.
326#ifdef PARALLEL
327 SELECT TYPE(bound1 => this%Boundary(south)%p)
328 TYPE IS (boundary_shearing)
329 ! Apply shearing boundaries after MPI communication. Compare with
330 ! code passage further below.
331 CLASS DEFAULT
332 IF (mesh%JNUM.GT.1) THEN
333 CALL this%Boundary(south)%p%SetBoundaryData(mesh,physics,time,pvar)
334 IF (this%err.GT.0) CALL this%Error("boundary_generic::CenterBoundary", &
335 "southern boundary condition failed")
336 END IF
337 END SELECT
338
339 SELECT TYPE(bound2 => this%Boundary(north)%p)
340 TYPE IS (boundary_shearing)
341 ! Apply shearing boundaries after MPI communication. Compare with
342 ! code passage further below.
343 CLASS DEFAULT
344 IF (mesh%JNUM.GT.1) THEN
345 CALL this%Boundary(north)%p%SetBoundaryData(mesh,physics,time,pvar)
346 IF (this%err.GT.0) CALL this%Error("boundary_generic::CenterBoundary", &
347 "northern boundary condition failed")
348 END IF
349 END SELECT
350#else
351 IF (mesh%JNUM.GT.1) THEN
352 CALL this%Boundary(south)%p%SetBoundaryData(mesh,physics,time,pvar)
353 IF (this%err.GT.0) CALL this%Error("boundary_generic::CenterBoundary", &
354 "southern boundary condition failed")
355 CALL this%Boundary(north)%p%SetBoundaryData(mesh,physics,time,pvar)
356 IF (this%err.GT.0) CALL this%Error("boundary_generic::CenterBoundary", &
357 "northern boundary condition failed")
358 END IF
359#endif
360
361 ! set physical boundary conditions at top and bottom boundaries
362 IF (mesh%KNUM.GT.1) THEN
363 CALL this%Boundary(bottom)%p%SetBoundaryData(mesh,physics,time,pvar)
364 IF (this%err.GT.0) CALL this%Error("boundary_generic::CenterBoundary", &
365 "bottom boundary condition failed")
366 CALL this%Boundary(top)%p%SetBoundaryData(mesh,physics,time,pvar)
367 IF (this%err.GT.0) CALL this%Error("boundary_generic::CenterBoundary", &
368 "top boundary condition failed")
369 END IF
370
371#ifdef PARALLEL
372 CALL mpiboundarycommunication(this,mesh,physics,pvar%data4d)
373#endif
374
375#ifdef PARALLEL
376 ! set physical boundary conditions at southern and northern boundaries
377 ! Here an extra case for parallel execution is handled. When shearing
378 ! boundaries are applied, periodic boundaries are assumed.
379 SELECT TYPE(bound1 => this%Boundary(south)%p)
380 TYPE IS (boundary_shearing)
381 IF (mesh%JNUM.GT.1) THEN
382 CALL this%Boundary(south)%p%SetBoundaryData(mesh,physics,time,pvar)
383 END IF
384 CLASS DEFAULT
385 ! do nothing
386 END SELECT
387 SELECT TYPE(bound2 => this%Boundary(north)%p)
388 TYPE IS (boundary_shearing)
389 IF (mesh%JNUM.GT.1) THEN
390 CALL this%Boundary(north)%p%SetBoundaryData(mesh,physics,time,pvar)
391 END IF
392 CLASS DEFAULT
393 ! do nothing
394 END SELECT
395#endif
396
397 CALL setcorneredges(this,mesh,physics,pvar%data4d)
398
399 ! convert primitive variables in ghost cells
400 IF (mesh%INUM.GT.1) THEN
401 CALL physics%Convert2Conservative(mesh%IGMIN,mesh%IMIN-mesh%IP1, &
402 mesh%JGMIN,mesh%JGMAX,mesh%KGMIN,mesh%KGMAX,pvar,cvar)
403 CALL physics%Convert2Conservative(mesh%IMAX+mesh%IP1,mesh%IGMAX, &
404 mesh%JGMIN,mesh%JGMAX,mesh%KGMIN,mesh%KGMAX,pvar,cvar)
405 END IF
406 IF (mesh%JNUM.GT.1) THEN
407 CALL physics%Convert2Conservative(mesh%IMIN,mesh%IMAX, &
408 mesh%JGMIN,mesh%JMIN-mesh%JP1,mesh%KGMIN,mesh%KGMAX,pvar,cvar)
409 CALL physics%Convert2Conservative(mesh%IMIN,mesh%IMAX, &
410 mesh%JMAX+mesh%JP1,mesh%JGMAX,mesh%KGMIN,mesh%KGMAX,pvar,cvar)
411 END IF
412 IF (mesh%KNUM.GT.1) THEN
413 CALL physics%Convert2Conservative(mesh%IMIN,mesh%IMAX, &
414 mesh%JMIN,mesh%JMAX,mesh%KGMIN,mesh%KMIN-mesh%KP1,pvar,cvar)
415 CALL physics%Convert2Conservative(mesh%IMIN,mesh%IMAX, &
416 mesh%JMIN,mesh%JMAX,mesh%KMAX+mesh%KP1,mesh%KGMAX,pvar,cvar)
417 END IF
418 END SUBROUTINE centerboundary
419
420
434 SUBROUTINE setcorneredges(this,Mesh,Physics,pvar)
435 IMPLICIT NONE
436 !------------------------------------------------------------------------!
437 CLASS(boundary_generic),INTENT(INOUT) :: this
438 CLASS(mesh_base), INTENT(IN) :: Mesh
439 CLASS(physics_base), INTENT(IN) :: Physics
440 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VNUM), &
441 INTENT(INOUT) :: pvar
442 !------------------------------------------------------------------------!
443 INTEGER :: i,j,k
444 !------------------------------------------------------------------------!
445 IF ((mesh%INUM.NE.1.AND.mesh%JNUM.NE.1.AND.mesh%KNUM.EQ.1).OR. &
446 (mesh%INUM.NE.1.AND.mesh%JNUM.NE.1.AND.mesh%KNUM.NE.1)) THEN
447
448 ! south-west
449#ifdef PARALLEL
450 IF(mesh%mycoords(1).EQ.0.AND.mesh%mycoords(2).EQ.0) THEN
451#endif
452!NEC$ SHORTLOOP
453 DO j=1,mesh%GNUM
454!NEC$ SHORTLOOP
455 DO i=j+1,mesh%GNUM
456 pvar(mesh%IMIN-i,mesh%JMIN-j,:,:) = pvar(mesh%IMIN-i,mesh%JMIN-j+1,:,:)
457 pvar(mesh%IMIN-j,mesh%JMIN-i,:,:) = pvar(mesh%IMIN-j+1,mesh%JMIN-i,:,:)
458 END DO
459 pvar(mesh%IMIN-j,mesh%JMIN-j,:,:) = 0.5 * (pvar(mesh%IMIN-j,mesh%JMIN,:,:) &
460 + pvar(mesh%IMIN,mesh%JMIN-j,:,:))
461 END DO
462#ifdef PARALLEL
463 END IF
464#endif
465 ! south-east
466#ifdef PARALLEL
467 IF(mesh%mycoords(1).EQ.mesh%dims(1)-1.AND.mesh%mycoords(2).EQ.0) THEN
468#endif
469!NEC$ SHORTLOOP
470 DO j=1,mesh%GNUM
471!NEC$ SHORTLOOP
472 DO i=j+1,mesh%GNUM
473 pvar(mesh%IMAX+i,mesh%JMIN-j,:,:) = pvar(mesh%IMAX+i,mesh%JMIN-j+1,:,:)
474 pvar(mesh%IMAX+j,mesh%JMIN-i,:,:) = pvar(mesh%IMAX+j-1,mesh%JMIN-i,:,:)
475 END DO
476 pvar(mesh%IMAX+j,mesh%JMIN-j,:,:) = 0.5 * (pvar(mesh%IMAX+j,mesh%JMIN,:,:) &
477 + pvar(mesh%IMAX,mesh%JMIN-j,:,:))
478 END DO
479#ifdef PARALLEL
480 END IF
481#endif
482 ! north-west
483#ifdef PARALLEL
484 IF(mesh%mycoords(1).EQ.0.AND.mesh%mycoords(2).EQ.mesh%dims(2)-1) THEN
485#endif
486!NEC$ SHORTLOOP
487 DO j=1,mesh%GNUM
488!NEC$ SHORTLOOP
489 DO i=j+1,mesh%GNUM
490 pvar(mesh%IMIN-i,mesh%JMAX+j,:,:) = pvar(mesh%IMIN-i,mesh%JMAX+j-1,:,:)
491 pvar(mesh%IMIN-j,mesh%JMAX+i,:,:) = pvar(mesh%IMIN-j+1,mesh%JMAX+i,:,:)
492 END DO
493 pvar(mesh%IMIN-j,mesh%JMAX+j,:,:) = 0.5 * (pvar(mesh%IMIN-j,mesh%JMAX,:,:) &
494 + pvar(mesh%IMIN,mesh%JMAX+j,:,:))
495 END DO
496#ifdef PARALLEL
497 END IF
498#endif
499
500 ! north-east
501#ifdef PARALLEL
502 IF(mesh%mycoords(1).EQ.mesh%dims(1)-1.AND.mesh%mycoords(2).EQ.mesh%dims(2)-1) THEN
503#endif
504!NEC$ SHORTLOOP
505 DO j=1,mesh%GNUM
506!NEC$ SHORTLOOP
507 DO i=j+1,mesh%GNUM
508 pvar(mesh%IMAX+i,mesh%JMAX+j,:,:) = pvar(mesh%IMAX+i,mesh%JMAX+j-1,:,:)
509 pvar(mesh%IMAX+j,mesh%JMAX+i,:,:) = pvar(mesh%IMAX+j-1,mesh%JMAX+i,:,:)
510 END DO
511 pvar(mesh%IMAX+j,mesh%JMAX+j,:,:) = 0.5 * (pvar(mesh%IMAX+j,mesh%JMAX,:,:) &
512 + pvar(mesh%IMAX,mesh%JMAX+j,:,:))
513 END DO
514#ifdef PARALLEL
515 END IF
516#endif
517 END IF
518
519 IF ((mesh%INUM.NE.1.AND.mesh%JNUM.EQ.1.AND.mesh%KNUM.NE.1).OR. &
520 (mesh%INUM.NE.1.AND.mesh%JNUM.NE.1.AND.mesh%KNUM.NE.1)) THEN
521
522 ! bottom-west
523#ifdef PARALLEL
524 IF(mesh%mycoords(1).EQ.0.AND.mesh%mycoords(3).EQ.0) THEN
525#endif
526!NEC$ SHORTLOOP
527 DO i=1,mesh%GNUM
528!NEC$ SHORTLOOP
529 DO k=i+1,mesh%GNUM
530 pvar(mesh%IMIN-i,:,mesh%KMIN-k,:) = pvar(mesh%IMIN-i+1,:,mesh%KMIN-k,:)
531 pvar(mesh%IMIN-k,:,mesh%KMIN-i,:) = pvar(mesh%IMIN-k,:,mesh%KMIN-i+1,:)
532 END DO
533 pvar(mesh%IMIN-i,:,mesh%KMIN-i,:) = 0.5 * (pvar(mesh%IMIN-i,:,mesh%KMIN,:) &
534 + pvar(mesh%IMIN,:,mesh%KMIN-i,:))
535 END DO
536#ifdef PARALLEL
537 END IF
538#endif
539 ! bottom-east
540#ifdef PARALLEL
541 IF(mesh%mycoords(1).EQ.mesh%dims(1)-1.AND.mesh%mycoords(3).EQ.0) THEN
542#endif
543!NEC$ SHORTLOOP
544 DO i=1,mesh%GNUM
545!NEC$ SHORTLOOP
546 DO k=i+1,mesh%GNUM
547 pvar(mesh%IMAX+i,:,mesh%KMIN-k,:) = pvar(mesh%IMAX+i-1,:,mesh%KMIN-k,:)
548 pvar(mesh%IMAX+k,:,mesh%KMIN-i,:) = pvar(mesh%IMAX+k,:,mesh%KMIN-i+1,:)
549 END DO
550 pvar(mesh%IMAX+i,:,mesh%KMIN-i,:) = 0.5 * (pvar(mesh%IMAX,:,mesh%KMIN-i,:) &
551 + pvar(mesh%IMAX+i,:,mesh%KMIN,:))
552 END DO
553#ifdef PARALLEL
554 END IF
555#endif
556
557 ! top-west
558#ifdef PARALLEL
559 IF(mesh%mycoords(1).EQ.0.AND.mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
560#endif
561!NEC$ SHORTLOOP
562 DO i=1,mesh%GNUM
563!NEC$ SHORTLOOP
564 DO k=i+1,mesh%GNUM
565 pvar(mesh%IMIN-i,:,mesh%KMAX+k,:) = pvar(mesh%IMIN-i+1,:,mesh%KMAX+k,:)
566 pvar(mesh%IMIN-k,:,mesh%KMAX+i,:) = pvar(mesh%IMIN-k,:,mesh%KMAX+i-1,:)
567 END DO
568 pvar(mesh%IMIN-i,:,mesh%KMAX+i,:) = 0.5 * (pvar(mesh%IMIN,:,mesh%KMAX+i,:) &
569 + pvar(mesh%IMIN-i,:,mesh%KMAX,:))
570 END DO
571#ifdef PARALLEL
572 END IF
573#endif
574
575 ! top-east
576#ifdef PARALLEL
577 IF(mesh%mycoords(1).EQ.mesh%dims(1)-1.AND.mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
578#endif
579!NEC$ SHORTLOOP
580 DO i=1,mesh%GNUM
581!NEC$ SHORTLOOP
582 DO k=i+1,mesh%GNUM
583 pvar(mesh%IMAX+i,:,mesh%KMAX+k,:) = pvar(mesh%IMAX+i-1,:,mesh%KMAX+k,:)
584 pvar(mesh%IMAX+k,:,mesh%KMAX+i,:) = pvar(mesh%IMAX+k,:,mesh%KMAX+i-1,:)
585 END DO
586 pvar(mesh%IMAX+i,:,mesh%KMAX+i,:) = 0.5 * (pvar(mesh%IMAX,:,mesh%KMAX+i,:) &
587 + pvar(mesh%IMAX+i,:,mesh%KMAX,:))
588 END DO
589#ifdef PARALLEL
590 END IF
591#endif
592 END IF
593
594 IF ((mesh%INUM.EQ.1.AND.mesh%JNUM.NE.1.AND.mesh%KNUM.NE.1).OR. &
595 (mesh%INUM.NE.1.AND.mesh%JNUM.NE.1.AND.mesh%KNUM.NE.1)) THEN
596
597 ! bottom-south
598#ifdef PARALLEL
599 IF(mesh%mycoords(2).EQ.0.AND.mesh%mycoords(3).EQ.0) THEN
600#endif
601!NEC$ SHORTLOOP
602 DO j=1,mesh%GNUM
603!NEC$ SHORTLOOP
604 DO k=j+1,mesh%GNUM
605 pvar(:,mesh%JMIN-j,mesh%KMIN-k,:) = pvar(:,mesh%JMIN-j+1,mesh%KMIN-k,:)
606 pvar(:,mesh%JMIN-k,mesh%KMIN-j,:) = pvar(:,mesh%JMIN-k,mesh%KMIN-j+1,:)
607 END DO
608 pvar(:,mesh%JMIN-j,mesh%KMIN-j,:) = 0.5 * (pvar(:,mesh%JMIN,mesh%KMIN-j,:) &
609 + pvar(:,mesh%JMIN-j,mesh%KMIN,:))
610 END DO
611#ifdef PARALLEL
612 END IF
613#endif
614
615 ! bottom-north
616#ifdef PARALLEL
617 IF(mesh%mycoords(2).EQ.mesh%dims(2)-1.AND.mesh%mycoords(3).EQ.0) THEN
618#endif
619!NEC$ SHORTLOOP
620 DO j=1,mesh%GNUM
621!NEC$ SHORTLOOP
622 DO k=j+1,mesh%GNUM
623 pvar(:,mesh%JMAX+j,mesh%KMIN-k,:) = pvar(:,mesh%JMAX+j-1,mesh%KMIN-k,:)
624 pvar(:,mesh%JMAX+k,mesh%KMIN-j,:) = pvar(:,mesh%JMAX+k,mesh%KMIN-j+1,:)
625 END DO
626!NEC$ IVDEP
627 pvar(:,mesh%JMAX+j,mesh%KMIN-j,:) = 0.5 * (pvar(:,mesh%JMAX,mesh%KMIN-j,:) &
628 + pvar(:,mesh%JMAX+j,mesh%KMIN,:))
629 END DO
630#ifdef PARALLEL
631 END IF
632#endif
633
634 ! top-south
635#ifdef PARALLEL
636 IF(mesh%mycoords(2).EQ.0.AND.mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
637#endif
638!NEC$ SHORTLOOP
639 DO j=1,mesh%GNUM
640!NEC$ SHORTLOOP
641 DO k=j+1,mesh%GNUM
642 pvar(:,mesh%JMIN-j,mesh%KMAX+k,:) = pvar(:,mesh%JMIN-j+1,mesh%KMAX+k,:)
643 pvar(:,mesh%JMIN-k,mesh%KMAX+j,:) = pvar(:,mesh%JMIN-k,mesh%KMAX+j-1,:)
644 END DO
645 pvar(:,mesh%JMIN-j,mesh%KMAX+j,:) = 0.5 * (pvar(:,mesh%JMIN,mesh%KMAX+j,:) &
646 + pvar(:,mesh%JMIN-j,mesh%KMAX,:))
647 END DO
648#ifdef PARALLEL
649 END IF
650#endif
651
652 ! top-north
653#ifdef PARALLEL
654 IF(mesh%mycoords(2).EQ.mesh%dims(2)-1.AND.mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
655#endif
656!NEC$ SHORTLOOP
657 DO j=1,mesh%GNUM
658!NEC$ SHORTLOOP
659 DO k=j+1,mesh%GNUM
660 pvar(:,mesh%JMAX+j,mesh%KMAX+k,:) = pvar(:,mesh%JMAX+j-1,mesh%KMAX+k,:)
661 pvar(:,mesh%JMAX+k,mesh%KMAX+j,:) = pvar(:,mesh%JMAX+k,mesh%KMAX+j-1,:)
662 END DO
663 pvar(:,mesh%JMAX+j,mesh%KMAX+j,:) = 0.5 * (pvar(:,mesh%JMAX,mesh%KMAX+j,:) &
664 + pvar(:,mesh%JMAX+j,mesh%KMAX,:))
665 END DO
666#ifdef PARALLEL
667 END IF
668#endif
669 END IF
670
671 ! Set corner cells (only 3D)
672 IF (mesh%INUM.NE.1.AND.mesh%JNUM.NE.1.AND.mesh%KNUM.NE.1) THEN
673#ifdef PARALLEL
674 IF(mesh%mycoords(1).EQ.0.AND.mesh%mycoords(2).EQ.0.AND.mesh%mycoords(3).EQ.0) THEN
675#endif
676!NEC$ SHORTLOOP
677 DO i=1,mesh%GNUM
678 pvar(mesh%IMIN-i,mesh%JMIN-i,mesh%KMIN-i,:) = (pvar(mesh%IMIN-i,mesh%JMIN,mesh%KMIN,:) &
679 + pvar(mesh%IMIN,mesh%JMIN-i,mesh%KMIN,:) + pvar(mesh%IMIN,mesh%JMIN,mesh%KMIN-i,:))/3.0
680 END DO
681#ifdef PARALLEL
682 END IF
683#endif
684#ifdef PARALLEL
685 IF(mesh%mycoords(1).EQ.mesh%dims(1)-1.AND.mesh%mycoords(2).EQ.0.AND.mesh%mycoords(3).EQ.0) THEN
686#endif
687!NEC$ SHORTLOOP
688 DO i=1,mesh%GNUM
689 pvar(mesh%IMAX+i,mesh%JMIN-i,mesh%KMIN-i,:) = (pvar(mesh%IMAX+i,mesh%JMIN,mesh%KMIN,:) &
690 + pvar(mesh%IMAX,mesh%JMIN-i,mesh%KMIN,:) + pvar(mesh%IMAX,mesh%JMIN,mesh%KMIN-i,:))/3.0
691 END DO
692#ifdef PARALLEL
693 END IF
694#endif
695#ifdef PARALLEL
696 IF(mesh%mycoords(1).EQ.0.AND.mesh%mycoords(2).EQ.mesh%dims(2)-1.AND.mesh%mycoords(3).EQ.0) THEN
697#endif
698!NEC$ SHORTLOOP
699 DO i=1,mesh%GNUM
700 pvar(mesh%IMIN-i,mesh%JMAX+i,mesh%KMIN-i,:) = (pvar(mesh%IMIN-i,mesh%JMAX,mesh%KMIN,:) &
701 + pvar(mesh%IMIN,mesh%JMAX+i,mesh%KMIN,:) + pvar(mesh%IMIN,mesh%JMAX,mesh%KMIN-i,:))/3.0
702 END DO
703#ifdef PARALLEL
704 END IF
705#endif
706#ifdef PARALLEL
707 IF(mesh%mycoords(1).EQ.mesh%dims(1)-1.AND.mesh%mycoords(2).EQ.mesh%dims(2)-1.AND.mesh%mycoords(3).EQ.0) THEN
708#endif
709!NEC$ SHORTLOOP
710 DO i=1,mesh%GNUM
711 pvar(mesh%IMAX+i,mesh%JMAX+i,mesh%KMIN-i,:) = (pvar(mesh%IMAX+i,mesh%JMAX,mesh%KMIN,:) &
712 + pvar(mesh%IMAX,mesh%JMAX+i,mesh%KMIN,:) + pvar(mesh%IMAX,mesh%JMAX,mesh%KMIN-i,:))/3.0
713 END DO
714#ifdef PARALLEL
715 END IF
716#endif
717#ifdef PARALLEL
718 IF(mesh%mycoords(1).EQ.0.AND.mesh%mycoords(2).EQ.0.AND.mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
719#endif
720!NEC$ SHORTLOOP
721 DO i=1,mesh%GNUM
722 pvar(mesh%IMIN-i,mesh%JMIN-i,mesh%KMAX+i,:) = (pvar(mesh%IMIN-i,mesh%JMIN,mesh%KMAX,:) &
723 + pvar(mesh%IMIN,mesh%JMIN-i,mesh%KMAX,:) + pvar(mesh%IMIN,mesh%JMIN,mesh%KMAX+i,:))/3.0
724 END DO
725#ifdef PARALLEL
726 END IF
727#endif
728#ifdef PARALLEL
729 IF(mesh%mycoords(1).EQ.mesh%dims(1)-1.AND.mesh%mycoords(2).EQ.0.AND.mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
730#endif
731!NEC$ SHORTLOOP
732 DO i=1,mesh%GNUM
733 pvar(mesh%IMAX+i,mesh%JMIN-i,mesh%KMAX+i,:) = (pvar(mesh%IMAX+i,mesh%JMIN,mesh%KMAX,:) &
734 + pvar(mesh%IMAX,mesh%JMIN-i,mesh%KMAX,:) + pvar(mesh%IMAX,mesh%JMIN,mesh%KMAX+i,:))/3.0
735 END DO
736#ifdef PARALLEL
737 END IF
738#endif
739#ifdef PARALLEL
740 IF(mesh%mycoords(1).EQ.0.AND.mesh%mycoords(2).EQ.mesh%dims(2)-1.AND.mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
741#endif
742!NEC$ SHORTLOOP
743 DO i=1,mesh%GNUM
744 pvar(mesh%IMIN-i,mesh%JMAX+i,mesh%KMAX+i,:) = (pvar(mesh%IMIN-i,mesh%JMAX,mesh%KMAX,:) &
745 + pvar(mesh%IMIN,mesh%JMAX+i,mesh%KMAX,:) + pvar(mesh%IMIN,mesh%JMAX,mesh%KMAX+i,:))/3.0
746 END DO
747#ifdef PARALLEL
748 END IF
749#endif
750#ifdef PARALLEL
751 IF(mesh%mycoords(1).EQ.mesh%dims(1)-1.AND.mesh%mycoords(2).EQ.mesh%dims(2)-1.AND.mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
752#endif
753!NEC$ SHORTLOOP
754 DO i=1,mesh%GNUM
755 pvar(mesh%IMAX+i,mesh%JMAX+i,mesh%KMAX+i,:) = (pvar(mesh%IMAX+i,mesh%JMAX,mesh%KMAX,:) &
756 + pvar(mesh%IMAX,mesh%JMAX+i,mesh%KMAX,:) + pvar(mesh%IMAX,mesh%JMAX,mesh%KMAX+i,:))/3.0
757 END DO
758#ifdef PARALLEL
759 END IF
760#endif
761 END IF
762
763 END SUBROUTINE
764
765
766#ifdef PARALLEL
767
768 SUBROUTINE initboundary_mpi(this,Mesh,Physics,periods)
769#ifdef HAVE_MPI_MOD
770 USE mpi
771#endif
772 IMPLICIT NONE
773#ifdef HAVE_MPIF_H
774 include 'mpif.h'
775#endif
776 !------------------------------------------------------------------------!
777 CLASS(boundary_generic),INTENT(INOUT) :: this
778 CLASS(mesh_base), INTENT(INOUT) :: Mesh
779 CLASS(physics_base), INTENT(IN) :: Physics
780 LOGICAL, DIMENSION(3), INTENT(IN) :: periods
781 !------------------------------------------------------------------------!
782 INTEGER :: comm_old
783 INTEGER :: ierr, dir
784 LOGICAL, DIMENSION(SIZE(Mesh%dims)) :: remain_dims = .false.
785 !------------------------------------------------------------------------!
786 ! create new cartesian communicator using Mesh%comm_cart
787 ! and account for the periodicity
788 ! IMPORTANT: disable reordering of nodes
789 comm_old = mesh%comm_cart
790 CALL mpi_cart_create(comm_old,SIZE(mesh%dims),mesh%dims,periods,.false.,mesh%comm_cart,ierr)
791
792 ! save ranks of neighbor processes
793 CALL mpi_cart_shift(mesh%comm_cart,0,1,mesh%neighbor(west),mesh%neighbor(east),ierr)
794 CALL mpi_cart_shift(mesh%comm_cart,1,1,mesh%neighbor(south),mesh%neighbor(north),ierr)
795 CALL mpi_cart_shift(mesh%comm_cart,2,1,mesh%neighbor(bottom),mesh%neighbor(top),ierr)
796
797 ! create communicators for every column and row of the cartesian
798 ! topology (used eg. for fargo shifts)
799 remain_dims = (/ .false., .true., .true. /)
800 CALL mpi_cart_sub(mesh%comm_cart,remain_dims,mesh%Icomm,ierr)
801 remain_dims = (/ .true., .false., .true. /)
802 CALL mpi_cart_sub(mesh%comm_cart,remain_dims,mesh%Jcomm,ierr)
803 remain_dims = (/ .true., .true., .false. /)
804 CALL mpi_cart_sub(mesh%comm_cart,remain_dims,mesh%Kcomm,ierr)
805
806 ! allocate memory for boundary data buffers
807 ALLOCATE(&
808 this%boundary(west)%p%sendbuf(mesh%GINUM,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
809 this%boundary(west)%p%recvbuf(mesh%GINUM,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
810 this%boundary(east)%p%sendbuf(mesh%GINUM,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
811 this%boundary(east)%p%recvbuf(mesh%GINUM,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
812 this%boundary(south)%p%sendbuf(mesh%IGMIN:mesh%IGMAX,mesh%GJNUM,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
813 this%boundary(south)%p%recvbuf(mesh%IGMIN:mesh%IGMAX,mesh%GJNUM,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
814 this%boundary(north)%p%sendbuf(mesh%IGMIN:mesh%IGMAX,mesh%GJNUM,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
815 this%boundary(north)%p%recvbuf(mesh%IGMIN:mesh%IGMAX,mesh%GJNUM,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
816 this%boundary(bottom)%p%sendbuf(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%GKNUM,physics%VNUM), &
817 this%boundary(bottom)%p%recvbuf(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%GKNUM,physics%VNUM), &
818 this%boundary(top)%p%sendbuf(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%GKNUM,physics%VNUM), &
819 this%boundary(top)%p%recvbuf(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%GKNUM,physics%VNUM), &
820 stat=ierr)
821 IF (ierr.NE.0) THEN
822 CALL this%boundary(west)%p%Error("boundary_generic::InitBoundary_MPI", &
823 "Unable to allocate memory for data buffers.")
824 END IF
825
826 ! initialize all buffers with 0
827 DO dir=west,top
828 this%boundary(dir)%p%recvbuf = 0.
829 this%boundary(dir)%p%sendbuf = 0.
830 END DO
831 END SUBROUTINE initboundary_mpi
832
834 SUBROUTINE mpiboundarycommunication(this,Mesh,Physics,pvar)
835#ifdef HAVE_MPI_MOD
836 USE mpi
837#endif
838 IMPLICIT NONE
839#ifdef HAVE_MPIF_H
840 include 'mpif.h'
841#endif
842 !------------------------------------------------------------------------!
843 CLASS(boundary_generic),INTENT(INOUT) :: this
844 CLASS(mesh_base), INTENT(IN) :: Mesh
845 CLASS(physics_base), INTENT(IN) :: Physics
846 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VNUM), &
847 INTENT(INOUT) :: pvar
848 !------------------------------------------------------------------------!
849 INTEGER :: ierr
850 INTEGER :: req(4)
851#ifdef MPI_USE_SENDRECV
852 INTEGER :: status(MPI_STATUS_SIZE)
853#else
854 INTEGER :: status(MPI_STATUS_SIZE,4)
855#endif
856 !------------------------------------------------------------------------!
857 ! NOTE: There are two different implementations of the MPI non-blocking
858 ! communication, one uses MPI_Sendrecv (default) the other MPI_Irecv and
859 ! MPI_Issend; to switch to the latter one has to remove -DMPI_USE_SENDRECV
860 ! from the compile command which can be done by setting -DMPI_USE_SENDRECV=OFF
861 ! when invoking cmake
862
863 ! western <-> eastern MPI communication
864#ifndef MPI_USE_SENDRECV
865 ! receive boundary data from eastern neighbor (non-blocking)
866 CALL mpi_irecv(this%Boundary(east)%p%recvbuf, &
867 mesh%GINUM*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
868 default_mpi_real,mesh%neighbor(east),10+west,mesh%comm_cart,req(1),ierr)
869#endif
870 ! fill send buffer if western neighbor exists
871 IF (mesh%neighbor(west).NE.mpi_proc_null) &
872 CALL this%MPIpvar2buffer(mesh,physics,west,pvar,this%Boundary(west)%p%sendbuf)
873#ifdef MPI_USE_SENDRECV
874 ! send boundary data to western and receive from eastern neighbor
875 CALL mpi_sendrecv(this%Boundary(west)%p%sendbuf, &
876 mesh%GINUM*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
877 default_mpi_real,mesh%neighbor(west),10+west,this%Boundary(east)%p%recvbuf, &
878 mesh%GINUM*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
879 default_mpi_real,mesh%neighbor(east),mpi_any_tag,mesh%comm_cart,status,ierr)
880#else
881 ! send boundary data to western neighbor
882 CALL mpi_issend(this%Boundary(west)%p%sendbuf, &
883 mesh%GINUM*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
884 default_mpi_real,mesh%neighbor(west),10+west,mesh%comm_cart,req(2),ierr)
885
886 ! receive boundary data from western neighbor
887 CALL mpi_irecv(this%Boundary(west)%p%recvbuf, &
888 mesh%GINUM*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
889 default_mpi_real,mesh%neighbor(west),10+east,mesh%comm_cart,req(3),ierr)
890#endif
891 ! fill send buffer if eastern neighbor exists
892 IF (mesh%neighbor(east).NE.mpi_proc_null) &
893 CALL this%MPIpvar2buffer(mesh,physics,east,pvar,this%Boundary(east)%p%sendbuf)
894#ifdef MPI_USE_SENDRECV
895 ! send boundary data to eastern and receive from western neighbor
896 CALL mpi_sendrecv(this%Boundary(east)%p%sendbuf, &
897 mesh%GINUM*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
898 default_mpi_real,mesh%neighbor(east),10+east,this%Boundary(west)%p%recvbuf, &
899 mesh%GINUM*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
900 default_mpi_real,mesh%neighbor(west),mpi_any_tag,mesh%comm_cart,status,ierr)
901#else
902 ! send boundary data to eastern neighbor
903 CALL mpi_issend(this%Boundary(east)%p%sendbuf, &
904 mesh%GINUM*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
905 default_mpi_real,mesh%neighbor(east),10+east,mesh%comm_cart,req(4),ierr)
906 ! wait for unfinished MPI communication
907 CALL mpi_waitall(4,req,status,ierr)
908#endif
909 ! copy data from receive buffers into ghosts cells
910 IF (mesh%neighbor(west).NE.mpi_proc_null) &
911 CALL this%MPIbuffer2pvar(mesh,physics,west,this%Boundary(west)%p%recvbuf,pvar)
912 IF (mesh%neighbor(east).NE.mpi_proc_null) &
913 CALL this%MPIbuffer2pvar(mesh,physics,east,this%Boundary(east)%p%recvbuf,pvar)
914
915 ! southern <-> northern MPI communication
916#ifndef MPI_USE_SENDRECV
917 ! receive boundary data from northern neighbor
918 CALL mpi_irecv(this%Boundary(north)%p%recvbuf, &
919 mesh%GJNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
920 default_mpi_real,mesh%neighbor(north),10+south,mesh%comm_cart,req(1),ierr)
921#endif
922 ! fill send buffer if southern neighbor exists
923 IF (mesh%neighbor(south).NE.mpi_proc_null) &
924 CALL this%MPIpvar2buffer(mesh,physics,south,pvar,this%Boundary(south)%p%sendbuf)
925#ifdef MPI_USE_SENDRECV
926 ! send boundary data to southern and receive from northern neighbor
927 CALL mpi_sendrecv(this%Boundary(south)%p%sendbuf, &
928 mesh%GJNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
929 default_mpi_real,mesh%neighbor(south),10+south,this%Boundary(north)%p%recvbuf, &
930 mesh%GJNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
931 default_mpi_real,mesh%neighbor(north),mpi_any_tag,mesh%comm_cart,status,ierr)
932#else
933 ! send boundary data to southern neighbor
934 CALL mpi_issend(this%Boundary(south)%p%sendbuf, &
935 mesh%GJNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
936 default_mpi_real,mesh%neighbor(south),10+south,mesh%comm_cart,req(2),ierr)
937
938 ! receive boundary data from southern neighbor
939 CALL mpi_irecv(this%Boundary(south)%p%recvbuf, &
940 mesh%GJNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
941 default_mpi_real,mesh%neighbor(south),10+north,mesh%comm_cart,req(3),ierr)
942#endif
943 ! fill send buffer if northern neighbor exists
944 IF (mesh%neighbor(north).NE.mpi_proc_null) &
945 CALL this%MPIpvar2buffer(mesh,physics,north,pvar,this%Boundary(north)%p%sendbuf)
946#ifdef MPI_USE_SENDRECV
947 ! send boundary data to northern and receive from southern neighbor
948 CALL mpi_sendrecv(this%Boundary(north)%p%sendbuf, &
949 mesh%GJNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
950 default_mpi_real,mesh%neighbor(north),10+north,this%Boundary(south)%p%recvbuf, &
951 mesh%GJNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
952 default_mpi_real,mesh%neighbor(south),mpi_any_tag,mesh%comm_cart,status,ierr)
953#else
954 ! send boundary data to northern neighbor
955 CALL mpi_issend(this%Boundary(north)%p%sendbuf, &
956 mesh%GJNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)*physics%VNUM, &
957 default_mpi_real,mesh%neighbor(north),10+north,mesh%comm_cart,req(4),ierr)
958 ! wait for unfinished MPI communication
959 CALL mpi_waitall(4,req,status,ierr)
960#endif
961 ! copy data from receive buffers into ghosts cells
962 IF (mesh%neighbor(south).NE.mpi_proc_null) &
963 CALL this%MPIbuffer2pvar(mesh,physics,south,this%Boundary(south)%p%recvbuf,pvar)
964 IF (mesh%neighbor(north).NE.mpi_proc_null) &
965 CALL this%MPIbuffer2pvar(mesh,physics,north,this%Boundary(north)%p%recvbuf,pvar)
966
967 ! bottom <-> top MPI communication
968#ifndef MPI_USE_SENDRECV
969 ! receive boundary data from top neighbor
970 CALL mpi_irecv(this%Boundary(top)%p%recvbuf, &
971 mesh%GKNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*physics%VNUM, &
972 default_mpi_real,mesh%neighbor(top),10+bottom,mesh%comm_cart,req(1),ierr)
973#endif
974 ! fill send buffer if bottomer neighbor exists
975 IF (mesh%neighbor(bottom).NE.mpi_proc_null) &
976 CALL this%MPIpvar2buffer(mesh,physics,bottom,pvar,this%Boundary(bottom)%p%sendbuf)
977#ifdef MPI_USE_SENDRECV
978 ! send boundary data to bottom and receive from top neighbor
979 CALL mpi_sendrecv(this%Boundary(bottom)%p%sendbuf,&
980 mesh%GKNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*physics%VNUM, &
981 default_mpi_real,mesh%neighbor(bottom),10+bottom,this%Boundary(top)%p%recvbuf, &
982 mesh%GKNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*physics%VNUM, &
983 default_mpi_real,mesh%neighbor(top),mpi_any_tag,mesh%comm_cart,status,ierr)
984#else
985 ! send boundary data to bottom neighbor
986 CALL mpi_issend(this%Boundary(bottom)%p%sendbuf, &
987 mesh%GKNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*physics%VNUM, &
988 default_mpi_real,mesh%neighbor(bottom),10+bottom,mesh%comm_cart,req(2),ierr)
989
990 ! receive boundary data from bottom neighbor
991 CALL mpi_irecv(this%Boundary(bottom)%p%recvbuf, &
992 mesh%GKNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*physics%VNUM, &
993 default_mpi_real,mesh%neighbor(bottom),10+top,mesh%comm_cart,req(3),ierr)
994#endif
995 ! send boundary data to top and receive from bottom neighbor
996 IF (mesh%neighbor(top).NE.mpi_proc_null) &
997 CALL this%MPIpvar2buffer(mesh,physics,top,pvar,this%Boundary(top)%p%sendbuf)
998#ifdef MPI_USE_SENDRECV
999 ! send boundary data to top and receive from bottom neighbor
1000 CALL mpi_sendrecv(this%Boundary(top)%p%sendbuf, &
1001 mesh%GKNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*physics%VNUM, &
1002 default_mpi_real,mesh%neighbor(top),10+top,this%Boundary(bottom)%p%recvbuf, &
1003 mesh%GKNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*physics%VNUM, &
1004 default_mpi_real,mesh%neighbor(bottom),mpi_any_tag,mesh%comm_cart,status,ierr)
1005#else
1006 ! send boundary data to top neighbor
1007 CALL mpi_issend(this%Boundary(top)%p%sendbuf, &
1008 mesh%GKNUM*(mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*physics%VNUM, &
1009 default_mpi_real,mesh%neighbor(top),10+top,mesh%comm_cart,req(4),ierr)
1010 ! wait for unfinished MPI communication
1011 CALL mpi_waitall(4,req,status,ierr)
1012#endif
1013 ! copy data from receive buffers into ghosts cells
1014 IF (mesh%neighbor(bottom).NE.mpi_proc_null) &
1015 CALL this%MPIbuffer2pvar(mesh,physics,bottom,this%Boundary(bottom)%p%recvbuf,pvar)
1016 IF (mesh%neighbor(top).NE.mpi_proc_null) &
1017 CALL this%MPIbuffer2pvar(mesh,physics,top,this%Boundary(top)%p%recvbuf,pvar)
1018 END SUBROUTINE mpiboundarycommunication
1019
1021 SUBROUTINE mpibuffer2pvar(this,Mesh,Physics,output_dir,buffer,pvar)
1022 IMPLICIT NONE
1023 !------------------------------------------------------------------------!
1024 CLASS(boundary_generic),INTENT(INOUT) :: this
1025 CLASS(mesh_base), INTENT(IN) :: Mesh
1026 CLASS(physics_base), INTENT(IN) :: Physics
1027 INTEGER, INTENT(IN) :: output_dir
1028 REAL, DIMENSION(:,:,:,:), INTENT(IN) :: buffer
1029 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VNUM), &
1030 INTENT(OUT) :: pvar
1031 !------------------------------------------------------------------------!
1032 INTEGER :: i,j,k,n
1033 !------------------------------------------------------------------------!
1034 SELECT CASE(output_dir)
1035 CASE(west)
1036!NEC$ shortloop
1037 DO n=1,physics%VNUM
1038 DO k=mesh%KGMIN,mesh%KGMAX
1039 DO j=mesh%JGMIN,mesh%JGMAX
1040!NEC$ shortloop
1041 DO i=1,mesh%GINUM
1042 pvar(mesh%IGMIN+i-1,j,k,n) = buffer(i,j-mesh%JGMIN+1,k-mesh%KGMIN+1,n)
1043 END DO
1044 END DO
1045 END DO
1046 END DO
1047 CASE(east)
1048!NEC$ shortloop
1049 DO n=1,physics%VNUM
1050 DO k=mesh%KGMIN,mesh%KGMAX
1051 DO j=mesh%JGMIN,mesh%JGMAX
1052!NEC$ shortloop
1053 DO i=1,mesh%GINUM
1054 pvar(mesh%IMAX+i,j,k,n) = buffer(i,j-mesh%JGMIN+1,k-mesh%KGMIN+1,n)
1055 END DO
1056 END DO
1057 END DO
1058 END DO
1059 CASE(south)
1060!NEC$ shortloop
1061 DO n=1,physics%VNUM
1062 DO k=mesh%KGMIN,mesh%KGMAX
1063!NEC$ shortloop
1064 DO j=1,mesh%GJNUM
1065 DO i=mesh%IGMIN,mesh%IGMAX
1066 pvar(i,mesh%JGMIN+j-1,k,n) = buffer(i-mesh%IGMIN+1,j,k-mesh%KGMIN+1,n)
1067 END DO
1068 END DO
1069 END DO
1070 END DO
1071 CASE(north)
1072!NEC$ shortloop
1073 DO n=1,physics%VNUM
1074 DO k=mesh%KGMIN,mesh%KGMAX
1075!NEC$ shortloop
1076 DO j=1,mesh%GJNUM
1077 DO i=mesh%IGMIN,mesh%IGMAX
1078 pvar(i,mesh%JMAX+j,k,n) = buffer(i-mesh%IGMIN+1,j,k-mesh%KGMIN+1,n)
1079 END DO
1080 END DO
1081 END DO
1082 END DO
1083 CASE(bottom)
1084!NEC$ shortloop
1085 DO n=1,physics%VNUM
1086!NEC$ shortloop
1087 DO k=1,mesh%GKNUM
1088 DO j=mesh%JGMIN,mesh%JGMAX
1089 DO i=mesh%IGMIN,mesh%IGMAX
1090 pvar(i,j,mesh%KGMIN+k-1,n) = buffer(i-mesh%IGMIN+1,j-mesh%JGMIN+1,k,n)
1091 END DO
1092 END DO
1093 END DO
1094 END DO
1095 CASE(top)
1096!NEC$ shortloop
1097 DO n=1,physics%VNUM
1098!NEC$ shortloop
1099 DO k=1,mesh%GKNUM
1100 DO j=mesh%JGMIN,mesh%JGMAX
1101 DO i=mesh%IGMIN,mesh%IGMAX
1102 pvar(i,j,mesh%KMAX+k,n) = buffer(i-mesh%IGMIN+1,j-mesh%JGMIN+1,k,n)
1103 END DO
1104 END DO
1105 END DO
1106 END DO
1107 END SELECT
1108 END SUBROUTINE mpibuffer2pvar
1109
1111 SUBROUTINE mpipvar2buffer(this,Mesh,Physics,input_dir,pvar,buffer)
1112 IMPLICIT NONE
1113 !------------------------------------------------------------------------!
1114 CLASS(boundary_generic),INTENT(INOUT) :: this
1115 CLASS(mesh_base), INTENT(IN) :: Mesh
1116 CLASS(physics_base), INTENT(IN) :: Physics
1117 INTEGER, INTENT(IN) :: input_dir
1118 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VNUM), &
1119 INTENT(IN) :: pvar
1120 REAL, DIMENSION(:,:,:,:), INTENT(OUT) :: buffer
1121 !------------------------------------------------------------------------!
1122 INTEGER :: i,j,k,n
1123 !------------------------------------------------------------------------!
1124 SELECT CASE(input_dir)
1125 CASE(west)
1126!NEC$ shortloop
1127 DO n=1,physics%VNUM
1128 DO k=mesh%KGMIN,mesh%KGMAX
1129 DO j=mesh%JGMIN,mesh%JGMAX
1130!NEC$ shortloop
1131 DO i=1,mesh%GINUM
1132 buffer(i,j-mesh%JGMIN+1,k-mesh%KGMIN+1,n) = pvar(mesh%IMIN+i-1,j,k,n)
1133 END DO
1134 END DO
1135 END DO
1136 END DO
1137 CASE(east)
1138!NEC$ shortloop
1139 DO n=1,physics%VNUM
1140 DO k=mesh%KGMIN,mesh%KGMAX
1141 DO j=mesh%JGMIN,mesh%JGMAX
1142!NEC$ shortloop
1143 DO i=1,mesh%GINUM
1144 buffer(i,j-mesh%JGMIN+1,k-mesh%KGMIN+1,n) = pvar(mesh%IMAX-mesh%GINUM+i,j,k,n)
1145 END DO
1146 END DO
1147 END DO
1148 END DO
1149 CASE(south)
1150!NEC$ shortloop
1151 DO n=1,physics%VNUM
1152 DO k=mesh%KGMIN,mesh%KGMAX
1153!NEC$ shortloop
1154 DO j=1,mesh%GJNUM
1155 DO i=mesh%IGMIN,mesh%IGMAX
1156 buffer(i-mesh%IGMIN+1,j,k-mesh%KGMIN+1,n) = pvar(i,mesh%JMIN+j-1,k,n)
1157 END DO
1158 END DO
1159 END DO
1160 END DO
1161 CASE(north)
1162!NEC$ shortloop
1163 DO n=1,physics%VNUM
1164 DO k=mesh%KGMIN,mesh%KGMAX
1165!NEC$ shortloop
1166 DO j=1,mesh%GJNUM
1167 DO i=mesh%IGMIN,mesh%IGMAX
1168 buffer(i-mesh%IGMIN+1,j,k-mesh%KGMIN+1,n) = pvar(i,mesh%JMAX-mesh%GJNUM+j,k,n)
1169 END DO
1170 END DO
1171 END DO
1172 END DO
1173 CASE(bottom)
1174!NEC$ shortloop
1175 DO n=1,physics%VNUM
1176!NEC$ shortloop
1177 DO k=1,mesh%GKNUM
1178 DO j=mesh%JGMIN,mesh%JGMAX
1179 DO i=mesh%IGMIN,mesh%IGMAX
1180 buffer(i-mesh%IGMIN+1,j-mesh%JGMIN+1,k,n) = pvar(i,j,mesh%KMIN+k-1,n)
1181 END DO
1182 END DO
1183 END DO
1184 END DO
1185 CASE(top)
1186!NEC$ shortloop
1187 DO n=1,physics%VNUM
1188!NEC$ shortloop
1189 DO k=1,mesh%GKNUM
1190 DO j=mesh%JGMIN,mesh%JGMAX
1191 DO i=mesh%IGMIN,mesh%IGMAX
1192 buffer(i-mesh%IGMIN+1,j-mesh%JGMIN+1,k,n) = pvar(i,j,mesh%KMAX-mesh%GKNUM+k,n)
1193 END DO
1194 END DO
1195 END DO
1196 END DO
1197 END SELECT
1198 END SUBROUTINE mpipvar2buffer
1199#endif
1200
1201 SUBROUTINE finalize(this)
1202 IMPLICIT NONE
1203 !------------------------------------------------------------------------!
1204 CLASS(boundary_generic), INTENT(INOUT) :: this
1205 !------------------------------------------------------------------------!
1206 INTEGER :: dir
1207 !------------------------------------------------------------------------!
1208 ! loop over all boundaries
1209 DO dir=1,6
1210 CALL this%boundary(dir)%p%Finalize()
1211#ifdef PARALLEL
1212 ! deallocate MPI send/recv buffers
1213 DEALLOCATE(this%boundary(dir)%p%sendbuf,this%boundary(dir)%p%recvbuf)
1214#endif
1215 END DO
1216 END SUBROUTINE finalize
1217
1218END MODULE boundary_generic_mod
Boundary module for absorbing (non-reflecting) conditions.
Boundary module for axis boundaries.
integer, parameter axis
axis
integer, parameter periodic
connects opposite boundaries
integer, parameter custom
user defined
integer, parameter none
integer, parameter no_gradients
copy data from last cell in comp. domain in ghost zones
integer, parameter farfield
uses far-field data and Riemann invariants
integer, parameter shearing
periodic with shear for shearing sheet/box
integer, parameter reflecting
reflecting, i.e. wall
integer, parameter absorbing
vanishing characteristic pseudo-variables for incomming waves
integer, parameter noslip
reflecting, but with moving wall
integer, parameter fixed
set fixed boundary data
subroutine initboundary(this, Mesh, Physics, bctype, bcname, dir, config)
Boundary module for custom conditions.
Boundary module for far field conditions.
Boundary module for fixed in/outflow conditions.
Generic boundary module.
subroutine initboundary_mpi(this, Mesh, Physics, periods)
initializes the MPI communication
subroutine, public new_boundary(Boundary, Mesh, Physics, config, IO)
subroutine, private centerboundary(this, Mesh, Physics, time, pvar, cvar)
Sets boundaries in all directions.
subroutine mpipvar2buffer(this, Mesh, Physics, input_dir, pvar, buffer)
Copys ghost cell data to buffer.
subroutine mpiboundarycommunication(this, Mesh, Physics, pvar)
Handles the MPI communication of the inner (and physical periodic) boundaries.
subroutine setcorneredges(this, Mesh, Physics, pvar)
Calculates the corner in 2D and the corners and edges in 3D.
subroutine mpibuffer2pvar(this, Mesh, Physics, output_dir, buffer, pvar)
Copys buffer data to ghost cells.
Boundary module for the inner boundaries (only necessary in parallel runs)
Boundary module for reflecting boundaries.
Boundary module for noslip conditions (see wikipedia )
Boundary module for periodic boundary conditions.
Boundary module for reflecting boundaries.
Boundary module for a shearingsheet/shearingbox. (see e.g. the standard run )
Dictionary for generic data types.
Definition: common_dict.f90:61
type(logging_base), save this
Basic fosite module.
subroutine finalize(this)
Destructor of logging_base class.
integer, save default_mpi_real
default real type for MPI
derived class for compound of mesh arrays
basic mesh module
Definition: mesh_base.f90:72
integer, parameter east
named constant for eastern boundary
Definition: mesh_base.f90:101
integer, parameter bottom
named constant for bottom boundary
Definition: mesh_base.f90:101
integer, parameter south
named constant for southern boundary
Definition: mesh_base.f90:101
integer, parameter top
named constant for top boundary
Definition: mesh_base.f90:101
integer, parameter north
named constant for northern boundary
Definition: mesh_base.f90:101
integer, parameter west
named constant for western boundary
Definition: mesh_base.f90:101
Basic physics module.
physics module for 1D,2D and 3D isothermal Euler equations
common data structure
mesh data structure
Definition: mesh_base.f90:122