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