physics_eulerisotherm.f90
Go to the documentation of this file.
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: physics_eulerisotherm.f90 #
5 !# #
6 !# Copyright (C) 2007-2018 #
7 !# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8 !# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
9 !# Manuel Jung <mjung@astrophysik.uni-kiel.de> #
10 !# Lars Bösch <lboesch@astrophysik.uni-kiel.de> #
11 !# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
12 !# #
13 !# This program is free software; you can redistribute it and/or modify #
14 !# it under the terms of the GNU General Public License as published by #
15 !# the Free Software Foundation; either version 2 of the License, or (at #
16 !# your option) any later version. #
17 !# #
18 !# This program is distributed in the hope that it will be useful, but #
19 !# WITHOUT ANY WARRANTY; without even the implied warranty of #
20 !# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
21 !# NON INFRINGEMENT. See the GNU General Public License for more #
22 !# details. #
23 !# #
24 !# You should have received a copy of the GNU General Public License #
25 !# along with this program; if not, write to the Free Software #
26 !# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
27 !# #
28 !#############################################################################
32 !----------------------------------------------------------------------------!
43 !----------------------------------------------------------------------------!
44 MODULE physics_eulerisotherm_mod
45  USE logging_base_mod
46  USE physics_base_mod
47  USE mesh_base_mod
48  USE marray_base_mod
49  USE marray_compound_mod
50  USE common_dict
51  IMPLICIT NONE
52  !--------------------------------------------------------------------------!
53  PRIVATE
54  INTEGER, PARAMETER :: num_var = 3 ! number of variables !
55  CHARACTER(LEN=32), PARAMETER :: problem_name = "Euler isotherm"
56  !--------------------------------------------------------------------------!
57  TYPE, EXTENDS(physics_base) :: physics_eulerisotherm
60  CLASS(marray_base), ALLOCATABLE &
61  :: bccsound, & !< at cell bary centers
62  fcsound
65  REAL :: csiso
66  CONTAINS
67  PROCEDURE :: initphysics_eulerisotherm
68  PROCEDURE :: printconfiguration_eulerisotherm
69  PROCEDURE :: enableoutput
70  PROCEDURE :: new_statevector
71  !------Convert2Primitive-------!
72  PROCEDURE :: convert2primitive_all
73  PROCEDURE :: convert2primitive_subset
74  !------Convert2Conservative----!
75  PROCEDURE :: convert2conservative_all
76  PROCEDURE :: convert2conservative_subset
77  !------Soundspeed Routines-----!
78  PROCEDURE :: setsoundspeeds_center
79  PROCEDURE :: setsoundspeeds_faces
80  generic :: setsoundspeeds => setsoundspeeds_center, setsoundspeeds_faces
81  !------Wavespeed Routines------!
82  PROCEDURE :: calcwavespeeds_center
83  PROCEDURE :: calcwavespeeds_faces
84  !------Flux Routines-----------!
85  PROCEDURE :: calcfluxesx
86  PROCEDURE :: calcfluxesy
87  PROCEDURE :: calcfluxesz
88  !------Fargo Routines----------!
89  PROCEDURE :: addbackgroundvelocityx
90  PROCEDURE :: addbackgroundvelocityy
91  PROCEDURE :: addbackgroundvelocityz
92  PROCEDURE :: subtractbackgroundvelocityx
93  PROCEDURE :: subtractbackgroundvelocityy
94  PROCEDURE :: subtractbackgroundvelocityz
95  PROCEDURE :: addfargosources
96 
97  PROCEDURE :: geometricalsources
98  PROCEDURE :: externalsources
99  PROCEDURE :: viscositysources
100 
101  PROCEDURE :: calcstresses
102 ! PROCEDURE :: CalcIntermediateStateX_eulerisotherm ! for HLLC
103 ! PROCEDURE :: CalcIntermediateStateY_eulerisotherm ! for HLLC
104  PROCEDURE :: calculatecharsystemx ! for absorbing boundaries
105  PROCEDURE :: calculatecharsystemy ! for absorbing boundaries
106  PROCEDURE :: calculatecharsystemz ! for absorbing boundaries
107  PROCEDURE :: calculateboundarydatax ! for absorbing boundaries
108  PROCEDURE :: calculateboundarydatay ! for absorbing boundaries
109  PROCEDURE :: calculateboundarydataz ! for absorbing boundaries
110 ! PROCEDURE :: CalcPrim2RiemannX_eulerisotherm ! for farfield boundaries
111 ! PROCEDURE :: CalcPrim2RiemannY_eulerisotherm ! for farfield boundaries
112 ! PROCEDURE :: CalcRiemann2PrimX_eulerisotherm ! for farfield boundaries
113 ! PROCEDURE :: CalcRiemann2PrimY_eulerisotherm ! for farfield boundaries
114 ! PROCEDURE :: CalcRoeAverages_eulerisotherm ! for advanced wavespeeds
115 ! PROCEDURE :: ExternalSources_eulerisotherm
116  PROCEDURE :: reflectionmasks ! for reflecting boundaries
117  PROCEDURE :: axismasks
118 
119  PROCEDURE :: finalize
120  END TYPE
121  TYPE, EXTENDS(marray_compound) :: statevector_eulerisotherm
122  INTEGER :: flavour = undefined
123  TYPE(marray_base), POINTER &
124  :: density => null(), &
125  velocity => null(), &
126  momentum => null()
127  CONTAINS
128  PROCEDURE :: assignmarray_0
129  END TYPE
130 ! INTERFACE statevector_eulerisotherm
131 ! MODULE PROCEDURE CreateStateVector
132 ! END INTERFACE
133  INTERFACE setflux
134  MODULE PROCEDURE setflux1d, setflux2d, setflux3d
135  END INTERFACE
136  INTERFACE cons2prim
137  MODULE PROCEDURE cons2prim1d, cons2prim2d, cons2prim3d
138  END INTERFACE
139  INTERFACE prim2cons
140  MODULE PROCEDURE prim2cons1d, prim2cons2d, prim2cons3d
141  END INTERFACE
142  !--------------------------------------------------------------------------!
143  PUBLIC :: &
144  ! types
145  physics_eulerisotherm, &
146  statevector_eulerisotherm, &
147  createstatevector_eulerisotherm
148  !--------------------------------------------------------------------------!
149 
150 CONTAINS
151 
152 !----------------------------------------------------------------------------!
154 
159  SUBROUTINE initphysics_eulerisotherm(this,Mesh,config,IO)
160  IMPLICIT NONE
161  !------------------------------------------------------------------------!
162  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
163  CLASS(mesh_base), INTENT(IN) :: Mesh
164  TYPE(Dict_TYP), POINTER, INTENT(IN) :: config, IO
165  !------------------------------------------------------------------------!
166  INTEGER :: next_idx,err
167  !------------------------------------------------------------------------!
168  CALL this%InitPhysics(mesh,config,io,euler_isotherm,problem_name)
169 
170  ! set the total number of variables in a state vector
171  this%VNUM = this%VDIM + 1
172 
173  ! get isothermal sound speed from configuration
174  CALL getattr(config, "cs", this%csiso, 0.0)
175 
176  ! allocate memory for arrays used in eulerisotherm
177  ALLOCATE(this%pvarname(this%VNUM),this%cvarname(this%VNUM),this%bccsound, &
178  this%fcsound, &
179  stat = err)
180  IF (err.NE.0) &
181  CALL this%Error("InitPhysics_eulerisotherm", "Unable to allocate memory.")
182 
186  this%DENSITY = 1 ! mass density !
187  this%pvarname(this%DENSITY) = "density"
188  this%cvarname(this%DENSITY) = "density"
189  this%XVELOCITY = 2 ! x-velocity !
190  this%XMOMENTUM = 2 ! x-momentum !
191  IF (this%VDIM.GE.2) THEN
192  this%YVELOCITY = 3 ! y-velocity !
193  this%YMOMENTUM = 3 ! y-momentum !
194  ELSE
195  this%YVELOCITY = 0 ! no y-velocity !
196  this%YMOMENTUM = 0 ! no y-momentum !
197  END IF
198  IF (this%VDIM.EQ.3) THEN
199  this%ZVELOCITY = 4 ! z-velocity !
200  this%ZMOMENTUM = 4 ! z-momentum !
201  ELSE
202  this%ZVELOCITY = 0 ! no z-velocity !
203  this%ZMOMENTUM = 0 ! no z-momentum !
204  END IF
205  this%PRESSURE = 0 ! no pressure !
206  this%ENERGY = 0 ! no total energy !
207 
208  ! check which vector components are available and
209  ! set names shown in the data file
210  next_idx = 2
211  IF (btest(mesh%VECTOR_COMPONENTS,0)) THEN
212  this%pvarname(next_idx) = "xvelocity"
213  this%cvarname(next_idx) = "xmomentum"
214  next_idx = next_idx + 1
215  END IF
216  IF (btest(mesh%VECTOR_COMPONENTS,1)) THEN
217  this%pvarname(next_idx) = "yvelocity"
218  this%cvarname(next_idx) = "ymomentum"
219  next_idx = next_idx + 1
220  END IF
221  IF (btest(mesh%VECTOR_COMPONENTS,2)) THEN
222  this%pvarname(next_idx) = "zvelocity"
223  this%cvarname(next_idx) = "zmomentum"
224  END IF
225 
226  ! create new mesh array bccsound
227  this%bccsound = marray_base()
228  this%fcsound = marray_base(mesh%NFACES)
229 
230  IF(this%csiso.GT.0.) THEN
231  this%bccsound%data1d(:) = this%csiso
232  this%fcsound%data1d(:) = this%csiso
233  ELSE
234  this%bccsound%data1d(:) = 0.
235  this%fcsound%data1d(:) = 0.
236  END IF
237 
238  ! enable support for absorbing boundary conditions
239  this%supports_absorbing = .true.
240 
241  CALL this%EnableOutput(mesh,config,io)
242  END SUBROUTINE initphysics_eulerisotherm
243 
244  SUBROUTINE printconfiguration_eulerisotherm(this)
245  IMPLICIT NONE
246  !------------------------------------------------------------------------!
247  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
248  !------------------------------------------------------------------------!
249  CALL this%PrintConfiguration()
250  END SUBROUTINE printconfiguration_eulerisotherm
251 
253  SUBROUTINE enableoutput(this,Mesh,config,IO)
254  IMPLICIT NONE
255  !------------------------------------------------------------------------!
256  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
257  CLASS(mesh_base), INTENT(IN) :: Mesh
258  TYPE(Dict_TYP), POINTER, INTENT(IN) :: config, IO
259  !------------------------------------------------------------------------!
260  INTEGER :: valwrite
261  !------------------------------------------------------------------------!
262  ! check if output of sound speeds is requested
263  CALL getattr(config, "output/bccsound", valwrite, 0)
264  IF (valwrite .EQ. 1) &
265  CALL setattr(io, "bccsound",&
266  this%bccsound%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
267 
268  CALL getattr(config, "output/fcsound", valwrite, 0)
269  IF (valwrite .EQ. 1) &
270  CALL setattr(io, "fcsound",&
271  this%fcsound%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,:))
272  END SUBROUTINE enableoutput
273 
275  SUBROUTINE new_statevector(this,new_sv,flavour,num)
276  IMPLICIT NONE
277  !------------------------------------------------------------------------!
278  CLASS(physics_eulerisotherm), INTENT(IN) :: this
279  CLASS(marray_compound), POINTER :: new_sv
280  INTEGER, OPTIONAL, INTENT(IN) :: flavour,num
281  !------------------------------------------------------------------------!
282  ALLOCATE(statevector_eulerisotherm::new_sv)
283  SELECT TYPE(sv => new_sv)
284  TYPE IS (statevector_eulerisotherm)
285 ! sv = statevector_eulerisotherm(this,flavour,num)
286  CALL createstatevector_eulerisotherm(this,sv,flavour,num)
287  END SELECT
288  END SUBROUTINE new_statevector
289 
291  PURE SUBROUTINE setsoundspeeds_center(this,Mesh,bccsound)
292  IMPLICIT NONE
293  !------------------------------------------------------------------------!
294  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
295  CLASS(mesh_base), INTENT(IN) :: mesh
296  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
297  INTENT(IN) :: bccsound
298  !------------------------------------------------------------------------!
299  this%bccsound = bccsound
300  END SUBROUTINE setsoundspeeds_center
301 
303  PURE SUBROUTINE setsoundspeeds_faces(this,Mesh,fcsound)
304  IMPLICIT NONE
305  !------------------------------------------------------------------------!
306  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
307  CLASS(mesh_base), INTENT(IN) :: mesh
308  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Mesh%NFACES), &
309  INTENT(IN) :: fcsound
310  !------------------------------------------------------------------------!
311  this%fcsound = fcsound
312  END SUBROUTINE setsoundspeeds_faces
313 
315  PURE SUBROUTINE convert2primitive_all(this,cvar,pvar)
316  IMPLICIT NONE
317  !------------------------------------------------------------------------!
318  CLASS(physics_eulerisotherm), INTENT(IN) :: this
319  CLASS(marray_compound), INTENT(INOUT) :: cvar,pvar
320  !------------------------------------------------------------------------!
321  SELECT TYPE(c => cvar)
322  TYPE IS (statevector_eulerisotherm)
323  SELECT TYPE(p => pvar)
324  TYPE IS (statevector_eulerisotherm)
325  IF (c%flavour.EQ.conservative.AND.p%flavour.EQ.primitive) THEN
326  ! perform the transformation depending on dimensionality
327  SELECT CASE(this%VDIM)
328  CASE(1) ! 1D velocity / momentum
329  CALL cons2prim(c%density%data1d(:),c%momentum%data1d(:), &
330  p%density%data1d(:),p%velocity%data1d(:))
331  CASE(2) ! 2D velocity / momentum
332  CALL cons2prim(c%density%data1d(:),c%momentum%data2d(:,1), &
333  c%momentum%data2d(:,2),p%density%data1d(:), &
334  p%velocity%data2d(:,1),p%velocity%data2d(:,2))
335  CASE(3) ! 3D velocity / momentum
336  CALL cons2prim(c%density%data1d(:),c%momentum%data2d(:,1), &
337  c%momentum%data2d(:,2),c%momentum%data2d(:,3), &
338  p%density%data1d(:),p%velocity%data2d(:,1), &
339  p%velocity%data2d(:,2),p%velocity%data2d(:,3))
340  END SELECT
341  ELSE
342  ! do nothing
343  END IF
344  END SELECT
345  END SELECT
346  END SUBROUTINE convert2primitive_all
347 
349  PURE SUBROUTINE convert2primitive_subset(this,i1,i2,j1,j2,k1,k2,cvar,pvar)
350  IMPLICIT NONE
351  !------------------------------------------------------------------------!
352  CLASS(physics_eulerisotherm), INTENT(IN) :: this
353  INTEGER, INTENT(IN) :: i1,i2,j1,j2,k1,k2
354  CLASS(marray_compound), INTENT(INOUT) :: cvar,pvar
355  !------------------------------------------------------------------------!
356  SELECT TYPE(c => cvar)
357  TYPE IS (statevector_eulerisotherm)
358  SELECT TYPE(p => pvar)
359  TYPE IS (statevector_eulerisotherm)
360  IF (c%flavour.EQ.conservative.AND.p%flavour.EQ.primitive) THEN
361  SELECT CASE (c%density%RANK)
362  CASE(0) ! state vector contains cell center values
363  ! perform the transformation depending on dimensionality
364  SELECT CASE(this%VDIM)
365  CASE(1) ! 1D velocity / momentum
366  CALL cons2prim(c%density%data3d(i1:i2,j1:j2,k1:k2), &
367  c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
368  p%density%data3d(i1:i2,j1:j2,k1:k2), &
369  p%velocity%data4d(i1:i2,j1:j2,k1:k2,1))
370  CASE(2) ! 2D velocity / momentum
371  CALL cons2prim(c%density%data3d(i1:i2,j1:j2,k1:k2), &
372  c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
373  c%momentum%data4d(i1:i2,j1:j2,k1:k2,2), &
374  p%density%data3d(i1:i2,j1:j2,k1:k2), &
375  p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
376  p%velocity%data4d(i1:i2,j1:j2,k1:k2,2))
377  CASE(3) ! 3D velocity / momentum
378  CALL cons2prim(c%density%data3d(i1:i2,j1:j2,k1:k2), &
379  c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
380  c%momentum%data4d(i1:i2,j1:j2,k1:k2,2), &
381  c%momentum%data4d(i1:i2,j1:j2,k1:k2,3), &
382  p%density%data3d(i1:i2,j1:j2,k1:k2), &
383  p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
384  p%velocity%data4d(i1:i2,j1:j2,k1:k2,2), &
385  p%velocity%data4d(i1:i2,j1:j2,k1:k2,3))
386  END SELECT
387  CASE(1) ! state vector contains cell face / corner values
388  ! perform the transformation depending on dimensionality
389  SELECT CASE(this%VDIM)
390  CASE(1) ! 1D velocity / momentum
391  CALL cons2prim(c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
392  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
393  p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
394  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1))
395  CASE(2) ! 2D velocity / momentum
396  CALL cons2prim(c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
397  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
398  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,2), &
399  p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
400  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
401  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,2))
402  CASE(3) ! 3D velocity / momentum
403  CALL cons2prim(c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
404  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
405  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,2), &
406  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,3), &
407  p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
408  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
409  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,2), &
410  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,3))
411  END SELECT
412  CASE DEFAULT
413  ! do nothing
414  END SELECT
415  ELSE
416  ! do nothing
417  END IF
418  END SELECT
419  END SELECT
420  END SUBROUTINE convert2primitive_subset
421 
423  PURE SUBROUTINE convert2conservative_all(this,pvar,cvar)
424  IMPLICIT NONE
425  !------------------------------------------------------------------------!
426  CLASS(physics_eulerisotherm), INTENT(IN) :: this
427  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
428  !------------------------------------------------------------------------!
429  SELECT TYPE(p => pvar)
430  TYPE IS (statevector_eulerisotherm)
431  SELECT TYPE(c => cvar)
432  TYPE IS (statevector_eulerisotherm)
433  IF (p%flavour.EQ.primitive.AND.c%flavour.EQ.conservative) THEN
434  ! perform the transformation depending on dimensionality
435  SELECT CASE(this%VDIM)
436  CASE(1) ! 1D velocity / momentum
437  CALL prim2cons(p%density%data1d(:),p%velocity%data1d(:), &
438  c%density%data1d(:),c%momentum%data1d(:))
439  CASE(2) ! 2D velocity / momentum
440  CALL prim2cons(p%density%data1d(:),p%velocity%data2d(:,1), &
441  p%velocity%data2d(:,2),c%density%data1d(:), &
442  c%momentum%data2d(:,1),c%momentum%data2d(:,2))
443  CASE(3) ! 3D velocity / momentum
444  CALL prim2cons(p%density%data1d(:),p%velocity%data2d(:,1), &
445  p%velocity%data2d(:,2),p%velocity%data2d(:,3), &
446  c%density%data1d(:),c%momentum%data2d(:,1), &
447  c%momentum%data2d(:,2),c%momentum%data2d(:,3))
448  END SELECT
449  ELSE
450  ! do nothing
451  END IF
452  END SELECT
453  END SELECT
454  END SUBROUTINE convert2conservative_all
455 
457  PURE SUBROUTINE convert2conservative_subset(this,i1,i2,j1,j2,k1,k2,pvar,cvar)
458  IMPLICIT NONE
459  !------------------------------------------------------------------------!
460  CLASS(physics_eulerisotherm), INTENT(IN) :: this
461  INTEGER, INTENT(IN) :: i1,i2,j1,j2,k1,k2
462  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
463  !------------------------------------------------------------------------!
464  SELECT TYPE(p => pvar)
465  TYPE IS (statevector_eulerisotherm)
466  SELECT TYPE(c => cvar)
467  TYPE IS (statevector_eulerisotherm)
468  IF (p%flavour.EQ.primitive.AND.c%flavour.EQ.conservative) THEN
469  SELECT CASE (p%density%RANK)
470  CASE(0) ! state vector contains cell center values
471  ! perform the transformation depending on dimensionality
472  SELECT CASE(this%VDIM)
473  CASE(1) ! 1D velocity / momentum
474  CALL prim2cons(p%density%data3d(i1:i2,j1:j2,k1:k2), &
475  p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
476  c%density%data3d(i1:i2,j1:j2,k1:k2), &
477  c%momentum%data4d(i1:i2,j1:j2,k1:k2,1))
478  CASE(2) ! 2D velocity / momentum
479  CALL prim2cons(p%density%data3d(i1:i2,j1:j2,k1:k2), &
480  p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
481  p%velocity%data4d(i1:i2,j1:j2,k1:k2,2), &
482  c%density%data3d(i1:i2,j1:j2,k1:k2), &
483  c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
484  c%momentum%data4d(i1:i2,j1:j2,k1:k2,2))
485  CASE(3) ! 3D velocity / momentum
486  CALL prim2cons(p%density%data3d(i1:i2,j1:j2,k1:k2), &
487  p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
488  p%velocity%data4d(i1:i2,j1:j2,k1:k2,2), &
489  p%velocity%data4d(i1:i2,j1:j2,k1:k2,3), &
490  c%density%data3d(i1:i2,j1:j2,k1:k2), &
491  c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
492  c%momentum%data4d(i1:i2,j1:j2,k1:k2,2), &
493  c%momentum%data4d(i1:i2,j1:j2,k1:k2,3))
494  END SELECT
495  CASE(1) ! state vector contains cell face / corner values
496  ! perform the transformation depending on dimensionality
497  SELECT CASE(this%VDIM)
498  CASE(1) ! 1D velocity / momentum
499  CALL prim2cons(p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
500  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
501  c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
502  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1))
503  CASE(2) ! 2D velocity / momentum
504  CALL prim2cons(p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
505  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
506  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,2), &
507  c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
508  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
509  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,2))
510  CASE(3) ! 3D velocity / momentum
511  CALL prim2cons(p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
512  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
513  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,2), &
514  p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,3), &
515  c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
516  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
517  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,2), &
518  c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,3))
519  END SELECT
520  CASE DEFAULT
521  ! do nothing
522  END SELECT
523  ELSE
524  ! do nothing
525  END IF
526  END SELECT
527  END SELECT
528  END SUBROUTINE convert2conservative_subset
529 
531  PURE SUBROUTINE calcwavespeeds_center(this,Mesh,pvar,minwav,maxwav)
532  IMPLICIT NONE
533  !------------------------------------------------------------------------!
534  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
535  CLASS(mesh_base), INTENT(IN) :: mesh
536  CLASS(marray_compound), INTENT(INOUT) :: pvar
537  TYPE(marray_base), INTENT(INOUT) :: minwav,maxwav
538  !------------------------------------------------------------------------!
539  INTEGER :: m,n
540  !------------------------------------------------------------------------!
541  ! compute minimal and maximal wave speeds at cell centers
542  SELECT TYPE(p => pvar)
543  TYPE IS(statevector_eulerisotherm)
544 !NEC$ SHORTLOOP
545  m = 1
546  DO n=1,mesh%NDIMS
547  IF (mesh%ROTSYM.EQ.n) m = m + 1 ! skip this velocity
548  CALL setwavespeeds(this%bccsound%data1d(:),p%velocity%data2d(:,m),&
549  minwav%data2d(:,n),maxwav%data2d(:,n))
550  m = m + 1
551  END DO
552  END SELECT
553  END SUBROUTINE calcwavespeeds_center
554 
556  PURE SUBROUTINE calcwavespeeds_faces(this,Mesh,prim,cons,minwav,maxwav)
557  IMPLICIT NONE
558  !------------------------------------------------------------------------!
559  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
560  CLASS(mesh_base), INTENT(IN) :: mesh
561  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Mesh%NFACES,this%VNUM), &
562  INTENT(IN) :: prim,cons
563  TYPE(marray_base), INTENT(INOUT) :: minwav,maxwav
564  !------------------------------------------------------------------------!
565  INTEGER :: i,j,k,m,n
566  !------------------------------------------------------------------------!
567  m = this%XVELOCITY
568  n = 1
569  IF (mesh%INUM.GT.1) THEN
570  ! compute minimal and maximal western/eastern wave speeds
571  CALL setwavespeeds(this%fcsound%data4d(:,:,:,2*n-1),prim(:,:,:,2*n-1,m), &
572  this%tmp(:,:,:),this%tmp1(:,:,:))
573  CALL setwavespeeds(this%fcsound%data4d(:,:,:,2*n),prim(:,:,:,2*n,m), &
574  minwav%data4d(:,:,:,n),maxwav%data4d(:,:,:,n))
575  ! determine the minimum and maximum at cell interfaces
576  DO k=mesh%KGMIN,mesh%KGMAX
577  DO j=mesh%JGMIN,mesh%JGMAX
578 !NEC$ IVDEP
579  DO i=mesh%IMIN+mesh%IM1,mesh%IMAX
580  ! western & eastern interfaces
581  minwav%data4d(i,j,k,n) = min(0.0,this%tmp(i+mesh%IP1,j,k) ,minwav%data4d(i,j,k,n))
582  maxwav%data4d(i,j,k,n) = max(0.0,this%tmp1(i+mesh%IP1,j,k),maxwav%data4d(i,j,k,n))
583  END DO
584  END DO
585  END DO
586  n = n + 1
587  m = m + 1
588  ELSE
589  IF (mesh%ROTSYM.EQ.1) m = m + 1 ! increases the velocity index
590  END IF
591 
592  IF (mesh%JNUM.GT.1) THEN
593  ! compute minimal and maximal southern/northern wave speeds
594  CALL setwavespeeds(this%fcsound%data4d(:,:,:,2*n-1),prim(:,:,:,2*n-1,m), &
595  this%tmp(:,:,:),this%tmp1(:,:,:))
596  CALL setwavespeeds(this%fcsound%data4d(:,:,:,2*n),prim(:,:,:,2*n,m), &
597  minwav%data4d(:,:,:,n),maxwav%data4d(:,:,:,n))
598  ! determine the minimum and maximum at cell interfaces
599  DO k=mesh%KGMIN,mesh%KGMAX
600  DO j=mesh%JMIN+mesh%JM1,mesh%JMAX
601 !NEC$ IVDEP
602  DO i=mesh%IGMIN,mesh%IGMAX
603  ! southern & northern interfaces
604  minwav%data4d(i,j,k,n) = min(0.0,this%tmp(i,j+mesh%JP1,k),minwav%data4d(i,j,k,n))
605  maxwav%data4d(i,j,k,n) = max(0.0,this%tmp1(i,j+mesh%JP1,k),maxwav%data4d(i,j,k,n))
606  END DO
607  END DO
608  END DO
609  n = n + 1
610  m = m + 1
611  ELSE
612  IF (mesh%ROTSYM.EQ.2) m = m + 1 ! increases the velocity index
613  END IF
614 
615  IF (mesh%KNUM.GT.1) THEN
616  ! compute minimal and maximal lower/upper wave speeds
617  CALL setwavespeeds(this%fcsound%data4d(:,:,:,2*n-1),prim(:,:,:,2*n-1,m), &
618  this%tmp(:,:,:),this%tmp1(:,:,:))
619  CALL setwavespeeds(this%fcsound%data4d(:,:,:,2*n),prim(:,:,:,2*n,m), &
620  minwav%data4d(:,:,:,n),maxwav%data4d(:,:,:,n))
621  ! determine the minimum and maximum at cell interfaces
622  DO k=mesh%KMIN+mesh%KM1,mesh%KMAX
623  DO j=mesh%JGMIN,mesh%JGMAX
624 !NEC$ IVDEP
625  DO i=mesh%IGMIN,mesh%IGMAX
626  ! bottom & top interfaces
627  minwav%data4d(i,j,k,n) = min(0.0,this%tmp(i,j,k+mesh%KP1),minwav%data4d(i,j,k,n))
628  maxwav%data4d(i,j,k,n) = max(0.0,this%tmp1(i,j,k+mesh%KP1),maxwav%data4d(i,j,k,n))
629  END DO
630  END DO
631  END DO
632  END IF
633 
634  !!! THIS IS MOST PROBABLY BROKEN !!!
635  ! set minimal and maximal wave speeds at cell interfaces of neighboring cells
636 ! IF (this%advanced_wave_speeds) THEN
637 ! DO k=Mesh%KGMIN,Mesh%KGMAX
638 ! DO j=Mesh%JGMIN,Mesh%JGMAX
639 !!NEC$ IVDEP
640 ! DO i=Mesh%IMIN-1,Mesh%IMAX
641 ! ! western & eastern interfaces
642 ! ! get Roe averaged x-velocity
643 ! CALL SetRoeAverages(prim(i,j,k,2,this%DENSITY),prim(i+1,j,k,1,this%DENSITY), &
644 ! prim(i,j,k,2,this%XVELOCITY),prim(i+1,j,k,1,this%XVELOCITY), &
645 ! uRoe)
646 ! ! compute Roe averaged wave speeds
647 ! CALL SetWaveSpeeds(this%fcsound(i,j,k,2),uRoe,aminRoe,amaxRoe)
648 ! minwav(i,j,k,1) = MIN(aminRoe,this%tmp(i+1,j,k),minwav(i,j,k,1))
649 ! maxwav(i,j,k,1) = MAX(amaxRoe,this%tmp1(i+1,j,k),maxwav(i,j,k,1))
650 ! END DO
651 ! END DO
652 ! END DO
653 ! DO k=Mesh%KGMIN,Mesh%KGMAX
654 ! DO j=Mesh%JMIN-1,Mesh%JMAX
655 !!NEC$ IVDEP
656 ! DO i=Mesh%IGMIN,Mesh%IGMAX
657 ! ! southern & northern interfaces
658 ! ! get Roe averaged y-velocity
659 ! CALL SetRoeAverages(prim(i,j,k,4,this%DENSITY),prim(i,j+1,k,3,this%DENSITY), &
660 ! prim(i,j,k,4,this%YVELOCITY),prim(i,j+1,k,3,this%YVELOCITY), &
661 ! uRoe)
662 ! ! compute Roe averaged wave speeds
663 ! CALL SetWaveSpeeds(this%fcsound(i,j,k,4),uRoe,aminRoe,amaxRoe)
664 ! minwav(i,j,k,2) = MIN(aminRoe,this%tmp2(i,j+1,k),minwav(i,j,k,2))
665 ! maxwav(i,j,k,2) = MAX(amaxRoe,this%tmp3(i,j+1,k),maxwav(i,j,k,2))
666 ! END DO
667 ! END DO
668 ! END DO
669 ! DO k=Mesh%KGMIN-1,Mesh%KGMAX
670 ! DO j=Mesh%JMIN,Mesh%JMAX
671 !!NEC$ IVDEP
672 ! DO i=Mesh%IGMIN,Mesh%IGMAX
673 ! ! topper & bottomer interfaces
674 ! ! get Roe averaged z-velocity
675 ! CALL SetRoeAverages(prim(i,j,k,6,this%DENSITY),prim(i,j,k+1,5,this%DENSITY), &
676 ! prim(i,j,k,6,this%ZVELOCITY),prim(i,j,k+1,5,this%ZVELOCITY), &
677 ! uRoe)
678 ! ! compute Roe averaged wave speeds
679 ! CALL SetWaveSpeeds(this%fcsound(i,j,k,6),uRoe,aminRoe,amaxRoe)
680 ! minwav(i,j,k,3) = MIN(aminRoe,this%tmp3(i,j,k+Mesh%kp1),minwav(i,j,k,2))
681 ! maxwav(i,j,k,3) = MAX(amaxRoe,this%tmp4(i,j,k+Mesh%kp1),maxwav(i,j,k,2))
682 ! END DO
683 ! END DO
684 ! END DO
685 ! ELSE
686 ! END IF
687  END SUBROUTINE calcwavespeeds_faces
688 
690  PURE SUBROUTINE geometricalsources(this,Mesh,pvar,cvar,sterm)
691  IMPLICIT NONE
692  !------------------------------------------------------------------------!
693  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
694  CLASS(mesh_base), INTENT(IN) :: mesh
695  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar,sterm
696  !------------------------------------------------------------------------!
697  ! compute geometrical source only for non-cartesian mesh
698  IF (mesh%Geometry%GetType().NE.cartesian) THEN
699  SELECT TYPE(p => pvar)
700  TYPE IS(statevector_eulerisotherm)
701  SELECT TYPE(c => cvar)
702  TYPE IS(statevector_eulerisotherm)
703  SELECT TYPE(s => sterm)
704  TYPE IS(statevector_eulerisotherm)
705  ! no source terms
706  s%density%data1d(:) = 0.0
707  SELECT CASE(mesh%VECTOR_COMPONENTS)
708  CASE(vector_x) ! 1D momentum in x-direction
709  ! vy = vz = my = mz = 0
710  s%momentum%data2d(:,1) = getgeometricalsourcex( &
711  mesh%cxyx%data2d(:,2),mesh%cxzx%data2d(:,2), &
712  mesh%cyxy%data2d(:,2),mesh%czxz%data2d(:,2), &
713  p%velocity%data2d(:,1),0.0,0.0, &
714  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
715  0.0,0.0)
716  CASE(vector_y) ! 1D momentum in y-direction
717  ! vx = vz = mx = mz = 0
718  s%momentum%data2d(:,1) = getgeometricalsourcey( &
719  mesh%cxyx%data2d(:,2),mesh%cyxy%data2d(:,2), &
720  mesh%cyzy%data2d(:,2),mesh%czyz%data2d(:,2), &
721  0.0,p%velocity%data2d(:,1),0.0, &
722  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
723  0.0,0.0)
724  CASE(vector_z) ! 1D momentum in z-direction
725  ! vx = vy = mx = my = 0
726  s%momentum%data2d(:,1) = getgeometricalsourcez( &
727  mesh%cxzx%data2d(:,2),mesh%cyzy%data2d(:,2), &
728  mesh%czxz%data2d(:,2),mesh%czyz%data2d(:,2), &
729  0.0,0.0,p%velocity%data2d(:,1), &
730  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
731  0.0,0.0)
732  CASE(ior(vector_x,vector_y)) ! 2D momentum in x-y-plane
733  ! vz = mz = 0
734  ! x-momentum
735  s%momentum%data2d(:,1) = getgeometricalsourcex( &
736  mesh%cxyx%data2d(:,2),mesh%cxzx%data2d(:,2), &
737  mesh%cyxy%data2d(:,2),mesh%czxz%data2d(:,2), &
738  p%velocity%data2d(:,1),p%velocity%data2d(:,2),0.0, &
739  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
740  c%momentum%data2d(:,2),0.0)
741  ! y-momentum
742  s%momentum%data2d(:,2) = getgeometricalsourcey( &
743  mesh%cxyx%data2d(:,2),mesh%cyxy%data2d(:,2), &
744  mesh%cyzy%data2d(:,2),mesh%czyz%data2d(:,2), &
745  p%velocity%data2d(:,1),p%velocity%data2d(:,2),0.0, &
746  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
747  c%momentum%data2d(:,1),0.0)
748  CASE(ior(vector_x,vector_z)) ! 2D momentum in x-z-plane
749  ! vy = my = 0
750  ! x-momentum
751  s%momentum%data2d(:,1) = getgeometricalsourcex( &
752  mesh%cxyx%data2d(:,2),mesh%cxzx%data2d(:,2), &
753  mesh%cyxy%data2d(:,2),mesh%czxz%data2d(:,2), &
754  p%velocity%data2d(:,1),0.0,p%velocity%data2d(:,2), &
755  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
756  0.0,c%momentum%data2d(:,2))
757  ! z-momentum
758  s%momentum%data2d(:,2) = getgeometricalsourcez( &
759  mesh%cxzx%data2d(:,2),mesh%cyzy%data2d(:,2), &
760  mesh%czxz%data2d(:,2),mesh%czyz%data2d(:,2), &
761  p%velocity%data2d(:,1),0.0,p%velocity%data2d(:,2), &
762  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
763  c%momentum%data2d(:,1),0.0)
764  CASE(ior(vector_y,vector_z)) ! 2D momentum in y-z-plane
765  ! vx = mx = 0
766  ! y-momentum
767  s%momentum%data2d(:,1) = getgeometricalsourcey( &
768  mesh%cxyx%data2d(:,2),mesh%cyxy%data2d(:,2), &
769  mesh%cyzy%data2d(:,2),mesh%czyz%data2d(:,2), &
770  0.0,p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
771  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
772  0.0,c%momentum%data2d(:,2))
773  ! z-momentum
774  s%momentum%data2d(:,2) = getgeometricalsourcez( &
775  mesh%cxzx%data2d(:,2),mesh%cyzy%data2d(:,2), &
776  mesh%czxz%data2d(:,2),mesh%czyz%data2d(:,2), &
777  0.0,p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
778  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
779  0.0,c%momentum%data2d(:,1))
780  CASE(ior(ior(vector_x,vector_y),vector_z)) ! 3D momentum
781  ! x-momentum
782  s%momentum%data2d(:,1) = getgeometricalsourcex( &
783  mesh%cxyx%data2d(:,2),mesh%cxzx%data2d(:,2), &
784  mesh%cyxy%data2d(:,2),mesh%czxz%data2d(:,2), &
785  p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
786  p%velocity%data2d(:,3), &
787  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
788  c%momentum%data2d(:,2),c%momentum%data2d(:,3))
789  ! y-momentum
790  s%momentum%data2d(:,2) = getgeometricalsourcey( &
791  mesh%cxyx%data2d(:,2),mesh%cyxy%data2d(:,2), &
792  mesh%cyzy%data2d(:,2),mesh%czyz%data2d(:,2), &
793  p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
794  p%velocity%data2d(:,3), &
795  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
796  c%momentum%data2d(:,1),c%momentum%data2d(:,3))
797  ! z-momentum
798  s%momentum%data2d(:,3) = getgeometricalsourcez( &
799  mesh%cxzx%data2d(:,2),mesh%cyzy%data2d(:,2), &
800  mesh%czxz%data2d(:,2),mesh%czyz%data2d(:,2), &
801  p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
802  p%velocity%data2d(:,3), &
803  p%density%data1d(:)*this%bccsound%data1d(:)**2, &
804  c%momentum%data2d(:,1),c%momentum%data2d(:,2))
805  CASE DEFAULT
806  ! return NaN
807  s%momentum%data1d(:) = nan_default_real
808  END SELECT
809  END SELECT
810  END SELECT
811  END SELECT
812  ! reset ghost cell data
813  ! reset ghost cell data
814  IF (mesh%INUM.GT.1) THEN
815  sterm%data4d(mesh%IGMIN:mesh%IMIN+mesh%IM1,:,:,:) = 0.0
816  sterm%data4d(mesh%IMAX+mesh%IP1:mesh%IGMAX,:,:,:) = 0.0
817  END IF
818  IF (mesh%JNUM.GT.1) THEN
819  sterm%data4d(:,mesh%JGMIN:mesh%JMIN+mesh%JM1,:,:) = 0.0
820  sterm%data4d(:,mesh%JMAX+mesh%JP1:mesh%JGMAX,:,:) = 0.0
821  END IF
822  IF (mesh%KNUM.GT.1) THEN
823  ! collapse the first 2 dimensions
824  sterm%data3d(:,mesh%KGMIN:mesh%KMIN+mesh%KM1,:) = 0.0
825  sterm%data3d(:,mesh%KMAX+mesh%KP1:mesh%KGMAX,:) = 0.0
826  END IF
827  END IF
828  END SUBROUTINE geometricalsources
829 
831  PURE SUBROUTINE calcfluxesx(this,Mesh,nmin,nmax,prim,cons,xfluxes)
832  IMPLICIT NONE
833  !------------------------------------------------------------------------!
834  CLASS(physics_eulerisotherm), INTENT(IN) :: this
835  CLASS(mesh_base), INTENT(IN) :: mesh
836  INTEGER, INTENT(IN) :: nmin,nmax
837  CLASS(marray_compound), INTENT(INOUT) :: prim,cons,xfluxes
838  !------------------------------------------------------------------------!
839  SELECT TYPE(p => prim)
840  TYPE IS(statevector_eulerisotherm)
841  SELECT TYPE(c => cons)
842  TYPE IS(statevector_eulerisotherm)
843  SELECT TYPE(f => xfluxes)
844  TYPE IS(statevector_eulerisotherm)
845  SELECT CASE(this%VDIM)
846  CASE(1) ! 1D flux
847  CALL setflux(this%fcsound%data2d(:,nmin:nmax), &
848  p%density%data2d(:,nmin:nmax), &
849  p%velocity%data3d(:,nmin:nmax,1), &
850  c%momentum%data3d(:,nmin:nmax,1), &
851  f%density%data2d(:,nmin:nmax), &
852  f%momentum%data3d(:,nmin:nmax,1))
853  CASE(2) ! 2D flux
854  CALL setflux(this%fcsound%data2d(:,nmin:nmax), &
855  p%density%data2d(:,nmin:nmax), &
856  p%velocity%data3d(:,nmin:nmax,1), &
857  c%momentum%data3d(:,nmin:nmax,1), &
858  c%momentum%data3d(:,nmin:nmax,2), &
859  f%density%data2d(:,nmin:nmax), &
860  f%momentum%data3d(:,nmin:nmax,1), &
861  f%momentum%data3d(:,nmin:nmax,2))
862  CASE(3) ! 3D flux
863  CALL setflux(this%fcsound%data2d(:,nmin:nmax), &
864  p%density%data2d(:,nmin:nmax), &
865  p%velocity%data3d(:,nmin:nmax,1), &
866  c%momentum%data3d(:,nmin:nmax,1), &
867  c%momentum%data3d(:,nmin:nmax,2), &
868  c%momentum%data3d(:,nmin:nmax,3), &
869  f%density%data2d(:,nmin:nmax), &
870  f%momentum%data3d(:,nmin:nmax,1), &
871  f%momentum%data3d(:,nmin:nmax,2), &
872  f%momentum%data3d(:,nmin:nmax,3))
873  END SELECT
874  END SELECT
875  END SELECT
876  END SELECT
877  END SUBROUTINE calcfluxesx
878 
880  PURE SUBROUTINE calcfluxesy(this,Mesh,nmin,nmax,prim,cons,yfluxes)
881  IMPLICIT NONE
882  !------------------------------------------------------------------------!
883  CLASS(physics_eulerisotherm), INTENT(IN) :: this
884  CLASS(mesh_base), INTENT(IN) :: mesh
885  INTEGER, INTENT(IN) :: nmin,nmax
886  CLASS(marray_compound), INTENT(INOUT) :: prim,cons,yfluxes
887  !------------------------------------------------------------------------!
888  SELECT TYPE(p => prim)
889  TYPE IS(statevector_eulerisotherm)
890  SELECT TYPE(c => cons)
891  TYPE IS(statevector_eulerisotherm)
892  SELECT TYPE(f => yfluxes)
893  TYPE IS(statevector_eulerisotherm)
894  SELECT CASE(this%VDIM)
895  CASE(1) ! 1D flux
896  CALL setflux(this%fcsound%data2d(:,nmin:nmax), &
897  p%density%data2d(:,nmin:nmax), &
898  p%velocity%data3d(:,nmin:nmax,1), &
899  c%momentum%data3d(:,nmin:nmax,1), &
900  f%density%data2d(:,nmin:nmax), &
901  f%momentum%data3d(:,nmin:nmax,1))
902  CASE(2) ! 2D flux
903  CALL setflux(this%fcsound%data2d(:,nmin:nmax), &
904  p%density%data2d(:,nmin:nmax), &
905  p%velocity%data3d(:,nmin:nmax,2), &
906  c%momentum%data3d(:,nmin:nmax,2), &
907  c%momentum%data3d(:,nmin:nmax,1), &
908  f%density%data2d(:,nmin:nmax), &
909  f%momentum%data3d(:,nmin:nmax,2), &
910  f%momentum%data3d(:,nmin:nmax,1))
911  CASE(3) ! 3D flux
912  CALL setflux(this%fcsound%data2d(:,nmin:nmax), &
913  p%density%data2d(:,nmin:nmax), &
914  p%velocity%data3d(:,nmin:nmax,2), &
915  c%momentum%data3d(:,nmin:nmax,2), &
916  c%momentum%data3d(:,nmin:nmax,1), &
917  c%momentum%data3d(:,nmin:nmax,3), &
918  f%density%data2d(:,nmin:nmax), &
919  f%momentum%data3d(:,nmin:nmax,2), &
920  f%momentum%data3d(:,nmin:nmax,1), &
921  f%momentum%data3d(:,nmin:nmax,3))
922  END SELECT
923  END SELECT
924  END SELECT
925  END SELECT
926  END SUBROUTINE calcfluxesy
927 
929  PURE SUBROUTINE calcfluxesz(this,Mesh,nmin,nmax,prim,cons,zfluxes)
930  IMPLICIT NONE
931  !------------------------------------------------------------------------!
932  CLASS(physics_eulerisotherm), INTENT(IN) :: this
933  CLASS(mesh_base), INTENT(IN) :: mesh
934  INTEGER, INTENT(IN) :: nmin,nmax
935  CLASS(marray_compound), INTENT(INOUT) :: prim,cons,zfluxes
936  !------------------------------------------------------------------------!
937  SELECT TYPE(p => prim)
938  TYPE IS(statevector_eulerisotherm)
939  SELECT TYPE(c => cons)
940  TYPE IS(statevector_eulerisotherm)
941  SELECT TYPE(f => zfluxes)
942  TYPE IS(statevector_eulerisotherm)
943  SELECT CASE(this%VDIM)
944  CASE(1) ! 1D flux
945  CALL setflux(this%fcsound%data2d(:,nmin:nmax), &
946  p%density%data2d(:,nmin:nmax), &
947  p%velocity%data3d(:,nmin:nmax,1), &
948  c%momentum%data3d(:,nmin:nmax,1), &
949  f%density%data2d(:,nmin:nmax), &
950  f%momentum%data3d(:,nmin:nmax,1))
951  CASE(2) ! 2D flux
952  CALL setflux(this%fcsound%data2d(:,nmin:nmax), &
953  p%density%data2d(:,nmin:nmax), &
954  p%velocity%data3d(:,nmin:nmax,2), &
955  c%momentum%data3d(:,nmin:nmax,2), &
956  c%momentum%data3d(:,nmin:nmax,1), &
957  f%density%data2d(:,nmin:nmax), &
958  f%momentum%data3d(:,nmin:nmax,2), &
959  f%momentum%data3d(:,nmin:nmax,1))
960  CASE(3) ! 3D flux
961  CALL setflux(this%fcsound%data2d(:,nmin:nmax), &
962  p%density%data2d(:,nmin:nmax), &
963  p%velocity%data3d(:,nmin:nmax,3), &
964  c%momentum%data3d(:,nmin:nmax,3), &
965  c%momentum%data3d(:,nmin:nmax,1), &
966  c%momentum%data3d(:,nmin:nmax,2), &
967  f%density%data2d(:,nmin:nmax), &
968  f%momentum%data3d(:,nmin:nmax,3), &
969  f%momentum%data3d(:,nmin:nmax,1), &
970  f%momentum%data3d(:,nmin:nmax,2))
971  END SELECT
972  END SELECT
973  END SELECT
974  END SELECT
975  END SUBROUTINE calcfluxesz
976 
978  PURE SUBROUTINE viscositysources(this,Mesh,pvar,btxx,btxy,btxz,btyy,btyz,btzz,sterm)
979  IMPLICIT NONE
980  !------------------------------------------------------------------------!
981  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
982  CLASS(mesh_base), INTENT(IN) :: mesh
983  CLASS(marray_compound), INTENT(INOUT) :: pvar,sterm
984  REAL, INTENT(IN), &
985  DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
986  :: btxx,btxy,btxz,btyy,btyz,btzz
987  !------------------------------------------------------------------------!
988  SELECT TYPE(p => pvar)
989  CLASS IS(statevector_eulerisotherm)
990  SELECT TYPE(s => sterm)
991  CLASS IS(statevector_eulerisotherm)
992  ! no viscous sources in continuity equation
993  s%density%data1d(:) = 0.0
994  ! viscous momentum sources
995  SELECT CASE(this%VDIM)
996 ! CASE(1) ! 1D velocities are currently not supported
997 ! CALL Mesh%Divergence(btxx,sterm(:,:,:,this%XMOMENTUM))
998  CASE(2)
999  ! divergence of stress tensor with symmetry btyx=btxy
1000  CALL mesh%Divergence(btxx,btxy,btxy,btyy,s%momentum%data4d(:,:,:,1), &
1001  s%momentum%data4d(:,:,:,2))
1002  CASE(3)
1003  ! divergence of stress tensor with symmetry btyx=btxy, btxz=btzx, btyz=btzy
1004  CALL mesh%Divergence(btxx,btxy,btxz,btxy,btyy,btyz,btxz,btyz,btzz, &
1005  s%momentum%data4d(:,:,:,1),s%momentum%data4d(:,:,:,2), &
1006  s%momentum%data4d(:,:,:,3))
1007  CASE DEFAULT
1008  ! return NaN
1009  s%data1d(:) = nan_default_real
1010  END SELECT
1011  END SELECT
1012  END SELECT
1013  END SUBROUTINE viscositysources
1014 
1019  PURE SUBROUTINE calcstresses(this,Mesh,pvar,dynvis,bulkvis, &
1020  btxx,btxy,btxz,btyy,btyz,btzz)
1021  IMPLICIT NONE
1022  !------------------------------------------------------------------------!
1023  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
1024  CLASS(mesh_base), INTENT(IN) :: mesh
1025  CLASS(marray_compound), INTENT(INOUT) :: pvar
1026  CLASS(marray_base), INTENT(INOUT) :: dynvis,bulkvis
1027  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: &
1028  btxx,btxy,btxz,btyy,btyz,btzz
1029  !------------------------------------------------------------------------!
1030  INTEGER :: i,j,k
1031  !------------------------------------------------------------------------!
1032  INTENT(OUT) :: btxx,btxy,btxz,btyy,btyz,btzz
1033  !------------------------------------------------------------------------!
1034  SELECT TYPE(p => pvar)
1035  CLASS IS(statevector_eulerisotherm)
1036  SELECT CASE(mesh%VECTOR_COMPONENTS)
1037 !!!! 1D velocities are currently not supported !!!!
1038 ! CASE(VECTOR_X) ! 1D velocity in x-direction
1039 ! CASE(VECTOR_Y) ! 1D velocity in y-direction
1040 ! CASE(VECTOR_Z) ! 1D velocity in z-direction
1041  CASE(ior(vector_x,vector_y)) ! 2D velocities in x-y-plane
1042  ! compute bulk viscosity first and store the result in this%tmp
1043  CALL mesh%Divergence(p%velocity%data4d(:,:,:,1),p%velocity%data4d(:,:,:,2),&
1044  this%tmp(:,:,:))
1045  this%tmp(:,:,:) = bulkvis%data3d(:,:,:)*this%tmp(:,:,:)
1046  DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
1047  DO j=mesh%JMIN-mesh%JP1,mesh%JMAX+mesh%JP1
1048 !NEC$ IVDEP
1049  DO i=mesh%IMIN-mesh%IP1,mesh%IMAX+mesh%IP1
1050  ! compute the diagonal elements of the stress tensor
1051  btxx(i,j,k) = dynvis%data3d(i,j,k) * &
1052  ((p%velocity%data4d(i+mesh%IP1,j,k,1) - p%velocity%data4d(i-mesh%IP1,j,k,1)) &
1053  / mesh%dlx%data3d(i,j,k) &
1054  + 2.0 * mesh%cxyx%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2)) + this%tmp(i,j,k)
1055  btyy(i,j,k) = dynvis%data3d(i,j,k) * &
1056  ((p%velocity%data4d(i,j+mesh%JP1,k,2) - p%velocity%data4d(i,j-mesh%JP1,k,2)) &
1057  / mesh%dly%data3d(i,j,k) &
1058  + 2.0 * mesh%cyxy%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) ) + this%tmp(i,j,k)
1059  ! compute the off-diagonal elements (no bulk viscosity)
1060  btxy(i,j,k) = dynvis%data3d(i,j,k) * ( 0.5 * &
1061  ((p%velocity%data4d(i+mesh%IP1,j,k,2) - p%velocity%data4d(i-mesh%IP1,j,k,2)) &
1062  / mesh%dlx%data3d(i,j,k) &
1063  + (p%velocity%data4d(i,j+mesh%JP1,k,1) - p%velocity%data4d(i,j-mesh%JP1,k,1)) &
1064  / mesh%dly%data3d(i,j,k) ) &
1065  - mesh%cxyx%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) &
1066  - mesh%cyxy%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2) )
1067  END DO
1068  END DO
1069  END DO
1070  CASE(ior(vector_x,vector_z)) ! 2D velocities in x-z-plane
1071  ! compute bulk viscosity first and store the result in this%tmp
1072  CALL mesh%Divergence(p%velocity%data4d(:,:,:,1),p%velocity%data4d(:,:,:,2),&
1073  this%tmp(:,:,:))
1074  this%tmp(:,:,:) = bulkvis%data3d(:,:,:)*this%tmp(:,:,:)
1075  DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
1076  DO j=mesh%JMIN-mesh%JP1,mesh%JMAX+mesh%JP1
1077 !NEC$ IVDEP
1078  DO i=mesh%IMIN-mesh%IP1,mesh%IMAX+mesh%IP1
1079  ! compute the diagonal elements btxx, btzz => btyy of the stress tensor
1080  btxx(i,j,k) = dynvis%data3d(i,j,k) * &
1081  ((p%velocity%data4d(i+mesh%IP1,j,k,1) - p%velocity%data4d(i-mesh%IP1,j,k,1)) &
1082  / mesh%dlx%data3d(i,j,k) &
1083  + 2.0 * mesh%cxzx%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2) ) & ! vz => vy
1084  + this%tmp(i,j,k)
1085  btyy(i,j,k) = dynvis%data3d(i,j,k) * & ! btzz => btyy
1086  ((p%velocity%data4d(i,j,k+mesh%KP1,2) - p%velocity%data4d(i,j,k-mesh%KP1,2)) & ! vz => vy
1087  / mesh%dlz%data3d(i,j,k) &
1088  + 2.0 * mesh%czxz%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) ) &
1089  + this%tmp(i,j,k)
1090  ! compute the off-diagonal elements btxz => btxy (no bulk viscosity)
1091  btxy(i,j,k) = dynvis%data3d(i,j,k) * ( 0.5 * & ! btxz => btxy
1092  ((p%velocity%data4d(i+mesh%IP1,j,k,2) - p%velocity%data4d(i-mesh%IP1,j,k,2)) & ! vz => vy
1093  / mesh%dlx%data3d(i,j,k) &
1094  + (p%velocity%data4d(i,j,k+mesh%KP1,1) - p%velocity%data4d(i,j,k-mesh%KP1,1)) &
1095  / mesh%dlz%data3d(i,j,k) ) &
1096  - mesh%czxz%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2) & ! vz => vy
1097  - mesh%cxzx%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) )
1098  END DO
1099  END DO
1100  END DO
1101  CASE(ior(vector_y,vector_z)) ! 2D velocities in y-z-plane
1102  CALL mesh%Divergence(p%velocity%data4d(:,:,:,1),p%velocity%data4d(:,:,:,2),&
1103  this%tmp(:,:,:))
1104  this%tmp(:,:,:) = bulkvis%data3d(:,:,:)*this%tmp(:,:,:)
1105  DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
1106  DO j=mesh%JMIN-mesh%JP1,mesh%JMAX+mesh%JP1
1107 !NEC$ IVDEP
1108  DO i=mesh%IMIN-mesh%IP1,mesh%IMAX+mesh%IP1
1109  ! compute the diagonal elements btyy => btxx, btzz => btyy of the stress tensor
1110  btxx(i,j,k) = dynvis%data3d(i,j,k) * &
1111  ((p%velocity%data4d(i,j+mesh%JP1,k,1) - p%velocity%data4d(i,j-mesh%JP1,k,1)) &
1112  / mesh%dly%data3d(i,j,k) &
1113  + 2.0 * mesh%cyzy%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2) ) & ! vz => vy
1114  + this%tmp(i,j,k)
1115  btyy(i,j,k) = dynvis%data3d(i,j,k) * &
1116  ((p%velocity%data4d(i,j,k+mesh%KP1,2) - p%velocity%data4d(i,j,k-mesh%KP1,2)) & ! vz => vy
1117  / mesh%dlz%data3d(i,j,k) &
1118  + 2.0 * mesh%czyz%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) ) &
1119  + this%tmp(i,j,k)
1120  ! compute the off-diagonal elements btyz => btxy (no bulk viscosity)
1121  btxy(i,j,k) = dynvis%data3d(i,j,k) * ( 0.5 * &
1122  ((p%velocity%data4d(i,j,k+mesh%KP1,1) - p%velocity%data4d(i,j,k-mesh%KP1,1)) & ! vy => vx
1123  / mesh%dlz%data3d(i,j,k) &
1124  + (p%velocity%data4d(i,j+mesh%JP1,k,2) - p%velocity%data4d(i,j-mesh%JP1,k,2)) & ! vz => vy
1125  / mesh%dly%data3d(i,j,k) ) &
1126  - mesh%czyz%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2) & ! vz => vy
1127  - mesh%cyzy%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) ) ! vy => vx
1128  END DO
1129  END DO
1130  END DO
1131  CASE(ior(ior(vector_x,vector_y),vector_z)) ! 3D velocities
1132  ! compute bulk viscosity first and store the result in this%tmp
1133  CALL mesh%Divergence(p%velocity%data4d(:,:,:,1),p%velocity%data4d(:,:,:,2),&
1134  p%velocity%data4d(:,:,:,3),this%tmp(:,:,:))
1135  this%tmp(:,:,:) = bulkvis%data3d(:,:,:)*this%tmp(:,:,:)
1136  DO k=mesh%KMIN-mesh%KP1,mesh%KMAX+mesh%KP1
1137  DO j=mesh%JMIN-mesh%JP1,mesh%JMAX+mesh%JP1
1138 !NEC$ IVDEP
1139  DO i=mesh%IMIN-mesh%IP1,mesh%IMAX+mesh%IP1
1140  ! compute the diagonal elements of the stress tensor
1141  btxx(i,j,k) = dynvis%data3d(i,j,k) * &
1142  ((p%velocity%data4d(i+mesh%IP1,j,k,1) - p%velocity%data4d(i-mesh%IP1,j,k,1)) &
1143  / mesh%dlx%data3d(i,j,k) &
1144  + 2.0 * mesh%cxyx%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2) &
1145  + 2.0 * mesh%cxzx%bcenter(i,j,k) * p%velocity%data4d(i,j,k,3) ) &
1146  + this%tmp(i,j,k)
1147  btyy(i,j,k) = dynvis%data3d(i,j,k) * &
1148  ((p%velocity%data4d(i,j+mesh%JP1,k,2) - p%velocity%data4d(i,j-mesh%JP1,k,2)) &
1149  / mesh%dly%data3d(i,j,k) &
1150  + 2.0 * mesh%cyxy%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) &
1151  + 2.0 * mesh%cyzy%bcenter(i,j,k) * p%velocity%data4d(i,j,k,3) ) &
1152  + this%tmp(i,j,k)
1153  btzz(i,j,k) = dynvis%data3d(i,j,k) * &
1154  ((p%velocity%data4d(i,j,k+mesh%KP1,3) - p%velocity%data4d(i,j,k-mesh%KP1,3)) &
1155  / mesh%dlz%data3d(i,j,k) &
1156  + 2.0 * mesh%czxz%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) &
1157  + 2.0 * mesh%czyz%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2) ) &
1158  + this%tmp(i,j,k)
1159  ! compute the off-diagonal elements (no bulk viscosity)
1160  btxy(i,j,k) = dynvis%data3d(i,j,k) * ( 0.5 * &
1161  ((p%velocity%data4d(i+mesh%IP1,j,k,2) - p%velocity%data4d(i-mesh%IP1,j,k,2)) &
1162  / mesh%dlx%data3d(i,j,k) &
1163  + (p%velocity%data4d(i,j+mesh%JP1,k,1) - p%velocity%data4d(i,j-mesh%JP1,k,1)) &
1164  / mesh%dly%data3d(i,j,k) ) &
1165  - mesh%cxyx%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) &
1166  - mesh%cyxy%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2) )
1167  btxz(i,j,k) = dynvis%data3d(i,j,k) * ( 0.5 * &
1168  ((p%velocity%data4d(i+mesh%IP1,j,k,3) - p%velocity%data4d(i-mesh%IP1,j,k,3)) &
1169  / mesh%dlx%data3d(i,j,k) &
1170  + (p%velocity%data4d(i,j,k+mesh%KP1,1) - p%velocity%data4d(i,j,k-mesh%KP1,1)) &
1171  / mesh%dlz%data3d(i,j,k) ) &
1172  - mesh%czxz%bcenter(i,j,k) * p%velocity%data4d(i,j,k,3) &
1173  - mesh%cxzx%bcenter(i,j,k) * p%velocity%data4d(i,j,k,1) )
1174  btyz(i,j,k) = dynvis%data3d(i,j,k) * ( 0.5 * &
1175  ((p%velocity%data4d(i,j,k+mesh%KP1,2) - p%velocity%data4d(i,j,k-mesh%KP1,2)) &
1176  / mesh%dlz%data3d(i,j,k) &
1177  + (p%velocity%data4d(i,j+mesh%JP1,k,3) - p%velocity%data4d(i,j-mesh%JP1,k,3)) &
1178  / mesh%dly%data3d(i,j,k) ) &
1179  - mesh%czyz%bcenter(i,j,k) * p%velocity%data4d(i,j,k,3) &
1180  - mesh%cyzy%bcenter(i,j,k) * p%velocity%data4d(i,j,k,2) )
1181  END DO
1182  END DO
1183  END DO
1184  CASE DEFAULT
1185  ! return NaN
1186  btxx(:,:,:) = nan_default_real
1187  btxy(:,:,:) = nan_default_real
1188  btxz(:,:,:) = nan_default_real
1189  btyy(:,:,:) = nan_default_real
1190  btyz(:,:,:) = nan_default_real
1191  btzz(:,:,:) = nan_default_real
1192  END SELECT
1193  END SELECT
1194  END SUBROUTINE calcstresses
1195 
1197  PURE SUBROUTINE externalsources(this,accel,pvar,cvar,sterm)
1198  !------------------------------------------------------------------------!
1199  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1200  CLASS(marray_base), INTENT(IN) :: accel
1201  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar,sterm
1202  !------------------------------------------------------------------------!
1203  INTEGER :: n
1204  !------------------------------------------------------------------------!
1205  SELECT TYPE(p => pvar)
1206  CLASS IS(statevector_eulerisotherm)
1207  SELECT TYPE(c => cvar)
1208  CLASS IS(statevector_eulerisotherm)
1209  SELECT TYPE(s => sterm)
1210  CLASS IS(statevector_eulerisotherm)
1211  s%density%data1d(:) = 0.0
1212 !NEC$ UNROLL(3)
1213  DO n=1,this%VDIM
1214  s%momentum%data2d(:,n) = c%density%data1d(:) * accel%data2d(:,n)
1215  END DO
1216  END SELECT
1217  END SELECT
1218  END SELECT
1219  END SUBROUTINE externalsources
1220 
1221 
1222  PURE SUBROUTINE addbackgroundvelocityx(this,Mesh,w,pvar,cvar)
1223  IMPLICIT NONE
1224  !------------------------------------------------------------------------!
1225  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
1226  CLASS(mesh_base), INTENT(IN) :: mesh
1227  !------------------------------------------------------------------------!
1228  REAL,DIMENSION(Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
1229  INTENT(IN) :: w
1230  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
1231  !------------------------------------------------------------------------!
1232  INTEGER :: i,j,k
1233  !------------------------------------------------------------------------!
1234  IF (this%transformed_xvelocity) THEN
1235  SELECT TYPE(p => pvar)
1236  TYPE IS(statevector_eulerisotherm)
1237  SELECT TYPE(c => cvar)
1238  TYPE IS(statevector_eulerisotherm)
1239  DO k=mesh%KGMIN,mesh%KGMAX
1240  DO j=mesh%JGMIN,mesh%JGMAX
1241  DO i=mesh%IGMIN,mesh%IGMAX
1242  p%velocity%data4d(i,j,k,1) = p%velocity%data4d(i,j,k,1) + w(j,k)
1243  c%momentum%data4d(i,j,k,1) = c%momentum%data4d(i,j,k,1) &
1244  + c%density%data3d(i,j,k)*w(j,k)
1245  END DO
1246  END DO
1247  END DO
1248  this%transformed_xvelocity = .false.
1249  END SELECT
1250  END SELECT
1251  this%transformed_xvelocity = .false.
1252  END IF
1253  END SUBROUTINE addbackgroundvelocityx
1254 
1255  PURE SUBROUTINE addbackgroundvelocityy(this,Mesh,w,pvar,cvar)
1256  IMPLICIT NONE
1257  !------------------------------------------------------------------------!
1258  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
1259  CLASS(mesh_base), INTENT(IN) :: mesh
1260  !------------------------------------------------------------------------!
1261  REAL,DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%KGMIN:Mesh%KGMAX), &
1262  INTENT(IN) :: w
1263  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
1264  !------------------------------------------------------------------------!
1265  INTEGER :: i,j,k
1266  !------------------------------------------------------------------------!
1267  IF (this%transformed_yvelocity) THEN
1268  SELECT TYPE(p => pvar)
1269  TYPE IS(statevector_eulerisotherm)
1270  SELECT TYPE(c => cvar)
1271  TYPE IS(statevector_eulerisotherm)
1272  DO k=mesh%KGMIN,mesh%KGMAX
1273  DO j=mesh%JGMIN,mesh%JGMAX
1274  DO i=mesh%IGMIN,mesh%IGMAX
1275  p%velocity%data4d(i,j,k,2) = p%velocity%data4d(i,j,k,2) + w(i,k)
1276  c%momentum%data4d(i,j,k,2) = c%momentum%data4d(i,j,k,2) &
1277  + c%density%data3d(i,j,k)*w(i,k)
1278  END DO
1279  END DO
1280  END DO
1281  this%transformed_yvelocity = .false.
1282  END SELECT
1283  END SELECT
1284  END IF
1285  END SUBROUTINE addbackgroundvelocityy
1286 
1287  PURE SUBROUTINE addbackgroundvelocityz(this,Mesh,w,pvar,cvar)
1288  IMPLICIT NONE
1289  !------------------------------------------------------------------------!
1290  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
1291  CLASS(mesh_base), INTENT(IN) :: mesh
1292  !------------------------------------------------------------------------!
1293  REAL,DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX), &
1294  INTENT(IN) :: w
1295  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
1296  !------------------------------------------------------------------------!
1297  INTEGER :: i,j,k
1298  !------------------------------------------------------------------------!
1299  IF (this%transformed_zvelocity) THEN
1300  SELECT TYPE(p => pvar)
1301  TYPE IS(statevector_eulerisotherm)
1302  SELECT TYPE(c => cvar)
1303  TYPE IS(statevector_eulerisotherm)
1304  DO k=mesh%KGMIN,mesh%KGMAX
1305  DO j=mesh%JGMIN,mesh%JGMAX
1306  DO i=mesh%IGMIN,mesh%IGMAX
1307  p%velocity%data4d(i,j,k,3) = p%velocity%data4d(i,j,k,3) + w(i,j)
1308  c%momentum%data4d(i,j,k,3) = c%momentum%data4d(i,j,k,3) &
1309  + c%density%data3d(i,j,k)*w(i,j)
1310  END DO
1311  END DO
1312  END DO
1313  this%transformed_zvelocity = .false.
1314  END SELECT
1315  END SELECT
1316  END IF
1317  END SUBROUTINE addbackgroundvelocityz
1318 
1319  PURE SUBROUTINE subtractbackgroundvelocityx(this,Mesh,w,pvar,cvar)
1320  IMPLICIT NONE
1321  !------------------------------------------------------------------------!
1322  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
1323  CLASS(mesh_base), INTENT(IN) :: mesh
1324  !------------------------------------------------------------------------!
1325  REAL,DIMENSION(Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
1326  INTENT(IN) :: w
1327  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
1328  !------------------------------------------------------------------------!
1329  INTEGER :: i,j,k
1330  !------------------------------------------------------------------------!
1331  IF (.NOT.this%transformed_xvelocity) THEN
1332  SELECT TYPE(p => pvar)
1333  TYPE IS(statevector_eulerisotherm)
1334  SELECT TYPE(c => cvar)
1335  TYPE IS(statevector_eulerisotherm)
1336  DO k=mesh%KGMIN,mesh%KGMAX
1337  DO j=mesh%JGMIN,mesh%JGMAX
1338  DO i=mesh%IGMIN,mesh%IGMAX
1339  p%velocity%data4d(i,j,k,1) = p%velocity%data4d(i,j,k,1) - w(j,k)
1340  c%momentum%data4d(i,j,k,1) = c%momentum%data4d(i,j,k,1) &
1341  - c%density%data3d(i,j,k)*w(j,k)
1342  END DO
1343  END DO
1344  END DO
1345  this%transformed_xvelocity = .true.
1346  END SELECT
1347  END SELECT
1348  END IF
1349  END SUBROUTINE subtractbackgroundvelocityx
1350 
1351  PURE SUBROUTINE subtractbackgroundvelocityy(this,Mesh,w,pvar,cvar)
1352  IMPLICIT NONE
1353  !------------------------------------------------------------------------!
1354  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
1355  CLASS(mesh_base), INTENT(IN) :: mesh
1356  !------------------------------------------------------------------------!
1357  REAL,DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%KGMIN:Mesh%KGMAX), &
1358  INTENT(IN) :: w
1359  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
1360  !------------------------------------------------------------------------!
1361  INTEGER :: i,j,k
1362  !------------------------------------------------------------------------!
1363  IF (.NOT.this%transformed_yvelocity) THEN
1364  SELECT TYPE(p => pvar)
1365  TYPE IS(statevector_eulerisotherm)
1366  SELECT TYPE(c => cvar)
1367  TYPE IS(statevector_eulerisotherm)
1368  DO k=mesh%KGMIN,mesh%KGMAX
1369  DO j=mesh%JGMIN,mesh%JGMAX
1370  DO i=mesh%IGMIN,mesh%IGMAX
1371  p%velocity%data4d(i,j,k,2) = p%velocity%data4d(i,j,k,2) - w(i,k)
1372  c%momentum%data4d(i,j,k,2) = c%momentum%data4d(i,j,k,2) &
1373  - c%density%data3d(i,j,k)*w(i,k)
1374  END DO
1375  END DO
1376  END DO
1377  this%transformed_yvelocity = .true.
1378  END SELECT
1379  END SELECT
1380  END IF
1381  END SUBROUTINE subtractbackgroundvelocityy
1382 
1383  PURE SUBROUTINE subtractbackgroundvelocityz(this,Mesh,w,pvar,cvar)
1384  IMPLICIT NONE
1385  !------------------------------------------------------------------------!
1386  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
1387  CLASS(mesh_base), INTENT(IN) :: mesh
1388  !------------------------------------------------------------------------!
1389  REAL,DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX), &
1390  INTENT(IN) :: w
1391  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
1392  !------------------------------------------------------------------------!
1393  INTEGER :: i,j,k
1394  !------------------------------------------------------------------------!
1395  IF (.NOT.this%transformed_zvelocity) THEN
1396  SELECT TYPE(p => pvar)
1397  TYPE IS(statevector_eulerisotherm)
1398  SELECT TYPE(c => cvar)
1399  TYPE IS(statevector_eulerisotherm)
1400  DO k=mesh%KGMIN,mesh%KGMAX
1401  DO j=mesh%JGMIN,mesh%JGMAX
1402  DO i=mesh%IGMIN,mesh%IGMAX
1403  p%velocity%data4d(i,j,k,3) = p%velocity%data4d(i,j,k,3) - w(i,j)
1404  c%momentum%data4d(i,j,k,3) = c%momentum%data4d(i,j,k,3) &
1405  - c%density%data3d(i,j,k)*w(i,j)
1406  END DO
1407  END DO
1408  END DO
1409  this%transformed_zvelocity = .true.
1410  END SELECT
1411  END SELECT
1412  END IF
1413  END SUBROUTINE subtractbackgroundvelocityz
1414 
1415 
1416 
1427  PURE SUBROUTINE addfargosources(this,Mesh,w,pvar,cvar,sterm)
1428  IMPLICIT NONE
1429  !------------------------------------------------------------------------!
1430  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1431  CLASS(mesh_base), INTENT(IN) :: mesh
1432  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%KGMIN:Mesh%KGMAX), INTENT(IN) &
1433  :: w
1434  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar,sterm
1435  !------------------------------------------------------------------------!
1436  INTEGER :: i,j,k
1437  !------------------------------------------------------------------------!
1438  SELECT TYPE(c => cvar)
1439  CLASS IS(statevector_eulerisotherm)
1440  SELECT TYPE(s => sterm)
1441  CLASS IS(statevector_eulerisotherm)
1442  DO k=mesh%KMIN,mesh%KMAX
1443  DO j=mesh%JMIN,mesh%JMAX
1444  DO i=mesh%IMIN,mesh%IMAX
1445  ! ATTENTION: fargo sources are added to the given sterm
1446  s%momentum%data4d(i,j,k,2) = s%momentum%data4d(i,j,k,2) &
1447  - c%momentum%data4d(i,j,k,1) * 0.5 * (w(i+1,k)-w(i-1,k)) &
1448  / mesh%dlx%data3d(i,j,k)
1449  END DO
1450  END DO
1451  END DO
1452  END SELECT
1453  END SELECT
1454  END SUBROUTINE addfargosources
1455 
1459  PURE SUBROUTINE reflectionmasks(this,Mesh,reflX,reflY,reflZ)
1460  IMPLICIT NONE
1461  !------------------------------------------------------------------------!
1462  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1463  CLASS(mesh_base), INTENT(IN) :: mesh
1464  LOGICAL, DIMENSION(this%VNUM), INTENT(OUT) :: reflx,refly,reflz
1465  !------------------------------------------------------------------------!
1466  reflx(:) = .false.
1467  refly(:) = .false.
1468  reflz(:) = .false.
1469  SELECT CASE(this%VDIM)
1470  CASE(1) ! 1D velocity
1471  IF (mesh%INUM.GT.1) reflx(2) = .true. ! reflect vx at east/west-boundaries
1472  IF (mesh%JNUM.GT.1) refly(2) = .true. ! reflect vy at south/north-boundaries
1473  IF (mesh%KNUM.GT.1) reflz(2) = .true. ! reflect vz at bottom/top-boundaries
1474  CASE(2) ! 2D velocity
1475  IF (mesh%KNUM.EQ.1.AND..NOT.mesh%ROTSYM.EQ.3) THEN
1476  ! transport in x-y-plane
1477  reflx(2) = .true. ! reflect vx at east/west-boundaries
1478  refly(3) = .true. ! reflect vy at south/north-boundaries
1479  ELSE IF (mesh%JNUM.EQ.1.AND..NOT.mesh%ROTSYM.EQ.2) THEN
1480  ! transport in x-z-plane
1481  reflx(2) = .true. ! reflect vx at east/west-boundaries
1482  reflz(3) = .true. ! reflect vz at bottom/top-boundaries
1483  ELSE IF (mesh%INUM.EQ.1.AND..NOT.mesh%ROTSYM.EQ.1) THEN
1484  ! transport in y-z-plane
1485  refly(2) = .true. ! reflect vy at south/north-boundaries
1486  reflz(3) = .true. ! reflect vz at bottom/top-boundaries
1487  END IF
1488  CASE(3) ! 3D velocity
1489  reflx(2) = .true. ! reflect vx at east/west-boundaries
1490  refly(3) = .true. ! reflect vy at south/north-boundaries
1491  reflz(4) = .true. ! reflect vz at bottom/top-boundaries
1492  END SELECT
1493  END SUBROUTINE reflectionmasks
1494 
1501  PURE SUBROUTINE axismasks(this,Mesh,reflX,reflY,reflZ)
1502  IMPLICIT NONE
1503  !------------------------------------------------------------------------!
1504  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1505  CLASS(mesh_base), INTENT(IN) :: mesh
1506  LOGICAL, DIMENSION(this%VNUM), INTENT(OUT) :: reflx,refly,reflz
1507  !------------------------------------------------------------------------!
1508  INTEGER :: aidx
1509  !------------------------------------------------------------------------!
1510  ! set masks according to reflecting boundaries, i. e. change sign
1511  ! of normal velocities
1512  CALL this%ReflectionMasks(mesh,reflx,refly,reflz)
1513  ! get coordinate index of azimuthal angle (=0 if there is no azimuthal angle)
1514  aidx = mesh%geometry%GetAzimuthIndex()
1515  IF (aidx.GT.0) THEN
1516  ! geometry has azimuthal angle -> determine which velocity changes sign
1517  SELECT CASE(this%VDIM)
1518  CASE(1) ! 1D velocity -> do nothing
1519  CASE(2) ! 2D velocity
1520  IF (mesh%KNUM.EQ.1.AND..NOT.mesh%ROTSYM.EQ.3) THEN
1521  ! transport in x-y-plane
1522  SELECT CASE(aidx)
1523  CASE(1) ! 1st coordinate is the azimuthal angle
1524  refly(2) = .true. ! vx is tangential at south/north-boundaries
1525  CASE(2) ! 2nd coordinate is the azimuthal angle
1526  reflx(3) = .true. ! vy is tangential at east/west-boundaries
1527  CASE(3) ! 3rd coordinate is the azimuthal angle
1528  ! do nothing
1529  END SELECT
1530  ELSE IF (mesh%JNUM.EQ.1.AND..NOT.mesh%ROTSYM.EQ.2) THEN
1531  ! transport in x-z-plane
1532  SELECT CASE(aidx)
1533  CASE(1) ! 1st coordinate is the azimuthal angle
1534  reflz(2) = .true. ! vx is tangential at bottom/top-boundaries
1535  CASE(2) ! 2nd coordinate is the azimuthal angle
1536  ! do nothing
1537  CASE(3) ! 3rd coordinate is the azimuthal angle
1538  reflx(3) = .true. ! vz is tangential at east/west-boundaries
1539  END SELECT
1540  ELSE IF (mesh%INUM.EQ.1.AND..NOT.mesh%ROTSYM.EQ.1) THEN
1541  ! transport in y-z-plane
1542  SELECT CASE(aidx)
1543  CASE(1) ! 1st coordinate is the azimuthal angle
1544  ! do nothing
1545  CASE(2) ! 2nd coordinate is the azimuthal angle
1546  reflz(2) = .true. ! vy is tangential at bottom/top-boundaries
1547  CASE(3) ! 3rd coordinate is the azimuthal angle
1548  refly(3) = .true. ! vz is tangential at south/north-boundaries
1549  END SELECT
1550  END IF
1551  CASE(3) ! 3D velocity
1552  SELECT CASE(aidx)
1553  CASE(1) ! 1st coordinate is the azimuthal angle
1554  reflz(2) = .true. ! vx is tangential at bottom/top-boundaries
1555  CASE(2) ! 2nd coordinate is the azimuthal angle
1556  reflx(3) = .true. ! vy is tangential at east/west-boundaries
1557  CASE(3) ! 3rd coordinate is the azimuthal angle
1558  refly(4) = .true. ! vz is tangential at south/north-boundaries
1559  END SELECT
1560  END SELECT
1561  END IF
1562  END SUBROUTINE axismasks
1563 
1564  PURE SUBROUTINE calculatecharsystemx(this,Mesh,i1,i2,pvar,lambda,xvar)
1565  IMPLICIT NONE
1566  !------------------------------------------------------------------------!
1567  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1568  CLASS(mesh_base), INTENT(IN) :: mesh
1569  INTEGER, INTENT(IN) :: i1,i2
1570  CLASS(marray_compound), INTENT(INOUT) :: pvar
1571  REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM), &
1572  INTENT(OUT) :: lambda,xvar
1573  !------------------------------------------------------------------------!
1574  INTEGER :: il,ir
1575  !------------------------------------------------------------------------!
1576  SELECT TYPE(p => pvar)
1577  TYPE IS(statevector_eulerisotherm)
1578  il = min(i1,i2)
1579  ir = max(i1,i2)
1580  SELECT CASE(this%VDIM)
1581  CASE(1) ! 1D
1582  ! compute eigenvalues at i1
1583  CALL seteigenvalues1d( &
1584  this%bccsound%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1585  p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1586  lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1587  lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2))
1588  ! compute characteristic variables using cell mean values of adjacent
1589  ! cells to calculate derivatives and the isothermal speed of sound
1590  ! at the intermediate cell face
1591  CALL setcharvars1d( &
1592  this%fcsound%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,west), &
1593  p%density%data3d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1594  p%density%data3d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1595  p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1596  p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1597  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1598  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2))
1599  CASE(2) ! 2D
1600  ! compute eigenvalues at i1
1601  CALL seteigenvalues2d( &
1602  this%bccsound%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1603  p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1604  lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1605  lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1606  lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3))
1607  ! compute characteristic variables using cell mean values of adjacent
1608  ! cells to calculate derivatives and the isothermal speed of sound
1609  ! at the intermediate cell face
1610  CALL setcharvars2d( &
1611  this%fcsound%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,west), &
1612  p%density%data3d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1613  p%density%data3d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1614  p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1615  p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1616  p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1617  p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1618  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1619  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1620  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3))
1621  CASE(3) ! 3D
1622  ! compute eigenvalues at i1
1623  CALL seteigenvalues3d( &
1624  this%bccsound%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1625  p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1626  lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1627  lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1628  lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1629  lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4))
1630  ! compute characteristic variables using cell mean values of adjacent
1631  ! cells to calculate derivatives and the isothermal speed of sound
1632  ! at the intermediate cell face
1633  CALL setcharvars3d( &
1634  this%fcsound%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,west), &
1635  p%density%data3d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1636  p%density%data3d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1637  p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1638  p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1639  p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1640  p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1641  p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1642  p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1643  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1644  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1645  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1646  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4))
1647  END SELECT
1648  END SELECT
1649  END SUBROUTINE calculatecharsystemx
1650 
1651 
1652  PURE SUBROUTINE calculatecharsystemy(this,Mesh,j1,j2,pvar,lambda,xvar)
1653  IMPLICIT NONE
1654  !------------------------------------------------------------------------!
1655  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1656  CLASS(mesh_base), INTENT(IN) :: mesh
1657  INTEGER, INTENT(IN) :: j1,j2
1658  CLASS(marray_compound), INTENT(INOUT) :: pvar
1659  REAL, DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM), &
1660  INTENT(OUT) :: lambda,xvar
1661  !------------------------------------------------------------------------!
1662  INTEGER :: jl,jr,vn,vt
1663  !------------------------------------------------------------------------!
1664  SELECT TYPE(p => pvar)
1665  TYPE IS(statevector_eulerisotherm)
1666  jl = min(j1,j2)
1667  jr = max(j1,j2)
1668  SELECT CASE(this%VDIM)
1669  CASE(1) ! 1D
1670  ! compute eigenvalues at j
1671  CALL seteigenvalues1d(this%bccsound%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1672  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,1), &
1673  lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1674  lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2))
1675  ! compute characteristic variables using cell mean values of adjacent
1676  ! cells to calculate derivatives and the isothermal speed of sound
1677  ! at the intermediate cell face
1678  CALL setcharvars1d(this%fcsound%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,south), &
1679  p%density%data3d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX), &
1680  p%density%data3d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX), &
1681  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,1), &
1682  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,1), &
1683  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1684  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2))
1685  CASE(2) ! 2D
1686  ! check which velocity component is along y-direction
1687  SELECT CASE(mesh%VECTOR_COMPONENTS)
1688  CASE(ior(vector_x,vector_y)) ! 2D velocities in x-y-plane
1689  vt = 1
1690  vn = 2
1691  CASE(ior(vector_y,vector_z)) ! 2D velocities in y-z-plane
1692  vt = 2
1693  vn = 1
1694  CASE DEFAULT
1695  ! this should not happen
1696  RETURN
1697  END SELECT
1698  ! compute eigenvalues at j
1699  CALL seteigenvalues2d(this%bccsound%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1700  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,vn), &
1701  lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1702  lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1703  lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3))
1704  ! compute characteristic variables using cell mean values of adjacent
1705  ! cells to calculate derivatives and the isothermal speed of sound
1706  ! at the intermediate cell face
1707  CALL setcharvars2d(this%fcsound%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,south), &
1708  p%density%data3d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX), &
1709  p%density%data3d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX), &
1710  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,vn), &
1711  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,vn), &
1712  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,vt), &
1713  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,vt), &
1714  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1715  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1716  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3))
1717  CASE(3) ! 3D
1718  ! compute eigenvalues at j
1719  CALL seteigenvalues3d(this%bccsound%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1720  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,2), &
1721  lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1722  lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1723  lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1724  lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4))
1725  ! compute characteristic variables using cell mean values of adjacent
1726  ! cells to calculate derivatives and the isothermal speed of sound
1727  ! at the intermediate cell face
1728  CALL setcharvars3d(this%fcsound%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,south), &
1729  p%density%data3d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX), &
1730  p%density%data3d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX), &
1731  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,2), &
1732  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,2), &
1733  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,1), &
1734  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,1), &
1735  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,3), &
1736  p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,3), &
1737  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1738  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1739  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1740  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4))
1741  END SELECT
1742  END SELECT
1743  END SUBROUTINE calculatecharsystemy
1744 
1745 
1746  PURE SUBROUTINE calculatecharsystemz(this,Mesh,k1,k2,pvar,lambda,xvar)
1747  IMPLICIT NONE
1748  !------------------------------------------------------------------------!
1749  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1750  CLASS(mesh_base), INTENT(IN) :: mesh
1751  INTEGER, INTENT(IN) :: k1,k2
1752  CLASS(marray_compound), INTENT(INOUT) :: pvar
1753  REAL, DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,this%VNUM), &
1754  INTENT(OUT) :: lambda,xvar
1755  !------------------------------------------------------------------------!
1756  INTEGER :: kl,kr
1757  !------------------------------------------------------------------------!
1758  SELECT TYPE(p => pvar)
1759  TYPE IS(statevector_eulerisotherm)
1760  kl = min(k1,k2)
1761  kr = max(k1,k2)
1762  SELECT CASE(this%VDIM)
1763  CASE(1) ! 1D
1764  ! compute eigenvalues at k
1765  CALL seteigenvalues1d(this%bccsound%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1766  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,1), &
1767  lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1768  lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2))
1769  ! compute characteristic variables using cell mean values of adjacent
1770  ! cells to calculate derivatives and the isothermal speed of sound
1771  ! at the intermediate cell face
1772  CALL setcharvars1d(this%fcsound%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,bottom), &
1773  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl), &
1774  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr), &
1775  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,1), &
1776  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,1), &
1777  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1778  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2))
1779  CASE(2) ! 2D
1780  ! compute eigenvalues at k
1781  CALL seteigenvalues2d(this%bccsound%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1782  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,2), & ! 2nd component is vz
1783  lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1784  lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1785  lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3))
1786  ! compute characteristic variables using cell mean values of adjacent
1787  ! cells to calculate derivatives and the isothermal speed of sound
1788  ! at the intermediate cell face
1789  CALL setcharvars2d(this%fcsound%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,bottom), &
1790  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl), &
1791  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr), &
1792  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,2), & ! 2nd component is vz
1793  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,2), &
1794  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,1), & ! 1st component: vx or vy
1795  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,1), &
1796  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1797  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1798  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3))
1799  CASE(3) ! 3D
1800  ! compute eigenvalues at k
1801  CALL seteigenvalues3d(this%bccsound%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1802  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,3), &
1803  lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1804  lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1805  lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1806  lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4))
1807  ! compute characteristic variables using cell mean values of adjacent
1808  ! cells to calculate derivatives and the isothermal speed of sound
1809  ! at the intermediate cell face
1810  CALL setcharvars3d(this%fcsound%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,bottom), &
1811  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl), &
1812  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr), &
1813  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,3), &
1814  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,3), &
1815  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,1), &
1816  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,1), &
1817  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,2), &
1818  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,2), &
1819  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1820  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1821  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1822  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4))
1823  END SELECT
1824  END SELECT
1825  END SUBROUTINE calculatecharsystemz
1826 
1827  PURE SUBROUTINE calculateboundarydatax(this,Mesh,i1,i2,xvar,pvar)
1828  IMPLICIT NONE
1829  !------------------------------------------------------------------------!
1830  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1831  CLASS(mesh_base), INTENT(IN) :: mesh
1832  INTEGER, INTENT(IN) :: i1,i2
1833  REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM), &
1834  INTENT(IN) :: xvar
1835  CLASS(marray_compound), INTENT(INOUT) :: pvar
1836  !------------------------------------------------------------------------!
1837  INTEGER :: fidx
1838  !------------------------------------------------------------------------!
1839  IF (i2.LT.i1) THEN
1840  fidx = west
1841  ELSE
1842  fidx = east
1843  END IF
1844  SELECT TYPE(p => pvar)
1845  TYPE IS(statevector_eulerisotherm)
1846  SELECT CASE(this%VDIM)
1847  CASE(1) ! 1D
1848  CALL setboundarydata1d(i2-i1, &
1849  this%fcsound%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,fidx), &
1850  p%density%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1851  p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1852  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1853  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1854  p%density%data3d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1855  p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1))
1856  CASE(2) ! 2D
1857  CALL setboundarydata2d(i2-i1, &
1858  this%fcsound%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,fidx), &
1859  p%density%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1860  p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1861  p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1862  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1863  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1864  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1865  p%density%data3d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1866  p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1867  p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2))
1868  CASE(3) ! 3D
1869  CALL setboundarydata3d(i2-i1, &
1870  this%fcsound%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,fidx), &
1871  p%density%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1872  p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1873  p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1874  p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1875  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1876  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1877  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1878  xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
1879  p%density%data3d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1880  p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1881  p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1882  p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3))
1883  END SELECT
1884  END SELECT
1885  END SUBROUTINE calculateboundarydatax
1886 
1887 
1888  PURE SUBROUTINE calculateboundarydatay(this,Mesh,j1,j2,xvar,pvar)
1889  IMPLICIT NONE
1890  !------------------------------------------------------------------------!
1891  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1892  CLASS(mesh_base), INTENT(IN) :: mesh
1893  INTEGER, INTENT(IN) :: j1,j2
1894  REAL, DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM), &
1895  INTENT(IN) :: xvar
1896  CLASS(marray_compound), INTENT(INOUT) :: pvar
1897  !------------------------------------------------------------------------!
1898  INTEGER :: fidx,vt,vn
1899  !------------------------------------------------------------------------!
1900  IF (j2.LT.j1) THEN
1901  fidx = south
1902  ELSE
1903  fidx = north
1904  END IF
1905  SELECT TYPE(p => pvar)
1906  TYPE IS(statevector_eulerisotherm)
1907  SELECT CASE(this%VDIM)
1908  CASE(1) ! 1D
1909  CALL setboundarydata1d(j2-j1, &
1910  this%fcsound%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,fidx), &
1911  p%density%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1912  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,1), &
1913  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1914  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1915  p%density%data3d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX), &
1916  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,1))
1917  CASE(2) ! 2D
1918  ! check which velocity component is along y-direction
1919  SELECT CASE(mesh%VECTOR_COMPONENTS)
1920  CASE(ior(vector_x,vector_y)) ! 2D velocities in x-y-plane
1921  vt = 1
1922  vn = 2
1923  CASE(ior(vector_y,vector_z)) ! 2D velocities in y-z-plane
1924  vt = 2
1925  vn = 1
1926  CASE DEFAULT
1927  ! this should not happen
1928  RETURN
1929  END SELECT
1930  CALL setboundarydata2d(j2-j1, &
1931  this%fcsound%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,fidx), &
1932  p%density%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1933  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,vn), &
1934  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,vt), &
1935  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1936  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1937  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1938  p%density%data3d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX), &
1939  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,vn), &
1940  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,vt))
1941  CASE(3) ! 3D
1942  CALL setboundarydata3d(j2-j1, &
1943  this%fcsound%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,fidx), &
1944  p%density%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1945  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,2), &
1946  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,1), &
1947  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,3), &
1948  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1949  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1950  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1951  xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
1952  p%density%data3d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX), &
1953  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,2), &
1954  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,1), &
1955  p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,3))
1956  END SELECT
1957  END SELECT
1958  END SUBROUTINE calculateboundarydatay
1959 
1960  PURE SUBROUTINE calculateboundarydataz(this,Mesh,k1,k2,xvar,pvar)
1961  IMPLICIT NONE
1962  !------------------------------------------------------------------------!
1963  CLASS(physics_eulerisotherm), INTENT(IN) :: this
1964  CLASS(mesh_base), INTENT(IN) :: mesh
1965  INTEGER, INTENT(IN) :: k1,k2
1966  REAL, DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,this%VNUM), &
1967  INTENT(IN) :: xvar
1968  CLASS(marray_compound), INTENT(INOUT) :: pvar
1969  !------------------------------------------------------------------------!
1970  INTEGER :: fidx
1971  !------------------------------------------------------------------------!
1972  IF (k2.LT.k1) THEN
1973  fidx = bottom
1974  ELSE
1975  fidx = top
1976  END IF
1977  SELECT TYPE(p => pvar)
1978  TYPE IS(statevector_eulerisotherm)
1979  SELECT CASE(this%VDIM)
1980  CASE(1) ! 1D
1981  CALL setboundarydata1d(k2-k1, &
1982  this%fcsound%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,fidx), &
1983  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1984  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,1), &
1985  xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1986  xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1987  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2), &
1988  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,1))
1989  CASE(2) ! 2D
1990  CALL setboundarydata2d(k2-k1, &
1991  this%fcsound%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,fidx), &
1992  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1993  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,2), &
1994  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,1), &
1995  xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1996  xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1997  xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1998  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2), &
1999  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,2), &
2000  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,1))
2001  CASE(3) ! 3D
2002  CALL setboundarydata3d(k2-k1, &
2003  this%fcsound%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,fidx), &
2004  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
2005  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,3), &
2006  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,1), &
2007  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,2), &
2008  xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
2009  xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
2010  xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
2011  xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
2012  p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2), &
2013  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,3), &
2014  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,1), &
2015  p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,2))
2016  END SELECT
2017  END SELECT
2018  END SUBROUTINE calculateboundarydataz
2019 
2021  SUBROUTINE finalize(this)
2022  IMPLICIT NONE
2023  !------------------------------------------------------------------------!
2024  CLASS(physics_eulerisotherm), INTENT(INOUT) :: this
2025  !------------------------------------------------------------------------!
2026  CALL this%bccsound%Destroy() ! call destructor of the mesh array
2027  CALL this%fcsound%Destroy()
2028  DEALLOCATE(this%bccsound,this%fcsound)
2029  CALL this%Finalize_base()
2030  END SUBROUTINE finalize
2031 
2032 !----------------------------------------------------------------------------!
2034 
2040  FUNCTION createstatevector(Physics,flavour,num) RESULT(new_sv)
2041  IMPLICIT NONE
2042  !-------------------------------------------------------------------!
2043  CLASS(physics_eulerisotherm), INTENT(IN) :: physics
2044  INTEGER, OPTIONAL, INTENT(IN) :: flavour,num
2045  TYPE(statevector_eulerisotherm) :: new_sv
2046  !-------------------------------------------------------------------!
2047  IF (.NOT.physics%Initialized()) &
2048  CALL physics%Error("physics_eulerisotherm::CreateStatevector", "Physics not initialized.")
2049 
2050  ! create a new empty compound of marrays
2051  new_sv = marray_compound(num)
2052  IF (PRESENT(flavour)) THEN
2053  SELECT CASE(flavour)
2054  CASE(primitive,conservative)
2055  new_sv%flavour = flavour
2056  CASE DEFAULT
2057  new_sv%flavour = undefined
2058  END SELECT
2059  END IF
2060  SELECT CASE(new_sv%flavour)
2061  CASE(primitive)
2062  ! allocate memory for density and velocity mesh arrays
2063  ALLOCATE(new_sv%density,new_sv%velocity)
2064  IF (PRESENT(num).AND.num.GT.0) THEN
2065  ! create a bunch of scalars and vectors
2066  new_sv%density = marray_base(num) ! num scalars
2067  new_sv%velocity = marray_base(num,physics%VDIM) ! num vectors
2068  ELSE
2069  new_sv%density = marray_base() ! one scalar
2070  new_sv%velocity = marray_base(physics%VDIM) ! one vector
2071  END IF
2072  ! append to compound
2073  CALL new_sv%AppendMArray(new_sv%density)
2074  CALL new_sv%AppendMArray(new_sv%velocity)
2075  CASE(conservative)
2076  ! allocate memory for density and momentum mesh arrays
2077  ALLOCATE(new_sv%density,new_sv%momentum)
2078  IF (PRESENT(num).AND.num.GT.0) THEN
2079  ! create a bunch of scalars and vectors
2080  new_sv%density = marray_base(num) ! num scalars
2081  new_sv%momentum = marray_base(num,physics%VDIM) ! num vectors
2082  ELSE
2083  new_sv%density = marray_base() ! one scalar
2084  new_sv%momentum = marray_base(physics%VDIM) ! one vector
2085  END IF
2086  ! append to compound
2087  CALL new_sv%AppendMArray(new_sv%density)
2088  CALL new_sv%AppendMArray(new_sv%momentum)
2089  CASE DEFAULT
2090  CALL physics%Warning("physics_eulerisotherm::CreateStateVector", "Empty state vector created.")
2091  END SELECT
2092  END FUNCTION createstatevector
2093 
2094  SUBROUTINE createstatevector_eulerisotherm(Physics,new_sv,flavour,num)
2095  IMPLICIT NONE
2096  !-------------------------------------------------------------------!
2097  CLASS(physics_eulerisotherm), INTENT(IN) :: physics
2098  TYPE(statevector_eulerisotherm), INTENT(INOUT) :: new_sv
2099  INTEGER, OPTIONAL, INTENT(IN) :: flavour,num
2100  !-------------------------------------------------------------------!
2101  IF (.NOT.physics%Initialized()) &
2102  CALL physics%Error("physics_eulerisotherm::CreateStatevector", "Physics not initialized.")
2103 
2104  ! create a new empty compound of marrays
2105  new_sv = marray_compound(num)
2106  IF (PRESENT(flavour)) THEN
2107  SELECT CASE(flavour)
2108  CASE(primitive,conservative)
2109  new_sv%flavour = flavour
2110  CASE DEFAULT
2111  new_sv%flavour = undefined
2112  END SELECT
2113  END IF
2114  SELECT CASE(new_sv%flavour)
2115  CASE(primitive)
2116  ! allocate memory for density and velocity mesh arrays
2117  ALLOCATE(new_sv%density,new_sv%velocity)
2118  ! create a bunch of scalars and vectors
2119  IF (PRESENT(num)) THEN
2120  new_sv%density = marray_base(num) ! num scalars
2121  new_sv%velocity = marray_base(num,physics%VDIM) ! num vectors
2122  ELSE
2123  new_sv%density = marray_base() ! one scalar
2124  new_sv%velocity = marray_base(physics%VDIM) ! one vector
2125  END IF
2126  ! append to compound
2127  CALL new_sv%AppendMArray(new_sv%density)
2128  CALL new_sv%AppendMArray(new_sv%velocity)
2129  CASE(conservative)
2130  ! allocate memory for density and momentum mesh arrays
2131  ALLOCATE(new_sv%density,new_sv%momentum)
2132  ! create a bunch of scalars and vectors
2133  IF (PRESENT(num)) THEN
2134  new_sv%density = marray_base(num) ! num scalars
2135  new_sv%momentum = marray_base(num,physics%VDIM) ! num vectors
2136  ELSE
2137  new_sv%density = marray_base() ! one scalar
2138  new_sv%momentum = marray_base(physics%VDIM) ! one vector
2139  END IF
2140  ! append to compound
2141  CALL new_sv%AppendMArray(new_sv%density)
2142  CALL new_sv%AppendMArray(new_sv%momentum)
2143  CASE DEFAULT
2144  CALL physics%Warning("physics_eulerisotherm::CreateStateVector_eulerisotherm", &
2145  "Empty state vector created.")
2146  END SELECT
2147  END SUBROUTINE createstatevector_eulerisotherm
2148 
2150  SUBROUTINE assignmarray_0(this,ma)
2151  IMPLICIT NONE
2152  !------------------------------------------------------------------------!
2153  CLASS(statevector_eulerisotherm),INTENT(INOUT) :: this
2154  CLASS(marray_base),INTENT(IN) :: ma
2155  !------------------------------------------------------------------------!
2156  CALL this%marray_compound%AssignMArray_0(ma)
2157  IF (SIZE(this%data1d).GT.0) THEN
2158  SELECT TYPE(src => ma)
2159  CLASS IS(statevector_eulerisotherm)
2160  ! copy flavour
2161  this%flavour = src%flavour
2162  ! set pointer to the data structures in the compound
2165  this%density => this%GetItem(this%FirstItem()) ! density is the first item
2166  SELECT CASE(this%flavour)
2167  CASE(primitive)
2168  ! velocity is the second item
2169  this%velocity => this%GetItem(this%NextItem(this%FirstItem()))
2170  CASE(conservative)
2171  ! momentum is the second item
2172  this%momentum => this%GetItem(this%NextItem(this%FirstItem()))
2173  CASE DEFAULT
2174  ! error, this should not happen
2175  END SELECT
2176  CLASS DEFAULT
2177  ! error, this should not happen
2178  END SELECT
2179  END IF
2180  END SUBROUTINE assignmarray_0
2181 
2182 
2183 !----------------------------------------------------------------------------!
2185 
2187  ELEMENTAL SUBROUTINE setwavespeeds(cs,v,minwav,maxwav)
2188  IMPLICIT NONE
2189  !------------------------------------------------------------------------!
2190  REAL, INTENT(IN) :: cs,v
2191  REAL, INTENT(OUT) :: minwav,maxwav
2192  !------------------------------------------------------------------------!
2193  minwav = min(0.,v-cs) ! minimal wave speed
2194  maxwav = max(0.,v+cs) ! maximal wave speed
2195  END SUBROUTINE setwavespeeds
2196 
2198  ELEMENTAL SUBROUTINE seteigenvalues1d(cs,v,l1,l2)
2199  IMPLICIT NONE
2200  !------------------------------------------------------------------------!
2201  REAL, INTENT(IN) :: cs,v
2202  REAL, INTENT(OUT) :: l1,l2
2203  !------------------------------------------------------------------------!
2204  l1 = v - cs
2205  l2 = v + cs
2206  END SUBROUTINE seteigenvalues1d
2207 
2208  ELEMENTAL SUBROUTINE seteigenvalues2d(cs,v,l1,l2,l3)
2209  IMPLICIT NONE
2210  !------------------------------------------------------------------------!
2211  REAL, INTENT(IN) :: cs,v
2212  REAL, INTENT(OUT) :: l1,l2,l3
2213  !------------------------------------------------------------------------!
2214  l1 = v - cs
2215  l2 = v
2216  l3 = v + cs
2217  END SUBROUTINE seteigenvalues2d
2218 
2219  ELEMENTAL SUBROUTINE seteigenvalues3d(cs,v,l1,l2,l3,l4)
2220  IMPLICIT NONE
2221  !------------------------------------------------------------------------!
2222  REAL, INTENT(IN) :: cs,v
2223  REAL, INTENT(OUT) :: l1,l2,l3,l4
2224  !------------------------------------------------------------------------!
2225  l1 = v - cs
2226  l2 = v
2227  l3 = v
2228  l4 = v + cs
2229  END SUBROUTINE seteigenvalues3d
2230 
2232  ELEMENTAL SUBROUTINE setflux1d(cs,rho,u,mu,f1,f2)
2233  IMPLICIT NONE
2234  !------------------------------------------------------------------------!
2235  REAL, INTENT(IN) :: cs,rho,u,mu
2236  REAL, INTENT(OUT) :: f1,f2
2237  !------------------------------------------------------------------------!
2238  f1 = rho*u
2239  f2 = mu*u + rho*cs*cs
2240  END SUBROUTINE setflux1d
2241 
2243  ELEMENTAL SUBROUTINE setflux2d(cs,rho,u,mu,mv,f1,f2,f3)
2244  IMPLICIT NONE
2245  !------------------------------------------------------------------------!
2246  REAL, INTENT(IN) :: cs,rho,u,mu,mv
2247  REAL, INTENT(OUT) :: f1,f2,f3
2248  !------------------------------------------------------------------------!
2249  CALL setflux1d(cs,rho,u,mu,f1,f2)
2250  f3 = mv*u
2251  END SUBROUTINE setflux2d
2252 
2254  ELEMENTAL SUBROUTINE setflux3d(cs,rho,u,mu,mv,mw,f1,f2,f3,f4)
2255  IMPLICIT NONE
2256  !------------------------------------------------------------------------!
2257  REAL, INTENT(IN) :: cs,rho,u,mu,mv,mw
2258  REAL, INTENT(OUT) :: f1,f2,f3,f4
2259  !------------------------------------------------------------------------!
2260  CALL setflux2d(cs,rho,u,mu,mv,f1,f2,f3)
2261  f4 = mw*u
2262  END SUBROUTINE setflux3d
2263 
2265  ELEMENTAL SUBROUTINE cons2prim1d(rho_in,mu,rho_out,u)
2266  IMPLICIT NONE
2267  !------------------------------------------------------------------------!
2268  REAL, INTENT(IN) :: rho_in,mu
2269  REAL, INTENT(OUT) :: rho_out,u
2270  !------------------------------------------------------------------------!
2271  rho_out = rho_in
2272  u = mu / rho_in
2273  END SUBROUTINE cons2prim1d
2274 
2276  ELEMENTAL SUBROUTINE cons2prim2d(rho_in,mu,mv,rho_out,u,v)
2277  IMPLICIT NONE
2278  !------------------------------------------------------------------------!
2279  REAL, INTENT(IN) :: rho_in,mu,mv
2280  REAL, INTENT(OUT) :: rho_out,u,v
2281  !------------------------------------------------------------------------!
2282  REAL :: inv_rho
2283  !------------------------------------------------------------------------!
2284  inv_rho = 1./rho_in
2285  rho_out = rho_in
2286  u = mu * inv_rho
2287  v = mv * inv_rho
2288  END SUBROUTINE cons2prim2d
2289 
2291  ELEMENTAL SUBROUTINE cons2prim3d(rho_in,mu,mv,mw,rho_out,u,v,w)
2292  IMPLICIT NONE
2293  !------------------------------------------------------------------------!
2294  REAL, INTENT(IN) :: rho_in,mu,mv,mw
2295  REAL, INTENT(OUT) :: rho_out,u,v,w
2296  !------------------------------------------------------------------------!
2297  REAL :: inv_rho
2298  !------------------------------------------------------------------------!
2299  inv_rho = 1./rho_in
2300  rho_out = rho_in
2301  u = mu * inv_rho
2302  v = mv * inv_rho
2303  w = mw * inv_rho
2304  END SUBROUTINE cons2prim3d
2305 
2307  ELEMENTAL SUBROUTINE prim2cons1d(rho_in,u,rho_out,mu)
2308  IMPLICIT NONE
2309  !------------------------------------------------------------------------!
2310  REAL, INTENT(IN) :: rho_in,u
2311  REAL, INTENT(OUT) :: rho_out,mu
2312  !------------------------------------------------------------------------!
2313  rho_out = rho_in
2314  mu = rho_in * u
2315  END SUBROUTINE prim2cons1d
2316 
2318  ELEMENTAL SUBROUTINE prim2cons2d(rho_in,u,v,rho_out,mu,mv)
2319  IMPLICIT NONE
2320  !------------------------------------------------------------------------!
2321  REAL, INTENT(IN) :: rho_in,u,v
2322  REAL, INTENT(OUT) :: rho_out,mu,mv
2323  !------------------------------------------------------------------------!
2324  CALL prim2cons1d(rho_in,u,rho_out,mu)
2325  mv = rho_in * v
2326  END SUBROUTINE prim2cons2d
2327 
2329  ELEMENTAL SUBROUTINE prim2cons3d(rho_in,u,v,w,rho_out,mu,mv,mw)
2330  IMPLICIT NONE
2331  !------------------------------------------------------------------------!
2332  REAL, INTENT(IN) :: rho_in,u,v,w
2333  REAL, INTENT(OUT) :: rho_out,mu,mv,mw
2334  !------------------------------------------------------------------------!
2335  CALL prim2cons2d(rho_in,u,v,rho_out,mu,mv)
2336  mw = rho_in * w
2337  END SUBROUTINE prim2cons3d
2338 
2340  ELEMENTAL SUBROUTINE setcharvars1d(cs,rho1,rho2,u1,u2,xvar1,xvar2)
2341  IMPLICIT NONE
2342  !------------------------------------------------------------------------!
2343  REAL, INTENT(IN) :: cs,rho1,rho2,u1,u2
2344  REAL, INTENT(OUT) :: xvar1,xvar2
2345  !------------------------------------------------------------------------!
2346  REAL :: dlnrho,ducs
2347  !------------------------------------------------------------------------!
2348  dlnrho = log(rho2/rho1)
2349  ducs = (u2-u1)/cs
2350  ! characteristic variables
2351  xvar1 = dlnrho - ducs
2352  xvar2 = dlnrho + ducs
2353  END SUBROUTINE setcharvars1d
2354 
2356  ELEMENTAL SUBROUTINE setcharvars2d(cs,rho1,rho2,u1,u2,v1,v2,xvar1,xvar2,xvar3)
2357  IMPLICIT NONE
2358  !------------------------------------------------------------------------!
2359  REAL, INTENT(IN) :: cs,rho1,rho2,u1,u2,v1,v2
2360  REAL, INTENT(OUT) :: xvar1,xvar2,xvar3
2361  !------------------------------------------------------------------------!
2362  CALL setcharvars1d(cs,rho1,rho2,u1,u2,xvar1,xvar3)
2363  xvar2 = v2-v1
2364  END SUBROUTINE setcharvars2d
2365 
2367  ELEMENTAL SUBROUTINE setcharvars3d(cs,rho1,rho2,u1,u2,v1,v2,w1,w2,&
2368  xvar1,xvar2,xvar3,xvar4)
2369  IMPLICIT NONE
2370  !------------------------------------------------------------------------!
2371  REAL, INTENT(IN) :: cs,rho1,rho2,u1,u2,v1,v2,w1,w2
2372  REAL, INTENT(OUT) :: xvar1,xvar2,xvar3,xvar4
2373  !------------------------------------------------------------------------!
2374  REAL :: dlnrho,du
2375  !------------------------------------------------------------------------!
2376  CALL setcharvars2d(cs,rho1,rho2,u1,u2,v1,v2,xvar1,xvar2,xvar4)
2377  xvar3 = w2-w1
2378  END SUBROUTINE setcharvars3d
2379 
2381  ELEMENTAL SUBROUTINE setboundarydata1d(delta,cs,rho1,u1,xvar1,xvar2,rho2,u2)
2382  IMPLICIT NONE
2383  !------------------------------------------------------------------------!
2384  INTEGER, INTENT(IN) :: delta
2385  REAL, INTENT(IN) :: cs,rho1,u1,xvar1,xvar2
2386  REAL, INTENT(OUT) :: rho2,u2
2387  !------------------------------------------------------------------------!
2388  rho2 = rho1 * exp(delta*0.5*(xvar2+xvar1))
2389  u2 = u1 + delta*0.5*cs*(xvar2-xvar1)
2390  END SUBROUTINE setboundarydata1d
2391 
2392  ELEMENTAL SUBROUTINE setboundarydata2d(delta,cs,rho1,u1,v1,xvar1,xvar2,xvar3, &
2393  rho2,u2,v2)
2394  IMPLICIT NONE
2395  !------------------------------------------------------------------------!
2396  INTEGER, INTENT(IN) :: delta
2397  REAL, INTENT(IN) :: cs,rho1,u1,v1,xvar1,xvar2,xvar3
2398  REAL, INTENT(OUT) :: rho2,u2,v2
2399  !------------------------------------------------------------------------!
2400  CALL setboundarydata1d(delta,cs,rho1,u1,xvar1,xvar3,rho2,u2)
2401  v2 = v1 + delta*xvar2
2402  END SUBROUTINE setboundarydata2d
2403 
2404  ELEMENTAL SUBROUTINE setboundarydata3d(delta,cs,rho1,u1,v1,w1,xvar1,xvar2,xvar3, &
2405  xvar4,rho2,u2,v2,w2)
2406  IMPLICIT NONE
2407  !------------------------------------------------------------------------!
2408  INTEGER, INTENT(IN) :: delta
2409  REAL, INTENT(IN) :: cs,rho1,u1,v1,w1,xvar1,xvar2,xvar3,xvar4
2410  REAL, INTENT(OUT) :: rho2,u2,v2,w2
2411  !------------------------------------------------------------------------!
2412  CALL setboundarydata2d(delta,cs,rho1,u1,v1,xvar1,xvar2,xvar4,rho2,u2,v2)
2413  w2 = w1 + delta*xvar3
2414  END SUBROUTINE setboundarydata3d
2415 
2425  ELEMENTAL FUNCTION getgeometricalsourcex(cxyx,cxzx,cyxy,czxz,vx,vy,vz,P,my,mz)
2426  IMPLICIT NONE
2427  !------------------------------------------------------------------------!
2428  REAL, INTENT(IN) :: cxyx,cxzx,cyxy,czxz,vx,vy,vz,p,my,mz
2429  REAL :: getgeometricalsourcex
2430  !------------------------------------------------------------------------!
2431  getgeometricalsourcex = my * (cyxy*vy - cxyx*vx) + mz * (czxz*vz - cxzx*vx) &
2432  + (cyxy + czxz) * p
2433  END FUNCTION getgeometricalsourcex
2434 
2436  ELEMENTAL FUNCTION getgeometricalsourcey(cxyx,cyxy,cyzy,czyz,vx,vy,vz,P,mx,mz)
2437  IMPLICIT NONE
2438  !------------------------------------------------------------------------!
2439  REAL, INTENT(IN) :: cxyx,cyxy,cyzy,czyz,vx,vy,vz,p,mx,mz
2440  REAL :: getgeometricalsourcey
2441  !------------------------------------------------------------------------!
2442  getgeometricalsourcey = mz * (czyz*vz - cyzy*vy) + mx * (cxyx*vx - cyxy*vy) &
2443  + (cxyx + czyz) * p
2444  END FUNCTION getgeometricalsourcey
2445 
2447  ELEMENTAL FUNCTION getgeometricalsourcez(cxzx,cyzy,czxz,czyz,vx,vy,vz,P,mx,my)
2448  IMPLICIT NONE
2449  !------------------------------------------------------------------------!
2450  REAL, INTENT(IN) :: cxzx,cyzy,czxz,czyz,vx,vy,vz,p,mx,my
2451  REAL :: getgeometricalsourcez
2452  !------------------------------------------------------------------------!
2453  getgeometricalsourcez = mx * (cxzx*vx - czxz*vz) + my * (cyzy*vy - czyz*vz) &
2454  + (cxzx + cyzy) * p
2455  END FUNCTION getgeometricalsourcez
2456 
2457 ! ! TODO: Not verified
2458 ! !!
2459 ! !! non-global elemental routine
2460 ! ELEMENTAL SUBROUTINE SetRoeAverages(rhoL,rhoR,vL,vR,v)
2461 ! IMPLICIT NONE
2462 ! !------------------------------------------------------------------------!
2463 ! REAL, INTENT(IN) :: rhoL,rhoR,vL,vR
2464 ! REAL, INTENT(OUT) :: v
2465 ! !------------------------------------------------------------------------!
2466 ! REAL :: sqrtrhoL,sqrtrhoR
2467 ! !------------------------------------------------------------------------!
2468 ! sqrtrhoL = SQRT(rhoL)
2469 ! sqrtrhoR = SQRT(rhoR)
2470 ! v = 0.5*(sqrtrhoL*vL + sqrtrhoR*vR) / (sqrtrhoL + sqrtrhoR)
2471 ! END SUBROUTINE SetRoeAverages
2472 
2473 END MODULE physics_eulerisotherm_mod
physics_eulerisotherm_mod::subtractbackgroundvelocityy
pure subroutine subtractbackgroundvelocityy(this, Mesh, w, pvar, cvar)
Definition: physics_eulerisotherm.f90:1352
mesh_base
logging_base_mod::finalize
subroutine finalize(this)
Destructor of common class.
Definition: logging_base.f90:155
physics_eulerisotherm_mod::calculateboundarydataz
pure subroutine calculateboundarydataz(this, Mesh, k1, k2, xvar, pvar)
Definition: physics_eulerisotherm.f90:1961
physics_eulerisotherm_mod::setwavespeeds
elemental subroutine setwavespeeds(cs, v, minwav, maxwav)
set minimal and maximal wave speeds
Definition: physics_eulerisotherm.f90:2188
physics_eulerisotherm_mod::calcwavespeeds_center
pure subroutine calcwavespeeds_center(this, Mesh, pvar, minwav, maxwav)
Calculates wave speeds at cell-centers.
Definition: physics_eulerisotherm.f90:532
physics_eulerisotherm_mod::printconfiguration_eulerisotherm
subroutine printconfiguration_eulerisotherm(this)
Definition: physics_eulerisotherm.f90:245
physics_eulerisotherm_mod::prim2cons1d
elemental subroutine prim2cons1d(rho_in, u, rho_out, mu)
Convert from 1D primitive to conservative variables.
Definition: physics_eulerisotherm.f90:2308
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
physics_eulerisotherm_mod::num_var
integer, parameter num_var
Definition: physics_eulerisotherm.f90:54
physics_eulerisotherm_mod::cons2prim2d
elemental subroutine cons2prim2d(rho_in, mu, mv, rho_out, u, v)
Convert from 2D conservative to primitive variables.
Definition: physics_eulerisotherm.f90:2277
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
physics_eulerisotherm_mod::convert2primitive_subset
pure subroutine convert2primitive_subset(this, i1, i2, j1, j2, k1, k2, cvar, pvar)
Converts to primitives at cell centers using state vectors.
Definition: physics_eulerisotherm.f90:350
physics_eulerisotherm_mod::setcharvars3d
elemental subroutine setcharvars3d(cs, rho1, rho2, u1, u2, v1, v2, w1, w2, xvar1, xvar2, xvar3, xvar4)
compute characteristic variables using adjacent primitve states
Definition: physics_eulerisotherm.f90:2369
logging_base_mod::nan_default_real
real, save nan_default_real
NaN real constant.
Definition: logging_base.f90:85
physics_eulerisotherm_mod::seteigenvalues2d
elemental subroutine seteigenvalues2d(cs, v, l1, l2, l3)
Definition: physics_eulerisotherm.f90:2209
physics_eulerisotherm_mod::prim2cons2d
elemental subroutine prim2cons2d(rho_in, u, v, rho_out, mu, mv)
Convert from 2D primitive to conservative variables.
Definition: physics_eulerisotherm.f90:2319
physics_eulerisotherm_mod::addbackgroundvelocityz
pure subroutine addbackgroundvelocityz(this, Mesh, w, pvar, cvar)
Definition: physics_eulerisotherm.f90:1288
physics_eulerisotherm_mod::physics_eulerisotherm
Definition: physics_eulerisotherm.f90:57
physics_eulerisotherm_mod::calcfluxesy
pure subroutine calcfluxesy(this, Mesh, nmin, nmax, prim, cons, yfluxes)
Calculate Fluxes in y-direction.
Definition: physics_eulerisotherm.f90:881
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
physics_eulerisotherm_mod::problem_name
character(len=32), parameter problem_name
Definition: physics_eulerisotherm.f90:55
physics_eulerisotherm_mod::calculatecharsystemz
pure subroutine calculatecharsystemz(this, Mesh, k1, k2, pvar, lambda, xvar)
Definition: physics_eulerisotherm.f90:1747
physics_eulerisotherm_mod::viscositysources
pure subroutine viscositysources(this, Mesh, pvar, btxx, btxy, btxz, btyy, btyz, btzz, sterm)
compute viscous source terms
Definition: physics_eulerisotherm.f90:979
physics_eulerisotherm_mod::addbackgroundvelocityx
pure subroutine addbackgroundvelocityx(this, Mesh, w, pvar, cvar)
Definition: physics_eulerisotherm.f90:1223
physics_eulerisotherm_mod::prim2cons
Definition: physics_eulerisotherm.f90:139
physics_eulerisotherm_mod::initphysics_eulerisotherm
subroutine initphysics_eulerisotherm(this, Mesh, config, IO)
Intialization of isothermal physics.
Definition: physics_eulerisotherm.f90:160
physics_eulerisotherm_mod::setcharvars1d
elemental subroutine setcharvars1d(cs, rho1, rho2, u1, u2, xvar1, xvar2)
compute characteristic variables using adjacent primitve states
Definition: physics_eulerisotherm.f90:2341
physics_eulerisotherm_mod::seteigenvalues1d
elemental subroutine seteigenvalues1d(cs, v, l1, l2)
compute all eigenvalues
Definition: physics_eulerisotherm.f90:2199
physics_eulerisotherm_mod::calcstresses
pure subroutine calcstresses(this, Mesh, pvar, dynvis, bulkvis, btxx, btxy, btxz, btyy, btyz, btzz)
calculate components of the stress tensor
Definition: physics_eulerisotherm.f90:1021
physics_eulerisotherm_mod::convert2conservative_subset
pure subroutine convert2conservative_subset(this, i1, i2, j1, j2, k1, k2, pvar, cvar)
Converts to primitive to conservative variables on a subset of the data.
Definition: physics_eulerisotherm.f90:458
physics_eulerisotherm_mod::setsoundspeeds_center
pure subroutine setsoundspeeds_center(this, Mesh, bccsound)
Sets soundspeeds at cell-centers.
Definition: physics_eulerisotherm.f90:292
physics_eulerisotherm_mod::externalsources
pure subroutine externalsources(this, accel, pvar, cvar, sterm)
compute momentum sources given an external force
Definition: physics_eulerisotherm.f90:1198
physics_eulerisotherm_mod::geometricalsources
pure subroutine geometricalsources(this, Mesh, pvar, cvar, sterm)
Calculates geometrical sources.
Definition: physics_eulerisotherm.f90:691
physics_eulerisotherm_mod::setsoundspeeds_faces
pure subroutine setsoundspeeds_faces(this, Mesh, fcsound)
Sets soundspeeds at cell-faces.
Definition: physics_eulerisotherm.f90:304
physics_eulerisotherm_mod::setflux1d
elemental subroutine setflux1d(cs, rho, u, mu, f1, f2)
set mass and 1D momentum flux for transport along the 1st dimension
Definition: physics_eulerisotherm.f90:2233
physics_eulerisotherm_mod::createstatevector_eulerisotherm
subroutine, public createstatevector_eulerisotherm(Physics, new_sv, flavour, num)
Definition: physics_eulerisotherm.f90:2095
physics_eulerisotherm_mod::cons2prim1d
elemental subroutine cons2prim1d(rho_in, mu, rho_out, u)
Convert from 1D conservative to primitive variables.
Definition: physics_eulerisotherm.f90:2266
physics_eulerisotherm_mod::convert2conservative_all
pure subroutine convert2conservative_all(this, pvar, cvar)
Converts primitive to conservative variables on the whole mesh.
Definition: physics_eulerisotherm.f90:424
physics_eulerisotherm_mod::setflux3d
elemental subroutine setflux3d(cs, rho, u, mu, mv, mw, f1, f2, f3, f4)
set mass and 3D momentum flux for transport along the 1st dimension
Definition: physics_eulerisotherm.f90:2255
physics_eulerisotherm_mod::axismasks
pure subroutine axismasks(this, Mesh, reflX, reflY, reflZ)
return masks for axis boundaries
Definition: physics_eulerisotherm.f90:1502
physics_eulerisotherm_mod::calcwavespeeds_faces
pure subroutine calcwavespeeds_faces(this, Mesh, prim, cons, minwav, maxwav)
Calculates wave speeds at cell-faces.
Definition: physics_eulerisotherm.f90:557
physics_eulerisotherm_mod::reflectionmasks
pure subroutine reflectionmasks(this, Mesh, reflX, reflY, reflZ)
return masks for reflecting boundaries
Definition: physics_eulerisotherm.f90:1460
physics_eulerisotherm_mod::subtractbackgroundvelocityx
pure subroutine subtractbackgroundvelocityx(this, Mesh, w, pvar, cvar)
Definition: physics_eulerisotherm.f90:1320
physics_eulerisotherm_mod::cons2prim
Definition: physics_eulerisotherm.f90:136
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:64
physics_eulerisotherm_mod
physics module for 1D,2D and 3D isothermal Euler equations
Definition: physics_eulerisotherm.f90:44
physics_eulerisotherm_mod::enableoutput
subroutine enableoutput(this, Mesh, config, IO)
Enables output of certain arrays defined in this class.
Definition: physics_eulerisotherm.f90:254
physics_eulerisotherm_mod::seteigenvalues3d
elemental subroutine seteigenvalues3d(cs, v, l1, l2, l3, l4)
Definition: physics_eulerisotherm.f90:2220
physics_base_mod::physics_base
named integer constants for flavour of state vectors
Definition: physics_base.f90:69
physics_eulerisotherm_mod::setboundarydata1d
elemental subroutine setboundarydata1d(delta, cs, rho1, u1, xvar1, xvar2, rho2, u2)
extrapolate boundary values using primitve and characteristic variables
Definition: physics_eulerisotherm.f90:2382
physics_eulerisotherm_mod::createstatevector
type(statevector_eulerisotherm) function createstatevector(Physics, flavour, num)
Constructor of statevector_eulerisotherm.
Definition: physics_eulerisotherm.f90:2041
physics_eulerisotherm_mod::calculatecharsystemx
pure subroutine calculatecharsystemx(this, Mesh, i1, i2, pvar, lambda, xvar)
Definition: physics_eulerisotherm.f90:1565
physics_eulerisotherm_mod::addfargosources
pure subroutine addfargosources(this, Mesh, w, pvar, cvar, sterm)
sources terms for fargo advection
Definition: physics_eulerisotherm.f90:1428
physics_base_mod::euler_isotherm
integer, parameter, public euler_isotherm
Definition: physics_base.f90:391
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
physics_eulerisotherm_mod::getgeometricalsourcex
elemental real function getgeometricalsourcex(cxyx, cxzx, cyxy, czxz, vx, vy, vz, P, my, mz)
geometrical momentum source terms P is the either isothermal pressure rho*cs**2 or the real pressure...
Definition: physics_eulerisotherm.f90:2426
physics_eulerisotherm_mod::calculateboundarydatax
pure subroutine calculateboundarydatax(this, Mesh, i1, i2, xvar, pvar)
Definition: physics_eulerisotherm.f90:1828
physics_eulerisotherm_mod::addbackgroundvelocityy
pure subroutine addbackgroundvelocityy(this, Mesh, w, pvar, cvar)
Definition: physics_eulerisotherm.f90:1256
physics_eulerisotherm_mod::calculateboundarydatay
pure subroutine calculateboundarydatay(this, Mesh, j1, j2, xvar, pvar)
Definition: physics_eulerisotherm.f90:1889
physics_eulerisotherm_mod::cons2prim3d
elemental subroutine cons2prim3d(rho_in, mu, mv, mw, rho_out, u, v, w)
Convert from 3D conservative to primitive variables.
Definition: physics_eulerisotherm.f90:2292
physics_eulerisotherm_mod::getgeometricalsourcey
elemental real function getgeometricalsourcey(cxyx, cyxy, cyzy, czyz, vx, vy, vz, P, mx, mz)
y-momentum geometrical source term
Definition: physics_eulerisotherm.f90:2437
physics_eulerisotherm_mod::calculatecharsystemy
pure subroutine calculatecharsystemy(this, Mesh, j1, j2, pvar, lambda, xvar)
Definition: physics_eulerisotherm.f90:1653
physics_eulerisotherm_mod::convert2primitive_all
pure subroutine convert2primitive_all(this, cvar, pvar)
Converts to primitives at cell centers using state vectors.
Definition: physics_eulerisotherm.f90:316
physics_eulerisotherm_mod::calcfluxesx
pure subroutine calcfluxesx(this, Mesh, nmin, nmax, prim, cons, xfluxes)
Calculate Fluxes in x-direction.
Definition: physics_eulerisotherm.f90:832
physics_eulerisotherm_mod::setboundarydata3d
elemental subroutine setboundarydata3d(delta, cs, rho1, u1, v1, w1, xvar1, xvar2, xvar3, xvar4, rho2, u2, v2, w2)
Definition: physics_eulerisotherm.f90:2406
physics_eulerisotherm_mod::getgeometricalsourcez
elemental real function getgeometricalsourcez(cxzx, cyzy, czxz, czyz, vx, vy, vz, P, mx, my)
z-momentum geometrical source term
Definition: physics_eulerisotherm.f90:2448
physics_eulerisotherm_mod::statevector_eulerisotherm
Definition: physics_eulerisotherm.f90:121
physics_eulerisotherm_mod::calcfluxesz
pure subroutine calcfluxesz(this, Mesh, nmin, nmax, prim, cons, zfluxes)
Calculate Fluxes in z-direction.
Definition: physics_eulerisotherm.f90:930
physics_eulerisotherm_mod::setboundarydata2d
elemental subroutine setboundarydata2d(delta, cs, rho1, u1, v1, xvar1, xvar2, xvar3, rho2, u2, v2)
Definition: physics_eulerisotherm.f90:2394
physics_eulerisotherm_mod::subtractbackgroundvelocityz
pure subroutine subtractbackgroundvelocityz(this, Mesh, w, pvar, cvar)
Definition: physics_eulerisotherm.f90:1384
physics_eulerisotherm_mod::setflux
Definition: physics_eulerisotherm.f90:133
physics_eulerisotherm_mod::setflux2d
elemental subroutine setflux2d(cs, rho, u, mu, mv, f1, f2, f3)
set mass and 2D momentum flux for transport along the 1st dimension
Definition: physics_eulerisotherm.f90:2244
marray_base_mod::assignmarray_0
subroutine assignmarray_0(this, ma)
assigns one mesh array to another mesh array
Definition: marray_base.f90:254
physics_eulerisotherm_mod::setcharvars2d
elemental subroutine setcharvars2d(cs, rho1, rho2, u1, u2, v1, v2, xvar1, xvar2, xvar3)
compute characteristic variables using adjacent primitve states
Definition: physics_eulerisotherm.f90:2357
physics_eulerisotherm_mod::prim2cons3d
elemental subroutine prim2cons3d(rho_in, u, v, w, rho_out, mu, mv, mw)
Convert from 3D primitive to conservative variables.
Definition: physics_eulerisotherm.f90:2330
physics_base_mod::new_statevector
Definition: physics_base.f90:188
Generated by   doxygen 1.8.14