mesh_base.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: mesh_generic.f90 #
5!# #
6!# Copyright (C) 2006-2024 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# #
9!# This program is free software; you can redistribute it and/or modify #
10!# it under the terms of the GNU General Public License as published by #
11!# the Free Software Foundation; either version 2 of the License, or (at #
12!# your option) any later version. #
13!# #
14!# This program is distributed in the hope that it will be useful, but #
15!# WITHOUT ANY WARRANTY; without even the implied warranty of #
16!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
17!# NON INFRINGEMENT. See the GNU General Public License for more #
18!# details. #
19!# #
20!# You should have received a copy of the GNU General Public License #
21!# along with this program; if not, write to the Free Software #
22!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
23!# #
24!#############################################################################
29!----------------------------------------------------------------------------!
64
65!----------------------------------------------------------------------------!
71!----------------------------------------------------------------------------!
72MODULE mesh_base_mod
73 USE logging_base_mod
74 USE marray_base_mod
75 USE marray_cellscalar_mod
76 USE marray_cellvector_mod
77 USE geometry_base_mod
78 USE geometry_generic_mod
79 USE common_dict
80#ifdef PARALLEL
81#ifdef HAVE_MPI_MOD
82 USE mpi
83#endif
84#endif
85 IMPLICIT NONE
86#ifdef PARALLEL
87#ifdef HAVE_MPIF_H
88 include 'mpif.h'
89#endif
90#endif
91!--------------------------------------------------------------------------!
92 PRIVATE
95 INTEGER, PARAMETER :: midpoint = 1
96! INTEGER, PARAMETER :: TRAPEZOIDAL = 2 !< use trapezoidal rule to approximate flux integrals
97 !! #### parameters depending on dimensionality
98 INTEGER, PARAMETER :: nfaces(3) = (/ 2, 4, 6 /)
99 INTEGER, PARAMETER :: ncorners(3) = (/ 2, 4, 8 /)
101 INTEGER, PARAMETER :: &
102 west = 1, & !< named constant for western boundary
103 east = 2, &
104 south = 3, &
105 north = 4, &
106 bottom = 5, &
107 top = 6
109 INTEGER, PARAMETER :: &
110 vector_x = int(b'001'), &
111 vector_y = int(b'010'), &
112 vector_z = int(b'100')
114 TYPE, EXTENDS(logging_base) :: fargo_base
115 INTEGER :: direction = 0
116 CONTAINS
117 PROCEDURE :: getdirection
118 PROCEDURE :: getdirectionname
119 END TYPE fargo_base
121 !PRIVATE
122 TYPE,ABSTRACT, EXTENDS(logging_base) :: mesh_base
124 CLASS(geometry_base),ALLOCATABLE :: geometry
125 TYPE(fargo_base) :: fargo
126 TYPE(selection_base) :: without_ghost_zones
127 INTEGER :: gnum
128 INTEGER :: ginum,gjnum,gknum
129 INTEGER :: inum,jnum,knum
130 INTEGER :: imin,imax
131 INTEGER :: jmin,jmax
132 INTEGER :: kmin,kmax
133 INTEGER :: igmin,igmax
134 INTEGER :: jgmin,jgmax
135 INTEGER :: kgmin,kgmax
136 INTEGER :: ndims
137 INTEGER :: nfaces
138 INTEGER :: ncorners
139 INTEGER :: ip1,ip2,im1,im2
140 INTEGER :: jp1,jp2,jm1,jm2
141 INTEGER :: kp1,kp2,km1,km2
142 REAL :: xmin, xmax
143 REAL :: ymin, ymax
144 REAL :: zmin, zmax
145 REAL :: dx,dy,dz
146 REAL :: invdx,invdy,invdz
147 REAL :: omega
148 REAL :: q
149 INTEGER :: shear_dir
150 INTEGER :: rotsym
151 INTEGER :: vector_components
152 REAL :: rotcent(3)
155 TYPE(marray_cellvector) :: curv, & !< curvilinear coordinates
156 cart
157 REAL, DIMENSION(:,:,:,:), POINTER :: &
158 center => null(),&
159 bcenter=> null(),&
160 bccart => null()
161 REAL, DIMENSION(:,:,:,:,:), POINTER :: &
162 fccart=> null(), &
163 ccart => null()
166 TYPE(marray_base) :: dlx,dly,dlz, & !< cell centered line elements
167 volume, & !< cell volumes
168 dxdydv,dydzdv,dzdxdv
169 REAL, DIMENSION(:,:,:,:), POINTER :: &
170 dax,day,daz, & !< cell surface area elements
171 daxdydz,daydzdx,dazdxdy
174 TYPE(marray_cellscalar) :: hx,hy,hz, & !< scale factors
175 sqrtg, & !< sqrt(det(metric))
176 cyxy,cyzy,cxzx,cxyx,czxz,czyz
179 TYPE(marray_cellscalar) :: radius
180 TYPE(marray_cellvector) :: posvec
183 REAL, DIMENSION(:,:,:), POINTER :: &
184 rotation => null()
185 REAL, DIMENSION(:,:,:,:,:,:), POINTER :: &
186 weights => null()
187#ifdef PARALLEL
188
189 INTEGER :: maxinum,maxjnum,maxknum
190 INTEGER :: mininum,minjnum,minknum
191 INTEGER :: comm_cart
192 INTEGER :: icomm,jcomm,kcomm
193 ! always setup a 3D cartesian process topology, even for 1D/2D simulations
194 INTEGER, DIMENSION(NFACES(3)) :: comm_boundaries
195 INTEGER, DIMENSION(NFACES(3)) :: rank0_boundaries
196 INTEGER, DIMENSION(NFACES(3)) :: neighbor
197 INTEGER, DIMENSION(3) :: dims
198 INTEGER, DIMENSION(3) :: mycoords
199#endif
200 CONTAINS
201 PROCEDURE :: remapbounds_1
202 PROCEDURE :: remapbounds_2
203 PROCEDURE :: remapbounds_3
204 PROCEDURE :: remapbounds_4
205 generic :: remapbounds => remapbounds_1, remapbounds_2, remapbounds_3, remapbounds_4
206 PROCEDURE :: initmesh
207 PROCEDURE :: finalize_base
208 procedure(finalize), DEFERRED :: finalize
209 PROCEDURE :: internalpoint
210 procedure(tensordivergence3d), DEFERRED :: tensordivergence3d
211 procedure(vectordivergence3d), DEFERRED :: vectordivergence3d
212 procedure(tensordivergence2d_1), DEFERRED :: tensordivergence2d_1
213 procedure(vectordivergence2d_1), DEFERRED :: vectordivergence2d_1
214 generic :: divergence => tensordivergence3d, vectordivergence3d, &
215 vectordivergence2d_1, &! VectorDivergence2D_2, &
216 tensordivergence2d_1 !, TensorDivergence2D_2, &
217 END TYPE mesh_base
218 abstract INTERFACE
219 PURE SUBROUTINE tensordivergence3d(this,Txx,Txy,Txz,Tyx,Tyy,Tyz,Tzx,Tzy,Tzz, &
220 divTx,divTy,divTz)
221 IMPORT mesh_base
222 IMPLICIT NONE
223 !------------------------------------------------------------------------!
224 CLASS(mesh_base), INTENT(IN) :: this
225 REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX), &
226 INTENT(IN) :: txx,txy,txz,tyx,tyy,tyz,tzx,tzy,tzz
227 REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX), &
228 INTENT(OUT) :: divtx,divty,divtz
229 END SUBROUTINE
230 PURE SUBROUTINE vectordivergence3d(this,vx,vy,vz,divv)
231 IMPORT mesh_base
232 IMPLICIT NONE
233 !------------------------------------------------------------------------!
234 CLASS(mesh_base),INTENT(IN) :: this
235 REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
236 :: vx,vy,vz,divv
237 !------------------------------------------------------------------------!
238 INTENT(IN) :: vx,vy,vz
239 INTENT(OUT) :: divv
240 END SUBROUTINE
241 PURE SUBROUTINE vectordivergence2d_1(this,vx,vy,divv)
242 IMPORT mesh_base
243 IMPLICIT NONE
244 !------------------------------------------------------------------------!
245 CLASS(mesh_base),INTENT(IN) :: this
246 REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
247 :: vx,vy,divv
248 !------------------------------------------------------------------------!
249 INTENT(IN) :: vx,vy
250 INTENT(OUT) :: divv
251 END SUBROUTINE
252 PURE SUBROUTINE tensordivergence2d_1(this,Txx,Txy,Tyx,Tyy,divTx,divTy)
253 IMPORT mesh_base
254 IMPLICIT NONE
255 !------------------------------------------------------------------------!
256 CLASS(mesh_base),INTENT(IN) :: this
257 REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
258 :: txx,txy,tyx,tyy,divtx,divty
259 !------------------------------------------------------------------------!
260 INTENT(IN) :: txx,txy,tyx,tyy
261 INTENT(OUT) :: divtx,divty
262 END SUBROUTINE
263 SUBROUTINE finalize(this)
264 IMPORT mesh_base
265 IMPLICIT NONE
266 !------------------------------------------------------------------------!
267 CLASS(mesh_base), INTENT(INOUT) :: this
268 END SUBROUTINE
269 END INTERFACE
270 !--------------------------------------------------------------------------!
271 PUBLIC :: &
272 ! types
273 mesh_base, fargo_base, &
274 ! constants
275 pi, &
276#ifdef PARALLEL
277 default_mpi_real, &
278#endif
279 midpoint, &
280 cartesian, &
281 cylindrical, logcylindrical, tancylindrical, &
282 spherical, logspherical, spherical_planet, &
283 vector_x, vector_y, vector_z, &
284 west, east, south, north, bottom, top
285 !--------------------------------------------------------------------------!
286
287CONTAINS
288
295 SUBROUTINE initmesh(this,config,IO,mtype,mname)
296 IMPLICIT NONE
297 !------------------------------------------------------------------------!
298 CLASS(mesh_base) :: this
299 TYPE(dict_typ),POINTER :: config
301 TYPE(dict_typ),POINTER :: IO
302 INTEGER :: mtype
303 CHARACTER(LEN=32) :: mname
304 !------------------------------------------------------------------------!
305 CHARACTER(LEN=32) :: xres,yres,zres,somega,fargo_method(4)
306 INTEGER :: fargo,fargo_dir
307 INTEGER :: i,j,k
308 REAL :: mesh_dx,mesh_dy,mesh_dz
309#ifdef PARALLEL
310 INTEGER :: inum,jnum,knum,err
311#endif
312 REAL, DIMENSION(6,3) :: cfaces
313 REAL, DIMENSION(8,3) :: ccorners
314 !------------------------------------------------------------------------!
315 INTENT(INOUT) :: this
316 !------------------------------------------------------------------------!
317 CALL this%logging_base%InitLogging(mtype,mname)
318
319 CALL new_geometry(this%geometry, config)
320
321 ! check whether shearing box simulation is enabled and along which
322 ! coordinate direction the shear should be applied
323 ! default: disable shearing box
324 CALL getattr(config,"shearingbox",this%shear_dir,0)
325 IF (this%shear_dir.LT.0.OR.this%shear_dir.GT.2) &
326 CALL this%Error("mesh_base::InitMesh","direction of the shear should be one of {1,2} or 0")
327
328 ! Here the angular velocity and the center of rotation are defined, if
329 ! the mesh is in a rotating frame of reference. The fictious forces must
330 ! be added through one of the following methods:
331 ! 1. rotframe external source module
332 ! 2. special physics module with angular momentum transport *IAMT, *IAMROT
333 ! 3. special rhstype with angular momentum conservation
334 ! Only the IAMROT physics module and rotframe source module allow
335 ! for a different center of rotation than (/0, 0, 0/)
336
337 ! angular velocity of the rotating reference frame
338 CALL getattr(config, "omega", this%omega, 0.0)
339
340 ! shearing force in shearingsheet (default: Keplerian rotation)
341 CALL getattr(config, "Qshear", this%Q, 1.5)
342
343 ! center of rotation
344 this%rotcent = (/ 0., 0., 0./)
345 CALL getattr(config, "rotation_center", this%rotcent, this%rotcent)
346
347
348 ! total resolution
349 ! IMPORTANT: The resolution is the key value in order to determine the
350 ! used dimensions below
351 CALL getattr(config, "inum", this%INUM)
352 CALL getattr(config, "jnum", this%JNUM)
353 CALL getattr(config, "knum", this%KNUM)
354
355 ! do some sanity checks
356 IF (any((/ this%INUM.LE.0, this%JNUM.LE.0, this%KNUM.LE.0 /))) &
357 CALL this%Error("InitMesh","zero or negative resolution not allowed.")
358 IF (all((/ this%INUM.EQ.1, this%JNUM.EQ.1, this%KNUM.EQ.1 /))) &
359 CALL this%Error("InitMesh","at least one of inum,jnum,knum must be > 1.")
360
361 ! coordinate domain
362 CALL getattr(config, "xmin", this%xmin)
363 CALL getattr(config, "xmax", this%xmax)
364 CALL getattr(config, "ymin", this%ymin)
365 CALL getattr(config, "ymax", this%ymax)
366 CALL getattr(config, "zmin", this%zmin)
367 CALL getattr(config, "zmax", this%zmax)
368
371
372 ! check for fargo transport
373 ! default: enable fargo in shearing sheet/box, i.e. if shear_dir is set to something >0
374 ! disable fargo otherwise
375 fargo = 0
376 IF(this%shear_dir.GT.0) fargo = 3
377 ! check if user request fargo transport
378 CALL getattr(config, "fargo/method", fargo, fargo)
379 SELECT CASE(fargo)
380 CASE(0,1,2,3)
381 ! 0 = disable fargo (use this for shearing sheet/box simulations without fargo)
382 ! 1 = calculated mean background velocity
383 ! 2 = fixed user supplied background velocity
384 ! 3 = shearingsheet/box fixed background velocity
385 fargo_method = [CHARACTER(LEN=32) :: "disabled", "dynamic velocity",&
386 "user supplied fixed, velocity","shearingsheet/box shear velocity" ]
387 CALL this%fargo%logging_base%InitLogging(fargo,fargo_method(fargo+1))
388 IF (fargo.GT.0) THEN
389 ! set/check fargo transport direction
390 fargo_dir = this%Geometry%GetAzimuthIndex() ! fargo transport direction corresponds to azimuthal direction if available
391 SELECT CASE(fargo_dir)
392 CASE(0) ! no azimuthal coordinate (probably cartesian geometry)
393 SELECT TYPE(mgeo => this%Geometry)
394 TYPE IS(geometry_cartesian) ! fargo transport could be along any cartesian direction
395 IF (this%fargo%GetType().EQ.3) THEN
396 ! in the shearing sheet the fargo direction is given by the shearing direction
397 ! the user supplied fargo direction from the dictionary is ignored
398 this%fargo%direction = this%shear_dir
399 ELSE
400 ! read user supplied direction from dictionary
401 CALL getattr(config, "fargo/direction", this%fargo%direction)
402 IF (this%fargo%direction.GT.3.OR.this%fargo%direction.LT.1) &
403 CALL this%Error("mesh_base::InitMesh","invalid fargo transport direction, should be one of {1,2,3}")
404 END IF
405 CLASS DEFAULT
406 CALL this%Error("mesh_base::InitMesh","fargo transport not supported for this geometry")
407 END SELECT
408 CASE(1,2,3)
409 ! curvilinear geometry with azimuthal (i.e. periodic) coordinate;
410 ! the user supplied fargo direction from the dictionary is ignored
411 this%fargo%direction = fargo_dir
412 CASE DEFAULT
413 ! this should not happen
414 CALL this%Error("mesh_base::InitMesh","wrong coordinate index returned from GetAzimuthIndex")
415 END SELECT
416 ELSE
417 ! fargo disabled -> set direction to 0
418 this%fargo%direction = 0
419 END IF
420 CASE DEFAULT
421 ! wrong input -> abort
422 CALL this%Error("mesh_base::InitMesh","fargo method should be one of {0,1,2,3}")
423 END SELECT
424
425 ! check if rotational symmetry is assumed, i.e., if the azimuthal
426 ! direction consists of only one cell and a full 2*PI range
427 ! this%ROTSYM is either set to the direction of the azimuthal angle {1,2,3}
428 ! or 0 if rotational symmetry is disabled
429 this%ROTSYM = 0 ! disable rotational symmetry by default
430 SELECT CASE(this%geometry%GetAzimuthIndex())
431 CASE(1)
432 IF (this%INUM.EQ.1.AND.(abs(this%xmax-this%xmin-2*pi).LE.epsilon(this%xmin))) &
433 this%ROTSYM = 1
434 CASE(2)
435 IF (this%JNUM.EQ.1.AND.(abs(this%ymax-this%ymin-2*pi).LE.epsilon(this%ymin))) &
436 this%ROTSYM = 2
437 CASE(3)
438 IF (this%KNUM.EQ.1.AND.(abs(this%zmax-this%zmin-2*pi).LE.epsilon(this%zmin))) &
439 this%ROTSYM = 3
440 END SELECT
441
442 ! These constants are used to access date from adjacent cells, i.e., with
443 ! indices i,j,k +/- 1 and +/- 2. This should always be used instead of direct
444 ! indexing, e.g. i-1 should become i+this%IM1. Some of the constants are
445 ! set to 0 below, if there is only one cell in the particular direction.
446 this%ip1 = 1
447 this%ip2 = 2
448 this%im1 = -1
449 this%im2 = -2
450 this%jp1 = 1
451 this%jp2 = 2
452 this%jm1 = -1
453 this%jm2 = -2
454 this%kp1 = 1
455 this%kp2 = 2
456 this%km1 = -1
457 this%km2 = -2
458
459 ! number of ghost cells; default is 2 at all boundaries
460 ! it may be set to 0, if there is only one cell in
461 ! a particular dimension (see below)
462 this%GNUM = 2
463 this%GINUM = this%GNUM
464 this%GJNUM = this%GNUM
465 this%GKNUM = this%GNUM
466
467 ! check dimensionality and suppress boundaries
468 ! if there is only one cell in that direction
469 this%NDIMS = 3 ! assume 3D simulation
470 this%VECTOR_COMPONENTS = int(b'111') ! enable all vector components
471 IF (this%INUM.EQ.1) THEN
472 IF (this%shear_dir.EQ.1) &
473 CALL this%Error("mesh_base:Initmesh","shearing box enabled with direction 1, but INUM set to 1")
474 IF (this%fargo%direction.EQ.1) &
475 CALL this%Error("mesh_base:Initmesh","fargo transport enabled with direction 1, but INUM set to 1")
476 this%NDIMS = this%NDIMS-1
477 IF (this%ROTSYM.NE.1) this%VECTOR_COMPONENTS = ieor(this%VECTOR_COMPONENTS,vector_x)
478 this%GINUM = 0
479 this%ip1 = 0
480 this%ip2 = 0
481 this%im1 = 0
482 this%im2 = 0
483 END IF
484 IF (this%JNUM.EQ.1) THEN
485 IF (this%shear_dir.EQ.2) &
486 CALL this%Error("mesh_base:Initmesh","shearing box enabled with direction 2, but JNUM set to 1")
487 IF (this%fargo%direction.EQ.2) &
488 CALL this%Error("mesh_base:Initmesh","fargo transport enabled with direction 2, but JNUM set to 1")
489 this%NDIMS = this%NDIMS-1
490 IF (this%ROTSYM.NE.2) this%VECTOR_COMPONENTS = ieor(this%VECTOR_COMPONENTS,vector_y)
491 this%GJNUM = 0
492 this%jp1 = 0
493 this%jp2 = 0
494 this%jm1 = 0
495 this%jm2 = 0
496 END IF
497 IF (this%KNUM.EQ.1) THEN
498 IF (this%shear_dir.EQ.3) &
499 CALL this%Error("mesh_base:Initmesh","shearing box enabled with direction 3, but KNUM set to 1")
500 IF (this%fargo%direction.EQ.3) &
501 CALL this%Error("mesh_base:Initmesh","fargo transport enabled with direction 3, but KNUM set to 1")
502 this%NDIMS = this%NDIMS-1
503 IF (this%ROTSYM.NE.3) this%VECTOR_COMPONENTS = ieor(this%VECTOR_COMPONENTS,vector_z)
504 this%GKNUM = 0
505 this%kp1 = 0
506 this%kp2 = 0
507 this%km1 = 0
508 this%km2 = 0
509 END IF
510 ! set number of faces and corners depending on dimensionality
511 this%NFACES = nfaces(this%NDIMS)
512 this%NCORNERS = ncorners(this%NDIMS)
513
514 ! set index ranges
515 ! ATTENTION: all variables are local on MPI processes, hence this%IMIN on
516 ! one process may differ from this%IMIN on another process
517#ifdef PARALLEL
518 ! compute the MPI partitioning, i.e. the decomposition of the computational domain,
519 ! and setup various data structures MPI uses for communication
520 CALL initmesh_parallel(this, config)
521 CALL mpi_barrier(mpi_comm_world,err)
522 ! determine the maximal amount of cells in 1st dimension
523 inum = this%IMAX - this%IMIN + 1
524 CALL mpi_allreduce(inum,this%MININUM,1,mpi_integer,mpi_min,mpi_comm_world,err)
525 CALL mpi_allreduce(inum,this%MAXINUM,1,mpi_integer,mpi_max,mpi_comm_world,err)
526 ! determine the maximal amount of cells in 2nd dimension
527 jnum = this%JMAX - this%JMIN + 1
528 CALL mpi_allreduce(jnum,this%MINJNUM,1,mpi_integer,mpi_min,mpi_comm_world,err)
529 CALL mpi_allreduce(jnum,this%MAXJNUM,1,mpi_integer,mpi_max,mpi_comm_world,err)
530 ! determine the maximal amount of cells in 3rd dimension
531 knum = this%KMAX - this%KMIN + 1
532 CALL mpi_allreduce(knum,this%MINKNUM,1,mpi_integer,mpi_min,mpi_comm_world,err)
533 CALL mpi_allreduce(knum,this%MAXKNUM,1,mpi_integer,mpi_max,mpi_comm_world,err)
534#else
535 this%IMIN = 1
536 this%IMAX = this%INUM
537 this%JMIN = 1
538 this%JMAX = this%JNUM
539 this%KMIN = 1
540 this%KMAX = this%KNUM
541#endif
542
543 ! index ranges with ghost cells
544 this%IGMIN = this%IMIN - this%GINUM
545 this%IGMAX = this%IMAX + this%GINUM
546 this%JGMIN = this%JMIN - this%GJNUM
547 this%JGMAX = this%JMAX + this%GJNUM
548 this%KGMIN = this%KMIN - this%GKNUM
549 this%KGMAX = this%KMAX + this%GKNUM
550
551 ! initialize mesh arrays using indices with ghost cells
552 CALL initmeshproperties(this%IGMIN,this%IGMAX,this%JGMIN,this%JGMAX,this%KGMIN,this%KGMAX)
553
554 ! create selection for the internal region
555 this%without_ghost_zones = selection_base((/this%IMIN,this%IMAX,this%JMIN,this%JMAX,this%KMIN,this%KMAX/))
556
557 ! coordinate differences in each direction
558 this%dx = (this%xmax - this%xmin) / this%INUM
559 mesh_dx = this%dx
560 this%dy = (this%ymax - this%ymin) / this%JNUM
561 mesh_dy = this%dy
562 this%dz = (this%zmax - this%zmin) / this%KNUM
563 mesh_dz = this%dz
564 ! reset dx,dy,dz if direction is suppressed
565 ! additionally set a mesh_dx,mesh_dy,mesh_dz to set corners, etc.
566 ! for output correctly (see below)
567 IF (this%dx.LT.2*epsilon(this%dx)) THEN
568 IF (this%INUM.EQ.1) THEN
569 this%dx = 1.0
570 mesh_dx = 0.0
571 ELSE
572 CALL this%Error("mesh_base::InitMesh","INUM > 1 conflicts with zero x-extent")
573 END IF
574 END IF
575 IF (this%dy.LT.2*epsilon(this%dy)) THEN
576 IF (this%JNUM.EQ.1) THEN
577 this%dy = 1.0
578 mesh_dy = 0.0
579 ELSE
580 CALL this%Error("mesh_base::InitMesh","JNUM > 1 conflicts with zero y-extent")
581 END IF
582 END IF
583 IF (this%dz.LT.2*epsilon(this%dz)) THEN
584 IF (this%KNUM.EQ.1) THEN
585 this%dz = 1.0
586 mesh_dz = 0.0
587 ELSE
588 CALL this%Error("mesh_base::InitMesh","KNUM > 1 conflicts with zero z-extent")
589 END IF
590 END IF
591
592 ! inverse coordinate differences
593 this%invdx = 1./this%dx
594 this%invdy = 1./this%dy
595 this%invdz = 1./this%dz
596
597 ! allocate memory for curvilinear positions
598 this%curv = marray_cellvector()
599 this%center => this%curv%RemapBounds(this%curv%center)
600 this%bcenter => this%curv%RemapBounds(this%curv%bcenter)
601
602 ! translation vectors for cell faces and cell corners
603 ! (with respect to geometrical cell center)
604 cfaces(1,1) = -0.5*mesh_dx ! western x coordinate
605 cfaces(1,2) = 0.0 ! western y coordinate
606 cfaces(1,3) = 0.0 ! western z coordinate
607 cfaces(2,1) = 0.5*mesh_dx ! eastern x coordinate
608 cfaces(2,2) = 0.0 ! eastern y coordinate
609 cfaces(2,3) = 0.0 ! eastern z coordinate
610 cfaces(3,1) = 0.0 ! southern x coordinate
611 cfaces(3,2) = -0.5*mesh_dy ! southern y coordinate
612 cfaces(3,3) = 0.0 ! southern z coordinate
613 cfaces(4,1) = 0.0 ! northern x coordinate
614 cfaces(4,2) = 0.5*mesh_dy ! northern y coordinate
615 cfaces(4,3) = 0.0 ! northern z coordinate
616 cfaces(5,1) = 0.0 ! bottom x coordinate
617 cfaces(5,2) = 0.0 ! bottom y coordinate
618 cfaces(5,3) = -0.5*mesh_dz ! bottom z coordinate
619 cfaces(6,1) = 0.0 ! top x coordinate
620 cfaces(6,2) = 0.0 ! top y coordinate
621 cfaces(6,3) = 0.5*mesh_dz ! top z coordinate
622
623 ccorners(1,1) = cfaces(1,1) ! bottom-south-west
624 ccorners(1,2) = cfaces(3,2)
625 ccorners(1,3) = cfaces(5,3)
626 ccorners(2,1) = cfaces(2,1) ! bottom-south-east
627 ccorners(2,2) = cfaces(3,2)
628 ccorners(2,3) = cfaces(5,3)
629 ccorners(3,1) = cfaces(1,1) ! bottom-north-west
630 ccorners(3,2) = cfaces(4,2)
631 ccorners(3,3) = cfaces(5,3)
632 ccorners(4,1) = cfaces(2,1) ! bottom-north-east
633 ccorners(4,2) = cfaces(4,2)
634 ccorners(4,3) = cfaces(5,3)
635 ccorners(5,1) = cfaces(1,1) ! top-south-west
636 ccorners(5,2) = cfaces(3,2)
637 ccorners(5,3) = cfaces(6,3)
638 ccorners(6,1) = cfaces(2,1) ! top-south-east
639 ccorners(6,2) = cfaces(3,2)
640 ccorners(6,3) = cfaces(6,3)
641 ccorners(7,1) = cfaces(1,1) ! top-north-west
642 ccorners(7,2) = cfaces(4,2)
643 ccorners(7,3) = cfaces(6,3)
644 ccorners(8,1) = cfaces(2,1) ! top-north-east
645 ccorners(8,2) = cfaces(4,2)
646 ccorners(8,3) = cfaces(6,3)
647 ! calculate coordinates (geometric centers)
648 DO concurrent(i=this%IGMIN:this%IGMAX,j=this%JGMIN:this%JGMAX,k=this%KGMIN:this%KGMAX)
649 ! geometrical cell centers
650 this%curv%center(i,j,k,1) = this%xmin + (2*i-1)*0.5*mesh_dx ! x coord
651 this%curv%center(i,j,k,2) = this%ymin + (2*j-1)*0.5*mesh_dy ! y coord
652 this%curv%center(i,j,k,3) = this%zmin + (2*k-1)*0.5*mesh_dz ! z coord
653 ! set bary center curvilinear coordinates to geometric centers
654 ! will be overwritten later
655 this%curv%bcenter(i,j,k,:) = this%curv%center(i,j,k,:)
656 ! cell face centered positions
657 this%curv%faces(i,j,k,1,:) = this%curv%center(i,j,k,:) + cfaces(1,:) ! western
658 this%curv%faces(i,j,k,2,:) = this%curv%center(i,j,k,:) + cfaces(2,:) ! eastern
659 this%curv%faces(i,j,k,3,:) = this%curv%center(i,j,k,:) + cfaces(3,:) ! southern
660 this%curv%faces(i,j,k,4,:) = this%curv%center(i,j,k,:) + cfaces(4,:) ! northern
661 this%curv%faces(i,j,k,5,:) = this%curv%center(i,j,k,:) + cfaces(5,:) ! bottom
662 this%curv%faces(i,j,k,6,:) = this%curv%center(i,j,k,:) + cfaces(6,:) ! top
663 ! cell corner positions
664 this%curv%corners(i,j,k,1,:) = this%curv%center(i,j,k,:) + ccorners(1,:) ! bottom-south-west
665 this%curv%corners(i,j,k,2,:) = this%curv%center(i,j,k,:) + ccorners(2,:) ! bottom-south-east
666 this%curv%corners(i,j,k,3,:) = this%curv%center(i,j,k,:) + ccorners(3,:) ! bottom-north-west
667 this%curv%corners(i,j,k,4,:) = this%curv%center(i,j,k,:) + ccorners(4,:) ! bottom-north-east
668 this%curv%corners(i,j,k,5,:) = this%curv%center(i,j,k,:) + ccorners(5,:) ! top-south-west
669 this%curv%corners(i,j,k,6,:) = this%curv%center(i,j,k,:) + ccorners(6,:) ! top-south-east
670 this%curv%corners(i,j,k,7,:) = this%curv%center(i,j,k,:) + ccorners(7,:) ! top-north-west
671 this%curv%corners(i,j,k,8,:) = this%curv%center(i,j,k,:) + ccorners(8,:) ! top-north-east
672 END DO
673
674 ! basic mesh initialization
675
676 ! create mesh arrays for scale factors
677 this%hx = marray_cellscalar()
678 this%hy = marray_cellscalar()
679 this%hz = marray_cellscalar()
680
681 ! get geometrical scale factors for all cell positions
682 ! bary center values are overwritten below
683 CALL this%Geometry%ScaleFactors(this%curv,this%hx,this%hy,this%hz)
684
685 ! get square root of determinant of the metric
686 ! bary center values are overwritten below
687 ! ATTENTION: it seems that the intel compiler does not
688 ! assign the result of marray operations correctly
689 ! leading to segfaults in some cases
690 this%sqrtg = marray_cellscalar()
691#if defined (__INTEL_COMPILER) || defined (__flang__)
692 this%sqrtg%data1d(:) = this%hx%data1d(:)*(this%hy%data1d(:)*this%hz%data1d(:))
693#else
694 this%sqrtg = this%hx*(this%hy*this%hz)
695#endif
696
697 ! create mesh arrays for commutator coefficients
698 this%cxyx = marray_cellscalar()
699 this%cxzx = marray_cellscalar()
700 this%cyxy = marray_cellscalar()
701 this%cyzy = marray_cellscalar()
702 this%czxz = marray_cellscalar()
703 this%czyz = marray_cellscalar()
704
705 ! create mesh arrays for line and volume elements and related arrays
706 this%volume = marray_base()
707 this%dxdydV = marray_base()
708 this%dydzdV = marray_base()
709 this%dzdxdV = marray_base()
710 this%dlx = marray_base()
711 this%dly = marray_base()
712 this%dlz = marray_base()
713
714 ! create mesh array for cartesian coordinates
715 this%cart = marray_cellvector()
716 this%bccart => this%cart%RemapBounds(this%cart%bcenter)
717 this%fccart => this%cart%RemapBounds(this%cart%faces)
718 this%ccart => this%cart%RemapBounds(this%cart%corners)
719
720 ! compute cartesian coordinates for all cell positions (center,bcenter,faces,corners)
721 CALL this%Geometry%Convert2Cartesian(this%curv,this%cart)
722
723 ! create mesh array for distance to the origin of the mesh
724 this%radius = marray_cellscalar()
725 CALL this%geometry%Radius(this%curv,this%radius)
726
727 ! create mesh array for position vector
728 this%posvec = marray_cellvector()
729 CALL this%geometry%PositionVector(this%curv,this%posvec)
730
731 ! nullify remaining (mesh) arrays
732 NULLIFY(this%dAx,this%dAy,this%dAz, &
733 this%dAxdydz,this%dAydzdx,this%dAzdxdy, &
734 this%rotation,this%weights)
735
736! ! This is a quick hack for VTK and XDMF output on bipolar mesh.
737! ! It maps infinite coordinates on the boundary in cart%corners to finite values.
738! ! The bipolar mesh has a degeneracy at sigma=0(or 2*PI) and tau=0. If we
739! ! we transform to cartesian coordinates, these points would be mapped onto
740! ! a circle at infinity.
741! !> \todo NOT VERIFIED here is most probably something wrong in the 3D version
742! !! in most recent fosite 2d version this part is declared obsolete
743! IF (this%Geometry%GetType().EQ.BIPOLAR) THEN
744! ! get the maximal cartesian y coordinate at multiply by some appropriate number > 1
745! cart_max = 1.5*MAXVAL(ABS(this%cart%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX,2)))
746! DO j=this%JMIN,this%JMAX+1
747! IF (this%IMIN.EQ.1) THEN ! this is important for MPI
748! i = 1 ! boundary at tau_min
749! k = 1
750! ! check if tau changes sign
751! IF ((this%curv%corners(i,j,k,1,2).GT.0.0).AND.(this%curv%corners(i,j-1,k,1,2).LE.0.0)) THEN
752! ! map the cartesian corner to a value (x,y) = (0,cart_max)
753! IF (ABS(this%cart%corners(i,j-1,k,1,1)).LT.ABS(this%cart%corners(i,j,k,1,1))) THEN
754! this%cart%corners(i,j,k,1,1) = 0.0
755! this%cart%corners(i,j,k,1,2) = cart_max
756! this%cart%corners(i,j-1,k,3,:) = this%cart%corners(i,j,k,1,:)
757! ELSE
758! this%cart%corners(i,j-1,k,1,1) = 0.0
759! this%cart%corners(i,j-1,k,1,2) = cart_max
760! this%cart%corners(i,j-2,k,3,:) = this%cart%corners(i,j-1,k,1,:)
761! END IF
762! END IF
763! END IF
764! IF (this%IMAX.EQ.this%INUM) THEN ! this is important for MPI
765! i = this%IMAX+1 ! boundary at tau_max
766! ! check if tau changes sign
767! IF ((this%curv%corners(i,j,k,1,2).GT.0.0).AND.(this%curv%corners(i,j-1,k,1,2).LE.0.0)) THEN
768! ! map the cartesian corner to a value (x,y) = (0,cart_max)
769! IF (ABS(this%cart%corners(i,j-1,k,1,1)).GT.ABS(this%cart%corners(i,j,k,1,1))) THEN
770! this%cart%corners(i,j-1,k,1,1) = 0.0
771! this%cart%corners(i,j-1,k,1,2) = -cart_max
772! this%cart%corners(i-1,j-1,k,2,:) = this%cart%corners(i,j-1,k,1,:)
773! ELSE
774! this%cart%corners(i,j,k,1,1) = 0.0
775! this%cart%corners(i,j,k,1,2) = -cart_max
776! this%cart%corners(i-1,j,k,2,:) = this%cart%corners(i,j,k,1,:)
777! END IF
778! EXIT
779! END IF
780! END IF
781! END DO
782! END IF
783
784 ! print some information
785 CALL this%Info(" MESH-----> quadrature rule: " // trim(this%GetName()))
786 WRITE (xres, '(I0)') this%INUM ! this is just for better looking output
787 WRITE (yres, '(I0)') this%JNUM
788 WRITE (zres, '(I0)') this%KNUM
789 CALL this%Info(" resolution: " // trim(xres) // " x " // trim(yres) // " y "// trim(zres) // " z ")
790 WRITE (xres, '(ES9.2,A,ES9.2)') this%xmin, " ..", this%xmax
791 WRITE (yres, '(ES9.2,A,ES9.2)') this%ymin, " ..", this%ymax
792 WRITE (zres, '(ES9.2,A,ES9.2)') this%zmin, " ..", this%zmax
793 CALL this%Info(" computat. domain: x=" // trim(xres))
794 CALL this%Info(" y=" // trim(yres))
795 CALL this%Info(" z=" // trim(zres))
796 IF(this%omega.NE.0.) THEN
797 WRITE (somega, '(ES9.2)') this%omega
798 CALL this%Info(" ang. velocity: " // trim(somega))
799 END IF
800 SELECT CASE(this%ROTSYM)
801 CASE(0)
802 somega = "disabled"
803 CASE(1)
804 somega = "x-coordinate"
805 CASE(2)
806 somega = "y-coordinate"
807 CASE(3)
808 somega = "z-coordinate"
809 END SELECT
810 CALL this%Info(" rot. symmetry: " // trim(somega))
811#ifdef PARALLEL
812 WRITE (xres, '(I0)') this%dims(1)
813 WRITE (yres, '(I0)') this%dims(2)
814 WRITE (zres, '(I0)') this%dims(3)
815 CALL this%Info(" MPI partition: " // trim(xres) // " x " // trim(yres) // " y "// trim(zres) // " z ")
816#endif
817
818 ! \todo implement a correct check of the
819 ! curvilinear range (e.g. xi > 0 in polar)
820
821
822 CALL setoutput(this,config,io)
823
824 END SUBROUTINE initmesh
825
827 SUBROUTINE setoutput(this,config,IO)
828 IMPLICIT NONE
829 !------------------------------------------------------------------------!
830 CLASS(mesh_base),INTENT(INOUT) :: this
831 TYPE(dict_typ),POINTER :: config,IO
832 !------------------------------------------------------------------------!
833 INTEGER :: writefields
834 !------------------------------------------------------------------------!
835 !Set OutputDict
836 CALL getattr(config, "output/corners", writefields, 0)
837 IF((writefields.EQ.1).AND.ASSOCIATED(this%ccart)) &
838 CALL setattr(io,"corners",this%ccart(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX,:,:))
839
840 CALL getattr(config, "output/grid", writefields, 1)
841 IF((writefields.EQ.1).AND.ASSOCIATED(this%ccart)) THEN
842 CALL setattr(io,"grid_x",this%cart%corners(this%IMIN:this%IMAX+this%ip1, &
843 this%JMIN:this%JMAX+this%jp1, &
844 this%KMIN:this%KMAX+this%kp1,1,1))
845 CALL setattr(io,"grid_y",this%cart%corners(this%IMIN:this%IMAX+this%ip1, &
846 this%JMIN:this%JMAX+this%jp1, &
847 this%KMIN:this%KMAX+this%kp1,1,2))
848 CALL setattr(io,"grid_z",this%cart%corners(this%IMIN:this%IMAX+this%ip1, &
849 this%JMIN:this%JMAX+this%jp1, &
850 this%KMIN:this%KMAX+this%kp1,1,3))
851 END IF
852
853 CALL getattr(config, "output/bary_curv", writefields, 1)
854 IF((writefields.EQ.1).AND.ASSOCIATED(this%bcenter)) &
855 CALL setattr(io,"bary_curv",this%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX,:))
856! Dict("name" / "coordinates:bary_curv"))
857 CALL getattr(config, "output/bary", writefields, 1)
858 IF((writefields.EQ.1).AND.ASSOCIATED(this%bccart)) &
859 CALL setattr(io,"bary_centers",this%bccart(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX,:))
860! Dict("name" / "coordinates:bary"))
861
862 CALL getattr(config, "output/rotation", writefields, 0)
863 IF(writefields.EQ.1) THEN
864 CALL calculaterotation(this)
865 CALL setattr(io,"rotation",this%rotation(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
866 END IF
867
868 CALL getattr(config, "output/dl", writefields, 0)
869 IF(writefields.EQ.1) THEN
870 IF (ASSOCIATED(this%dlx%data3d)) &
871 CALL setattr(io,"dlx",this%dlx%data3d(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
872 IF (ASSOCIATED(this%dly%data3d)) &
873 CALL setattr(io,"dly",this%dly%data3d(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
874 IF (ASSOCIATED(this%dlz%data3d)) &
875 CALL setattr(io,"dlz",this%dlz%data3d(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
876 END IF
877
878 CALL getattr(config, "output/bh", writefields, 0)
879 IF(writefields.EQ.1) THEN
880 IF (ASSOCIATED(this%hx%bcenter)) &
881 CALL setattr(io,"bhx",this%hx%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
882 IF (ASSOCIATED(this%hy%bcenter)) &
883 CALL setattr(io,"bhy",this%hy%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
884 IF (ASSOCIATED(this%hz%bcenter)) &
885 CALL setattr(io,"bhz",this%hz%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
886 END IF
887
888 CALL getattr(config, "output/commutator", writefields, 0)
889 IF(writefields.EQ.1) THEN
890 IF (ASSOCIATED(this%cxyx%bcenter)) &
891 CALL setattr(io,"cxyx",this%cxyx%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
892 IF (ASSOCIATED(this%cxzx%bcenter)) &
893 CALL setattr(io,"cxzx",this%cxzx%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
894 IF (ASSOCIATED(this%cyxy%bcenter)) &
895 CALL setattr(io,"cyxy",this%cyxy%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
896 IF (ASSOCIATED(this%cyzy%bcenter)) &
897 CALL setattr(io,"cyzy",this%cyzy%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
898 IF (ASSOCIATED(this%czxz%bcenter)) &
899 CALL setattr(io,"czxz",this%czxz%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
900 IF (ASSOCIATED(this%czyz%bcenter)) &
901 CALL setattr(io,"czyz",this%czyz%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
902 END IF
903
904 CALL getattr(config, "output/volume", writefields, 0)
905 IF((writefields.EQ.1).AND.ASSOCIATED(this%volume%data3d)) &
906 CALL setattr(io,"volume",this%volume%data3d(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
907
908 CALL getattr(config, "output/dA", writefields, 0)
909 IF(writefields.EQ.1) THEN
910 IF (ASSOCIATED(this%dAx)) &
911 CALL setattr(io,"dAx",this%dAx(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX,:))
912 IF (ASSOCIATED(this%dAy)) &
913 CALL setattr(io,"dAy",this%dAy(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX,:))
914 IF (ASSOCIATED(this%dAz)) &
915 CALL setattr(io,"dAz",this%dAz(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX,:))
916 END IF
917
918 CALL getattr(config, "output/radius", writefields, 0)
919 IF((writefields.EQ.1).AND.ASSOCIATED(this%radius%bcenter)) &
920 CALL setattr(io,"radius",this%radius%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
921
922 CALL getattr(config, "output/position_vector", writefields, 0)
923 IF((writefields.EQ.1).AND.ASSOCIATED(this%posvec%bcenter)) &
924 CALL setattr(io,"position_vector",this%posvec%bcenter(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX,:))
925
926 END SUBROUTINE setoutput
927
928
935 SUBROUTINE calculaterotation(this)
936 IMPLICIT NONE
937 !------------------------------------------------------------------------!
938 CLASS(mesh_base) :: this
939 !------------------------------------------------------------------------!
940 INTEGER :: status
941 REAL,DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,3) :: cart,curv
942 !------------------------------------------------------------------------!
943 ALLOCATE(this%rotation(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX), &
944 stat=status)
945 IF(status.NE.0) &
946 CALL this%Error("CalculateRotation", "Couldn't allocate memory")
947 ! initialize with vectors pointing in x-direction
948 cart(:,:,:,1) = 1.
949 cart(:,:,:,2) = 0.
950 cart(:,:,:,3) = 0.
951 CALL this%geometry%Convert2Curvilinear(this%bcenter, cart, curv)
952 this%rotation(:,:,:) = atan2(curv(:,:,:,2),curv(:,:,:,1)) !!/todo{3D}
953 END SUBROUTINE calculaterotation
954
961 PURE FUNCTION internalpoint(this,x,y,z,mask) RESULT(ip)
962 IMPLICIT NONE
963 !------------------------------------------------------------------------!
964 CLASS(mesh_base),INTENT(IN) :: this
965 INTEGER, DIMENSION(6), OPTIONAL :: mask
966 REAL :: x,y,z
967 LOGICAL :: ip
968 !------------------------------------------------------------------------!
969 INTEGER :: imin,imax,jmin,jmax,kmin,kmax
970 !------------------------------------------------------------------------!
971 INTENT(IN) :: x,y,z,mask
972 !------------------------------------------------------------------------!
973 ip = .false.
974 ! compare with the curvilinear extend (curvilinear domain is always rectangular)
975 !
976 ! the function behaves in 2 different ways: global and local
977 !
978 ! 1. if mask is given, we do a local comparison, i.e.
979 ! we check if the point lies inside the specified rectangular
980 ! domain given by its cell indices: imin=mask(1),imax=mask(2),...
981 ! 2. otherwise we do a global check, i.e. with the whole computational
982 ! domain
983 IF (PRESENT(mask)) THEN
984 ! restrict indices to the actual domain,
985 ! which is different for each MPI process in parallel mode
986 imin = max(this%IGMIN,mask(1))
987 imax = min(this%IGMAX,mask(2))
988 jmin = max(this%JGMIN,mask(3))
989 jmax = min(this%JGMAX,mask(4))
990 kmin = max(this%KGMIN,mask(5))
991 kmax = min(this%KGMAX,mask(6))
992 ! first check if the masked region is at least a subdomain of the actual domain
993 IF (((imin.LE.this%IGMAX).AND.(imax.GE.imin)).AND. &
994 ((jmin.LE.this%JGMAX).AND.(jmax.GE.jmin)).AND. &
995 ((kmin.LE.this%KGMAX).AND.(kmax.GE.kmin))) THEN
996 ! compare the curvilinear coordinates at the boundaries of the masked region
997 ! with the transformed coordinates of the given point
998 IF ((x.GE.this%curv%faces(imin,jmin,kmin,1,1).AND.x.LE.this%curv%faces(imax,jmax,kmax,2,1)).AND. &
999 (y.GE.this%curv%faces(imin,jmin,kmin,1,2).AND.y.LE.this%curv%faces(imax,jmax,kmax,2,2)).AND. &
1000 (z.GE.this%curv%faces(imin,jmin,kmin,1,3).AND.z.LE.this%curv%faces(imax,jmax,kmax,2,3))) THEN
1001 ip = .true.
1002 END IF
1003 END IF
1004 ELSE
1005 ! do the global check
1006 IF (((x.GE.this%xmin).AND.(x.LE.this%xmax)).AND. &
1007 ((y.GE.this%ymin).AND.(y.LE.this%ymax)).AND. &
1008 ((z.GE.this%zmin).AND.(z.LE.this%zmax))) THEN
1009 ip = .true.
1010 END IF
1011 END IF
1012 END FUNCTION internalpoint
1013
1014
1015#ifdef PARALLEL
1016
1022 SUBROUTINE initmesh_parallel(this, config)
1023 IMPLICIT NONE
1024 !------------------------------------------------------------------------!
1025 CLASS(mesh_base), INTENT(INOUT) :: this
1026 TYPE(dict_typ),POINTER :: config
1027 !------------------------------------------------------------------------!
1028 CHARACTER(LEN=64) :: decomp_str
1029 LOGICAL, DIMENSION(3) :: remain_dims, periods
1030 INTEGER :: num,rem
1031 INTEGER :: ierror
1032 INTEGER :: i,j,k
1033 INTEGER :: worldgroup,newgroup
1034 INTEGER, DIMENSION(1) :: rank0in, rank0out
1035 INTEGER, DIMENSION(3) :: coords,ncells,dims
1036 INTEGER, ALLOCATABLE, DIMENSION(:) :: ranks
1037 !------------------------------------------------------------------------!
1038 ! 1. Check if the user has requests a particular decomposition.
1039 ! There are two ways a user may control automatic decomposition:
1040 ! (a) explicitly set the number of processes along a certain direction,
1041 ! e.g., decomposition = (/ 5,4,1 /) generates a cartesian communicator
1042 ! with 5 processes along the x-direction and 4 processes along the y-direction
1043 ! (b) if a number specified for certain dimension is negative let the algorithm
1044 ! find an optimal number of processes along that direction, e.g.,
1045 ! decomposition = (/ -1, 4, 1 /) generates a cartesian communicator
1046 ! with NumProcs / 4 processes along the x-direction and 4 processes
1047 ! along the y-direction with NumProcs beeing the total number of MPI
1048 ! processes given at the command line. If 4 is not a prime factor of
1049 ! NumProcs the program aborts with an error message.
1050 ! The default is automatic domain decomposition, i.e. finding the optimal
1051 ! distribution of processes along all directions considering the number of
1052 ! cells and the hardware vector length.
1053 ! This is done in CalculateDecomposition() (see below).
1054
1055 ! defaults to -1, i.e., fully automatic decomposition
1056 dims(:)= -1
1057 CALL getattr(config, "decomposition", dims, dims(1:3))
1058
1059 ! perform the decomposition on rank 0 and broadcast the result
1060 ! to other processes (see below)
1061 IF (this%GetRank().EQ.0) THEN
1062 ! for more elaborate error output
1063 WRITE (decomp_str,'(3(A,I0),A)') "[ ", dims(1), ", ", dims(2), ", ",dims(3), " ]"
1064
1065 ! perform some sanity checks
1066 IF (all(dims(:).GE.1).AND.product(dims(:)).NE.this%GetNumProcs()) &
1067 CALL this%Error("InitMesh_parallel","total number of processes in domain decomposition " &
1068 // trim(decomp_str) // achar(10) // repeat(' ',7) // &
1069 "does not match the number passed to mpirun")
1070
1071 IF (any(dims(:).EQ.0)) &
1072 CALL this%Error("InitMesh_parallel","numbers in decomposition should not be 0")
1073
1074 IF (mod(this%GetNumProcs(),product(dims(:),dims(:).GT.1)).NE.0) &
1075 CALL this%Error("InitMesh_parallel","numbers in decomposition " &
1076 // trim(decomp_str) // achar(10) // repeat(' ',7) // &
1077 "are not devisors of the total number of processes passed to mpirun")
1078
1079 ! set number of cells along each direction again
1080 ncells(1) = this%INUM
1081 ncells(2) = this%JNUM
1082 ncells(3) = this%KNUM
1083
1084 ! no decomposition for suppressed dimensions
1085 WHERE (ncells(:).EQ.1)
1086 dims(:) = 1
1087 END WHERE
1088
1089 ! balance number of processes if requested
1090 IF (all(dims(:).GT.0)) THEN
1091 ! (a) all dims user supplied -> do nothing
1092 ELSE IF (product(dims(:)).GT.0) THEN
1093 ! (b) two dims < 0 one dim > 0 -> find optimal decomposition fixing one of the dims
1094 ! using the user supplied number
1095 k = maxloc(dims(:),1) ! get the one index k with dims(k) > 0
1096 ! suppress decomposition along the k-direction
1097 ncells(k) = 1
1098 i = dims(k) ! remember dims(k)
1099 IF (k.EQ.1) THEN
1100 ! switch entries 1 and 3 to make sure the first entry is not the one we are not decomposing
1101 dims(3) = this%GetNumProcs() / dims(k) ! reduced total number of processes
1102 dims(2) = 1
1103 dims(1) = 1
1104 CALL calculatedecomposition(ncells(3),ncells(2),ncells(1),this%GKNUM, &
1105 dims(3),dims(2),dims(1))
1106 ELSE
1107 dims(1) = this%GetNumProcs() / dims(k) ! reduced total number of processes
1108 dims(2) = 1
1109 dims(3) = 1
1110 CALL calculatedecomposition(ncells(1),ncells(2),ncells(3),this%GKNUM, &
1111 dims(1),dims(2),dims(3))
1112 END IF
1113 dims(k) = i
1114 ELSE IF (all(dims(:).LT.0)) THEN
1115 ! (c) all dims < 0 -> find optimal decomposition using all dimensions
1116 dims(1) = this%GetNumProcs()
1117 dims(2) = 1
1118 dims(3) = 1
1119 CALL calculatedecomposition(ncells(1),ncells(2),ncells(3),this%GINUM, &
1120 dims(1),dims(2),dims(3))
1121 ELSE
1122 ! (d) one dim < 0 two dims > 0 -> fix the two dimensions and set the 3rd
1123 ! using the given number of processes
1124 k = minloc(dims(:),1) ! get the one index k with dims(k) < 0
1125 dims(k) = this%GetNumProcs() / product(dims(:),dims(:).GT.1)
1126 END IF
1127
1128 WRITE (decomp_str,'(3(A,I0),A)') "[ ", dims(1), ", ", dims(2), ", ",dims(3), " ]"
1129
1130 ! set number of cells along each direction again
1131 ncells(1) = this%INUM
1132 ncells(2) = this%JNUM
1133 ncells(3) = this%KNUM
1134 DO i=1,3
1135 IF (dims(i).GT.ncells(i)) &
1136 CALL this%Error("InitMesh_parallel","number of processes exceeds number of cells " &
1137 // achar(10) // repeat(' ',7) // "in dimension " // achar(48+i) // ", check decomposition " &
1138 // trim(decomp_str))
1139 END DO
1140
1141 IF (dims(3).LE.0) THEN
1142 CALL this%Error("InitMesh_parallel","automatic domain decomposition failed.")
1143 END IF
1144
1145 this%dims(:) = dims(:)
1146
1147 END IF
1148
1149 CALL mpi_bcast(this%dims,3,mpi_integer,0,mpi_comm_world,ierror)
1150
1151! PRINT *,this%dims(:)
1152! CALL this%Error("InitMesh_parallel","debug breakpoint")
1153
1154 ! 2. create the cartesian communicator
1155 ! IMPORTANT: disable reordering of nodes
1156 periods = .false.
1157 CALL mpi_cart_create(mpi_comm_world,3,this%dims,periods,.true., &
1158 this%comm_cart,ierror)
1159
1160 ! 3. inquire and save the own position
1161 this%mycoords(:) = 0
1162 CALL mpi_cart_get(this%comm_cart,3,this%dims,periods,this%mycoords,ierror)
1163
1164 ! 4. create communicators for every column and row of the cartesian topology
1165 remain_dims = (/ .true., .false., .false. /)
1166 CALL mpi_cart_sub(this%comm_cart,remain_dims,this%Icomm,ierror)
1167 remain_dims = (/ .false., .true., .false. /)
1168 CALL mpi_cart_sub(this%comm_cart,remain_dims,this%Jcomm,ierror)
1169 remain_dims = (/ .false., .false., .true. /)
1170 CALL mpi_cart_sub(this%comm_cart,remain_dims,this%Kcomm,ierror)
1171
1172 ! subdivide the mesh and set mesh indices
1173 ! x-direction
1174 rem = mod(this%INUM,this%dims(1)) ! remainder
1175 num = (this%INUM-rem) / this%dims(1) ! fraction
1176 IF (this%mycoords(1).LT.rem) THEN
1177 ! the first (rem-1) nodes get one more to account for the remainder
1178 this%IMIN = this%mycoords(1) * (num + 1) + 1
1179 this%IMAX = this%IMIN + num
1180 ELSE
1181 this%IMIN = rem * (num+1) + (this%mycoords(1) - rem) * num + 1
1182 this%IMAX = this%IMIN + num - 1
1183 END IF
1184 ! y-direction
1185 rem = mod(this%JNUM,this%dims(2)) ! remainder
1186 num = (this%JNUM-rem) / this%dims(2) ! fraction
1187 IF (this%mycoords(2).LT.rem) THEN
1188 ! the first (rem-1) nodes get one more to account for the remainder
1189 this%JMIN = this%mycoords(2) * (num + 1) + 1
1190 this%JMAX = this%JMIN + num
1191 ELSE
1192 this%JMIN = rem * (num+1) + (this%mycoords(2) - rem) * num + 1
1193 this%JMAX = this%JMIN + num - 1
1194 END IF
1195 ! z-direction
1196 rem = mod(this%KNUM,this%dims(3)) ! remainder
1197 num = (this%KNUM-rem) / this%dims(3) ! fraction
1198 IF (this%mycoords(3).LT.rem) THEN
1199 ! the first (rem-1) nodes get one more to account for the remainder
1200 this%KMIN = this%mycoords(3) * (num + 1) + 1
1201 this%KMAX = this%KMIN + num
1202 ELSE
1203 this%KMIN = rem * (num+1) + (this%mycoords(3) - rem) * num + 1
1204 this%KMAX = this%KMIN + num - 1
1205 END IF
1206
1207 ! create communicators for all boundaries
1208 ALLOCATE(ranks(max(this%dims(1)*this%dims(2),this%dims(1)*this%dims(3),this%dims(2)*this%dims(3))))
1209 ranks = 0
1210 CALL mpi_comm_group(mpi_comm_world,worldgroup,ierror)
1211 ! western boundary
1212 coords(1) = 0
1213 DO k=0,this%dims(3)-1
1214 DO j=0,this%dims(2)-1
1215 coords(2) = j
1216 coords(3) = k
1217 CALL mpi_cart_rank(this%comm_cart,coords,i,ierror)
1218 ranks(j+1+k*this%dims(2))=i
1219 END DO
1220 END DO
1221 CALL mpi_group_incl(worldgroup,this%dims(2)*this%dims(3),ranks,newgroup,ierror)
1222 CALL mpi_comm_create(this%comm_cart,newgroup,this%comm_boundaries(1),ierror)
1223 rank0in(1) = 0
1224 CALL mpi_group_translate_ranks(newgroup,1,rank0in,worldgroup,rank0out,ierror)
1225 this%rank0_boundaries(1) = rank0out(1)
1226 CALL mpi_group_free(newgroup,ierror)
1227 ! eastern boundary
1228 coords(1) = this%dims(1)-1
1229 DO k=0,this%dims(3)-1
1230 DO j=0,this%dims(2)-1
1231 coords(2) = j
1232 coords(3) = k
1233 CALL mpi_cart_rank(this%comm_cart,coords,i,ierror)
1234 ranks(j+1+k*this%dims(2))=i
1235 END DO
1236 END DO
1237 CALL mpi_group_incl(worldgroup,this%dims(2)*this%dims(3),ranks,newgroup,ierror)
1238 CALL mpi_comm_create(this%comm_cart,newgroup,this%comm_boundaries(2),ierror)
1239 rank0in(1) = 0
1240 CALL mpi_group_translate_ranks(newgroup,1,rank0in,worldgroup,rank0out,ierror)
1241 this%rank0_boundaries(2) = rank0out(1)
1242 CALL mpi_group_free(newgroup,ierror)
1243 ! southern boundary
1244 coords(2) = 0
1245 DO k=0,this%dims(3)-1
1246 DO i=0,this%dims(1)-1
1247 coords(1) = i
1248 coords(3) = k
1249 CALL mpi_cart_rank(this%comm_cart,coords,j,ierror)
1250 ranks(i+1+k*this%dims(1))=j
1251 END DO
1252 END DO
1253 CALL mpi_group_incl(worldgroup,this%dims(1)*this%dims(3),ranks,newgroup,ierror)
1254 CALL mpi_comm_create(this%comm_cart,newgroup,this%comm_boundaries(3),ierror)
1255 rank0in(1) = 0
1256 CALL mpi_group_translate_ranks(newgroup,1,rank0in,worldgroup,rank0out,ierror)
1257 this%rank0_boundaries(3) = rank0out(1)
1258 CALL mpi_group_free(newgroup,ierror)
1259 ! northern boundary
1260 coords(2) = this%dims(2)-1
1261 DO k=0,this%dims(3)-1
1262 DO i=0,this%dims(1)-1
1263 coords(1) = i
1264 coords(3) = k
1265 CALL mpi_cart_rank(this%comm_cart,coords,j,ierror)
1266 ranks(i+1+k*this%dims(1))=j
1267 END DO
1268 END DO
1269 CALL mpi_group_incl(worldgroup,this%dims(1)*this%dims(3),ranks,newgroup,ierror)
1270 CALL mpi_comm_create(this%comm_cart,newgroup,this%comm_boundaries(4),ierror)
1271 rank0in(1) = 0
1272 CALL mpi_group_translate_ranks(newgroup,1,rank0in,worldgroup,rank0out,ierror)
1273 this%rank0_boundaries(4) = rank0out(1)
1274 CALL mpi_group_free(newgroup,ierror)
1275 ! bottomer boundaries
1276 coords(3) = 0
1277 DO j=0,this%dims(2)-1
1278 DO i=0,this%dims(1)-1
1279 coords(1) = i
1280 coords(2) = j
1281 CALL mpi_cart_rank(this%comm_cart,coords,k,ierror)
1282 ranks(i+1+j*this%dims(1))=k
1283 END DO
1284 END DO
1285 CALL mpi_group_incl(worldgroup,this%dims(1)*this%dims(2),ranks,newgroup,ierror)
1286 CALL mpi_comm_create(this%comm_cart,newgroup,this%comm_boundaries(5),ierror)
1287 rank0in(1) = 0
1288 CALL mpi_group_translate_ranks(newgroup,1,rank0in,worldgroup,rank0out,ierror)
1289 this%rank0_boundaries(5) = rank0out(1)
1290 CALL mpi_group_free(newgroup,ierror)
1291 ! topper boundaries
1292 coords(3) = this%dims(3)-1
1293 DO j=0,this%dims(2)-1
1294 DO i=0,this%dims(1)-1
1295 coords(1) = i
1296 coords(2) = j
1297 CALL mpi_cart_rank(this%comm_cart,coords,k,ierror)
1298 ranks(i+1+j*this%dims(1))=k
1299 END DO
1300 END DO
1301 CALL mpi_group_incl(worldgroup,this%dims(1)*this%dims(2),ranks,newgroup,ierror)
1302 CALL mpi_comm_create(this%comm_cart,newgroup,this%comm_boundaries(6),ierror)
1303 rank0in(1) = 0
1304 CALL mpi_group_translate_ranks(newgroup,1,rank0in,worldgroup,rank0out,ierror)
1305 this%rank0_boundaries(6) = rank0out(1)
1306 CALL mpi_group_free(newgroup,ierror)
1307
1308 CALL mpi_group_free(worldgroup,ierror)
1309 DEALLOCATE(ranks)
1310 END SUBROUTINE initmesh_parallel
1311
1312
1320 SUBROUTINE calculatedecomposition(ni,nj,nk,ginum,pi,pj,pk)
1321 USE factors
1322 IMPLICIT NONE
1323 !------------------------------------------------------------------------!
1324 INTEGER, INTENT(IN) :: ni,nj,nk,ginum
1325 INTEGER, INTENT(INOUT) :: pi,pj,pk
1326 !------------------------------------------------------------------------!
1327 CHARACTER(LEN=8) :: svl_char
1328 INTEGER :: svl,err
1329 !------------------------------------------------------------------------!
1330 ! get the system vector length from preprocessor macro (string)
1331 svl_char = vector_length
1332 READ (svl_char,'(I8)',iostat=err) svl
1333 ! return immediatly for malformed input
1334 IF (any((/pi.LT.1,pi.GT.maxnum,pj.NE.1,pk.NE.1,err.NE.0, &
1335 min(ni,nj,nk).LT.1/))) THEN
1336 pk = 0
1337 RETURN
1338 END IF
1339 CALL decompose(pi,pj,pk)
1340
1341 CONTAINS
1342
1355 RECURSIVE SUBROUTINE decompose(pi,pj,pk)
1356 IMPLICIT NONE
1357 !-------------------------------------------------------------------!
1358 INTEGER, INTENT(INOUT) :: pi,pj,pk
1359 !-------------------------------------------------------------------!
1360 INTEGER :: pp,ptot
1361 INTEGER :: p1,p2,p3,pinew,pjnew,pknew,piold,pjold,pkold,pinew_,pjnew_,pknew_
1362 INTEGER :: pfmin,pfnew,pfold
1363 INTEGER :: bl,vl,blnew,vlnew,blnew_,vlnew_
1364 REAL :: bl_gain,vl_gain
1365 !-------------------------------------------------------------------!
1366 ! return immediatly if the number of processes exceeds the number of
1367 ! grid cells, i.e. no subdivision possible
1368 IF (pj.GT.nj.OR.pk.GT.nk) RETURN
1369 ! measure the costs of the given configuration
1370 CALL getcosts(ni,nj,nk,pi,pj,pk,bl,vl)
1371! PRINT '(3(A,I4),A,I7,A,I4)'," pi=",pi," pj=",pj," pk=",pk," boundary length=",bl," vector length=",vl
1372 ! save the configuration
1373 piold=pi
1374 pjold=pj
1375 pkold=pk
1376 p1 = pi
1377 p2 = pj
1378 p3 = pk
1379 ! compute the total number of processes
1380 ptot = pi*pj*pk
1381 pfmin = getfactor(pj) ! get smallest prime factor of pj
1382 pfold = 1
1383 pp = pi
1384 DO ! loop over all prime factors of pi which are larger than
1385 ! the smallest prime factor of pj
1386 IF (pp.LE.1) EXIT ! if pj has been reduced to 1
1387 ! get smallest prime factor of pp, i.e. pj
1388 pfnew = getfactor(pp)
1389 IF ((pfnew.GT.pfmin).AND.(pfmin.NE.1)) EXIT
1390 pp = pp/pfnew
1391 ! skip multiple prime factors
1392 IF (pfnew.NE.pfold) THEN
1393 ! create new configuration
1394 p1 = p1*pfold/pfnew
1395 p2 = p2/pfold*pfnew
1396 pfold = pfnew
1397 ! get the best configuration possible with p1 x p2 x p3 processes
1398 ! by reducing p1 => recursion
1399 pinew = p1
1400 pjnew = p2
1401 pknew = p3
1402 CALL decompose(pinew,pjnew,pknew)
1403 CALL getcosts(ni,nj,nk,pinew,pjnew,pknew,blnew,vlnew)
1404 ! do the same with p2 <-> p3 interchanged
1405 pinew_= p1
1406 pjnew_= p3
1407 pknew_= p2
1408 CALL decompose(pinew_,pjnew_,pknew_)
1409 CALL getcosts(ni,nj,nk,pinew_,pjnew_,pknew_,blnew_,vlnew_)
1410 ! check which of the 2 configurations is better
1411 bl_gain = blnew*(1.0/blnew_) ! smaller is better
1412 vl_gain = vlnew_*(1.0/vlnew) ! larger is better
1413!!$PRINT '(4I7)',pjold,pkold,bl,vl
1414!!$PRINT '(4I7)',pjnew,pknew,blnew,vlnew
1415!!$PRINT *,"--------------------------------------------"
1416!!$PRINT '(A,3F7.1)'," ",bl_gain,vl_gain,bl_gain*vl_gain
1417 IF (vl_gain*bl_gain.GT.1.0) THEN
1418 ! 2nd configuration is better
1419 pinew = pinew_
1420 pjnew = pjnew_
1421 pknew = pknew_
1422 blnew = blnew_
1423 vlnew = vlnew_
1424 END IF
1425 ! check if the old configuration is better
1426 bl_gain = bl*(1.0/blnew) ! smaller is better
1427 vl_gain = vlnew*(1.0/vl) ! larger is better
1428 IF (vl_gain*bl_gain.GT.1.0) THEN
1429 ! new configuration is better
1430 piold = pinew
1431 pjold = pjnew
1432 pkold = pknew
1433 bl = blnew
1434 vl = vlnew
1435 END IF
1436 END IF
1437 END DO
1438 ! return optimized number of processes in the first dimension
1439 ! (for the initial configuration pj x pk)
1440 pj=pjold
1441 pk=pkold
1442 pi=ptot/pj/pk
1443 END SUBROUTINE decompose
1444
1445 ! computes the sum of the length of all internal boundaries (bl)
1446 ! and the mean vector length (vl)
1447 PURE SUBROUTINE getcosts(n1,n2,n3,p1,p2,p3,bl,vl)
1448 IMPLICIT NONE
1449 !------------------------------------------------------------------------!
1450 INTEGER, INTENT(IN) :: n1,n2,n3,p1,p2,p3
1451 INTEGER, INTENT(OUT) :: bl,vl
1452 !------------------------------------------------------------------------!
1453 INTEGER :: num
1454 !------------------------------------------------------------------------!
1455 ! length of internal boundaries
1456 bl = n3*n2*(p1-1) + n1*n3*(p2-1) + n1*n2*(p3-1)
1457 ! estimate the average vector length if the system vector length is large than 1
1458 IF (svl.GT.1) THEN
1459 ! get the maximal number of cells along the first dimension
1460 ! including ghost cells on both ends
1461 num = n1 / p1 + 2*ginum
1462 IF (mod(n1,p1).NE.0) num = num + 1 ! if the remainder is not zero add 1
1463 IF (num.LE.svl) THEN
1464 ! num fits into one vector
1465 vl = num
1466 ELSE
1467 ! we need more the one vector which may be filled completely or only partly
1468 ! return the closest integer to the average filling
1469 vl = num / ((num / svl) + 1)
1470 END IF
1471 ELSE
1472 vl = 1
1473 END IF
1474 END SUBROUTINE getcosts
1475
1476 END SUBROUTINE calculatedecomposition
1477#endif
1478!---------------------END PARALLEL-------------------------------------------!
1479
1480
1488 FUNCTION remapbounds_1(this,array) RESULT(ptr)
1489 IMPLICIT NONE
1490 !------------------------------------------------------------------------!
1491 CLASS(mesh_base) :: this
1492 REAL, DIMENSION(this%IGMIN:,this%JGMIN:,this%KGMIN:), TARGET :: array
1493 !------------------------------------------------------------------------!
1494 REAL, DIMENSION(:,:,:), POINTER :: ptr
1495 !------------------------------------------------------------------------!
1496 INTENT(IN) :: array
1497 !------------------------------------------------------------------------!
1498 ptr => array
1499 END FUNCTION remapbounds_1
1500
1502 FUNCTION remapbounds_2(this,array) RESULT(ptr)
1503 IMPLICIT NONE
1504 !------------------------------------------------------------------------!
1505 CLASS(mesh_base) :: this
1506 REAL, DIMENSION(this%IGMIN:,this%JGMIN:,this%KGMIN:,:), TARGET &
1507 :: array
1508 !------------------------------------------------------------------------!
1509 REAL, DIMENSION(:,:,:,:), POINTER :: ptr
1510 !------------------------------------------------------------------------!
1511 INTENT(IN) :: array
1512 !------------------------------------------------------------------------!
1513 ptr => array
1514 END FUNCTION remapbounds_2
1515
1517 FUNCTION remapbounds_3(this,array) RESULT(ptr)
1518 IMPLICIT NONE
1519 !------------------------------------------------------------------------!
1520 CLASS(mesh_base) :: this
1521 REAL, DIMENSION(this%IGMIN:,this%JGMIN:,this%KGMIN:,:,:), TARGET &
1522 :: array
1523 !------------------------------------------------------------------------!
1524 REAL, DIMENSION(:,:,:,:,:), POINTER :: ptr
1525 !------------------------------------------------------------------------!
1526 INTENT(IN) :: array
1527 !------------------------------------------------------------------------!
1528 ptr => array
1529 END FUNCTION remapbounds_3
1530
1532 FUNCTION remapbounds_4(this,array) RESULT(ptr)
1533 IMPLICIT NONE
1534 !------------------------------------------------------------------------!
1535 CLASS(mesh_base) :: this
1536 REAL, DIMENSION(this%IGMIN:,this%JGMIN:,this%KGMIN:,:,:,:), TARGET &
1537 :: array
1538 !------------------------------------------------------------------------!
1539 REAL, DIMENSION(:,:,:,:,:,:), POINTER :: ptr
1540 !------------------------------------------------------------------------!
1541 INTENT(IN) :: array
1542 !------------------------------------------------------------------------!
1543 ptr => array
1544 END FUNCTION remapbounds_4
1545
1546
1548 SUBROUTINE finalize_base(this)
1549 IMPLICIT NONE
1550 !------------------------------------------------------------------------!
1551 CLASS(mesh_base),INTENT(INOUT) :: this
1552 !------------------------------------------------------------------------!
1553#ifdef PARALLEL
1554 INTEGER :: i,ierror
1555 !------------------------------------------------------------------------!
1556 IF (.NOT.this%Initialized()) &
1557 CALL this%Error("CloseMesh","not initialized")
1558
1559 DO i=1,6
1560 IF (this%comm_boundaries(i).NE.mpi_comm_null) &
1561 CALL mpi_comm_free(this%comm_boundaries(i),ierror)
1562 END DO
1563#endif
1564
1565 IF (ASSOCIATED(this%rotation)) DEALLOCATE(this%rotation)
1566
1567 CALL this%without_ghost_zones%Destroy()
1568
1569 CALL closemeshproperties
1570
1571 CALL this%Geometry%Finalize()
1572 DEALLOCATE(this%Geometry)
1573 END SUBROUTINE finalize_base
1574
1577 PURE FUNCTION getdirection(this) RESULT(dir)
1578 IMPLICIT NONE
1579 !------------------------------------------------------------------------!
1580 CLASS(fargo_base), INTENT(IN) :: this
1581 INTEGER :: dir
1582 !------------------------------------------------------------------------!
1583 dir = this%direction
1584 END FUNCTION getdirection
1585
1588 PURE FUNCTION getdirectionname(this) RESULT(dir)
1589 IMPLICIT NONE
1590 !------------------------------------------------------------------------!
1591 CLASS(fargo_base), INTENT(IN) :: this
1592 CHARACTER(LEN=4) :: dir, all_dir(3)
1593 !------------------------------------------------------------------------!
1594 all_dir = [CHARACTER(LEN=1) :: "x","y","z" ]
1595 SELECT CASE(this%direction)
1596 CASE(1,2,3)
1597 dir = all_dir(this%direction)
1598 CASE DEFAULT
1599 dir = "none"
1600 END SELECT
1601 END FUNCTION getdirectionname
1602
1603
1604END MODULE mesh_base_mod
getcosts
pure subroutine getcosts(n1, n2, n3, p1, p2, p3, bl, vl)
Definition: mesh_base.f90:1448
decompose
recursive subroutine decompose(pi, pj, pk)
searches for the best domain decomposition
Definition: mesh_base.f90:1356
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
factors
module for prime factorization
Definition: factors.f90:29
factors::getfactor
pure integer function, public getfactor(number)
Definition: factors.f90:63
factors::maxnum
integer, parameter, public maxnum
Definition: factors.f90:49
geometry_base_mod
base class for geometrical properties
Definition: geometry_base.f90:36
geometry_base_mod::pi
real, parameter, public pi
Definition: geometry_base.f90:43
geometry_base_mod::spherical
integer, parameter, public spherical
Definition: geometry_base.f90:359
geometry_base_mod::tancylindrical
integer, parameter, public tancylindrical
Definition: geometry_base.f90:357
geometry_base_mod::cartesian
integer, parameter, public cartesian
Definition: geometry_base.f90:345
geometry_base_mod::cylindrical
integer, parameter, public cylindrical
Definition: geometry_base.f90:355
geometry_base_mod::finalize_base
subroutine finalize_base(this)
Destructor of generic geometry module.
Definition: geometry_base.f90:623
geometry_base_mod::logcylindrical
integer, parameter, public logcylindrical
Definition: geometry_base.f90:356
geometry_base_mod::spherical_planet
integer, parameter, public spherical_planet
Definition: geometry_base.f90:361
geometry_base_mod::logspherical
integer, parameter, public logspherical
Definition: geometry_base.f90:360
geometry_generic_mod
constructor for geometry class
Definition: geometry_generic.f90:35
geometry_generic_mod::new_geometry
subroutine new_geometry(Geometry, config)
Definition: geometry_generic.f90:53
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_base_mod::closemeshproperties
subroutine, public closemeshproperties
unsets global mesh properties
Definition: marray_base.f90:247
marray_base_mod::igmin
integer, save igmin
Definition: marray_base.f90:42
marray_base_mod::knum
integer, save knum
array sizes
Definition: marray_base.f90:45
marray_base_mod::remapbounds_1
real function, dimension(:,:,:,:), pointer remapbounds_1(this, array)
remap lower bounds in the first 3 dimensions of rank 1 mesh arrays
Definition: marray_base.f90:316
marray_base_mod::initmeshproperties
subroutine, public initmeshproperties(igmin_, igmax_, jgmin_, jgmax_, kgmin_, kgmax_)
sets global mesh properties
Definition: marray_base.f90:225
marray_base_mod::inum
integer, save inum
Definition: marray_base.f90:45
marray_base_mod::jnum
integer, save jnum
Definition: marray_base.f90:45
marray_base_mod::igmax
integer, save igmax
1st dim
Definition: marray_base.f90:42
marray_base_mod::kgmin
integer, save kgmin
Definition: marray_base.f90:44
marray_base_mod::jgmax
integer, save jgmax
2nd dim
Definition: marray_base.f90:43
marray_base_mod::remapbounds_2
real function, dimension(:,:,:,:,:), pointer remapbounds_2(this, array)
remap lower bounds in the first 3 dimensions of rank 2 mesh arrays
Definition: marray_base.f90:331
marray_base_mod::kgmax
integer, save kgmax
3rd dim
Definition: marray_base.f90:44
marray_base_mod::jgmin
integer, save jgmin
Definition: marray_base.f90:43
marray_cellscalar_mod
derived mesh array class for scalar cell data
Definition: marray_cellscalar.f90:33
marray_cellvector_mod
derived mesh array class for vector cell data
Definition: marray_cellvector.f90:33
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
mesh_base_mod::ncorners
integer, dimension(3), parameter ncorners
number of corners
Definition: mesh_base.f90:99
mesh_base_mod::east
integer, parameter east
named constant for eastern boundary
Definition: mesh_base.f90:101
mesh_base_mod::bottom
integer, parameter bottom
named constant for bottom boundary
Definition: mesh_base.f90:101
mesh_base_mod::initmesh
subroutine initmesh(this, config, IO, mtype, mname)
Constructor of generic mesh module.
Definition: mesh_base.f90:296
mesh_base_mod::calculatedecomposition
subroutine calculatedecomposition(ni, nj, nk, ginum, pi, pj, pk)
return the best partitioning of processes
Definition: mesh_base.f90:1321
mesh_base_mod::vector_y
integer, parameter vector_y
Definition: mesh_base.f90:109
mesh_base_mod::initmesh_parallel
subroutine initmesh_parallel(this, config)
Initialize MPI (parallel mode only)
Definition: mesh_base.f90:1023
mesh_base_mod::getdirection
pure integer function getdirection(this)
Get the fargo transport direction.
Definition: mesh_base.f90:1578
mesh_base_mod::nfaces
integer, dimension(3), parameter nfaces
number of faces
Definition: mesh_base.f90:98
mesh_base_mod::vector_z
integer, parameter vector_z
Definition: mesh_base.f90:109
mesh_base_mod::internalpoint
pure logical function internalpoint(this, x, y, z, mask)
Check if a given coordinate pair represents an internal point.
Definition: mesh_base.f90:962
mesh_base_mod::south
integer, parameter south
named constant for southern boundary
Definition: mesh_base.f90:101
mesh_base_mod::midpoint
integer, parameter midpoint
use midpoint rule to approximate flux integrals
Definition: mesh_base.f90:95
mesh_base_mod::remapbounds_4
real function, dimension(:,:,:,:,:,:), pointer remapbounds_4(this, array)
remap lower bounds in the first 2 dimensions of rank 5 subarrays
Definition: mesh_base.f90:1533
mesh_base_mod::getdirectionname
pure character(len=4) function getdirectionname(this)
Get the fargo transport direction as string.
Definition: mesh_base.f90:1589
mesh_base_mod::vector_x
integer, parameter vector_x
flags to check which vector components are enabled
Definition: mesh_base.f90:109
mesh_base_mod::setoutput
subroutine setoutput(this, config, IO)
Setup mesh fields for i/o.
Definition: mesh_base.f90:828
mesh_base_mod::remapbounds_3
real function, dimension(:,:,:,:,:), pointer remapbounds_3(this, array)
remap lower bounds in the first 2 dimensions of rank 4 subarrays
Definition: mesh_base.f90:1518
mesh_base_mod::calculaterotation
subroutine calculaterotation(this)
initialize array for rotation angle
Definition: mesh_base.f90:936
mesh_base_mod::top
integer, parameter top
named constant for top boundary
Definition: mesh_base.f90:101
mesh_base_mod::north
integer, parameter north
named constant for northern boundary
Definition: mesh_base.f90:101
mesh_base_mod::west
integer, parameter west
named constant for western boundary
Definition: mesh_base.f90:101
omega
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
Definition: sources_viscosity.f90:432
common_dict::dict_typ
Definition: common_dict.f90:123
geometry_base_mod::geometry_base
Definition: geometry_base.f90:51
logging_base_mod::logging_base
common data structure
Definition: logging_base.f90:56
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:69
marray_base_mod::selection_base
type for selecting parts of an marray
Definition: marray_base.f90:48
marray_cellscalar_mod::marray_cellscalar
data types and methods
Definition: marray_cellscalar.f90:39
marray_cellvector_mod::marray_cellvector
data types and methods
Definition: marray_cellvector.f90:39
mesh_base_mod::fargo_base
base class for fargo transport properties and methods
Definition: mesh_base.f90:114
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
Generated by   doxygen 1.9.4