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-2017 #
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 !----------------------------------------------------------------------------!
51 MODULE 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 
91 CONTAINS
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,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  this%dx%data4d(:,:,:,m) = mesh%curv%faces(:,:,:,1,1) - mesh%curv%bcenter(:,:,:,1)
160  this%dx%data4d(:,:,:,m+1) = mesh%curv%faces(:,:,:,2,1) - mesh%curv%bcenter(:,:,:,1)
161  m = m + 2
162  END IF
163  IF (mesh%JNUM.GT.1) THEN
164  this%dx%data4d(:,:,:,m) = mesh%curv%faces(:,:,:,3,2) - mesh%curv%bcenter(:,:,:,2)
165  this%dx%data4d(:,:,:,m+1) = mesh%curv%faces(:,:,:,4,2) - mesh%curv%bcenter(:,:,:,2)
166  m = m + 2
167  END IF
168  IF (mesh%KNUM.GT.1) THEN
169  this%dx%data4d(:,:,:,m) = mesh%curv%faces(:,:,:,5,3) - mesh%curv%bcenter(:,:,:,3)
170  this%dx%data4d(:,:,:,m+1) = mesh%curv%faces(:,:,:,6,3) - mesh%curv%bcenter(:,:,:,3)
171  END IF
172 
173  CALL getattr(config, "output/slopes", valwrite, 0)
174  IF(valwrite.EQ.1) THEN
175  m = 1
176  IF (mesh%INUM.GT.1) THEN
177 !NEC$ NOVECTOR
178  DO l=1,physics%VNUM
179  key = trim(physics%pvarname(l)) // "_xslope"
180  CALL setattr(io, trim(key), &
181  this%slopes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m,l))
182  END DO
183  m = m + 1
184  END IF
185  IF (mesh%JNUM.GT.1) THEN
186 !NEC$ NOVECTOR
187  DO l=1,physics%VNUM
188  key = trim(physics%pvarname(l)) // "_yslope"
189  CALL setattr(io, trim(key), &
190  this%slopes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m,l))
191  END DO
192  m = m + 1
193  END IF
194  IF (mesh%KNUM.GT.1) THEN
195 !NEC$ NOVECTOR
196  DO l=1,physics%VNUM
197  key = trim(physics%pvarname(l)) // "_zslope"
198  CALL setattr(io, trim(key), &
199  this%slopes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m,l))
200  END DO
201  END IF
202  END IF
203 
204  CALL this%Info(" limiter: " // trim(this%GetName()))
205  END SUBROUTINE initreconstruction_linear
206 
207 
209  PURE SUBROUTINE calculatestates(this,Mesh,Physics,rvar,rstates)
210  IMPLICIT NONE
211  !------------------------------------------------------------------------!
212  CLASS(reconstruction_linear), INTENT(INOUT) :: this
213  CLASS(mesh_base), INTENT(IN) :: mesh
214  CLASS(physics_base), INTENT(IN) :: physics
215  CLASS(marray_compound), INTENT(INOUT) :: rvar
216  CLASS(marray_compound), INTENT(INOUT) :: rstates
217  !------------------------------------------------------------------------!
218  INTEGER :: l,n
219  !------------------------------------------------------------------------!
220  ! calculate slopes first
221  CALL this%CalculateSlopes(mesh,physics,rvar)
222 
223  ! reconstruct cell face values
224 !NEC$ SHORTLOOP
225  DO l=1,physics%VNUM
226 !NEC$ SHORTLOOP
227  DO n=1,mesh%NFACES
228  rstates%data3d(:,n,l) = rvar%data2d(:,l) + &
229  this%slopes%data3d(:,((n+1)/2),l)*this%dx%data2d(:,n)
230  END DO
231  END DO
232  END SUBROUTINE calculatestates
233 
235  PURE SUBROUTINE calculateslopes(this,Mesh,Physics,rvar)
236  IMPLICIT NONE
237  !------------------------------------------------------------------------!
238  CLASS(reconstruction_linear), INTENT(INOUT) :: this
239  CLASS(mesh_base), INTENT(IN) :: mesh
240  CLASS(physics_base), INTENT(IN) :: physics
241  CLASS(marray_compound), INTENT(INOUT) :: rvar
242  !------------------------------------------------------------------------!
243  INTEGER :: i,j,k,l,m
244  !------------------------------------------------------------------------!
245  m = 1 ! this counts the dimensions used for computation (last index in slopes)
246  ! choose limiter & check for cases where dimension needs to be neclected
247  SELECT CASE(this%GetType())
248  CASE(minmod)
249  ! calculate slopes in x-direction
250  IF (mesh%INUM.GT.1) THEN
251 !NEC$ SHORTLOOP
252  DO l=1,physics%VNUM
253  DO k=mesh%KGMIN,mesh%KGMAX
254  DO j=mesh%JGMIN,mesh%JGMAX
255 !NEC$ IVDEP
256  DO i=mesh%IGMIN+mesh%IP1,mesh%IGMAX+mesh%IM1
257  this%slopes%data5d(i,j,k,m,l) = mesh%invdx * minmod2_limiter(&
258  rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
259  rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
260  END DO
261  END DO
262  END DO
263  END DO
264  m = m + 1
265  END IF
266  ! calculate slopes in y-direction
267  IF (mesh%JNUM.GT.1) THEN
268  i = mesh%IGMAX-mesh%IGMIN+1
269 !NEC$ SHORTLOOP
270  DO l=1,physics%VNUM
271 !NEC$ NOVECTOR
272  DO k=mesh%KGMIN,mesh%KGMAX
273 ! use collapsed arrays for better vectorization
274 !NEC$ NOVECTOR
275  DO j=i+1,SIZE(rvar%data3d,1)-i
276  this%slopes%data4d(j,k,m,l) = mesh%invdy * minmod2_limiter(&
277  rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l), &
278  rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
279  END DO
280  END DO
281  END DO
282  m = m + 1
283  END IF
284  ! calculate slopes in z-direction
285  IF (mesh%KNUM.GT.1) THEN
286  i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
287 !NEC$ SHORTLOOP
288  DO l=1,physics%VNUM
289 ! use collapsed arrays for better vectorization
290 !NEC$ IVDEP
291  DO j=i+1,SIZE(rvar%data2d,1)-i
292  this%slopes%data3d(j,m,l) = mesh%invdz * minmod2_limiter(&
293  rvar%data2d(j,l) - rvar%data2d(j-i,l), &
294  rvar%data2d(j+i,l) - rvar%data2d(j,l))
295  END DO
296  END DO
297  END IF
298  CASE(monocent)
299  ! calculate slopes in x-direction
300  IF (mesh%INUM.GT.1) THEN
301 !NEC$ SHORTLOOP
302  DO l=1,physics%VNUM
303  DO k=mesh%KGMIN,mesh%KGMAX
304  DO j=mesh%JGMIN,mesh%JGMAX
305 !NEC$ IVDEP
306  DO i=mesh%IGMIN+mesh%IP1,mesh%IGMAX+mesh%IM1
307  this%slopes%data5d(i,j,k,m,l) = mesh%invdx * minmod3_limiter(&
308  this%limiter_param*(rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l)),&
309  this%limiter_param*(rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l)),&
310  0.5*(rvar%data4d(i+1,j,k,l) - rvar%data4d(i-1,j,k,l)))
311  END DO
312  END DO
313  END DO
314  END DO
315  m = m + 1
316  END IF
317  ! calculate slopes in y-direction
318  IF (mesh%JNUM.GT.1) THEN
319  i = mesh%IGMAX-mesh%IGMIN+1
320 !NEC$ SHORTLOOP
321  DO l=1,physics%VNUM
322  DO k=mesh%KGMIN,mesh%KGMAX
323 ! use collapsed arrays for better vectorization
324 !NEC$ IVDEP
325  DO j=i+1,SIZE(rvar%data3d,1)-i
326  this%slopes%data4d(j,k,m,l) = mesh%invdy * minmod3_limiter(&
327  this%limiter_param*(rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l)),&
328  this%limiter_param*(rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l)),&
329  0.5*(rvar%data3d(j+i,k,l) - rvar%data3d(j-i,k,l)))
330  END DO
331  END DO
332  END DO
333  m = m + 1
334  END IF
335  ! calculate slopes in z-direction
336  IF (mesh%KNUM.GT.1) THEN
337  i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
338 !NEC$ SHORTLOOP
339  DO l=1,physics%VNUM
340 ! use collapsed arrays for better vectorization
341 !NEC$ IVDEP
342  DO j=i+1,SIZE(rvar%data2d,1)-i
343  this%slopes%data3d(j,m,l) = mesh%invdz * minmod3_limiter(&
344  this%limiter_param*(rvar%data2d(j,l) - rvar%data2d(j-i,l)),&
345  this%limiter_param*(rvar%data2d(j+i,l) - rvar%data2d(j,l)),&
346  0.5*(rvar%data2d(j+i,l) - rvar%data2d(j-i,l)))
347  END DO
348  END DO
349  END IF
350  CASE(sweby)
351  ! calculate slopes in x-direction
352  IF (mesh%INUM.GT.1) THEN
353 !NEC$ SHORTLOOP
354  DO l=1,physics%VNUM
355  DO k=mesh%KGMIN,mesh%KGMAX
356  DO j=mesh%JGMIN,mesh%JGMAX
357 !NEC$ IVDEP
358  DO i=mesh%IGMIN+mesh%IP1,mesh%IGMAX+mesh%IM1
359  this%slopes%data5d(i,j,k,m,l) = mesh%invdx * sweby_limiter(&
360  rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
361  rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l),&
362  this%limiter_param)
363  END DO
364  END DO
365  END DO
366  END DO
367  m = m + 1
368  END IF
369  ! calculate slopes in y-direction
370  IF (mesh%JNUM.GT.1) THEN
371  i = mesh%IGMAX-mesh%IGMIN+1
372 !NEC$ SHORTLOOP
373  DO l=1,physics%VNUM
374  DO k=mesh%KGMIN,mesh%KGMAX
375 ! use collapsed arrays for better vectorization
376 !NEC$ IVDEP
377  DO j=i+1,SIZE(rvar%data3d,1)-i
378  this%slopes%data4d(j,k,m,l) = mesh%invdy * sweby_limiter(&
379  rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
380  rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l),&
381  this%limiter_param)
382  END DO
383  END DO
384  END DO
385  m = m + 1
386  END IF
387  ! calculate slopes in z-direction
388  IF (mesh%KNUM.GT.1) THEN
389  i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
390 !NEC$ SHORTLOOP
391  DO l=1,physics%VNUM
392 ! use collapsed arrays for better vectorization
393 !NEC$ IVDEP
394  DO j=i+1,SIZE(rvar%data2d,1)-i
395  this%slopes%data3d(j,m,l) = mesh%invdz * sweby_limiter(&
396  rvar%data2d(j,l) - rvar%data2d(j-i,l),&
397  rvar%data2d(j+i,l) - rvar%data2d(j,l),&
398  this%limiter_param)
399  END DO
400  END DO
401  END IF
402  CASE(superbee)
403  ! calculate slopes in x-direction
404  IF (mesh%INUM.GT.1) THEN
405 !NEC$ SHORTLOOP
406  DO l=1,physics%VNUM
407  DO k=mesh%KGMIN,mesh%KGMAX
408  DO j=mesh%JGMIN,mesh%JGMAX
409 !NEC$ IVDEP
410  DO i=mesh%IGMIN+mesh%IP1,mesh%IGMAX+mesh%IM1
411  this%slopes%data5d(i,j,k,m,l) = mesh%invdx * sweby_limiter(&
412  rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
413  rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l), 2.0)
414  END DO
415  END DO
416  END DO
417  END DO
418  m = m + 1
419  END IF
420  ! calculate slopes in y-direction
421  IF (mesh%JNUM.GT.1) THEN
422  i = mesh%IGMAX-mesh%IGMIN+1
423 !NEC$ SHORTLOOP
424  DO l=1,physics%VNUM
425  DO k=mesh%KGMIN,mesh%KGMAX
426 ! use collapsed arrays for better vectorization
427 !NEC$ IVDEP
428  DO j=i+1,SIZE(rvar%data3d,1)-i
429  this%slopes%data4d(j,k,m,l) = mesh%invdy * sweby_limiter(&
430  rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
431  rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l), 2.0)
432  END DO
433  END DO
434  END DO
435  m = m + 1
436  END IF
437  ! calculate slopes in z-direction
438  IF (mesh%KNUM.GT.1) THEN
439  i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
440 !NEC$ SHORTLOOP
441  DO l=1,physics%VNUM
442 ! use collapsed arrays for better vectorization
443 !NEC$ IVDEP
444  DO j=i+1,SIZE(rvar%data2d,1)-i
445  this%slopes%data3d(j,m,l) = mesh%invdz * sweby_limiter(&
446  rvar%data2d(j,l) - rvar%data2d(j-i,l),&
447  rvar%data2d(j+i,l) - rvar%data2d(j,l), 2.0)
448  END DO
449  END DO
450  END IF
451  CASE(ospre)
452  ! calculate slopes in x-direction
453  IF (mesh%INUM.GT.1) THEN
454 !NEC$ SHORTLOOP
455  DO l=1,physics%VNUM
456  DO k=mesh%KGMIN,mesh%KGMAX
457  DO j=mesh%JGMIN,mesh%JGMAX
458 !NEC$ IVDEP
459  DO i=mesh%IGMIN+mesh%IP1,mesh%IGMAX+mesh%IM1
460  this%slopes%data5d(i,j,k,m,l) = mesh%invdx * ospre_limiter(&
461  rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
462  rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
463  END DO
464  END DO
465  END DO
466  END DO
467  m = m + 1
468  END IF
469  ! calculate slopes in y-direction
470  IF (mesh%JNUM.GT.1) THEN
471  i = mesh%IGMAX-mesh%IGMIN+1
472 !NEC$ SHORTLOOP
473  DO l=1,physics%VNUM
474  DO k=mesh%KGMIN,mesh%KGMAX
475 ! use collapsed arrays for better vectorization
476 !NEC$ IVDEP
477  DO j=i+1,SIZE(rvar%data3d,1)-i
478  this%slopes%data4d(j,k,m,l) = mesh%invdy * ospre_limiter(&
479  rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
480  rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
481  END DO
482  END DO
483  END DO
484  m = m + 1
485  END IF
486  ! calculate slopes in z-direction
487  IF (mesh%KNUM.GT.1) THEN
488  i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
489 !NEC$ SHORTLOOP
490  DO l=1,physics%VNUM
491 ! use collapsed arrays for better vectorization
492 !NEC$ IVDEP
493  DO j=i+1,SIZE(rvar%data2d,1)-i
494  this%slopes%data3d(j,m,l) = mesh%invdz * ospre_limiter(&
495  rvar%data2d(j,l) - rvar%data2d(j-i,l),&
496  rvar%data2d(j+i,l) - rvar%data2d(j,l))
497  END DO
498  END DO
499  END IF
500 ! CASE(PP)
501 ! ! calculate slopes in both-directions
502 ! !NEC$ SHORTLOOP
503 ! DO l=1,Physics%VNUM
504 ! DO k=Mesh%KGMIN,Mesh%KGMAX
505 ! DO j=Mesh%JGMIN+Mesh%JP1,Mesh%JGMAX+Mesh%JM1
506 !!NEC$ IVDEP
507 ! DO i=Mesh%IGMIN+Mesh%IP1,Mesh%IGMAX+Mesh%IM1
508 ! CALL pp_limiter(this%slopes%data5d(i,j,k,m,l),&
509 ! this%slopes%data5d(i,j,k,m,l),&
510 ! rvar%data4d(i-1,j+1,k)-rvar%data4d(i,j,k),&
511 ! rvar%data4d(i ,j+1,k)-rvar%data4d(i,j,k),&
512 ! rvar%data4d(i+1,j+1,k)-rvar%data4d(i,j,k),&
513 ! rvar%data4d(i-1,j ,k)-rvar%data4d(i,j,k),&
514 ! this%limiter_param,&
515 ! rvar%data4d(i+1,j ,k)-rvar%data4d(i,j,k),&
516 ! rvar%data4d(i-1,j-1,k)-rvar%data4d(i,j,k),&
517 ! rvar%data4d(i ,j-1,k)-rvar%data4d(i,j,k),&
518 ! rvar%data4d(i+1,j-1,k)-rvar%data4d(i,j,k))
519 ! END DO
520 ! END DO
521 ! END DO
522 ! END DO
523  CASE(vanleer)
524  ! calculate slopes in x-direction
525  IF (mesh%INUM.GT.1) THEN
526 !NEC$ SHORTLOOP
527  DO l=1,physics%VNUM
528  DO k=mesh%KGMIN,mesh%KGMAX
529  DO j=mesh%JGMIN,mesh%JGMAX
530 !NEC$ IVDEP
531  DO i=mesh%IGMIN+mesh%IP1,mesh%IGMAX+mesh%IM1
532  this%slopes%data5d(i,j,k,m,l) = mesh%invdx * vanleer_limiter(&
533  rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
534  rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
535  END DO
536  END DO
537  END DO
538  END DO
539  m = m + 1
540  END IF
541  ! calculate slopes in y-direction
542  IF (mesh%JNUM.GT.1) THEN
543  i = mesh%IGMAX-mesh%IGMIN+1
544 !NEC$ SHORTLOOP
545  DO l=1,physics%VNUM
546  DO k=mesh%KGMIN,mesh%KGMAX
547 ! use collapsed arrays for better vectorization
548 !NEC$ IVDEP
549  DO j=i+1,SIZE(rvar%data3d,1)-i
550  this%slopes%data4d(j,k,m,l) = mesh%invdy * vanleer_limiter(&
551  rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
552  rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
553  END DO
554  END DO
555  END DO
556  m = m + 1
557  END IF
558  ! calculate slopes in z-direction
559  IF (mesh%KNUM.GT.1) THEN
560  i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
561 !NEC$ SHORTLOOP
562  DO l=1,physics%VNUM
563 ! use collapsed arrays for better vectorization
564 !NEC$ IVDEP
565  DO j=i+1,SIZE(rvar%data2d,1)-i
566  this%slopes%data3d(j,m,l) = mesh%invdz * vanleer_limiter(&
567  rvar%data2d(j,l) - rvar%data2d(j-i,l),&
568  rvar%data2d(j+i,l) - rvar%data2d(j,l))
569  END DO
570  END DO
571  END IF
572  CASE(nolimit)
573  ! calculate slopes in x-direction
574  IF (mesh%INUM.GT.1) THEN
575 !NEC$ SHORTLOOP
576  DO l=1,physics%VNUM
577  DO k=mesh%KGMIN,mesh%KGMAX
578  DO j=mesh%JGMIN,mesh%JGMAX
579 !NEC$ IVDEP
580  DO i=mesh%IGMIN+mesh%IP1,mesh%IGMAX+mesh%IM1
581  this%slopes%data5d(i,j,k,m,l) = mesh%invdx * nolimit_limiter(&
582  rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
583  rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
584  END DO
585  END DO
586  END DO
587  END DO
588  m = m + 1
589  END IF
590  ! calculate slopes in y-direction
591  IF (mesh%JNUM.GT.1) THEN
592  i = mesh%IGMAX-mesh%IGMIN+1
593 !NEC$ SHORTLOOP
594  DO l=1,physics%VNUM
595  DO k=mesh%KGMIN,mesh%KGMAX
596 ! use collapsed arrays for better vectorization
597 !NEC$ IVDEP
598  DO j=i+1,SIZE(rvar%data3d,1)-i
599  this%slopes%data4d(j,k,m,l) = mesh%invdy * nolimit_limiter(&
600  rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
601  rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
602  END DO
603  END DO
604  END DO
605  m = m + 1
606  END IF
607  ! calculate slopes in z-direction
608  IF (mesh%KNUM.GT.1) THEN
609  i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
610 !NEC$ SHORTLOOP
611  DO l=1,physics%VNUM
612 ! use collapsed arrays for better vectorization
613 !NEC$ IVDEP
614  DO j=i+1,SIZE(rvar%data2d,1)-i
615  this%slopes%data3d(j,m,l) = mesh%invdz * nolimit_limiter(&
616  rvar%data2d(j,l) - rvar%data2d(j-i,l),&
617  rvar%data2d(j+i,l) - rvar%data2d(j,l))
618  END DO
619  END DO
620  END IF
621  END SELECT
622  END SUBROUTINE calculateslopes
623 
624  SUBROUTINE finalize(this)
625  IMPLICIT NONE
626  !------------------------------------------------------------------------!
627  CLASS(reconstruction_linear), INTENT(INOUT) :: this
628  !------------------------------------------------------------------------!
629  CALL this%slopes%Destroy()
630  CALL this%dx%Destroy()
631  DEALLOCATE(this%slopes,this%dx)
632 
633  CALL this%Finalize_base()
634  END SUBROUTINE finalize
635 
636 !----------------------------------------------------------------------------!
638 
639  ELEMENTAL FUNCTION minmod2_limiter(arg1,arg2) RESULT(limarg)
640  IMPLICIT NONE
641  !--------------------------------------------------------------------!
642  REAL :: limarg
643  REAL, INTENT(IN) :: arg1, arg2
644  !--------------------------------------------------------------------!
645  IF (sign(1.0,arg1)*sign(1.0,arg2).GT.0) THEN
646  limarg = sign(min(abs(arg1),abs(arg2)),arg1)
647  ELSE
648  limarg = 0.
649  END IF
650  END FUNCTION minmod2_limiter
651 
652  ELEMENTAL FUNCTION minmod3_limiter(arg1,arg2,arg3) RESULT(limarg)
653  IMPLICIT NONE
654  !--------------------------------------------------------------------!
655  REAL :: limarg
656  REAL, INTENT(IN) :: arg1, arg2, arg3
657  !--------------------------------------------------------------------!
658  IF (((sign(1.0,arg1)*sign(1.0,arg2)).GT.0).AND.&
659  ((sign(1.0,arg2)*sign(1.0,arg3)).GT.0)) THEN
660  limarg = sign(min(abs(arg1),abs(arg2),abs(arg3)),arg1)
661  ELSE
662  limarg = 0.
663  END IF
664  END FUNCTION minmod3_limiter
665 
666  ELEMENTAL FUNCTION sweby_limiter(arg1,arg2,param) RESULT(limarg)
667  IMPLICIT NONE
668  !--------------------------------------------------------------------!
669  REAL :: limarg
670  REAL, INTENT(IN) :: arg1, arg2, param
671  !--------------------------------------------------------------------!
672  IF (sign(1.0,arg1)*sign(1.0,arg2).GT.0) THEN
673  limarg = sign(max(min(param*abs(arg1),abs(arg2)), &
674  min(abs(arg1),param*abs(arg2))),arg1)
675  ELSE
676  limarg = 0.
677  END IF
678  END FUNCTION sweby_limiter
679 
680  ELEMENTAL FUNCTION ospre_limiter(arg1,arg2) RESULT(limarg)
681  IMPLICIT NONE
682  !--------------------------------------------------------------------!
683  REAL :: limarg
684  REAL, INTENT(IN) :: arg1, arg2
685  !--------------------------------------------------------------------!
686  limarg = 1.5*arg1*arg2*(arg1 + arg2) / (arg1*(arg1 + 0.5*arg2) &
687  + arg2*(arg2 + 0.5*arg1) + tiny(1.0))
688  END FUNCTION ospre_limiter
689 
690 ! ELEMENTAL SUBROUTINE pp_limiter(xslope,yslope,arg1,arg2,arg3,arg4,param,arg6,arg7,arg8,arg9) ! TODO: 3D
691 ! IMPLICIT NONE
692 ! !--------------------------------------------------------------------!
693 ! REAL, INTENT(OUT):: xslope,yslope,zslope
694 ! REAL, INTENT(IN) :: arg1,arg2,arg3,arg4,param,arg6,arg7,arg8,arg9
695 ! !--------------------------------------------------------------------!
696 ! REAL :: Vmin,Vmax,V
697 ! !--------------------------------------------------------------------!
698 ! xslope = (arg6-arg4)*0.5
699 ! yslope = (arg2-arg8)*0.5 ! TODO: 3D add zslope
700 ! Vmin = MIN(arg1,arg2,arg3,arg4,-param,arg6,arg7,arg8,arg9)
701 ! Vmax = MAX(arg1,arg2,arg3,arg4,+param,arg6,arg7,arg8,arg9)
702 ! V = MIN(1.0,2.0*MIN(ABS(Vmin),ABS(Vmax))/(ABS(xslope)+ABS(yslope)))
703 ! xslope = V*xslope*Mesh%invdx
704 ! yslope = V*yslope*Mesh%invdy ! TODO: 3D add zslope
705 ! END SUBROUTINE pp_limiter
706 
707  ELEMENTAL FUNCTION vanleer_limiter(arg1,arg2) RESULT(limarg)
708  IMPLICIT NONE
709  !--------------------------------------------------------------------!
710  REAL :: limarg
711  REAL, INTENT(IN) :: arg1, arg2
712  !--------------------------------------------------------------------!
713  REAL :: a1,a2
714  !--------------------------------------------------------------------!
715  a1 = abs(arg1)
716  a2 = abs(arg2)
717  limarg = (arg1 * a2 + arg2 * a1)/(a1 + a2 + tiny(a1))
718  END FUNCTION vanleer_limiter
719 
720  ELEMENTAL FUNCTION nolimit_limiter(arg1,arg2) RESULT(limarg)
721  IMPLICIT NONE
722  !--------------------------------------------------------------------!
723  REAL :: limarg
724  REAL, INTENT(IN) :: arg1, arg2
725  !--------------------------------------------------------------------!
726  limarg = 0.5*(arg1+arg2)
727  END FUNCTION nolimit_limiter
728 
729 END MODULE reconstruction_linear_mod
mesh_base
logging_base_mod::finalize
subroutine finalize(this)
Destructor of common class.
Definition: logging_base.f90:155
reconstruction_linear_mod::nolimit
integer, parameter, public nolimit
Definition: reconstruction_linear.f90:78
reconstruction_linear_mod::limitertype_name
character(len=32), dimension(8), parameter limitertype_name
Definition: reconstruction_linear.f90:80
reconstruction_base_mod::linear
integer, parameter, public linear
Definition: reconstruction_base.f90:79
reconstruction_linear_mod::superbee
integer, parameter, public superbee
Definition: reconstruction_linear.f90:74
reconstruction_linear_mod::minmod2_limiter
elemental real function minmod2_limiter(arg1, arg2)
Definition: reconstruction_linear.f90:640
reconstruction_linear_mod::sweby_limiter
elemental real function sweby_limiter(arg1, arg2, param)
Definition: reconstruction_linear.f90:667
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
reconstruction_base_mod::reconstruction_base
Definition: reconstruction_base.f90:44
reconstruction_linear_mod::reconstruction_linear
Definition: reconstruction_linear.f90:62
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
reconstruction_linear_mod::vanleer_limiter
elemental real function vanleer_limiter(arg1, arg2)
Definition: reconstruction_linear.f90:708
reconstruction_linear_mod::calculatestates
pure subroutine calculatestates(this, Mesh, Physics, rvar, rstates)
Reconstructes states at cell boundaries.
Definition: reconstruction_linear.f90:210
reconstruction_linear_mod::recontype_name
character(len=32), parameter recontype_name
Definition: reconstruction_linear.f90:79
reconstruction_base_mod
base module for reconstruction process
Definition: reconstruction_base.f90:35
reconstruction_linear_mod::minmod
integer, parameter, public minmod
Definition: reconstruction_linear.f90:71
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:64
reconstruction_linear_mod::pp
integer, parameter, public pp
Definition: reconstruction_linear.f90:76
physics_base_mod::physics_base
named integer constants for flavour of state vectors
Definition: physics_base.f90:69
reconstruction_linear_mod::monocent
integer, parameter, public monocent
Definition: reconstruction_linear.f90:72
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
reconstruction_linear_mod::vanleer
integer, parameter, public vanleer
Definition: reconstruction_linear.f90:77
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::sweby
integer, parameter, public sweby
Definition: reconstruction_linear.f90:73
reconstruction_linear_mod
module for linear (first order) TVD reconstruction using slope limiters
Definition: reconstruction_linear.f90:51
reconstruction_linear_mod::ospre
integer, parameter, public ospre
Definition: reconstruction_linear.f90:75
reconstruction_linear_mod::ospre_limiter
elemental real function ospre_limiter(arg1, arg2)
Definition: reconstruction_linear.f90:681
reconstruction_linear_mod::nolimit_limiter
elemental real function nolimit_limiter(arg1, arg2)
Definition: reconstruction_linear.f90:721
reconstruction_linear_mod::calculateslopes
pure subroutine calculateslopes(this, Mesh, Physics, rvar)
computes the limited slopes
Definition: reconstruction_linear.f90:236
reconstruction_linear_mod::minmod3_limiter
elemental real function minmod3_limiter(arg1, arg2, arg3)
Definition: reconstruction_linear.f90:653
Generated by   doxygen 1.8.14