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!----------------------------------------------------------------------------!
58 USE common_dict
59 IMPLICIT NONE
60 !--------------------------------------------------------------------------!
61 PRIVATE
63 CLASS(marray_base), ALLOCATABLE :: slopes,dx
64 CONTAINS
66 PROCEDURE :: calculatestates
67 PROCEDURE :: calculateslopes
68 PROCEDURE :: finalize
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
87 ! constants
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)
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
Dictionary for generic data types.
Definition: common_dict.f90:61
Basic fosite module.
subroutine finalize(this)
Destructor of logging_base class.
base class for mesh arrays
Definition: marray_base.f90:36
derived class for compound of mesh arrays
basic mesh module
Definition: mesh_base.f90:72
Basic physics module.
base module for reconstruction process
integer, parameter, public linear
module for linear (first order) TVD reconstruction using slope limiters
subroutine calculateslopes(this, Mesh, Physics, rvar)
computes the limited slopes
integer, parameter, public nolimit
integer, parameter, public sweby
integer, parameter, public pp
elemental real function nolimit_limiter(arg1, arg2)
integer, parameter, public vanleer
integer, parameter, public superbee
integer, parameter, public minmod
subroutine initreconstruction_linear(this, Mesh, Physics, config, IO)
Constructor of linear reconstruction module.
character(len=32), parameter recontype_name
elemental real function sweby_limiter(arg1, arg2, param)
elemental real function minmod3_limiter(arg1, arg2, arg3)
elemental real function ospre_limiter(arg1, arg2)
character(len=32), dimension(8), parameter limitertype_name
subroutine calculatestates(this, Mesh, Physics, rvar, rstates)
Reconstructes states at cell boundaries.
elemental real function minmod2_limiter(arg1, arg2)
integer, parameter, public ospre
elemental real function vanleer_limiter(arg1, arg2)
integer, parameter, public monocent
basic mesh array class
Definition: marray_base.f90:69
mesh data structure
Definition: mesh_base.f90:122