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
109 INTEGER,
PARAMETER :: &
115 INTEGER :: direction = 0
128 INTEGER :: ginum,gjnum,gknum
139 INTEGER :: ip1,ip2,im1,im2
140 INTEGER :: jp1,jp2,jm1,jm2
141 INTEGER :: kp1,kp2,km1,km2
146 REAL :: invdx,invdy,invdz
151 INTEGER :: vector_components
157 REAL,
DIMENSION(:,:,:,:),
POINTER :: &
161 REAL,
DIMENSION(:,:,:,:,:),
POINTER :: &
166 TYPE(
marray_base) :: dlx,dly,dlz, & !< cell centered line elements
167 volume, & !< cell volumes
169 REAL,
DIMENSION(:,:,:,:),
POINTER :: &
170 dax,day,daz, & !< cell surface area elements
171 daxdydz,daydzdx,dazdxdy
175 sqrtg, & !< sqrt(det(metric))
176 cyxy,cyzy,cxzx,cxyx,czxz,czyz
183 REAL,
DIMENSION(:,:,:),
POINTER :: &
185 REAL,
DIMENSION(:,:,:,:,:,:),
POINTER :: &
189 INTEGER :: maxinum,maxjnum,maxknum
190 INTEGER :: mininum,minjnum,minknum
192 INTEGER :: icomm,jcomm,kcomm
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
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, &
219 PURE SUBROUTINE tensordivergence3d(this,Txx,Txy,Txz,Tyx,Tyy,Tyz,Tzx,Tzy,Tzz, &
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
230 PURE SUBROUTINE vectordivergence3d(this,vx,vy,vz,divv)
235 REAL,
DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
238 INTENT(IN) :: vx,vy,vz
241 PURE SUBROUTINE vectordivergence2d_1(this,vx,vy,divv)
246 REAL,
DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
252 PURE SUBROUTINE tensordivergence2d_1(this,Txx,Txy,Tyx,Tyy,divTx,divTy)
257 REAL,
DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
258 :: txx,txy,tyx,tyy,divtx,divty
260 INTENT(IN) :: txx,txy,tyx,tyy
261 INTENT(OUT) :: divtx,divty
303 CHARACTER(LEN=32) :: mname
305 CHARACTER(LEN=32) :: xres,yres,zres,somega,fargo_method(4)
306 INTEGER :: fargo,fargo_dir
308 REAL :: mesh_dx,mesh_dy,mesh_dz
310 INTEGER :: inum,jnum,knum,err
312 REAL,
DIMENSION(6,3) :: cfaces
313 REAL,
DIMENSION(8,3) :: ccorners
315 INTENT(INOUT) :: this
317 CALL this%logging_base%InitLogging(mtype,mname)
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")
338 CALL getattr(config,
"omega", this%omega, 0.0)
341 CALL getattr(config,
"Qshear", this%Q, 1.5)
344 this%rotcent = (/ 0., 0., 0./)
345 CALL getattr(config,
"rotation_center", this%rotcent, this%rotcent)
351 CALL getattr(config,
"inum", this%INUM)
352 CALL getattr(config,
"jnum", this%JNUM)
353 CALL getattr(config,
"knum", this%KNUM)
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.")
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)
376 IF(this%shear_dir.GT.0) fargo = 3
378 CALL getattr(config,
"fargo/method", fargo, fargo)
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))
390 fargo_dir = this%Geometry%GetAzimuthIndex()
391 SELECT CASE(fargo_dir)
393 SELECT TYPE(mgeo => this%Geometry)
394 TYPE IS(geometry_cartesian)
395 IF (this%fargo%GetType().EQ.3)
THEN
398 this%fargo%direction = this%shear_dir
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}")
406 CALL this%Error(
"mesh_base::InitMesh",
"fargo transport not supported for this geometry")
411 this%fargo%direction = fargo_dir
414 CALL this%Error(
"mesh_base::InitMesh",
"wrong coordinate index returned from GetAzimuthIndex")
418 this%fargo%direction = 0
422 CALL this%Error(
"mesh_base::InitMesh",
"fargo method should be one of {0,1,2,3}")
430 SELECT CASE(this%geometry%GetAzimuthIndex())
432 IF (this%INUM.EQ.1.AND.(abs(this%xmax-this%xmin-2*
pi).LE.epsilon(this%xmin))) &
435 IF (this%JNUM.EQ.1.AND.(abs(this%ymax-this%ymin-2*
pi).LE.epsilon(this%ymin))) &
438 IF (this%KNUM.EQ.1.AND.(abs(this%zmax-this%zmin-2*
pi).LE.epsilon(this%zmin))) &
463 this%GINUM = this%GNUM
464 this%GJNUM = this%GNUM
465 this%GKNUM = this%GNUM
470 this%VECTOR_COMPONENTS = int(b
'111')
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)
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)
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)
511 this%NFACES =
nfaces(this%NDIMS)
512 this%NCORNERS =
ncorners(this%NDIMS)
521 CALL mpi_barrier(mpi_comm_world,err)
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)
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)
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)
536 this%IMAX = this%INUM
538 this%JMAX = this%JNUM
540 this%KMAX = this%KNUM
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
552 CALL initmeshproperties(this%IGMIN,this%IGMAX,this%JGMIN,this%JGMAX,this%KGMIN,this%KGMAX)
555 this%without_ghost_zones =
selection_base((/this%IMIN,this%IMAX,this%JMIN,this%JMAX,this%KMIN,this%KMAX/))
558 this%dx = (this%xmax - this%xmin) / this%INUM
560 this%dy = (this%ymax - this%ymin) / this%JNUM
562 this%dz = (this%zmax - this%zmin) / this%KNUM
567 IF (this%dx.LT.2*epsilon(this%dx))
THEN
568 IF (this%INUM.EQ.1)
THEN
572 CALL this%Error(
"mesh_base::InitMesh",
"INUM > 1 conflicts with zero x-extent")
575 IF (this%dy.LT.2*epsilon(this%dy))
THEN
576 IF (this%JNUM.EQ.1)
THEN
580 CALL this%Error(
"mesh_base::InitMesh",
"JNUM > 1 conflicts with zero y-extent")
583 IF (this%dz.LT.2*epsilon(this%dz))
THEN
584 IF (this%KNUM.EQ.1)
THEN
588 CALL this%Error(
"mesh_base::InitMesh",
"KNUM > 1 conflicts with zero z-extent")
593 this%invdx = 1./this%dx
594 this%invdy = 1./this%dy
595 this%invdz = 1./this%dz
599 this%center => this%curv%RemapBounds(this%curv%center)
600 this%bcenter => this%curv%RemapBounds(this%curv%bcenter)
604 cfaces(1,1) = -0.5*mesh_dx
607 cfaces(2,1) = 0.5*mesh_dx
611 cfaces(3,2) = -0.5*mesh_dy
614 cfaces(4,2) = 0.5*mesh_dy
618 cfaces(5,3) = -0.5*mesh_dz
621 cfaces(6,3) = 0.5*mesh_dz
623 ccorners(1,1) = cfaces(1,1)
624 ccorners(1,2) = cfaces(3,2)
625 ccorners(1,3) = cfaces(5,3)
626 ccorners(2,1) = cfaces(2,1)
627 ccorners(2,2) = cfaces(3,2)
628 ccorners(2,3) = cfaces(5,3)
629 ccorners(3,1) = cfaces(1,1)
630 ccorners(3,2) = cfaces(4,2)
631 ccorners(3,3) = cfaces(5,3)
632 ccorners(4,1) = cfaces(2,1)
633 ccorners(4,2) = cfaces(4,2)
634 ccorners(4,3) = cfaces(5,3)
635 ccorners(5,1) = cfaces(1,1)
636 ccorners(5,2) = cfaces(3,2)
637 ccorners(5,3) = cfaces(6,3)
638 ccorners(6,1) = cfaces(2,1)
639 ccorners(6,2) = cfaces(3,2)
640 ccorners(6,3) = cfaces(6,3)
641 ccorners(7,1) = cfaces(1,1)
642 ccorners(7,2) = cfaces(4,2)
643 ccorners(7,3) = cfaces(6,3)
644 ccorners(8,1) = cfaces(2,1)
645 ccorners(8,2) = cfaces(4,2)
646 ccorners(8,3) = cfaces(6,3)
648 DO concurrent(i=this%IGMIN:this%IGMAX,j=this%JGMIN:this%JGMAX,k=this%KGMIN:this%KGMAX)
650 this%curv%center(i,j,k,1) = this%xmin + (2*i-1)*0.5*mesh_dx
651 this%curv%center(i,j,k,2) = this%ymin + (2*j-1)*0.5*mesh_dy
652 this%curv%center(i,j,k,3) = this%zmin + (2*k-1)*0.5*mesh_dz
655 this%curv%bcenter(i,j,k,:) = this%curv%center(i,j,k,:)
657 this%curv%faces(i,j,k,1,:) = this%curv%center(i,j,k,:) + cfaces(1,:)
658 this%curv%faces(i,j,k,2,:) = this%curv%center(i,j,k,:) + cfaces(2,:)
659 this%curv%faces(i,j,k,3,:) = this%curv%center(i,j,k,:) + cfaces(3,:)
660 this%curv%faces(i,j,k,4,:) = this%curv%center(i,j,k,:) + cfaces(4,:)
661 this%curv%faces(i,j,k,5,:) = this%curv%center(i,j,k,:) + cfaces(5,:)
662 this%curv%faces(i,j,k,6,:) = this%curv%center(i,j,k,:) + cfaces(6,:)
664 this%curv%corners(i,j,k,1,:) = this%curv%center(i,j,k,:) + ccorners(1,:)
665 this%curv%corners(i,j,k,2,:) = this%curv%center(i,j,k,:) + ccorners(2,:)
666 this%curv%corners(i,j,k,3,:) = this%curv%center(i,j,k,:) + ccorners(3,:)
667 this%curv%corners(i,j,k,4,:) = this%curv%center(i,j,k,:) + ccorners(4,:)
668 this%curv%corners(i,j,k,5,:) = this%curv%center(i,j,k,:) + ccorners(5,:)
669 this%curv%corners(i,j,k,6,:) = this%curv%center(i,j,k,:) + ccorners(6,:)
670 this%curv%corners(i,j,k,7,:) = this%curv%center(i,j,k,:) + ccorners(7,:)
671 this%curv%corners(i,j,k,8,:) = this%curv%center(i,j,k,:) + ccorners(8,:)
683 CALL this%Geometry%ScaleFactors(this%curv,this%hx,this%hy,this%hz)
691#if defined (__INTEL_COMPILER) || defined (__flang__)
692 this%sqrtg%data1d(:) = this%hx%data1d(:)*(this%hy%data1d(:)*this%hz%data1d(:))
694 this%sqrtg = this%hx*(this%hy*this%hz)
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)
721 CALL this%Geometry%Convert2Cartesian(this%curv,this%cart)
725 CALL this%geometry%Radius(this%curv,this%radius)
729 CALL this%geometry%PositionVector(this%curv,this%posvec)
732 NULLIFY(this%dAx,this%dAy,this%dAz, &
733 this%dAxdydz,this%dAydzdx,this%dAzdxdy, &
734 this%rotation,this%weights)
785 CALL this%Info(
" MESH-----> quadrature rule: " // trim(this%GetName()))
786 WRITE (xres,
'(I0)') this%INUM
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))
800 SELECT CASE(this%ROTSYM)
804 somega =
"x-coordinate"
806 somega =
"y-coordinate"
808 somega =
"z-coordinate"
810 CALL this%Info(
" rot. symmetry: " // trim(somega))
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 ")
833 INTEGER :: writefields
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,:,:))
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))
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,:))
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,:))
862 CALL getattr(config,
"output/rotation", writefields, 0)
863 IF(writefields.EQ.1)
THEN
865 CALL setattr(io,
"rotation",this%rotation(this%IMIN:this%IMAX,this%JMIN:this%JMAX,this%KMIN:this%KMAX))
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))
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))
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))
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))
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,:))
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))
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,:))
941 REAL,
DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,3) :: cart,curv
943 ALLOCATE(this%rotation(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX), &
946 CALL this%Error(
"CalculateRotation",
"Couldn't allocate memory")
951 CALL this%geometry%Convert2Curvilinear(this%bcenter, cart, curv)
952 this%rotation(:,:,:) = atan2(curv(:,:,:,2),curv(:,:,:,1))
965 INTEGER,
DIMENSION(6),
OPTIONAL :: mask
969 INTEGER :: imin,imax,jmin,jmax,kmin,kmax
971 INTENT(IN) :: x,y,z,mask
983 IF (
PRESENT(mask))
THEN
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))
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
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
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
1028 CHARACTER(LEN=64) :: decomp_str
1029 LOGICAL,
DIMENSION(3) :: remain_dims, periods
1033 INTEGER :: worldgroup,newgroup
1034 INTEGER,
DIMENSION(1) :: rank0in, rank0out
1035 INTEGER,
DIMENSION(3) :: coords,ncells,dims
1036 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ranks
1057 CALL getattr(config,
"decomposition", dims, dims(1:3))
1061 IF (this%GetRank().EQ.0)
THEN
1063 WRITE (decomp_str,
'(3(A,I0),A)')
"[ ", dims(1),
", ", dims(2),
", ",dims(3),
" ]"
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")
1071 IF (any(dims(:).EQ.0)) &
1072 CALL this%Error(
"InitMesh_parallel",
"numbers in decomposition should not be 0")
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")
1080 ncells(1) = this%INUM
1081 ncells(2) = this%JNUM
1082 ncells(3) = this%KNUM
1085 WHERE (ncells(:).EQ.1)
1090 IF (all(dims(:).GT.0))
THEN
1092 ELSE IF (product(dims(:)).GT.0)
THEN
1095 k = maxloc(dims(:),1)
1101 dims(3) = this%GetNumProcs() / dims(k)
1105 dims(3),dims(2),dims(1))
1107 dims(1) = this%GetNumProcs() / dims(k)
1111 dims(1),dims(2),dims(3))
1114 ELSE IF (all(dims(:).LT.0))
THEN
1116 dims(1) = this%GetNumProcs()
1120 dims(1),dims(2),dims(3))
1124 k = minloc(dims(:),1)
1125 dims(k) = this%GetNumProcs() / product(dims(:),dims(:).GT.1)
1128 WRITE (decomp_str,
'(3(A,I0),A)')
"[ ", dims(1),
", ", dims(2),
", ",dims(3),
" ]"
1131 ncells(1) = this%INUM
1132 ncells(2) = this%JNUM
1133 ncells(3) = this%KNUM
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))
1141 IF (dims(3).LE.0)
THEN
1142 CALL this%Error(
"InitMesh_parallel",
"automatic domain decomposition failed.")
1145 this%dims(:) = dims(:)
1149 CALL mpi_bcast(this%dims,3,mpi_integer,0,mpi_comm_world,ierror)
1157 CALL mpi_cart_create(mpi_comm_world,3,this%dims,periods,.true., &
1158 this%comm_cart,ierror)
1161 this%mycoords(:) = 0
1162 CALL mpi_cart_get(this%comm_cart,3,this%dims,periods,this%mycoords,ierror)
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)
1174 rem = mod(this%INUM,this%dims(1))
1175 num = (this%INUM-rem) / this%dims(1)
1176 IF (this%mycoords(1).LT.rem)
THEN
1178 this%IMIN = this%mycoords(1) * (num + 1) + 1
1179 this%IMAX = this%IMIN + num
1181 this%IMIN = rem * (num+1) + (this%mycoords(1) - rem) * num + 1
1182 this%IMAX = this%IMIN + num - 1
1185 rem = mod(this%JNUM,this%dims(2))
1186 num = (this%JNUM-rem) / this%dims(2)
1187 IF (this%mycoords(2).LT.rem)
THEN
1189 this%JMIN = this%mycoords(2) * (num + 1) + 1
1190 this%JMAX = this%JMIN + num
1192 this%JMIN = rem * (num+1) + (this%mycoords(2) - rem) * num + 1
1193 this%JMAX = this%JMIN + num - 1
1196 rem = mod(this%KNUM,this%dims(3))
1197 num = (this%KNUM-rem) / this%dims(3)
1198 IF (this%mycoords(3).LT.rem)
THEN
1200 this%KMIN = this%mycoords(3) * (num + 1) + 1
1201 this%KMAX = this%KMIN + num
1203 this%KMIN = rem * (num+1) + (this%mycoords(3) - rem) * num + 1
1204 this%KMAX = this%KMIN + num - 1
1208 ALLOCATE(ranks(max(this%dims(1)*this%dims(2),this%dims(1)*this%dims(3),this%dims(2)*this%dims(3))))
1210 CALL mpi_comm_group(mpi_comm_world,worldgroup,ierror)
1213 DO k=0,this%dims(3)-1
1214 DO j=0,this%dims(2)-1
1217 CALL mpi_cart_rank(this%comm_cart,coords,i,ierror)
1218 ranks(j+1+k*this%dims(2))=i
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)
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)
1228 coords(1) = this%dims(1)-1
1229 DO k=0,this%dims(3)-1
1230 DO j=0,this%dims(2)-1
1233 CALL mpi_cart_rank(this%comm_cart,coords,i,ierror)
1234 ranks(j+1+k*this%dims(2))=i
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)
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)
1245 DO k=0,this%dims(3)-1
1246 DO i=0,this%dims(1)-1
1249 CALL mpi_cart_rank(this%comm_cart,coords,j,ierror)
1250 ranks(i+1+k*this%dims(1))=j
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)
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)
1260 coords(2) = this%dims(2)-1
1261 DO k=0,this%dims(3)-1
1262 DO i=0,this%dims(1)-1
1265 CALL mpi_cart_rank(this%comm_cart,coords,j,ierror)
1266 ranks(i+1+k*this%dims(1))=j
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)
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)
1277 DO j=0,this%dims(2)-1
1278 DO i=0,this%dims(1)-1
1281 CALL mpi_cart_rank(this%comm_cart,coords,k,ierror)
1282 ranks(i+1+j*this%dims(1))=k
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)
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)
1292 coords(3) = this%dims(3)-1
1293 DO j=0,this%dims(2)-1
1294 DO i=0,this%dims(1)-1
1297 CALL mpi_cart_rank(this%comm_cart,coords,k,ierror)
1298 ranks(i+1+j*this%dims(1))=k
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)
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)
1308 CALL mpi_group_free(worldgroup,ierror)
1324 INTEGER,
INTENT(IN) :: ni,nj,nk,ginum
1325 INTEGER,
INTENT(INOUT) :: pi,pj,pk
1327 CHARACTER(LEN=8) :: svl_char
1331 svl_char = vector_length
1332 READ (svl_char,
'(I8)',iostat=err) svl
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
1358 INTEGER,
INTENT(INOUT) :: pi,pj,pk
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
1368 IF (pj.GT.nj.OR.pk.GT.nk)
RETURN
1370 CALL getcosts(ni,nj,nk,pi,pj,pk,bl,vl)
1389 IF ((pfnew.GT.pfmin).AND.(pfmin.NE.1))
EXIT
1392 IF (pfnew.NE.pfold)
THEN
1403 CALL getcosts(ni,nj,nk,pinew,pjnew,pknew,blnew,vlnew)
1409 CALL getcosts(ni,nj,nk,pinew_,pjnew_,pknew_,blnew_,vlnew_)
1411 bl_gain = blnew*(1.0/blnew_)
1412 vl_gain = vlnew_*(1.0/vlnew)
1417 IF (vl_gain*bl_gain.GT.1.0)
THEN
1426 bl_gain = bl*(1.0/blnew)
1427 vl_gain = vlnew*(1.0/vl)
1428 IF (vl_gain*bl_gain.GT.1.0)
THEN
1450 INTEGER,
INTENT(IN) :: n1,n2,n3,p1,p2,p3
1451 INTEGER,
INTENT(OUT) :: bl,vl
1456 bl = n3*n2*(p1-1) + n1*n3*(p2-1) + n1*n2*(p3-1)
1461 num = n1 / p1 + 2*ginum
1462 IF (mod(n1,p1).NE.0) num = num + 1
1463 IF (num.LE.svl)
THEN
1469 vl = num / ((num / svl) + 1)
1492 REAL,
DIMENSION(this%IGMIN:,this%JGMIN:,this%KGMIN:),
TARGET :: array
1494 REAL,
DIMENSION(:,:,:),
POINTER :: ptr
1506 REAL,
DIMENSION(this%IGMIN:,this%JGMIN:,this%KGMIN:,:),
TARGET &
1509 REAL,
DIMENSION(:,:,:,:),
POINTER :: ptr
1521 REAL,
DIMENSION(this%IGMIN:,this%JGMIN:,this%KGMIN:,:,:),
TARGET &
1524 REAL,
DIMENSION(:,:,:,:,:),
POINTER :: ptr
1536 REAL,
DIMENSION(this%IGMIN:,this%JGMIN:,this%KGMIN:,:,:,:),
TARGET &
1539 REAL,
DIMENSION(:,:,:,:,:,:),
POINTER :: ptr
1556 IF (.NOT.this%Initialized()) &
1557 CALL this%Error(
"CloseMesh",
"not initialized")
1560 IF (this%comm_boundaries(i).NE.mpi_comm_null) &
1561 CALL mpi_comm_free(this%comm_boundaries(i),ierror)
1565 IF (
ASSOCIATED(this%rotation))
DEALLOCATE(this%rotation)
1567 CALL this%without_ghost_zones%Destroy()
1571 CALL this%Geometry%Finalize()
1572 DEALLOCATE(this%Geometry)
1583 dir =
this%direction
1592 CHARACTER(LEN=4) :: dir, all_dir(3)
1594 all_dir = [
CHARACTER(LEN=1) ::
"x",
"y",
"z" ]
1595 SELECT CASE(
this%direction)
1597 dir = all_dir(
this%direction)
pure subroutine getcosts(n1, n2, n3, p1, p2, p3, bl, vl)
recursive subroutine decompose(pi, pj, pk)
searches for the best domain decomposition
Dictionary for generic data types.
type(logging_base), save this
module for prime factorization
pure integer function, public getfactor(number)
integer, parameter, public maxnum
base class for geometrical properties
real, parameter, public pi
integer, parameter, public spherical
integer, parameter, public tancylindrical
integer, parameter, public cartesian
integer, parameter, public cylindrical
subroutine finalize_base(this)
Destructor of generic geometry module.
integer, parameter, public logcylindrical
integer, parameter, public spherical_planet
integer, parameter, public logspherical
constructor for geometry class
subroutine new_geometry(Geometry, config)
subroutine finalize(this)
Destructor of logging_base class.
integer, save default_mpi_real
default real type for MPI
base class for mesh arrays
subroutine, public closemeshproperties
unsets global mesh properties
integer, save knum
array sizes
real function, dimension(:,:,:,:), pointer remapbounds_1(this, array)
remap lower bounds in the first 3 dimensions of rank 1 mesh arrays
subroutine, public initmeshproperties(igmin_, igmax_, jgmin_, jgmax_, kgmin_, kgmax_)
sets global mesh properties
integer, save igmax
1st dim
integer, save jgmax
2nd dim
real function, dimension(:,:,:,:,:), pointer remapbounds_2(this, array)
remap lower bounds in the first 3 dimensions of rank 2 mesh arrays
integer, save kgmax
3rd dim
derived mesh array class for scalar cell data
derived mesh array class for vector cell data
integer, dimension(3), parameter ncorners
number of corners
integer, parameter east
named constant for eastern boundary
integer, parameter bottom
named constant for bottom boundary
subroutine initmesh(this, config, IO, mtype, mname)
Constructor of generic mesh module.
subroutine calculatedecomposition(ni, nj, nk, ginum, pi, pj, pk)
return the best partitioning of processes
integer, parameter vector_y
subroutine initmesh_parallel(this, config)
Initialize MPI (parallel mode only)
pure integer function getdirection(this)
Get the fargo transport direction.
integer, dimension(3), parameter nfaces
number of faces
integer, parameter vector_z
pure logical function internalpoint(this, x, y, z, mask)
Check if a given coordinate pair represents an internal point.
integer, parameter south
named constant for southern boundary
integer, parameter midpoint
use midpoint rule to approximate flux integrals
real function, dimension(:,:,:,:,:,:), pointer remapbounds_4(this, array)
remap lower bounds in the first 2 dimensions of rank 5 subarrays
pure character(len=4) function getdirectionname(this)
Get the fargo transport direction as string.
integer, parameter vector_x
flags to check which vector components are enabled
subroutine setoutput(this, config, IO)
Setup mesh fields for i/o.
real function, dimension(:,:,:,:,:), pointer remapbounds_3(this, array)
remap lower bounds in the first 2 dimensions of rank 4 subarrays
subroutine calculaterotation(this)
initialize array for rotation angle
integer, parameter top
named constant for top boundary
integer, parameter north
named constant for northern boundary
integer, parameter west
named constant for western boundary
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
type for selecting parts of an marray
base class for fargo transport properties and methods