mesh_midpoint.f90
Go to the documentation of this file.
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: mesh_midpoint.f90 #
5 !# #
6 !# Copyright (C) 2006-2019 #
7 !# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8 !# Manuel Jung #
9 !# #
10 !# This program is free software; you can redistribute it and/or modify #
11 !# it under the terms of the GNU General Public License as published by #
12 !# the Free Software Foundation; either version 2 of the License, or (at #
13 !# your option) any later version. #
14 !# #
15 !# This program is distributed in the hope that it will be useful, but #
16 !# WITHOUT ANY WARRANTY; without even the implied warranty of #
17 !# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
18 !# NON INFRINGEMENT. See the GNU General Public License for more #
19 !# details. #
20 !# #
21 !# You should have received a copy of the GNU General Public License #
22 !# along with this program; if not, write to the Free Software #
23 !# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
24 !# #
25 !#############################################################################
26 
27 !----------------------------------------------------------------------------!
35 !----------------------------------------------------------------------------!
36 MODULE mesh_midpoint_mod
37  USE logging_base_mod
38  USE mesh_base_mod
39  USE common_dict
40  IMPLICIT NONE
41  !--------------------------------------------------------------------------!
42  TYPE, EXTENDS(mesh_base) :: mesh_midpoint
43  PRIVATE
44  CONTAINS
45  PRIVATE
46  PROCEDURE, PUBLIC :: vectordivergence2d_1
47 ! PROCEDURE :: VectorDivergence2D_2
48  PROCEDURE, PUBLIC :: tensordivergence2d_1
49 ! PROCEDURE :: TensorDivergence2D_2
50  PROCEDURE, PUBLIC:: tensordivergence3d
51  PROCEDURE, PUBLIC:: vectordivergence3d
52  PROCEDURE, PUBLIC :: initmesh_midpoint
53  PROCEDURE, PUBLIC :: finalize
54  END TYPE mesh_midpoint
55  ! exclude interface block from doxygen processing
58  !--------------------------------------------------------------------------!
59  PRIVATE
60  CHARACTER(LEN=32), PARAMETER :: mesh_name = "midpoint"
61  !--------------------------------------------------------------------------!
62  PUBLIC :: &
63  ! constants
64 #ifdef PARALLEL
65  default_mpi_real, &
66 #endif
67  ! types
68  mesh_midpoint
69  !--------------------------------------------------------------------------!
70 
71 CONTAINS
72 
73  SUBROUTINE initmesh_midpoint(this,config,IO)
74  IMPLICIT NONE
75  !------------------------------------------------------------------------!
76  CLASS(mesh_midpoint),INTENT(INOUT) :: this
77  TYPE(Dict_TYP),POINTER :: config,IO
78  !------------------------------------------------------------------------!
79  INTEGER :: err
80  !------------------------------------------------------------------------!
81 
82  ! basic mesh and geometry initialization
83  CALL this%InitMesh(config, io, midpoint, mesh_name)
84 
85  ! allocate memory for pointers that are specific for midpoint fluxes
86  ALLOCATE( &
87  this%dAx(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
88  this%dAy(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
89  this%dAz(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
90  this%dAxdydz(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
91  this%dAydzdx(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
92  this%dAzdxdy(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX,2), &
93  stat=err)
94  IF (err.NE.0) THEN
95  CALL this%Error("InitMesh_midpoint", "Unable to allocate memory.")
96  END IF
97 
98  ! surface elements divided by dxdy or dydz or dzdx
99  ! perpendicular to x-direction
100  this%dAxdydz(:,:,:,1:2) = this%hy%faces(:,:,:,1:2)*this%hz%faces(:,:,:,1:2)
101  ! perpendicular to y-direction
102  this%dAydzdx(:,:,:,1:2) = this%hz%faces(:,:,:,3:4)*this%hx%faces(:,:,:,3:4)
103  ! perpendicular to z-direction
104  this%dAzdxdy(:,:,:,1:2) = this%hx%faces(:,:,:,5:6)*this%hy%faces(:,:,:,5:6)
105 
106  ! surface elements
107  this%dAx(:,:,:,:) = this%dAxdydz(:,:,:,:)*this%dy*this%dz ! perpendicular to x-direction
108  this%dAy(:,:,:,:) = this%dAydzdx(:,:,:,:)*this%dz*this%dx ! perpendicular to y-direction
109  this%dAz(:,:,:,:) = this%dAzdxdy(:,:,:,:)*this%dx*this%dy ! perpendicular to z-direction
110 
111  ! volume elements = hx*hy*hz*dx*dy*dz
112  this%volume%data3d(:,:,:) = this%sqrtg%center(:,:,:)*this%dx*this%dy*this%dz
113 
114  ! inverse volume elements multiplied by dxdy or dydz or dzdx
115  this%dxdydV%data1d(:) = this%dx*this%dy/(this%volume%data1d(:)+tiny(this%dx))
116  this%dydzdV%data1d(:) = this%dy*this%dz/(this%volume%data1d(:)+tiny(this%dy))
117  this%dzdxdV%data1d(:) = this%dz*this%dx/(this%volume%data1d(:)+tiny(this%dz))
118 
119  ! commutator coefficients (geometric center values)
120  this%cyxy%center(:,:,:) = 0.5*(this%hz%faces(:,:,:,2)+this%hz%faces(:,:,:,1)) &
121  * (this%hy%faces(:,:,:,2)-this%hy%faces(:,:,:,1)) * this%dydzdV%data3d(:,:,:)
122  this%cyzy%center(:,:,:) = 0.5*(this%hx%faces(:,:,:,6)+this%hx%faces(:,:,:,5)) &
123  * (this%hy%faces(:,:,:,6)-this%hy%faces(:,:,:,5)) * this%dxdydV%data3d(:,:,:)
124  this%cxyx%center(:,:,:) = 0.5*(this%hz%faces(:,:,:,4)+this%hz%faces(:,:,:,3)) &
125  * (this%hx%faces(:,:,:,4)-this%hx%faces(:,:,:,3)) * this%dzdxdV%data3d(:,:,:)
126  this%cxzx%center(:,:,:) = 0.5*(this%hy%faces(:,:,:,6)+this%hy%faces(:,:,:,5)) &
127  * (this%hx%faces(:,:,:,6)-this%hx%faces(:,:,:,5)) * this%dxdydV%data3d(:,:,:)
128  this%czxz%center(:,:,:) = 0.5*(this%hy%faces(:,:,:,2)+this%hy%faces(:,:,:,1)) &
129  * (this%hz%faces(:,:,:,2)-this%hz%faces(:,:,:,1)) * this%dydzdV%data3d(:,:,:)
130  this%czyz%center(:,:,:) = 0.5*(this%hx%faces(:,:,:,4)+this%hx%faces(:,:,:,3)) &
131  * (this%hz%faces(:,:,:,4)-this%hz%faces(:,:,:,3)) * this%dzdxdV%data3d(:,:,:)
132 
133  ! set bary center values to geometric center values
134  this%cyxy%bcenter(:,:,:) = this%cyxy%center(:,:,:)
135  this%cyzy%bcenter(:,:,:) = this%cyzy%center(:,:,:)
136  this%cxyx%bcenter(:,:,:) = this%cxyx%center(:,:,:)
137  this%cxzx%bcenter(:,:,:) = this%cxzx%center(:,:,:)
138  this%czxz%bcenter(:,:,:) = this%czxz%center(:,:,:)
139  this%czyz%bcenter(:,:,:) = this%czyz%center(:,:,:)
140 
141  ! center line elements
142  this%dlx%data3d(:,:,:) = this%hx%center(:,:,:)*this%dx
143  this%dly%data3d(:,:,:) = this%hy%center(:,:,:)*this%dy
144  this%dlz%data3d(:,:,:) = this%hz%center(:,:,:)*this%dz
145  END SUBROUTINE initmesh_midpoint
146 
147 
155  PURE SUBROUTINE vectordivergence2d_1(this,vx,vy,divv)
156  IMPLICIT NONE
157  !------------------------------------------------------------------------!
158  CLASS(mesh_midpoint),INTENT(IN) :: this
159  REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
160  :: vx,vy,divv
161  !------------------------------------------------------------------------!
162  INTEGER :: i,j,k
163  !------------------------------------------------------------------------!
164  INTENT(IN) :: vx,vy
165  INTENT(OUT) :: divv
166  !------------------------------------------------------------------------!
167  ! determine which velocity components are given
168  SELECT CASE(this%VECTOR_COMPONENTS)
169  CASE(ior(vector_x,vector_y))
170  DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
171  DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
172 !NEC$ IVDEP
173  DO i=this%IGMIN+this%IP1,this%IGMAX-this%IP1
174  divv(i,j,k) = divergence3d(&
175  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
176  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
177  0.0,0.0, &
178  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),0.0, &
179  0.0,0.0,0.0,0.0, & ! no commutator coefficients
180  0.5*(vx(i-this%IP1,j,k)+vx(i,j,k)),0.5*(vx(i+this%IP1,j,k)+vx(i,j,k)), &
181  0.5*(vy(i,j-this%JP1,k)+vy(i,j,k)),0.5*(vy(i,j+this%JP1,k)+vy(i,j,k)), &
182  0.0, 0.0, & ! no vz component
183  0.0,0.0,0.0,0.0) ! no off-diagonal tensor components
184  END DO
185  END DO
186  END DO
187  CASE(ior(vector_x,vector_z))
188  DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
189  DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
190 !NEC$ IVDEP
191  DO i=this%IGMIN+this%IP1,this%IGMAX-this%IP1
192  ! we simply use the 3D tensor divergence and set the commutator coefficients
193  ! and the off-diagonal tensor components to 0
194  divv(i,j,k) = divergence3d(&
195  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
196  0.0,0.0, &
197  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
198  this%dydzdV%data3d(i,j,k),0.0,this%dxdydV%data3d(i,j,k), &
199  0.0,0.0,0.0,0.0, & ! no commutator coefficients
200  0.5*(vx(i-this%IP1,j,k)+vx(i,j,k)),0.5*(vx(i+this%IP1,j,k)+vx(i,j,k)), &
201  0.0, 0.0, & ! no vy component
202  0.5*(vy(i,j,k-this%KP1)+vy(i,j,k)),0.5*(vy(i,j,k)+vy(i,j,k+this%KP1)), &
203  0.0,0.0,0.0,0.0) ! no off-diagonal tensor components
204  END DO
205  END DO
206  END DO
207  CASE(ior(vector_y,vector_z))
208  DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
209  DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
210 !NEC$ IVDEP
211  DO i=this%IGMIN+this%IP1,this%IGMAX-this%IP1
212  ! we simply use the 3D tensor divergence and set the commutator coefficients
213  ! and the off-diagonal tensor components to 0
214  divv(i,j,k) = divergence3d(&
215  0.0,0.0, &
216  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
217  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
218  0.0,this%dzdxdV%data3d(i,j,k),this%dxdydV%data3d(i,j,k), &
219  0.0,0.0,0.0,0.0, & ! no commutator coefficients
220  0.0, 0.0, & ! no vx component
221  0.5*(vx(i,j-this%JP1,k)+vx(i,j,k)),0.5*(vx(i,j+this%JP1,k)+vx(i,j,k)), &
222  0.5*(vy(i,j,k-this%KP1)+vy(i,j,k)),0.5*(vy(i,j,k)+vy(i,j,k+this%KP1)), &
223  0.0,0.0,0.0,0.0) ! no off-diagonal tensor components
224  END DO
225  END DO
226  END DO
227  CASE DEFAULT
228  ! return NaN
229  divv(:,:,:) = nan_default_real
230  END SELECT
231  END SUBROUTINE vectordivergence2d_1
232 
233 
234 ! ! computes the cell centered curvilinear vector divergence
235 ! ! for 2D vector v given on the 4 face centered positions
236 ! PURE SUBROUTINE VectorDivergence2D_2(this,v,divv)
237 ! IMPLICIT NONE
238 ! !------------------------------------------------------------------------!
239 ! CLASS(mesh_midpoint),INTENT(IN) :: this
240 ! REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,4,2) &
241 ! :: v
242 ! REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX) &
243 ! :: divv
244 ! !------------------------------------------------------------------------!
245 ! INTEGER :: i,j
246 ! !------------------------------------------------------------------------!
247 ! INTENT(IN) :: v
248 ! INTENT(OUT) :: divv
249 ! !------------------------------------------------------------------------!
250 ! ! we simply use the 3D tensor divergence and set all commutator coefficients
251 ! ! and the tensor components Tyx, Tyy, Tzz to zero
252 !!NEC$ IVDEP
253 ! DO j=this%JGMIN,this%JGMAX
254 !!NEC$ IVDEP
255 ! DO i=this%IGMIN,this%IGMAX
256 ! divv(i,j) = Divergence3D(this%dAxdy(i,j,1),this%dAxdy(i,j,2), &
257 ! this%dAydx(i,j,1),this%dAydx(i,j,2), &
258 ! this%dxdV(i,j),this%dydV(i,j), &
259 ! 0.0,0.0,0.0, & ! vanishing commutator coefficients
260 ! v(i,j,WEST,1),v(i,j,EAST,1), &
261 ! v(i,j,SOUTH,2),v(i,j,NORTH,2), &
262 ! 0.0,0.0,0.0) ! vanishing tensor components
263 ! END DO
264 ! END DO
265 ! END SUBROUTINE VectorDivergence2D_2
266 !
267 
278  PURE SUBROUTINE tensordivergence2d_1(this,Txx,Txy,Tyx,Tyy,divTx,divTy)
279  IMPLICIT NONE
280  !------------------------------------------------------------------------!
281  CLASS(mesh_midpoint),INTENT(IN) :: this
282  REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
283  :: txx,txy,tyx,tyy,divtx,divty
284  !------------------------------------------------------------------------!
285  INTEGER :: i,j,k
286  !------------------------------------------------------------------------!
287  INTENT(IN) :: txx,txy,tyx,tyy
288  INTENT(OUT) :: divtx,divty
289  !------------------------------------------------------------------------!
290  ! determine which velocity components are given
291  SELECT CASE(this%VECTOR_COMPONENTS)
292  CASE(ior(vector_x,vector_y))
293  DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
294  DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
295 !NEC$ IVDEP
296  DO i=this%IGMIN+this%IP1,this%IGMAX-this%IP1
297  ! x component of tensor divergence
298  divtx(i,j,k) = divergence3d(&
299  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
300  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
301  0.0,0.0,&
302  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),0.0, &
303  this%cxyx%center(i,j,k),0.0, &
304  this%cyxy%center(i,j,k),0.0, &
305  0.5*(txx(i-this%IP1,j,k)+txx(i,j,k)),0.5*(txx(i+this%IP1,j,k)+txx(i,j,k)), &
306  0.5*(txy(i,j-this%JP1,k)+txy(i,j,k)),0.5*(txy(i,j+this%JP1,k)+txy(i,j,k)), &
307  0.0,0.0, &
308  tyx(i,j,k),0.0,tyy(i,j,k),0.0)
309  ! y component of tensor divergence
310  ! change input of Divergence3D according to the rules given in
311  ! the comments (see below)
312  divty(i,j,k) = divergence3d(&
313  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
314  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
315  0.0,0.0, &
316  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),0.0, &
317  -this%cxyx%center(i,j,k),0.0, &
318  -this%cyxy%center(i,j,k),0.0, &
319  0.5*(tyx(i-this%IP1,j,k)+tyx(i,j,k)),0.5*(tyx(i+this%IP1,j,k)+tyx(i,j,k)), &
320  0.5*(tyy(i,j-this%JP1,k)+tyy(i,j,k)),0.5*(tyy(i,j+this%JP1,k)+tyy(i,j,k)), &
321  0.0,0.0,&
322  txx(i,j,k),0.0,txy(i,j,k),0.0)
323  END DO
324  END DO
325  END DO
326  CASE(ior(vector_x,vector_z))
327  ! renaming scheme
328  ! Txz => Txy
329  ! Tzx => Tyx
330  ! Tzz => Tyy
331  ! divTz => divTy
332  DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
333  DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
334  !NEC$ IVDEP
335  DO i=this%IGMIN+this%IP1,this%IGMAX-this%IP1
336  ! x component of tensor divergence
337  divtx(i,j,k) = divergence3d(&
338  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
339  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
340  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
341  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),this%dxdydV%data3d(i,j,k), &
342  this%cxyx%center(i,j,k),this%cxzx%center(i,j,k), &
343  this%cyxy%center(i,j,k),this%czxz%center(i,j,k), &
344  0.5*(txx(i-this%IP1,j,k)+txx(i,j,k)),0.5*(txx(i+this%IP1,j,k)+txx(i,j,k)), &
345  0.0,0.0, & ! Txy = 0
346  0.5*(txy(i,j,k-this%KP1)+txy(i,j,k)),0.5*(txy(i,j,k+this%KP1)+txy(i,j,k)), & ! Txz => Txy
347  0.0,tyx(i,j,k),0.0,tyy(i,j,k)) ! Tyx = 0, Tzx => Tyx, Tyy = 0, Tzz => Tyy
348  ! z component of tensor divergence
349  ! change input of Divergence3D according to the rules given in
350  ! the comments (see below)
351  divty(i,j,k) = divergence3d(& ! divTz => divTy
352  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
353  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
354  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
355  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),this%dxdydV%data3d(i,j,k), &
356  this%czyz%center(i,j,k),-this%cxzx%center(i,j,k), &
357  this%cyzy%center(i,j,k),-this%czxz%center(i,j,k), &
358  0.5*(tyx(i-this%IP1,j,k)+tyx(i,j,k)),0.5*(tyx(i+this%IP1,j,k)+tyx(i,j,k)), & ! Tzx => Tyx
359  0.0,0.0, & ! Tzy = 0
360  0.5*(tyy(i,j,k-this%KP1)+tyy(i,j,k)),0.5*(tyy(i,j,k+this%KP1)+tyy(i,j,k)), & ! Tzz => Tyy
361  0.0,txx(i,j,k),0.0,txy(i,j,k)) ! Tyz = 0, Txz => Txy, Tyy = 0
362  END DO
363  END DO
364  END DO
365  CASE(ior(vector_y,vector_z))
366  ! renaming scheme
367  ! Tyz => Txy
368  ! Tyy => Txx
369  ! Tzy => Tyx
370  ! Tzz => Tyy
371  ! divTy => divTx
372  ! divTz => divTy
373  DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
374  DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
375  !NEC$ IVDEP
376  DO i=this%IGMIN+this%IP1,this%IGMAX-this%IP1
377  ! y component of tensor divergence
378  ! change input of Divergence3D according to the rules given in
379  ! the comments (see below)
380  divtx(i,j,k) = divergence3d(& ! divTy => divTx
381  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
382  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
383  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
384  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),this%dxdydV%data3d(i,j,k), &
385  -this%cxyx%center(i,j,k),this%cyzy%center(i,j,k), &
386  -this%cyxy%center(i,j,k),this%czyz%center(i,j,k), &
387  0.0,0.0, & ! Tyx = 0
388  0.5*(txx(i,j-this%JP1,k)+txx(i,j,k)),0.5*(txx(i,j+this%JP1,k)+tyy(i,j,k)), & ! Tyy => Txx
389  0.5*(txy(i,j,k-this%KP1)+txy(i,j,k)),0.5*(txy(i,j,k+this%KP1)+txy(i,j,k)), & ! Tyz => Txy
390  0.0,tyx(i,j,k),0.0,tyy(i,j,k)) ! Txx = 0, Tzy => Tyx, Txy = 0, Tzz => Tyy
391  ! z component of tensor divergence
392  ! change input of Divergence3D according to the rules given in
393  ! the comments (see below)
394  divty(i,j,k) = divergence3d(& ! divTz => divTy
395  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
396  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
397  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
398  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),this%dxdydV%data3d(i,j,k), &
399  this%czyz%center(i,j,k),-this%cxzx%center(i,j,k), &
400  this%cyzy%center(i,j,k),-this%czxz%center(i,j,k), &
401  0.0,0.0, & ! Tzx = 0
402  0.5*(tyx(i,j-this%JP1,k)+tyx(i,j,k)),0.5*(tyx(i,j+this%JP1,k)+tyx(i,j,k)), & ! Tzy => Tyx
403  0.5*(tyy(i,j,k-this%KP1)+tyy(i,j,k)),0.5*(tyy(i,j,k+this%KP1)+tyy(i,j,k)), & ! Tzz => Tyy
404  txy(i,j,k),0.0,txx(i,j,k),0.0) ! Tyz => Txy, Txx = 0, Tyy => Txx, Txz = 0
405  END DO
406  END DO
407  END DO
408  CASE DEFAULT
409  ! return NaN
410  divtx(:,:,:) = nan_default_real
411  divty(:,:,:) = nan_default_real
412  END SELECT
413  END SUBROUTINE tensordivergence2d_1
414 
415 
416 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
417 !! ATTENTION: TensorDivergence2D_2 is untested, use with care
418 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
419 ! ! computes the cell centered curvilinear tensor divergence
420 ! ! input: 2D rank 2 tensor T given on the 4 face centered positions
421 ! ! output: 2D vector divT
422 ! PURE SUBROUTINE TensorDivergence2D_2(this,T,divT)
423 ! IMPLICIT NONE
424 ! !------------------------------------------------------------------------!
425 ! CLASS(mesh_midpoint),INTENT(IN) :: this
426 ! REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,4,2,2) :: T
427 ! REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,2) :: divT
428 ! !------------------------------------------------------------------------!
429 ! REAL :: Txx,Txy,Tyx,Tyy
430 ! INTEGER :: i,j
431 ! !------------------------------------------------------------------------!
432 ! INTENT(IN) :: T
433 ! INTENT(OUT) :: divT
434 ! !------------------------------------------------------------------------!
435 ! ! we simply use the 3D tensor divergence and set the commutator coefficient
436 ! ! and the tensor component related to the z-direction to zero
437 !!NEC$ IVDEP
438 ! DO j=this%JGMIN,this%JGMAX
439 !!NEC$ IVDEP
440 ! DO i=this%IGMIN,this%IGMAX
441 ! ! compute mean value of all four face values to obtain cell centered values
442 ! Txx = 0.25*SUM(T(i,j,:,1,1))
443 ! Txy = 0.25*SUM(T(i,j,:,1,2))
444 ! Tyx = 0.25*SUM(T(i,j,:,2,1))
445 ! Tyy = 0.25*SUM(T(i,j,:,2,2))
446 ! ! x component of tensor divergence
447 ! divT(i,j,1) = Divergence3D(this%dAxdy(i,j,1),this%dAxdy(i,j,2), &
448 ! this%dAydx(i,j,1),this%dAydx(i,j,2), &
449 ! this%dxdV(i,j),this%dydV(i,j), &
450 ! this%cxyx%center(i,j),this%cyxy%center(i,j), &
451 ! 0.0, & ! czxz = 0 because of 2D
452 ! T(i,j,WEST,1,1),T(i,j,EAST,1,1), &
453 ! T(i,j,SOUTH,1,2),T(i,j,NORTH,1,2), &
454 ! Tyx,Tyy,0.0) ! Tzz = 0 because of 2D
455 ! ! y component of tensor divergence
456 ! divT(i,j,2) = Divergence3D(this%dAxdy(i,j,1),this%dAxdy(i,j,2), &
457 ! this%dAydx(i,j,1),this%dAydx(i,j,2), &
458 ! this%dxdV(i,j),this%dydV(i,j), &
459 ! -this%cxyx%center(i,j),-this%cyxy%center(i,j), &
460 ! 0.0, & ! czyz = 0 because of 2D
461 ! T(i,j,WEST,2,1),T(i,j,EAST,2,1), &
462 ! T(i,j,SOUTH,2,2),T(i,j,NORTH,2,2), &
463 ! Txx,Txy,0.0) ! Tzz = 0 because of 2D
464 ! END DO
465 ! END DO
466 ! END SUBROUTINE TensorDivergence2D_2
467 
468 
473  PURE SUBROUTINE vectordivergence3d(this,vx,vy,vz,divv)
474  IMPLICIT NONE
475  !------------------------------------------------------------------------!
476  CLASS(mesh_midpoint),INTENT(IN) :: this
477  REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
478  :: vx,vy,vz,divv
479  !------------------------------------------------------------------------!
480  INTEGER :: i,j,k
481  !------------------------------------------------------------------------!
482  INTENT(IN) :: vx,vy,vz
483  INTENT(OUT) :: divv
484  !------------------------------------------------------------------------!
485  ! we simply use the 3D tensor divergence and set the commutator coefficients
486  ! and the off-diagonal tensor components to 0
487  DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
488  DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
489 !NEC$ IVDEP
490  DO i=this%IGMIN+this%IP1,this%IGMAX-this%IP1
491  divv(i,j,k) = divergence3d(&
492  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
493  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
494  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
495  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),this%dxdydV%data3d(i,j,k), &
496  0.0,0.0,0.0,0.0, & ! no commutator coefficients
497  0.5*(vx(i-this%IP1,j,k)+vx(i,j,k)),0.5*(vx(i+this%IP1,j,k)+vx(i,j,k)), &
498  0.5*(vy(i,j-this%JP1,k)+vy(i,j,k)),0.5*(vy(i,j+this%JP1,k)+vy(i,j,k)), &
499  0.5*(vz(i,j,k-this%KP1)+vz(i,j,k)),0.5*(vz(i,j,k+this%KP1)+vz(i,j,k)), &
500  0.0,0.0,0.0,0.0) ! no off-diagonal tensor components
501  END DO
502  END DO
503  END DO
504  END SUBROUTINE vectordivergence3d
505 
506 
507 
516  PURE SUBROUTINE tensordivergence3d(this,Txx,Txy,Txz,Tyx,Tyy,Tyz,Tzx,Tzy,Tzz, &
517  divTx,divTy,divTz)
518  IMPLICIT NONE
519  !------------------------------------------------------------------------!
520  CLASS(mesh_midpoint), INTENT(IN) :: this
521  REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
522  :: txx,txy,txz,tyx,tyy,tyz,tzx,tzy,tzz,divtx,divty,divtz
523  !------------------------------------------------------------------------!
524  INTEGER :: i,j,k
525  !------------------------------------------------------------------------!
526  INTENT(IN) :: txx,txy,txz,tyx,tyy,tyz,tzx,tzy,tzz
527  INTENT(OUT) :: divtx,divty,divtz
528  !------------------------------------------------------------------------!
529  DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
530  DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
531 !NEC$ IVDEP
532  DO i=this%IGMIN+this%IP1,this%IGMAX-this%IP1
533  ! x component of tensor divergence
534  divtx(i,j,k) = divergence3d(&
535  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
536  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
537  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
538  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),this%dxdydV%data3d(i,j,k), &
539  this%cxyx%center(i,j,k),this%cxzx%center(i,j,k), &
540  this%cyxy%center(i,j,k),this%czxz%center(i,j,k), &
541  0.5*(txx(i-this%IP1,j,k)+txx(i,j,k)),0.5*(txx(i+this%IP1,j,k)+txx(i,j,k)), &
542  0.5*(txy(i,j-this%JP1,k)+txy(i,j,k)),0.5*(txy(i,j+this%JP1,k)+txy(i,j,k)), &
543  0.5*(txz(i,j,k-this%KP1)+txz(i,j,k)),0.5*(txz(i,j,k+this%KP1)+txz(i,j,k)), &
544  tyx(i,j,k),tzx(i,j,k),tyy(i,j,k),tzz(i,j,k))
545  ! y component of tensor divergence
546  ! change input of Divergence3D according to the rules given in
547  ! the comments (see below)
548  divty(i,j,k) = divergence3d(&
549  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
550  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
551  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
552  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),this%dxdydV%data3d(i,j,k), &
553  -this%cxyx%center(i,j,k),this%cyzy%center(i,j,k), &
554  -this%cyxy%center(i,j,k),this%czyz%center(i,j,k), &
555  0.5*(tyx(i-this%IP1,j,k)+tyx(i,j,k)),0.5*(tyx(i+this%IP1,j,k)+tyx(i,j,k)), &
556  0.5*(tyy(i,j-this%JP1,k)+tyy(i,j,k)),0.5*(tyy(i,j+this%JP1,k)+tyy(i,j,k)), &
557  0.5*(tyz(i,j,k-this%KP1)+tyz(i,j,k)),0.5*(tyz(i,j,k+this%KP1)+tyz(i,j,k)), &
558  txx(i,j,k),tzy(i,j,k),txy(i,j,k),tzz(i,j,k))
559  ! z component of tensor divergence
560  ! change input of Divergence3D according to the rules given in
561  ! the comments (see below)
562  divtz(i,j,k) = divergence3d(&
563  this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
564  this%dAydzdx(i,j,k,1),this%dAydzdx(i,j,k,2), &
565  this%dAzdxdy(i,j,k,1),this%dAzdxdy(i,j,k,2), &
566  this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),this%dxdydV%data3d(i,j,k), &
567  this%czyz%center(i,j,k),-this%cxzx%center(i,j,k), &
568  this%cyzy%center(i,j,k),-this%czxz%center(i,j,k), &
569  0.5*(tzx(i-this%IP1,j,k)+tzx(i,j,k)),0.5*(tzx(i+this%IP1,j,k)+tzx(i,j,k)), &
570  0.5*(tzy(i,j-this%JP1,k)+tzy(i,j,k)),0.5*(tzy(i,j+this%JP1,k)+tzy(i,j,k)), &
571  0.5*(tzz(i,j,k-this%KP1)+tzz(i,j,k)),0.5*(tzz(i,j,k+this%KP1)+tzz(i,j,k)), &
572  tyz(i,j,k),txx(i,j,k),tyy(i,j,k),txz(i,j,k))
573  END DO
574  END DO
575  END DO
576  END SUBROUTINE tensordivergence3d
577 
578 
585  ELEMENTAL FUNCTION divergence3d(dAxdydzWest,dAxdydzEast,dAydxdzSouth,dAydxdzNorth,dAzdxdyBottom,dAzdxdyTop, &
586  dydzdV,dxdzdV,dxdydV,cxyx,cxzx,cyxy,czxz, &
587  TxxWest,TxxEast,TxySouth,TxyNorth,TxzBottom,TxzTop, &
588  TyxCent,TzxCent,TyyCent,TzzCent) RESULT(div)
589  IMPLICIT NONE
590  !------------------------------------------------------------------------!
591  REAL,INTENT(IN) :: daxdydzwest,daxdydzeast, & ! surface elements
592  daydxdzsouth,daydxdznorth, & ! dAx/dy, dAy/dx
593  dazdxdybottom,dazdxdytop, &
594  dxdydv,dxdzdv,dydzdv, & ! dx/dV and dy/dV
595  cxyx,cyxy,czxz,cxzx, & ! commutator coeffs
596  txxeast,txxwest, & ! tensor components
597  txynorth,txysouth, & ! on cell faces
598  txzbottom,txztop, &
599  tyxcent,tyycent,tzzcent,tzxcent ! central
600  ! to compute the y-component of the tensor divergence,
601  ! one has to modify the input according to
602  ! Txx <-> Tyx
603  ! Txy <-> Tyy
604  ! Txz -> Tyz
605  ! Tzx -> Tzy
606  ! Tzz -> Tzz
607  ! cxyx -> -cxyx
608  ! cyxy -> -cyxy
609  ! cxzx -> cyzy
610  ! czxz -> czyz
611  ! to compute the z-component of the tensor divergence,
612  ! one has to modify the input according to
613  ! Txx <-> Tzx
614  ! Txy -> Tzy
615  ! Txz <-> Tzz
616  ! Tyx -> Tyz
617  ! Tyy -> Tyy
618  ! cxyx -> czyz
619  ! cyxy -> cyzy
620  ! cxzx -> -cxzx
621  ! czxz -> -czxz
622  !------------------------------------------------------------------------!
623  REAL :: div
624  !------------------------------------------------------------------------!
625  div = dydzdv*(daxdydzeast*txxeast-daxdydzwest*txxwest) &
626  + dxdzdv*(daydxdznorth*txynorth-daydxdzsouth*txysouth) &
627  + dxdydv*(dazdxdytop*txztop-dazdxdybottom*txzbottom) &
628  + cxyx*tyxcent + cxzx*tzxcent - cyxy*tyycent - czxz*tzzcent
629  END FUNCTION divergence3d
630 
631 
632 
633  SUBROUTINE finalize(this)
634  IMPLICIT NONE
635  !------------------------------------------------------------------------!
636  CLASS(mesh_midpoint),INTENT(INOUT) :: this
637  !------------------------------------------------------------------------!
638  DEALLOCATE(this%dAx,this%dAy,this%dAz,this%dAxdydz,this%dAydzdx,this%dAzdxdy)
639  CALL this%Finalize_base()
640  END SUBROUTINE finalize
641 
642 END MODULE mesh_midpoint_mod
mesh_base
logging_base_mod::finalize
subroutine finalize(this)
Destructor of common class.
Definition: logging_base.f90:155
mesh_midpoint_mod
mesh module for midpoint quadrature rule
Definition: mesh_midpoint.f90:36
mesh_midpoint_mod::vectordivergence2d_1
pure subroutine vectordivergence2d_1(this, vx, vy, divv)
compute the cell centered 2D vector divergence
Definition: mesh_midpoint.f90:156
logging_base_mod::default_mpi_real
integer, save default_mpi_real
default real type for MPI
Definition: logging_base.f90:80
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
logging_base_mod::nan_default_real
real, save nan_default_real
NaN real constant.
Definition: logging_base.f90:85
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
mesh_midpoint_mod::mesh_name
character(len=32), parameter mesh_name
Definition: mesh_midpoint.f90:60
mesh_midpoint_mod::tensordivergence2d_1
pure subroutine tensordivergence2d_1(this, Txx, Txy, Tyx, Tyy, divTx, divTy)
compute the cell centered 2D rank 2 tensor divergence
Definition: mesh_midpoint.f90:279
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
mesh_midpoint_mod::initmesh_midpoint
subroutine initmesh_midpoint(this, config, IO)
Definition: mesh_midpoint.f90:74
mesh_midpoint_mod::divergence3d
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
Definition: mesh_midpoint.f90:589
mesh_midpoint_mod::vectordivergence3d
pure subroutine vectordivergence3d(this, vx, vy, vz, divv)
compute the cell centered 3D vector divergence
Definition: mesh_midpoint.f90:474
mesh_midpoint_mod::mesh_midpoint
Definition: mesh_midpoint.f90:42
mesh_midpoint_mod::tensordivergence3d
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
Definition: mesh_midpoint.f90:518
Generated by   doxygen 1.8.14