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