60 CHARACTER(LEN=32),
PARAMETER ::
mesh_name =
"midpoint"
86 this%dAx(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
87 this%dAy(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
88 this%dAz(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
89 this%dAxdydz(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
90 this%dAydzdx(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
91 this%dAzdxdy(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
94 CALL this%Error(
"InitMesh_midpoint",
"Unable to allocate memory.")
97 DO concurrent(i=this%IGMIN:this%IGMAX,j=this%JGMIN:this%JGMAX,k=this%KGMIN:this%KGMAX)
100 this%dAxdydz(i,j,k,1:2) = this%hy%faces(i,j,k,1:2)*this%hz%faces(i,j,k,1:2)
102 this%dAydzdx(i,j,k,1:2) = this%hz%faces(i,j,k,3:4)*this%hx%faces(i,j,k,3:4)
104 this%dAzdxdy(i,j,k,1:2) = this%hx%faces(i,j,k,5:6)*this%hy%faces(i,j,k,5:6)
107 this%dAx(i,j,k,:) = this%dAxdydz(i,j,k,:)*this%dy*this%dz
108 this%dAy(i,j,k,:) = this%dAydzdx(i,j,k,:)*this%dz*this%dx
109 this%dAz(i,j,k,:) = this%dAzdxdy(i,j,k,:)*this%dx*this%dy
112 this%volume%data3d(i,j,k) = this%sqrtg%center(i,j,k)*this%dx*this%dy*this%dz
115 this%dxdydV%data3d(i,j,k) = this%dx*this%dy/(this%volume%data3d(i,j,k)+tiny(this%dx))
116 this%dydzdV%data3d(i,j,k) = this%dy*this%dz/(this%volume%data3d(i,j,k)+tiny(this%dy))
117 this%dzdxdV%data3d(i,j,k) = this%dz*this%dx/(this%volume%data3d(i,j,k)+tiny(this%dz))
120 this%cyxy%center(i,j,k) = 0.5*(this%hz%faces(i,j,k,2)+this%hz%faces(i,j,k,1)) &
121 * (this%hy%faces(i,j,k,2)-this%hy%faces(i,j,k,1)) * this%dydzdV%data3d(i,j,k)
122 this%cyzy%center(i,j,k) = 0.5*(this%hx%faces(i,j,k,6)+this%hx%faces(i,j,k,5)) &
123 * (this%hy%faces(i,j,k,6)-this%hy%faces(i,j,k,5)) * this%dxdydV%data3d(i,j,k)
124 this%cxyx%center(i,j,k) = 0.5*(this%hz%faces(i,j,k,4)+this%hz%faces(i,j,k,3)) &
125 * (this%hx%faces(i,j,k,4)-this%hx%faces(i,j,k,3)) * this%dzdxdV%data3d(i,j,k)
126 this%cxzx%center(i,j,k) = 0.5*(this%hy%faces(i,j,k,6)+this%hy%faces(i,j,k,5)) &
127 * (this%hx%faces(i,j,k,6)-this%hx%faces(i,j,k,5)) * this%dxdydV%data3d(i,j,k)
128 this%czxz%center(i,j,k) = 0.5*(this%hy%faces(i,j,k,2)+this%hy%faces(i,j,k,1)) &
129 * (this%hz%faces(i,j,k,2)-this%hz%faces(i,j,k,1)) * this%dydzdV%data3d(i,j,k)
130 this%czyz%center(i,j,k) = 0.5*(this%hx%faces(i,j,k,4)+this%hx%faces(i,j,k,3)) &
131 * (this%hz%faces(i,j,k,4)-this%hz%faces(i,j,k,3)) * this%dzdxdV%data3d(i,j,k)
134 this%cyxy%bcenter(i,j,k) = this%cyxy%center(i,j,k)
135 this%cyzy%bcenter(i,j,k) = this%cyzy%center(i,j,k)
136 this%cxyx%bcenter(i,j,k) = this%cxyx%center(i,j,k)
137 this%cxzx%bcenter(i,j,k) = this%cxzx%center(i,j,k)
138 this%czxz%bcenter(i,j,k) = this%czxz%center(i,j,k)
139 this%czyz%bcenter(i,j,k) = this%czyz%center(i,j,k)
142 this%dlx%data3d(i,j,k) = this%hx%center(i,j,k)*this%dx
143 this%dly%data3d(i,j,k) = this%hy%center(i,j,k)*this%dy
144 this%dlz%data3d(i,j,k) = this%hz%center(i,j,k)*this%dz
160 REAL,
DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
169 SELECT CASE(
this%VECTOR_COMPONENTS)
176 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
177 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
179 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),0.0, &
181 0.5*(vx(i-
this%IP1,j,k)+vx(i,j,k)),0.5*(vx(i+
this%IP1,j,k)+vx(i,j,k)), &
182 0.5*(vy(i,j-
this%JP1,k)+vy(i,j,k)),0.5*(vy(i,j+
this%JP1,k)+vy(i,j,k)), &
196 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
198 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
199 this%dydzdV%data3d(i,j,k),0.0,
this%dxdydV%data3d(i,j,k), &
201 0.5*(vx(i-
this%IP1,j,k)+vx(i,j,k)),0.5*(vx(i+
this%IP1,j,k)+vx(i,j,k)), &
203 0.5*(vy(i,j,k-
this%KP1)+vy(i,j,k)),0.5*(vy(i,j,k)+vy(i,j,k+
this%KP1)), &
217 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
218 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
219 0.0,
this%dzdxdV%data3d(i,j,k),
this%dxdydV%data3d(i,j,k), &
222 0.5*(vx(i,j-
this%JP1,k)+vx(i,j,k)),0.5*(vx(i,j+
this%JP1,k)+vx(i,j,k)), &
223 0.5*(vy(i,j,k-
this%KP1)+vy(i,j,k)),0.5*(vy(i,j,k)+vy(i,j,k+
this%KP1)), &
283 REAL,
DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
284 :: txx,txy,tyx,tyy,divtx,divty
288 INTENT(IN) :: txx,txy,tyx,tyy
289 INTENT(OUT) :: divtx,divty
292 SELECT CASE(
this%VECTOR_COMPONENTS)
300 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
301 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
303 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),0.0, &
304 this%cxyx%center(i,j,k),0.0, &
305 this%cyxy%center(i,j,k),0.0, &
306 0.5*(txx(i-
this%IP1,j,k)+txx(i,j,k)),0.5*(txx(i+
this%IP1,j,k)+txx(i,j,k)), &
307 0.5*(txy(i,j-
this%JP1,k)+txy(i,j,k)),0.5*(txy(i,j+
this%JP1,k)+txy(i,j,k)), &
309 tyx(i,j,k),0.0,tyy(i,j,k),0.0)
314 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
315 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
317 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),0.0, &
318 -
this%cxyx%center(i,j,k),0.0, &
319 -
this%cyxy%center(i,j,k),0.0, &
320 0.5*(tyx(i-
this%IP1,j,k)+tyx(i,j,k)),0.5*(tyx(i+
this%IP1,j,k)+tyx(i,j,k)), &
321 0.5*(tyy(i,j-
this%JP1,k)+tyy(i,j,k)),0.5*(tyy(i,j+
this%JP1,k)+tyy(i,j,k)), &
323 txx(i,j,k),0.0,txy(i,j,k),0.0)
339 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
340 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
341 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
342 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),
this%dxdydV%data3d(i,j,k), &
343 this%cxyx%center(i,j,k),
this%cxzx%center(i,j,k), &
344 this%cyxy%center(i,j,k),
this%czxz%center(i,j,k), &
345 0.5*(txx(i-
this%IP1,j,k)+txx(i,j,k)),0.5*(txx(i+
this%IP1,j,k)+txx(i,j,k)), &
347 0.5*(txy(i,j,k-
this%KP1)+txy(i,j,k)),0.5*(txy(i,j,k+
this%KP1)+txy(i,j,k)), &
348 0.0,tyx(i,j,k),0.0,tyy(i,j,k))
353 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
354 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
355 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
356 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),
this%dxdydV%data3d(i,j,k), &
357 this%czyz%center(i,j,k),-
this%cxzx%center(i,j,k), &
358 this%cyzy%center(i,j,k),-
this%czxz%center(i,j,k), &
359 0.5*(tyx(i-
this%IP1,j,k)+tyx(i,j,k)),0.5*(tyx(i+
this%IP1,j,k)+tyx(i,j,k)), &
361 0.5*(tyy(i,j,k-
this%KP1)+tyy(i,j,k)),0.5*(tyy(i,j,k+
this%KP1)+tyy(i,j,k)), &
362 0.0,txx(i,j,k),0.0,txy(i,j,k))
382 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
383 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
384 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
385 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),
this%dxdydV%data3d(i,j,k), &
386 -
this%cxyx%center(i,j,k),
this%cyzy%center(i,j,k), &
387 -
this%cyxy%center(i,j,k),
this%czyz%center(i,j,k), &
389 0.5*(txx(i,j-
this%JP1,k)+txx(i,j,k)),0.5*(txx(i,j+
this%JP1,k)+tyy(i,j,k)), &
390 0.5*(txy(i,j,k-
this%KP1)+txy(i,j,k)),0.5*(txy(i,j,k+
this%KP1)+txy(i,j,k)), &
391 0.0,tyx(i,j,k),0.0,tyy(i,j,k))
396 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
397 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
398 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
399 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),
this%dxdydV%data3d(i,j,k), &
400 this%czyz%center(i,j,k),-
this%cxzx%center(i,j,k), &
401 this%cyzy%center(i,j,k),-
this%czxz%center(i,j,k), &
403 0.5*(tyx(i,j-
this%JP1,k)+tyx(i,j,k)),0.5*(tyx(i,j+
this%JP1,k)+tyx(i,j,k)), &
404 0.5*(tyy(i,j,k-
this%KP1)+tyy(i,j,k)),0.5*(tyy(i,j,k+
this%KP1)+tyy(i,j,k)), &
405 txy(i,j,k),0.0,txx(i,j,k),0.0)
478 REAL,
DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
483 INTENT(IN) :: vx,vy,vz
493 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
494 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
495 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
496 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),
this%dxdydV%data3d(i,j,k), &
498 0.5*(vx(i-
this%IP1,j,k)+vx(i,j,k)),0.5*(vx(i+
this%IP1,j,k)+vx(i,j,k)), &
499 0.5*(vy(i,j-
this%JP1,k)+vy(i,j,k)),0.5*(vy(i,j+
this%JP1,k)+vy(i,j,k)), &
500 0.5*(vz(i,j,k-
this%KP1)+vz(i,j,k)),0.5*(vz(i,j,k+
this%KP1)+vz(i,j,k)), &
517 PURE SUBROUTINE tensordivergence3d(this,Txx,Txy,Txz,Tyx,Tyy,Tyz,Tzx,Tzy,Tzz, &
522 REAL,
DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
523 :: txx,txy,txz,tyx,tyy,tyz,tzx,tzy,tzz,divtx,divty,divtz
527 INTENT(IN) :: txx,txy,txz,tyx,tyy,tyz,tzx,tzy,tzz
528 INTENT(OUT) :: divtx,divty,divtz
536 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
537 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
538 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
539 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),
this%dxdydV%data3d(i,j,k), &
540 this%cxyx%center(i,j,k),
this%cxzx%center(i,j,k), &
541 this%cyxy%center(i,j,k),
this%czxz%center(i,j,k), &
542 0.5*(txx(i-
this%IP1,j,k)+txx(i,j,k)),0.5*(txx(i+
this%IP1,j,k)+txx(i,j,k)), &
543 0.5*(txy(i,j-
this%JP1,k)+txy(i,j,k)),0.5*(txy(i,j+
this%JP1,k)+txy(i,j,k)), &
544 0.5*(txz(i,j,k-
this%KP1)+txz(i,j,k)),0.5*(txz(i,j,k+
this%KP1)+txz(i,j,k)), &
545 tyx(i,j,k),tzx(i,j,k),tyy(i,j,k),tzz(i,j,k))
550 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
551 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
552 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
553 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),
this%dxdydV%data3d(i,j,k), &
554 -
this%cxyx%center(i,j,k),
this%cyzy%center(i,j,k), &
555 -
this%cyxy%center(i,j,k),
this%czyz%center(i,j,k), &
556 0.5*(tyx(i-
this%IP1,j,k)+tyx(i,j,k)),0.5*(tyx(i+
this%IP1,j,k)+tyx(i,j,k)), &
557 0.5*(tyy(i,j-
this%JP1,k)+tyy(i,j,k)),0.5*(tyy(i,j+
this%JP1,k)+tyy(i,j,k)), &
558 0.5*(tyz(i,j,k-
this%KP1)+tyz(i,j,k)),0.5*(tyz(i,j,k+
this%KP1)+tyz(i,j,k)), &
559 txx(i,j,k),tzy(i,j,k),txy(i,j,k),tzz(i,j,k))
564 this%dAxdydz(i,j,k,1),
this%dAxdydz(i,j,k,2), &
565 this%dAydzdx(i,j,k,1),
this%dAydzdx(i,j,k,2), &
566 this%dAzdxdy(i,j,k,1),
this%dAzdxdy(i,j,k,2), &
567 this%dydzdV%data3d(i,j,k),
this%dzdxdV%data3d(i,j,k),
this%dxdydV%data3d(i,j,k), &
568 this%czyz%center(i,j,k),-
this%cxzx%center(i,j,k), &
569 this%cyzy%center(i,j,k),-
this%czxz%center(i,j,k), &
570 0.5*(tzx(i-
this%IP1,j,k)+tzx(i,j,k)),0.5*(tzx(i+
this%IP1,j,k)+tzx(i,j,k)), &
571 0.5*(tzy(i,j-
this%JP1,k)+tzy(i,j,k)),0.5*(tzy(i,j+
this%JP1,k)+tzy(i,j,k)), &
572 0.5*(tzz(i,j,k-
this%KP1)+tzz(i,j,k)),0.5*(tzz(i,j,k+
this%KP1)+tzz(i,j,k)), &
573 tyz(i,j,k),txx(i,j,k),tyy(i,j,k),txz(i,j,k))
586 ELEMENTAL FUNCTION divergence3d(dAxdydzWest,dAxdydzEast,dAydxdzSouth,dAydxdzNorth,dAzdxdyBottom,dAzdxdyTop, &
587 dydzdV,dxdzdV,dxdydV,cxyx,cxzx,cyxy,czxz, &
588 TxxWest,TxxEast,TxySouth,TxyNorth,TxzBottom,TxzTop, &
589 TyxCent,TzxCent,TyyCent,TzzCent)
RESULT(div)
592 REAL,
INTENT(IN) :: daxdydzwest,daxdydzeast, &
593 daydxdzsouth,daydxdznorth, &
594 dazdxdybottom,dazdxdytop, &
595 dxdydv,dxdzdv,dydzdv, &
596 cxyx,cyxy,czxz,cxzx, &
600 tyxcent,tyycent,tzzcent,tzxcent
626 div = dydzdv*(daxdydzeast*txxeast-daxdydzwest*txxwest) &
627 + dxdzdv*(daydxdznorth*txynorth-daydxdzsouth*txysouth) &
628 + dxdydv*(dazdxdytop*txztop-dazdxdybottom*txzbottom) &
629 + cxyx*tyxcent + cxzx*tzxcent - cyxy*tyycent - czxz*tzzcent
639 DEALLOCATE(this%dAx,this%dAy,this%dAz,this%dAxdydz,this%dAydzdx,this%dAzdxdy)
640 CALL this%Finalize_base()
Dictionary for generic data types.
type(logging_base), save this
real, save nan_default_real
NaN real constant.
subroutine finalize(this)
Destructor of logging_base class.
integer, save default_mpi_real
default real type for MPI
integer, parameter vector_y
integer, parameter vector_z
integer, parameter midpoint
use midpoint rule to approximate flux integrals
integer, parameter vector_x
flags to check which vector components are enabled
mesh module for midpoint quadrature rule
pure subroutine vectordivergence2d_1(this, vx, vy, divv)
compute the cell centered 2D vector divergence
subroutine initmesh_midpoint(this, config, IO)
pure subroutine vectordivergence3d(this, vx, vy, vz, divv)
compute the cell centered 3D vector divergence
elemental real function divergence3d(dAxdydzWest, dAxdydzEast, dAydxdzSouth, dAydxdzNorth, dAzdxdyBottom, dAzdxdyTop, dydzdV, dxdzdV, dxdydV, cxyx, cxzx, cyxy, czxz, TxxWest, TxxEast, TxySouth, TxyNorth, TxzBottom, TxzTop, TyxCent, TzxCent, TyyCent, TzzCent)
elemental workhorse to compute divergence
pure subroutine tensordivergence2d_1(this, Txx, Txy, Tyx, Tyy, divTx, divTy)
compute the cell centered 2D rank 2 tensor divergence
pure subroutine tensordivergence3d(this, Txx, Txy, Txz, Tyx, Tyy, Tyz, Tzx, Tzy, Tzz, divTx, divTy, divTz)
compute the cell centered 3D rank 2 tensor divergence
character(len=32), parameter mesh_name