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!----------------------------------------------------------------------------!
36MODULE 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
71CONTAINS
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 :: i,j,k,err
80 !------------------------------------------------------------------------!
81
82 ! basic mesh and geometry initialization
83 CALL this%InitMesh(config, io, midpoint, mesh_name)
84 ! allocate memory for pointers that are specific for midpoint fluxes
85 ALLOCATE( &
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), &
92 stat=err)
93 IF (err.NE.0) THEN
94 CALL this%Error("InitMesh_midpoint", "Unable to allocate memory.")
95 END IF
96
97 DO concurrent(i=this%IGMIN:this%IGMAX,j=this%JGMIN:this%JGMAX,k=this%KGMIN:this%KGMAX)
98 ! surface elements divided by dxdy or dydz or dzdx
99 ! perpendicular to x-direction
100 this%dAxdydz(i,j,k,1:2) = this%hy%faces(i,j,k,1:2)*this%hz%faces(i,j,k,1:2)
101 ! perpendicular to y-direction
102 this%dAydzdx(i,j,k,1:2) = this%hz%faces(i,j,k,3:4)*this%hx%faces(i,j,k,3:4)
103 ! perpendicular to z-direction
104 this%dAzdxdy(i,j,k,1:2) = this%hx%faces(i,j,k,5:6)*this%hy%faces(i,j,k,5:6)
105
106 ! surface elements
107 this%dAx(i,j,k,:) = this%dAxdydz(i,j,k,:)*this%dy*this%dz ! perpendicular to x-direction
108 this%dAy(i,j,k,:) = this%dAydzdx(i,j,k,:)*this%dz*this%dx ! perpendicular to y-direction
109 this%dAz(i,j,k,:) = this%dAzdxdy(i,j,k,:)*this%dx*this%dy ! perpendicular to z-direction
110
111 ! volume elements = hx*hy*hz*dx*dy*dz
112 this%volume%data3d(i,j,k) = this%sqrtg%center(i,j,k)*this%dx*this%dy*this%dz
113
114 ! inverse volume elements multiplied by dxdy or dydz or dzdx
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))
118
119 ! commutator coefficients (geometric center values)
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)
132
133 ! set bary center values to geometric center values
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)
140
141 ! center line elements
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
145 END DO
146 END SUBROUTINE initmesh_midpoint
147
148
156 PURE SUBROUTINE vectordivergence2d_1(this,vx,vy,divv)
157 IMPLICIT NONE
158 !------------------------------------------------------------------------!
159 CLASS(mesh_midpoint),INTENT(IN) :: this
160 REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
161 :: vx,vy,divv
162 !------------------------------------------------------------------------!
163 INTEGER :: i,j,k
164 !------------------------------------------------------------------------!
165 INTENT(IN) :: vx,vy
166 INTENT(OUT) :: divv
167 !------------------------------------------------------------------------!
168 ! determine which velocity components are given
169 SELECT CASE(this%VECTOR_COMPONENTS)
170 CASE(ior(vector_x,vector_y))
171 DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
172 DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
173!NEC$ IVDEP
174 DO concurrent(i=this%IGMIN+this%IP1:this%IGMAX-this%IP1)
175 divv(i,j,k) = divergence3d(&
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), &
178 0.0,0.0, &
179 this%dydzdV%data3d(i,j,k),this%dzdxdV%data3d(i,j,k),0.0, &
180 0.0,0.0,0.0,0.0, & ! no commutator coefficients
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)), &
183 0.0, 0.0, & ! no vz component
184 0.0,0.0,0.0,0.0) ! no off-diagonal tensor components
185 END DO
186 END DO
187 END DO
188 CASE(ior(vector_x,vector_z))
189 DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
190 DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
191!NEC$ IVDEP
192 DO concurrent(i=this%IGMIN+this%IP1:this%IGMAX-this%IP1)
193 ! we simply use the 3D tensor divergence and set the commutator coefficients
194 ! and the off-diagonal tensor components to 0
195 divv(i,j,k) = divergence3d(&
196 this%dAxdydz(i,j,k,1),this%dAxdydz(i,j,k,2), &
197 0.0,0.0, &
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), &
200 0.0,0.0,0.0,0.0, & ! no commutator coefficients
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)), &
202 0.0, 0.0, & ! no vy component
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)), &
204 0.0,0.0,0.0,0.0) ! no off-diagonal tensor components
205 END DO
206 END DO
207 END DO
208 CASE(ior(vector_y,vector_z))
209 DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
210 DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
211!NEC$ IVDEP
212 DO concurrent(i=this%IGMIN+this%IP1:this%IGMAX-this%IP1)
213 ! we simply use the 3D tensor divergence and set the commutator coefficients
214 ! and the off-diagonal tensor components to 0
215 divv(i,j,k) = divergence3d(&
216 0.0,0.0, &
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), &
220 0.0,0.0,0.0,0.0, & ! no commutator coefficients
221 0.0, 0.0, & ! no vx component
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)), &
224 0.0,0.0,0.0,0.0) ! no off-diagonal tensor components
225 END DO
226 END DO
227 END DO
228 CASE DEFAULT
229 ! return NaN
230 divv(:,:,:) = nan_default_real
231 END SELECT
232 END SUBROUTINE vectordivergence2d_1
233
234
235! ! computes the cell centered curvilinear vector divergence
236! ! for 2D vector v given on the 4 face centered positions
237! PURE SUBROUTINE VectorDivergence2D_2(this,v,divv)
238! IMPLICIT NONE
239! !------------------------------------------------------------------------!
240! CLASS(mesh_midpoint),INTENT(IN) :: this
241! REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,4,2) &
242! :: v
243! REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX) &
244! :: divv
245! !------------------------------------------------------------------------!
246! INTEGER :: i,j
247! !------------------------------------------------------------------------!
248! INTENT(IN) :: v
249! INTENT(OUT) :: divv
250! !------------------------------------------------------------------------!
251! ! we simply use the 3D tensor divergence and set all commutator coefficients
252! ! and the tensor components Tyx, Tyy, Tzz to zero
253!!NEC$ IVDEP
254! DO j=this%JGMIN,this%JGMAX
255!!NEC$ IVDEP
256! DO i=this%IGMIN,this%IGMAX
257! divv(i,j) = Divergence3D(this%dAxdy(i,j,1),this%dAxdy(i,j,2), &
258! this%dAydx(i,j,1),this%dAydx(i,j,2), &
259! this%dxdV(i,j),this%dydV(i,j), &
260! 0.0,0.0,0.0, & ! vanishing commutator coefficients
261! v(i,j,WEST,1),v(i,j,EAST,1), &
262! v(i,j,SOUTH,2),v(i,j,NORTH,2), &
263! 0.0,0.0,0.0) ! vanishing tensor components
264! END DO
265! END DO
266! END SUBROUTINE VectorDivergence2D_2
267!
268
279 PURE SUBROUTINE tensordivergence2d_1(this,Txx,Txy,Tyx,Tyy,divTx,divTy)
280 IMPLICIT NONE
281 !------------------------------------------------------------------------!
282 CLASS(mesh_midpoint),INTENT(IN) :: this
283 REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
284 :: txx,txy,tyx,tyy,divtx,divty
285 !------------------------------------------------------------------------!
286 INTEGER :: i,j,k
287 !------------------------------------------------------------------------!
288 INTENT(IN) :: txx,txy,tyx,tyy
289 INTENT(OUT) :: divtx,divty
290 !------------------------------------------------------------------------!
291 ! determine which velocity components are given
292 SELECT CASE(this%VECTOR_COMPONENTS)
293 CASE(ior(vector_x,vector_y))
294 DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
295 DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
296!NEC$ IVDEP
297 DO concurrent(i=this%IGMIN+this%IP1:this%IGMAX-this%IP1)
298 ! x component of tensor divergence
299 divtx(i,j,k) = divergence3d(&
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), &
302 0.0,0.0,&
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)), &
308 0.0,0.0, &
309 tyx(i,j,k),0.0,tyy(i,j,k),0.0)
310 ! y component of tensor divergence
311 ! change input of Divergence3D according to the rules given in
312 ! the comments (see below)
313 divty(i,j,k) = divergence3d(&
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), &
316 0.0,0.0, &
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)), &
322 0.0,0.0,&
323 txx(i,j,k),0.0,txy(i,j,k),0.0)
324 END DO
325 END DO
326 END DO
327 CASE(ior(vector_x,vector_z))
328 ! renaming scheme
329 ! Txz => Txy
330 ! Tzx => Tyx
331 ! Tzz => Tyy
332 ! divTz => divTy
333 DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
334 DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
335 !NEC$ IVDEP
336 DO concurrent(i=this%IGMIN+this%IP1:this%IGMAX-this%IP1)
337 ! x component of tensor divergence
338 divtx(i,j,k) = divergence3d(&
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)), &
346 0.0,0.0, & ! Txy = 0
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)), & ! Txz => Txy
348 0.0,tyx(i,j,k),0.0,tyy(i,j,k)) ! Tyx = 0, Tzx => Tyx, Tyy = 0, Tzz => Tyy
349 ! z component of tensor divergence
350 ! change input of Divergence3D according to the rules given in
351 ! the comments (see below)
352 divty(i,j,k) = divergence3d(& ! divTz => divTy
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)), & ! Tzx => Tyx
360 0.0,0.0, & ! Tzy = 0
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)), & ! Tzz => Tyy
362 0.0,txx(i,j,k),0.0,txy(i,j,k)) ! Tyz = 0, Txz => Txy, Tyy = 0
363 END DO
364 END DO
365 END DO
366 CASE(ior(vector_y,vector_z))
367 ! renaming scheme
368 ! Tyz => Txy
369 ! Tyy => Txx
370 ! Tzy => Tyx
371 ! Tzz => Tyy
372 ! divTy => divTx
373 ! divTz => divTy
374 DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
375 DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
376 !NEC$ IVDEP
377 DO concurrent(i=this%IGMIN+this%IP1:this%IGMAX-this%IP1)
378 ! y component of tensor divergence
379 ! change input of Divergence3D according to the rules given in
380 ! the comments (see below)
381 divtx(i,j,k) = divergence3d(& ! divTy => divTx
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), &
388 0.0,0.0, & ! Tyx = 0
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)), & ! Tyy => Txx
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)), & ! Tyz => Txy
391 0.0,tyx(i,j,k),0.0,tyy(i,j,k)) ! Txx = 0, Tzy => Tyx, Txy = 0, Tzz => Tyy
392 ! z component of tensor divergence
393 ! change input of Divergence3D according to the rules given in
394 ! the comments (see below)
395 divty(i,j,k) = divergence3d(& ! divTz => divTy
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), &
402 0.0,0.0, & ! Tzx = 0
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)), & ! Tzy => Tyx
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)), & ! Tzz => Tyy
405 txy(i,j,k),0.0,txx(i,j,k),0.0) ! Tyz => Txy, Txx = 0, Tyy => Txx, Txz = 0
406 END DO
407 END DO
408 END DO
409 CASE DEFAULT
410 ! return NaN
411 divtx(:,:,:) = nan_default_real
412 divty(:,:,:) = nan_default_real
413 END SELECT
414 END SUBROUTINE tensordivergence2d_1
415
416
417!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
418!! ATTENTION: TensorDivergence2D_2 is untested, use with care
419!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
420! ! computes the cell centered curvilinear tensor divergence
421! ! input: 2D rank 2 tensor T given on the 4 face centered positions
422! ! output: 2D vector divT
423! PURE SUBROUTINE TensorDivergence2D_2(this,T,divT)
424! IMPLICIT NONE
425! !------------------------------------------------------------------------!
426! CLASS(mesh_midpoint),INTENT(IN) :: this
427! REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,4,2,2) :: T
428! REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,2) :: divT
429! !------------------------------------------------------------------------!
430! REAL :: Txx,Txy,Tyx,Tyy
431! INTEGER :: i,j
432! !------------------------------------------------------------------------!
433! INTENT(IN) :: T
434! INTENT(OUT) :: divT
435! !------------------------------------------------------------------------!
436! ! we simply use the 3D tensor divergence and set the commutator coefficient
437! ! and the tensor component related to the z-direction to zero
438!!NEC$ IVDEP
439! DO j=this%JGMIN,this%JGMAX
440!!NEC$ IVDEP
441! DO i=this%IGMIN,this%IGMAX
442! ! compute mean value of all four face values to obtain cell centered values
443! Txx = 0.25*SUM(T(i,j,:,1,1))
444! Txy = 0.25*SUM(T(i,j,:,1,2))
445! Tyx = 0.25*SUM(T(i,j,:,2,1))
446! Tyy = 0.25*SUM(T(i,j,:,2,2))
447! ! x component of tensor divergence
448! divT(i,j,1) = Divergence3D(this%dAxdy(i,j,1),this%dAxdy(i,j,2), &
449! this%dAydx(i,j,1),this%dAydx(i,j,2), &
450! this%dxdV(i,j),this%dydV(i,j), &
451! this%cxyx%center(i,j),this%cyxy%center(i,j), &
452! 0.0, & ! czxz = 0 because of 2D
453! T(i,j,WEST,1,1),T(i,j,EAST,1,1), &
454! T(i,j,SOUTH,1,2),T(i,j,NORTH,1,2), &
455! Tyx,Tyy,0.0) ! Tzz = 0 because of 2D
456! ! y component of tensor divergence
457! divT(i,j,2) = Divergence3D(this%dAxdy(i,j,1),this%dAxdy(i,j,2), &
458! this%dAydx(i,j,1),this%dAydx(i,j,2), &
459! this%dxdV(i,j),this%dydV(i,j), &
460! -this%cxyx%center(i,j),-this%cyxy%center(i,j), &
461! 0.0, & ! czyz = 0 because of 2D
462! T(i,j,WEST,2,1),T(i,j,EAST,2,1), &
463! T(i,j,SOUTH,2,2),T(i,j,NORTH,2,2), &
464! Txx,Txy,0.0) ! Tzz = 0 because of 2D
465! END DO
466! END DO
467! END SUBROUTINE TensorDivergence2D_2
468
469
474 PURE SUBROUTINE vectordivergence3d(this,vx,vy,vz,divv)
475 IMPLICIT NONE
476 !------------------------------------------------------------------------!
477 CLASS(mesh_midpoint),INTENT(IN) :: this
478 REAL, DIMENSION(this%IGMIN:this%IGMAX,this%JGMIN:this%JGMAX,this%KGMIN:this%KGMAX) &
479 :: vx,vy,vz,divv
480 !------------------------------------------------------------------------!
481 INTEGER :: i,j,k
482 !------------------------------------------------------------------------!
483 INTENT(IN) :: vx,vy,vz
484 INTENT(OUT) :: divv
485 !------------------------------------------------------------------------!
486 ! we simply use the 3D tensor divergence and set the commutator coefficients
487 ! and the off-diagonal tensor components to 0
488 DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
489 DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
490!NEC$ IVDEP
491 DO concurrent(i=this%IGMIN+this%IP1:this%IGMAX-this%IP1)
492 divv(i,j,k) = divergence3d(&
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), &
497 0.0,0.0,0.0,0.0, & ! no commutator coefficients
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)), &
501 0.0,0.0,0.0,0.0) ! no off-diagonal tensor components
502 END DO
503 END DO
504 END DO
505 END SUBROUTINE vectordivergence3d
506
507
508
517 PURE SUBROUTINE tensordivergence3d(this,Txx,Txy,Txz,Tyx,Tyy,Tyz,Tzx,Tzy,Tzz, &
518 divTx,divTy,divTz)
519 IMPLICIT NONE
520 !------------------------------------------------------------------------!
521 CLASS(mesh_midpoint), INTENT(IN) :: this
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
524 !------------------------------------------------------------------------!
525 INTEGER :: i,j,k
526 !------------------------------------------------------------------------!
527 INTENT(IN) :: txx,txy,txz,tyx,tyy,tyz,tzx,tzy,tzz
528 INTENT(OUT) :: divtx,divty,divtz
529 !------------------------------------------------------------------------!
530 DO k=this%KGMIN+this%KP1,this%KGMAX-this%KP1
531 DO j=this%JGMIN+this%JP1,this%JGMAX-this%JP1
532!NEC$ IVDEP
533 DO concurrent(i=this%IGMIN+this%IP1:this%IGMAX-this%IP1)
534 ! x component of tensor divergence
535 divtx(i,j,k) = divergence3d(&
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))
546 ! y component of tensor divergence
547 ! change input of Divergence3D according to the rules given in
548 ! the comments (see below)
549 divty(i,j,k) = divergence3d(&
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))
560 ! z component of tensor divergence
561 ! change input of Divergence3D according to the rules given in
562 ! the comments (see below)
563 divtz(i,j,k) = divergence3d(&
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))
574 END DO
575 END DO
576 END DO
577 END SUBROUTINE tensordivergence3d
578
579
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)
590 IMPLICIT NONE
591 !------------------------------------------------------------------------!
592 REAL,INTENT(IN) :: daxdydzwest,daxdydzeast, & ! surface elements
593 daydxdzsouth,daydxdznorth, & ! dAx/dy, dAy/dx
594 dazdxdybottom,dazdxdytop, &
595 dxdydv,dxdzdv,dydzdv, & ! dx/dV and dy/dV
596 cxyx,cyxy,czxz,cxzx, & ! commutator coeffs
597 txxeast,txxwest, & ! tensor components
598 txynorth,txysouth, & ! on cell faces
599 txzbottom,txztop, &
600 tyxcent,tyycent,tzzcent,tzxcent ! central
601 ! to compute the y-component of the tensor divergence,
602 ! one has to modify the input according to
603 ! Txx <-> Tyx
604 ! Txy <-> Tyy
605 ! Txz -> Tyz
606 ! Tzx -> Tzy
607 ! Tzz -> Tzz
608 ! cxyx -> -cxyx
609 ! cyxy -> -cyxy
610 ! cxzx -> cyzy
611 ! czxz -> czyz
612 ! to compute the z-component of the tensor divergence,
613 ! one has to modify the input according to
614 ! Txx <-> Tzx
615 ! Txy -> Tzy
616 ! Txz <-> Tzz
617 ! Tyx -> Tyz
618 ! Tyy -> Tyy
619 ! cxyx -> czyz
620 ! cyxy -> cyzy
621 ! cxzx -> -cxzx
622 ! czxz -> -czxz
623 !------------------------------------------------------------------------!
624 REAL :: div
625 !------------------------------------------------------------------------!
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
630 END FUNCTION divergence3d
631
632
633
634 SUBROUTINE finalize(this)
635 IMPLICIT NONE
636 !------------------------------------------------------------------------!
637 CLASS(mesh_midpoint),INTENT(INOUT) :: this
638 !------------------------------------------------------------------------!
639 DEALLOCATE(this%dAx,this%dAy,this%dAz,this%dAxdydz,this%dAydzdx,this%dAzdxdy)
640 CALL this%Finalize_base()
641 END SUBROUTINE finalize
642
643END MODULE mesh_midpoint_mod
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
logging_base_mod::nan_default_real
real, save nan_default_real
NaN real constant.
Definition: logging_base.f90:85
logging_base_mod::finalize
subroutine finalize(this)
Destructor of logging_base class.
Definition: logging_base.f90:160
logging_base_mod::default_mpi_real
integer, save default_mpi_real
default real type for MPI
Definition: logging_base.f90:80
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
mesh_base_mod::vector_y
integer, parameter vector_y
Definition: mesh_base.f90:109
mesh_base_mod::vector_z
integer, parameter vector_z
Definition: mesh_base.f90:109
mesh_base_mod::midpoint
integer, parameter midpoint
use midpoint rule to approximate flux integrals
Definition: mesh_base.f90:95
mesh_base_mod::vector_x
integer, parameter vector_x
flags to check which vector components are enabled
Definition: mesh_base.f90:109
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:157
mesh_midpoint_mod::initmesh_midpoint
subroutine initmesh_midpoint(this, config, IO)
Definition: mesh_midpoint.f90:74
mesh_midpoint_mod::vectordivergence3d
pure subroutine vectordivergence3d(this, vx, vy, vz, divv)
compute the cell centered 3D vector divergence
Definition: mesh_midpoint.f90:475
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:590
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:280
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:519
mesh_midpoint_mod::mesh_name
character(len=32), parameter mesh_name
Definition: mesh_midpoint.f90:60
common_dict::dict_typ
Definition: common_dict.f90:123
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
mesh_midpoint_mod::mesh_midpoint
Definition: mesh_midpoint.f90:42
Generated by   doxygen 1.9.4