reconstruction_linear.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: reconstruction_linear.f90 #
5!# #
6!# Copyright (C) 2007-2021 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
9!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
10!# #
11!# This program is free software; you can redistribute it and/or modify #
12!# it under the terms of the GNU General Public License as published by #
13!# the Free Software Foundation; either version 2 of the License, or (at #
14!# your option) any later version. #
15!# #
16!# This program is distributed in the hope that it will be useful, but #
17!# WITHOUT ANY WARRANTY; without even the implied warranty of #
18!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19!# NON INFRINGEMENT. See the GNU General Public License for more #
20!# details. #
21!# #
22!# You should have received a copy of the GNU General Public License #
23!# along with this program; if not, write to the Free Software #
24!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25!# #
26!#############################################################################
27
28!----------------------------------------------------------------------------!
50!----------------------------------------------------------------------------!
51MODULE reconstruction_linear_mod
52 USE reconstruction_base_mod
53 USE logging_base_mod
54 USE mesh_base_mod
55 USE marray_base_mod
56 USE marray_compound_mod
57 USE physics_base_mod
58 USE common_dict
59 IMPLICIT NONE
60 !--------------------------------------------------------------------------!
61 PRIVATE
62 TYPE, EXTENDS (reconstruction_base) :: reconstruction_linear
63 CLASS(marray_base), ALLOCATABLE :: slopes,dx
64 CONTAINS
65 PROCEDURE :: initreconstruction_linear
66 PROCEDURE :: calculatestates
67 PROCEDURE :: calculateslopes
68 PROCEDURE :: finalize
69 END TYPE reconstruction_linear
70 !--------------------------------------------------------------------------!
71 INTEGER, PARAMETER :: minmod = 1
72 INTEGER, PARAMETER :: monocent = 2
73 INTEGER, PARAMETER :: sweby = 3
74 INTEGER, PARAMETER :: superbee = 4
75 INTEGER, PARAMETER :: ospre = 5
76 INTEGER, PARAMETER :: pp = 6
77 INTEGER, PARAMETER :: vanleer = 7
78 INTEGER, PARAMETER :: nolimit = 8
79 CHARACTER(LEN=32), PARAMETER :: recontype_name = "linear"
80 CHARACTER(LEN=32), DIMENSION(8), PARAMETER :: limitertype_name = (/ &
81 "minmod ", "mc ", "sweby ", "superbee ", "ospre ", &
82 "pp ", "van leer ", "no limit "/)
83 !--------------------------------------------------------------------------!
84 PUBLIC :: &
85 ! types
86 reconstruction_linear, &
87 ! constants
88 minmod, monocent, sweby, superbee, ospre, pp, vanleer, nolimit
89 !--------------------------------------------------------------------------!
90
91CONTAINS
92
93
95 SUBROUTINE initreconstruction_linear(this,Mesh,Physics,config,IO)
96 IMPLICIT NONE
97 !------------------------------------------------------------------------!
98 CLASS(reconstruction_linear), INTENT(INOUT) :: this
99 CLASS(mesh_base), INTENT(IN) :: Mesh
100 CLASS(physics_base), INTENT(IN) :: Physics
101 TYPE(dict_typ), POINTER :: config
102 TYPE(dict_typ), POINTER :: IO
103 !------------------------------------------------------------------------!
104 CHARACTER(LEN=60) :: key
105 INTEGER :: err,limiter,valwrite,i,j,k,l,m
106 REAL :: theta
107 !------------------------------------------------------------------------!
108 ! allocate memory for all arrays used in reconstruction_linear
109 ALLOCATE(this%slopes,this%dx, &
110 stat = err)
111 IF (err.NE.0) THEN
112 CALL this%Error("InitReconstruction_linear", "Unable to allocate memory.")
113 END IF
114
115
116 ! initialize parent
117 CALL this%InitReconstruction(mesh,physics,config,io,linear,recontype_name)
118
119 ! set defaults for the limiter
120 CALL getattr(config, "limiter", limiter, minmod)
121 CALL getattr(config, "theta", theta, 1.0)
122
123 ! limiter settings
124 SELECT CASE(limiter)
125 CASE(minmod,monocent,sweby,superbee,ospre,pp,vanleer,nolimit)
126 CALL this%InitLogging(limiter,limitertype_name(limiter))
127 SELECT CASE(limiter)
128 CASE(sweby,monocent)
129 this%limiter_param = theta
130 END SELECT
131 CASE DEFAULT
132 CALL this%Error("InitReconstruction_linear", "Unknown limiter")
133 END SELECT
134
135 ! check parameter settings for PP limiter
136 IF (limiter.EQ.pp) THEN
137 CALL this%Error("InitReconstruction_linear","PP limiter currently not supported.")
138! IF (this%limiter_param.LT.0.0) &
139! CALL this%Error("InitReconstruction_linear", &
140! "Parameter of PP limiter must be greater than 0")
141! IF (this%limiter_param.GT.EPSILON(this%limiter_param)) &
142! CALL this%Warning("InitReconstruction_linear", &
143! "Parameter of PP limiter should be less than machine precision.")
144 END IF
145
146 ! create new rank 2 mesh array for slopes
147 this%slopes = marray_base(mesh%NDIMS,physics%VNUM)
148
149 ! create new rank 1 mesh array for reconstruction points
150 this%dx = marray_base(mesh%NFACES)
151
152 ! zero the slopes
153 this%slopes%data1d(:) = 0.0
154
155 ! initialize coordinate differences for reconstruction
156 m = 1
157 IF (mesh%INUM.GT.1) THEN
158 ! reconstruct along x-direction
159 DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
160 this%dx%data4d(i,j,k,m) = mesh%curv%faces(i,j,k,1,1) - mesh%curv%bcenter(i,j,k,1)
161 this%dx%data4d(i,j,k,m+1) = mesh%curv%faces(i,j,k,2,1) - mesh%curv%bcenter(i,j,k,1)
162 END DO
163 m = m + 2
164 END IF
165 IF (mesh%JNUM.GT.1) THEN
166 DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
167 this%dx%data4d(i,j,k,m) = mesh%curv%faces(i,j,k,3,2) - mesh%curv%bcenter(i,j,k,2)
168 this%dx%data4d(i,j,k,m+1) = mesh%curv%faces(i,j,k,4,2) - mesh%curv%bcenter(i,j,k,2)
169 END DO
170 m = m + 2
171 END IF
172 IF (mesh%KNUM.GT.1) THEN
173 DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
174 this%dx%data4d(i,j,k,m) = mesh%curv%faces(i,j,k,5,3) - mesh%curv%bcenter(i,j,k,3)
175 this%dx%data4d(i,j,k,m+1) = mesh%curv%faces(i,j,k,6,3) - mesh%curv%bcenter(i,j,k,3)
176 END DO
177 END IF
178
179 CALL getattr(config, "output/slopes", valwrite, 0)
180 IF(valwrite.EQ.1) THEN
181 m = 1
182 IF (mesh%INUM.GT.1) THEN
183!NEC$ NOVECTOR
184 DO l=1,physics%VNUM
185 key = trim(physics%pvarname(l)) // "_xslope"
186 CALL setattr(io, trim(key), &
187 this%slopes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m,l))
188 END DO
189 m = m + 1
190 END IF
191 IF (mesh%JNUM.GT.1) THEN
192!NEC$ NOVECTOR
193 DO l=1,physics%VNUM
194 key = trim(physics%pvarname(l)) // "_yslope"
195 CALL setattr(io, trim(key), &
196 this%slopes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m,l))
197 END DO
198 m = m + 1
199 END IF
200 IF (mesh%KNUM.GT.1) THEN
201!NEC$ NOVECTOR
202 DO l=1,physics%VNUM
203 key = trim(physics%pvarname(l)) // "_zslope"
204 CALL setattr(io, trim(key), &
205 this%slopes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m,l))
206 END DO
207 END IF
208 END IF
209
210 CALL this%Info(" limiter: " // trim(this%GetName()))
211 END SUBROUTINE initreconstruction_linear
212
213
215 SUBROUTINE calculatestates(this,Mesh,Physics,rvar,rstates)
216 IMPLICIT NONE
217 !------------------------------------------------------------------------!
218 CLASS(reconstruction_linear), INTENT(INOUT) :: this
219 CLASS(mesh_base), INTENT(IN) :: Mesh
220 CLASS(physics_base), INTENT(IN) :: Physics
221 CLASS(marray_compound), INTENT(INOUT) :: rvar
222 CLASS(marray_compound), INTENT(INOUT) :: rstates
223 !------------------------------------------------------------------------!
224 INTEGER :: i,l,n
225 !------------------------------------------------------------------------!
226 ! calculate slopes first
227 CALL this%CalculateSlopes(mesh,physics,rvar)
228
229 ! reconstruct cell face values
230!NEC$ SHORTLOOP
231 DO l=1,physics%VNUM
232!NEC$ SHORTLOOP
233 DO n=1,mesh%NFACES
234 DO concurrent(i=1:SIZE(rvar%data2d,dim=1))
235 rstates%data3d(i,n,l) = rvar%data2d(i,l) + &
236 this%slopes%data3d(i,((n+1)/2),l)*this%dx%data2d(i,n)
237 END DO
238 END DO
239 END DO
240 END SUBROUTINE calculatestates
241
243 SUBROUTINE calculateslopes(this,Mesh,Physics,rvar)
244 IMPLICIT NONE
245 !------------------------------------------------------------------------!
246 CLASS(reconstruction_linear), INTENT(INOUT) :: this
247 CLASS(mesh_base), INTENT(IN) :: Mesh
248 CLASS(physics_base), INTENT(IN) :: Physics
249 CLASS(marray_compound), INTENT(INOUT) :: rvar
250 !------------------------------------------------------------------------!
251 INTEGER :: i,j,k,l,m
252 !------------------------------------------------------------------------!
253 m = 1 ! this counts the dimensions used for computation (last index in slopes)
254 ! choose limiter & check for cases where dimension needs to be neclected
255 SELECT CASE(this%GetType())
256 CASE(minmod)
257 ! calculate slopes in x-direction
258 IF (mesh%INUM.GT.1) THEN
259!NEC$ SHORTLOOP
260 DO l=1,physics%VNUM
261 DO k=mesh%KGMIN,mesh%KGMAX
262 DO j=mesh%JGMIN,mesh%JGMAX
263!NEC$ IVDEP
264 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
265 this%slopes%data5d(i,j,k,m,l) = mesh%invdx * minmod2_limiter(&
266 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
267 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
268 END DO
269 END DO
270 END DO
271 END DO
272 m = m + 1
273 END IF
274 ! calculate slopes in y-direction
275 IF (mesh%JNUM.GT.1) THEN
276 i = mesh%IGMAX-mesh%IGMIN+1
277!NEC$ SHORTLOOP
278 DO l=1,physics%VNUM
279!NEC$ NOVECTOR
280 DO k=mesh%KGMIN,mesh%KGMAX
281! use collapsed arrays for better vectorization
282!NEC$ NOVECTOR
283 DO j=i+1,SIZE(rvar%data3d,1)-i
284 this%slopes%data4d(j,k,m,l) = mesh%invdy * minmod2_limiter(&
285 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l), &
286 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
287 END DO
288 END DO
289 END DO
290 m = m + 1
291 END IF
292 ! calculate slopes in z-direction
293 IF (mesh%KNUM.GT.1) THEN
294 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
295!NEC$ SHORTLOOP
296 DO l=1,physics%VNUM
297! use collapsed arrays for better vectorization
298!NEC$ IVDEP
299 DO j=i+1,SIZE(rvar%data2d,1)-i
300 this%slopes%data3d(j,m,l) = mesh%invdz * minmod2_limiter(&
301 rvar%data2d(j,l) - rvar%data2d(j-i,l), &
302 rvar%data2d(j+i,l) - rvar%data2d(j,l))
303 END DO
304 END DO
305 END IF
306 CASE(monocent)
307 ! calculate slopes in x-direction
308 IF (mesh%INUM.GT.1) THEN
309!NEC$ SHORTLOOP
310 DO l=1,physics%VNUM
311 DO k=mesh%KGMIN,mesh%KGMAX
312 DO j=mesh%JGMIN,mesh%JGMAX
313!NEC$ IVDEP
314 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
315 this%slopes%data5d(i,j,k,m,l) = mesh%invdx * minmod3_limiter(&
316 this%limiter_param*(rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l)),&
317 this%limiter_param*(rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l)),&
318 0.5*(rvar%data4d(i+1,j,k,l) - rvar%data4d(i-1,j,k,l)))
319 END DO
320 END DO
321 END DO
322 END DO
323 m = m + 1
324 END IF
325 ! calculate slopes in y-direction
326 IF (mesh%JNUM.GT.1) THEN
327 i = mesh%IGMAX-mesh%IGMIN+1
328!NEC$ SHORTLOOP
329 DO l=1,physics%VNUM
330 DO k=mesh%KGMIN,mesh%KGMAX
331! use collapsed arrays for better vectorization
332!NEC$ IVDEP
333 DO j=i+1,SIZE(rvar%data3d,1)-i
334 this%slopes%data4d(j,k,m,l) = mesh%invdy * minmod3_limiter(&
335 this%limiter_param*(rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l)),&
336 this%limiter_param*(rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l)),&
337 0.5*(rvar%data3d(j+i,k,l) - rvar%data3d(j-i,k,l)))
338 END DO
339 END DO
340 END DO
341 m = m + 1
342 END IF
343 ! calculate slopes in z-direction
344 IF (mesh%KNUM.GT.1) THEN
345 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
346!NEC$ SHORTLOOP
347 DO l=1,physics%VNUM
348! use collapsed arrays for better vectorization
349!NEC$ IVDEP
350 DO j=i+1,SIZE(rvar%data2d,1)-i
351 this%slopes%data3d(j,m,l) = mesh%invdz * minmod3_limiter(&
352 this%limiter_param*(rvar%data2d(j,l) - rvar%data2d(j-i,l)),&
353 this%limiter_param*(rvar%data2d(j+i,l) - rvar%data2d(j,l)),&
354 0.5*(rvar%data2d(j+i,l) - rvar%data2d(j-i,l)))
355 END DO
356 END DO
357 END IF
358 CASE(sweby)
359 ! calculate slopes in x-direction
360 IF (mesh%INUM.GT.1) THEN
361!NEC$ SHORTLOOP
362 DO l=1,physics%VNUM
363 DO k=mesh%KGMIN,mesh%KGMAX
364 DO j=mesh%JGMIN,mesh%JGMAX
365!NEC$ IVDEP
366 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
367 this%slopes%data5d(i,j,k,m,l) = mesh%invdx * sweby_limiter(&
368 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
369 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l),&
370 this%limiter_param)
371 END DO
372 END DO
373 END DO
374 END DO
375 m = m + 1
376 END IF
377 ! calculate slopes in y-direction
378 IF (mesh%JNUM.GT.1) THEN
379 i = mesh%IGMAX-mesh%IGMIN+1
380!NEC$ SHORTLOOP
381 DO l=1,physics%VNUM
382 DO k=mesh%KGMIN,mesh%KGMAX
383! use collapsed arrays for better vectorization
384!NEC$ IVDEP
385 DO j=i+1,SIZE(rvar%data3d,1)-i
386 this%slopes%data4d(j,k,m,l) = mesh%invdy * sweby_limiter(&
387 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
388 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l),&
389 this%limiter_param)
390 END DO
391 END DO
392 END DO
393 m = m + 1
394 END IF
395 ! calculate slopes in z-direction
396 IF (mesh%KNUM.GT.1) THEN
397 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
398!NEC$ SHORTLOOP
399 DO l=1,physics%VNUM
400! use collapsed arrays for better vectorization
401!NEC$ IVDEP
402 DO j=i+1,SIZE(rvar%data2d,1)-i
403 this%slopes%data3d(j,m,l) = mesh%invdz * sweby_limiter(&
404 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
405 rvar%data2d(j+i,l) - rvar%data2d(j,l),&
406 this%limiter_param)
407 END DO
408 END DO
409 END IF
410 CASE(superbee)
411 ! calculate slopes in x-direction
412 IF (mesh%INUM.GT.1) THEN
413!NEC$ SHORTLOOP
414 DO l=1,physics%VNUM
415 DO k=mesh%KGMIN,mesh%KGMAX
416 DO j=mesh%JGMIN,mesh%JGMAX
417!NEC$ IVDEP
418 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
419 this%slopes%data5d(i,j,k,m,l) = mesh%invdx * sweby_limiter(&
420 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
421 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l), 2.0)
422 END DO
423 END DO
424 END DO
425 END DO
426 m = m + 1
427 END IF
428 ! calculate slopes in y-direction
429 IF (mesh%JNUM.GT.1) THEN
430 i = mesh%IGMAX-mesh%IGMIN+1
431!NEC$ SHORTLOOP
432 DO l=1,physics%VNUM
433 DO k=mesh%KGMIN,mesh%KGMAX
434! use collapsed arrays for better vectorization
435!NEC$ IVDEP
436 DO j=i+1,SIZE(rvar%data3d,1)-i
437 this%slopes%data4d(j,k,m,l) = mesh%invdy * sweby_limiter(&
438 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
439 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l), 2.0)
440 END DO
441 END DO
442 END DO
443 m = m + 1
444 END IF
445 ! calculate slopes in z-direction
446 IF (mesh%KNUM.GT.1) THEN
447 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
448!NEC$ SHORTLOOP
449 DO l=1,physics%VNUM
450! use collapsed arrays for better vectorization
451!NEC$ IVDEP
452 DO j=i+1,SIZE(rvar%data2d,1)-i
453 this%slopes%data3d(j,m,l) = mesh%invdz * sweby_limiter(&
454 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
455 rvar%data2d(j+i,l) - rvar%data2d(j,l), 2.0)
456 END DO
457 END DO
458 END IF
459 CASE(ospre)
460 ! calculate slopes in x-direction
461 IF (mesh%INUM.GT.1) THEN
462!NEC$ SHORTLOOP
463 DO l=1,physics%VNUM
464 DO k=mesh%KGMIN,mesh%KGMAX
465 DO j=mesh%JGMIN,mesh%JGMAX
466!NEC$ IVDEP
467 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
468 this%slopes%data5d(i,j,k,m,l) = mesh%invdx * ospre_limiter(&
469 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
470 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
471 END DO
472 END DO
473 END DO
474 END DO
475 m = m + 1
476 END IF
477 ! calculate slopes in y-direction
478 IF (mesh%JNUM.GT.1) THEN
479 i = mesh%IGMAX-mesh%IGMIN+1
480!NEC$ SHORTLOOP
481 DO l=1,physics%VNUM
482 DO k=mesh%KGMIN,mesh%KGMAX
483! use collapsed arrays for better vectorization
484!NEC$ IVDEP
485 DO j=i+1,SIZE(rvar%data3d,1)-i
486 this%slopes%data4d(j,k,m,l) = mesh%invdy * ospre_limiter(&
487 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
488 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
489 END DO
490 END DO
491 END DO
492 m = m + 1
493 END IF
494 ! calculate slopes in z-direction
495 IF (mesh%KNUM.GT.1) THEN
496 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
497!NEC$ SHORTLOOP
498 DO l=1,physics%VNUM
499! use collapsed arrays for better vectorization
500!NEC$ IVDEP
501 DO j=i+1,SIZE(rvar%data2d,1)-i
502 this%slopes%data3d(j,m,l) = mesh%invdz * ospre_limiter(&
503 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
504 rvar%data2d(j+i,l) - rvar%data2d(j,l))
505 END DO
506 END DO
507 END IF
508! CASE(PP)
509! ! calculate slopes in both-directions
510! !NEC$ SHORTLOOP
511! DO l=1,Physics%VNUM
512! DO k=Mesh%KGMIN,Mesh%KGMAX
513! DO j=Mesh%JGMIN+Mesh%JP1,Mesh%JGMAX+Mesh%JM1
514!!NEC$ IVDEP
515! DO i=Mesh%IGMIN+Mesh%IP1,Mesh%IGMAX+Mesh%IM1
516! CALL pp_limiter(this%slopes%data5d(i,j,k,m,l),&
517! this%slopes%data5d(i,j,k,m,l),&
518! rvar%data4d(i-1,j+1,k)-rvar%data4d(i,j,k),&
519! rvar%data4d(i ,j+1,k)-rvar%data4d(i,j,k),&
520! rvar%data4d(i+1,j+1,k)-rvar%data4d(i,j,k),&
521! rvar%data4d(i-1,j ,k)-rvar%data4d(i,j,k),&
522! this%limiter_param,&
523! rvar%data4d(i+1,j ,k)-rvar%data4d(i,j,k),&
524! rvar%data4d(i-1,j-1,k)-rvar%data4d(i,j,k),&
525! rvar%data4d(i ,j-1,k)-rvar%data4d(i,j,k),&
526! rvar%data4d(i+1,j-1,k)-rvar%data4d(i,j,k))
527! END DO
528! END DO
529! END DO
530! END DO
531 CASE(vanleer)
532 ! calculate slopes in x-direction
533 IF (mesh%INUM.GT.1) THEN
534!NEC$ SHORTLOOP
535 DO l=1,physics%VNUM
536 DO k=mesh%KGMIN,mesh%KGMAX
537 DO j=mesh%JGMIN,mesh%JGMAX
538!NEC$ IVDEP
539 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
540 this%slopes%data5d(i,j,k,m,l) = mesh%invdx * vanleer_limiter(&
541 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
542 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
543 END DO
544 END DO
545 END DO
546 END DO
547 m = m + 1
548 END IF
549 ! calculate slopes in y-direction
550 IF (mesh%JNUM.GT.1) THEN
551 i = mesh%IGMAX-mesh%IGMIN+1
552!NEC$ SHORTLOOP
553 DO l=1,physics%VNUM
554 DO k=mesh%KGMIN,mesh%KGMAX
555! use collapsed arrays for better vectorization
556!NEC$ IVDEP
557 DO j=i+1,SIZE(rvar%data3d,1)-i
558 this%slopes%data4d(j,k,m,l) = mesh%invdy * vanleer_limiter(&
559 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
560 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
561 END DO
562 END DO
563 END DO
564 m = m + 1
565 END IF
566 ! calculate slopes in z-direction
567 IF (mesh%KNUM.GT.1) THEN
568 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
569!NEC$ SHORTLOOP
570 DO l=1,physics%VNUM
571! use collapsed arrays for better vectorization
572!NEC$ IVDEP
573 DO j=i+1,SIZE(rvar%data2d,1)-i
574 this%slopes%data3d(j,m,l) = mesh%invdz * vanleer_limiter(&
575 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
576 rvar%data2d(j+i,l) - rvar%data2d(j,l))
577 END DO
578 END DO
579 END IF
580 CASE(nolimit)
581 ! calculate slopes in x-direction
582 IF (mesh%INUM.GT.1) THEN
583!NEC$ SHORTLOOP
584 DO l=1,physics%VNUM
585 DO k=mesh%KGMIN,mesh%KGMAX
586 DO j=mesh%JGMIN,mesh%JGMAX
587!NEC$ IVDEP
588 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
589 this%slopes%data5d(i,j,k,m,l) = mesh%invdx * nolimit_limiter(&
590 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
591 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
592 END DO
593 END DO
594 END DO
595 END DO
596 m = m + 1
597 END IF
598 ! calculate slopes in y-direction
599 IF (mesh%JNUM.GT.1) THEN
600 i = mesh%IGMAX-mesh%IGMIN+1
601!NEC$ SHORTLOOP
602 DO l=1,physics%VNUM
603 DO k=mesh%KGMIN,mesh%KGMAX
604! use collapsed arrays for better vectorization
605!NEC$ IVDEP
606 DO j=i+1,SIZE(rvar%data3d,1)-i
607 this%slopes%data4d(j,k,m,l) = mesh%invdy * nolimit_limiter(&
608 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
609 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
610 END DO
611 END DO
612 END DO
613 m = m + 1
614 END IF
615 ! calculate slopes in z-direction
616 IF (mesh%KNUM.GT.1) THEN
617 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
618!NEC$ SHORTLOOP
619 DO l=1,physics%VNUM
620! use collapsed arrays for better vectorization
621!NEC$ IVDEP
622 DO j=i+1,SIZE(rvar%data2d,1)-i
623 this%slopes%data3d(j,m,l) = mesh%invdz * nolimit_limiter(&
624 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
625 rvar%data2d(j+i,l) - rvar%data2d(j,l))
626 END DO
627 END DO
628 END IF
629 END SELECT
630 END SUBROUTINE calculateslopes
631
632 SUBROUTINE finalize(this)
633 IMPLICIT NONE
634 !------------------------------------------------------------------------!
635 CLASS(reconstruction_linear), INTENT(INOUT) :: this
636 !------------------------------------------------------------------------!
637 DEALLOCATE(this%slopes,this%dx)
638
639 CALL this%Finalize_base()
640 END SUBROUTINE finalize
641
642!----------------------------------------------------------------------------!
644
645 ELEMENTAL FUNCTION minmod2_limiter(arg1,arg2) RESULT(limarg)
646 IMPLICIT NONE
647 !--------------------------------------------------------------------!
648 REAL :: limarg
649 REAL, INTENT(IN) :: arg1, arg2
650 !--------------------------------------------------------------------!
651 IF (sign(1.0,arg1)*sign(1.0,arg2).GT.0) THEN
652 limarg = sign(min(abs(arg1),abs(arg2)),arg1)
653 ELSE
654 limarg = 0.
655 END IF
656 END FUNCTION minmod2_limiter
657
658 ELEMENTAL FUNCTION minmod3_limiter(arg1,arg2,arg3) RESULT(limarg)
659 IMPLICIT NONE
660 !--------------------------------------------------------------------!
661 REAL :: limarg
662 REAL, INTENT(IN) :: arg1, arg2, arg3
663 !--------------------------------------------------------------------!
664 IF (((sign(1.0,arg1)*sign(1.0,arg2)).GT.0).AND.&
665 ((sign(1.0,arg2)*sign(1.0,arg3)).GT.0)) THEN
666 limarg = sign(min(abs(arg1),abs(arg2),abs(arg3)),arg1)
667 ELSE
668 limarg = 0.
669 END IF
670 END FUNCTION minmod3_limiter
671
672 ELEMENTAL FUNCTION sweby_limiter(arg1,arg2,param) RESULT(limarg)
673 IMPLICIT NONE
674 !--------------------------------------------------------------------!
675 REAL :: limarg
676 REAL, INTENT(IN) :: arg1, arg2, param
677 !--------------------------------------------------------------------!
678 IF (sign(1.0,arg1)*sign(1.0,arg2).GT.0) THEN
679 limarg = sign(max(min(param*abs(arg1),abs(arg2)), &
680 min(abs(arg1),param*abs(arg2))),arg1)
681 ELSE
682 limarg = 0.
683 END IF
684 END FUNCTION sweby_limiter
685
686 ELEMENTAL FUNCTION ospre_limiter(arg1,arg2) RESULT(limarg)
687 IMPLICIT NONE
688 !--------------------------------------------------------------------!
689 REAL :: limarg
690 REAL, INTENT(IN) :: arg1, arg2
691 !--------------------------------------------------------------------!
692 limarg = 1.5*arg1*arg2*(arg1 + arg2) / (arg1*(arg1 + 0.5*arg2) &
693 + arg2*(arg2 + 0.5*arg1) + tiny(1.0))
694 END FUNCTION ospre_limiter
695
696! ELEMENTAL SUBROUTINE pp_limiter(xslope,yslope,arg1,arg2,arg3,arg4,param,arg6,arg7,arg8,arg9) ! TODO: 3D
697! IMPLICIT NONE
698! !--------------------------------------------------------------------!
699! REAL, INTENT(OUT):: xslope,yslope,zslope
700! REAL, INTENT(IN) :: arg1,arg2,arg3,arg4,param,arg6,arg7,arg8,arg9
701! !--------------------------------------------------------------------!
702! REAL :: Vmin,Vmax,V
703! !--------------------------------------------------------------------!
704! xslope = (arg6-arg4)*0.5
705! yslope = (arg2-arg8)*0.5 ! TODO: 3D add zslope
706! Vmin = MIN(arg1,arg2,arg3,arg4,-param,arg6,arg7,arg8,arg9)
707! Vmax = MAX(arg1,arg2,arg3,arg4,+param,arg6,arg7,arg8,arg9)
708! V = MIN(1.0,2.0*MIN(ABS(Vmin),ABS(Vmax))/(ABS(xslope)+ABS(yslope)))
709! xslope = V*xslope*Mesh%invdx
710! yslope = V*yslope*Mesh%invdy ! TODO: 3D add zslope
711! END SUBROUTINE pp_limiter
712
713 ELEMENTAL FUNCTION vanleer_limiter(arg1,arg2) RESULT(limarg)
714 IMPLICIT NONE
715 !--------------------------------------------------------------------!
716 REAL :: limarg
717 REAL, INTENT(IN) :: arg1, arg2
718 !--------------------------------------------------------------------!
719 REAL :: a1,a2
720 !--------------------------------------------------------------------!
721 a1 = abs(arg1)
722 a2 = abs(arg2)
723 limarg = (arg1 * a2 + arg2 * a1)/(a1 + a2 + tiny(a1))
724 END FUNCTION vanleer_limiter
725
726 ELEMENTAL FUNCTION nolimit_limiter(arg1,arg2) RESULT(limarg)
727 IMPLICIT NONE
728 !--------------------------------------------------------------------!
729 REAL :: limarg
730 REAL, INTENT(IN) :: arg1, arg2
731 !--------------------------------------------------------------------!
732 limarg = 0.5*(arg1+arg2)
733 END FUNCTION nolimit_limiter
734
735END MODULE reconstruction_linear_mod
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
logging_base_mod::finalize
subroutine finalize(this)
Destructor of logging_base class.
Definition: logging_base.f90:160
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
reconstruction_base_mod
base module for reconstruction process
Definition: reconstruction_base.f90:35
reconstruction_base_mod::linear
integer, parameter, public linear
Definition: reconstruction_base.f90:79
reconstruction_linear_mod
module for linear (first order) TVD reconstruction using slope limiters
Definition: reconstruction_linear.f90:51
reconstruction_linear_mod::calculateslopes
subroutine calculateslopes(this, Mesh, Physics, rvar)
computes the limited slopes
Definition: reconstruction_linear.f90:244
reconstruction_linear_mod::nolimit
integer, parameter, public nolimit
Definition: reconstruction_linear.f90:78
reconstruction_linear_mod::sweby
integer, parameter, public sweby
Definition: reconstruction_linear.f90:73
reconstruction_linear_mod::pp
integer, parameter, public pp
Definition: reconstruction_linear.f90:76
reconstruction_linear_mod::nolimit_limiter
elemental real function nolimit_limiter(arg1, arg2)
Definition: reconstruction_linear.f90:727
reconstruction_linear_mod::vanleer
integer, parameter, public vanleer
Definition: reconstruction_linear.f90:77
reconstruction_linear_mod::superbee
integer, parameter, public superbee
Definition: reconstruction_linear.f90:74
reconstruction_linear_mod::minmod
integer, parameter, public minmod
Definition: reconstruction_linear.f90:71
reconstruction_linear_mod::initreconstruction_linear
subroutine initreconstruction_linear(this, Mesh, Physics, config, IO)
Constructor of linear reconstruction module.
Definition: reconstruction_linear.f90:96
reconstruction_linear_mod::recontype_name
character(len=32), parameter recontype_name
Definition: reconstruction_linear.f90:79
reconstruction_linear_mod::sweby_limiter
elemental real function sweby_limiter(arg1, arg2, param)
Definition: reconstruction_linear.f90:673
reconstruction_linear_mod::minmod3_limiter
elemental real function minmod3_limiter(arg1, arg2, arg3)
Definition: reconstruction_linear.f90:659
reconstruction_linear_mod::ospre_limiter
elemental real function ospre_limiter(arg1, arg2)
Definition: reconstruction_linear.f90:687
reconstruction_linear_mod::limitertype_name
character(len=32), dimension(8), parameter limitertype_name
Definition: reconstruction_linear.f90:80
reconstruction_linear_mod::calculatestates
subroutine calculatestates(this, Mesh, Physics, rvar, rstates)
Reconstructes states at cell boundaries.
Definition: reconstruction_linear.f90:216
reconstruction_linear_mod::minmod2_limiter
elemental real function minmod2_limiter(arg1, arg2)
Definition: reconstruction_linear.f90:646
reconstruction_linear_mod::ospre
integer, parameter, public ospre
Definition: reconstruction_linear.f90:75
reconstruction_linear_mod::vanleer_limiter
elemental real function vanleer_limiter(arg1, arg2)
Definition: reconstruction_linear.f90:714
reconstruction_linear_mod::monocent
integer, parameter, public monocent
Definition: reconstruction_linear.f90:72
common_dict::dict_typ
Definition: common_dict.f90:123
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:69
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
physics_base_mod::physics_base
Definition: physics_base.f90:69
reconstruction_base_mod::reconstruction_base
Definition: reconstruction_base.f90:44
reconstruction_linear_mod::reconstruction_linear
Definition: reconstruction_linear.f90:62
Generated by   doxygen 1.9.4