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-2024 #
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!----------------------------------------------------------------------------!
41MODULE 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
58 PROCEDURE :: new_statevector
59 !------Convert2Primitve--------!
60 PROCEDURE :: convert2primitive_all
61 PROCEDURE :: convert2primitive_subset
62 !------Convert2Conservative----!
63 PROCEDURE :: convert2conservative_all
64 PROCEDURE :: convert2conservative_subset
65 !------soundspeed routines-----!
66 PROCEDURE :: updatesoundspeed
67 !------flux routines-----------!
68 PROCEDURE :: calcfluxesx
69 PROCEDURE :: calcfluxesy
70 PROCEDURE :: calcfluxesz
71 !------fargo routines----------!
72 PROCEDURE :: addbackgroundvelocityx
73 PROCEDURE :: addbackgroundvelocityy
74 PROCEDURE :: addbackgroundvelocityz
75 PROCEDURE :: subtractbackgroundvelocityx
76 PROCEDURE :: subtractbackgroundvelocityy
77 PROCEDURE :: subtractbackgroundvelocityz
78 PROCEDURE :: addfargosourcesx
79 PROCEDURE :: addfargosourcesy
80 PROCEDURE :: addfargosourcesz
81 !------HLLC routines-----------!
82! PROCEDURE :: CalcIntermediateStateX
83! PROCEDURE :: CalcIntermediateStateY
84
85 PROCEDURE :: externalsources
86 PROCEDURE :: geometricalsources
87 PROCEDURE :: viscositysources
88
89 ! boundarie routines
90 PROCEDURE :: calculatecharsystemx ! for absorbing boundaries
91 PROCEDURE :: calculatecharsystemy ! for absorbing boundaries
92 PROCEDURE :: calculatecharsystemz ! for absorbing boundaries
93 PROCEDURE :: calculateboundarydatax ! for absorbing boundaries
94 PROCEDURE :: calculateboundarydatay ! for absorbing boundaries
95 PROCEDURE :: calculateboundarydataz ! for absorbing boundaries
96 PROCEDURE :: calculateprim2riemannx ! for farfield boundaries
97 PROCEDURE :: calculateprim2riemanny ! for farfield boundaries
98 PROCEDURE :: calculateprim2riemannz ! for farfield boundaries
99 PROCEDURE :: calculateriemann2primx ! for farfield boundaries
100 PROCEDURE :: calculateriemann2primy ! for farfield boundaries
101 PROCEDURE :: calculateriemann2primz ! for farfield boundaries
102
103 PROCEDURE :: finalize
104 END TYPE
105 TYPE, EXTENDS(statevector_eulerisotherm) :: statevector_euler
106 TYPE(marray_base), POINTER :: &
107 pressure => null(), &
108 energy => null()
109 CONTAINS
110 PROCEDURE :: assignmarray_0
111 final :: finalize_statevector
112 END TYPE
113 INTERFACE statevector_euler
114 MODULE PROCEDURE createstatevector
115 END INTERFACE
116 INTERFACE setflux
117 MODULE PROCEDURE setflux1d, setflux2d, setflux3d
118 END INTERFACE
119 INTERFACE cons2prim
120 MODULE PROCEDURE cons2prim1d, cons2prim2d, cons2prim3d
121 END INTERFACE
122 INTERFACE prim2cons
123 MODULE PROCEDURE prim2cons1d, prim2cons2d, prim2cons3d
124 END INTERFACE
125 !--------------------------------------------------------------------------!
126 PUBLIC :: &
127 ! types
128 physics_euler, &
129 statevector_euler
130 !--------------------------------------------------------------------------!
131
132CONTAINS
133
135 SUBROUTINE initphysics(this,Mesh,config,IO)
136 IMPLICIT NONE
137 !------------------------------------------------------------------------!
138 CLASS(physics_euler), INTENT(INOUT) :: this
139 CLASS(mesh_base), INTENT(IN) :: Mesh
140 TYPE(dict_typ), POINTER :: config, IO
141 !------------------------------------------------------------------------!
142 INTEGER :: err,next_idx
143 CHARACTER(LEN=64) :: info_str
144 !------------------------------------------------------------------------!
145 ! call InitPhysics from base class
146 CALL this%InitPhysics_base(mesh,config,io,euler,problem_name)
147
148 ! set the total number of variables in a state vector
149 this%VNUM = this%VDIM + 2
150
151 ! ratio of specific heats
152 CALL getattr(config, "gamma", this%gamma, 1.4)
153
154 ! allocate memory for arrays used in euler
155 ALLOCATE(this%pvarname(this%VNUM),this%cvarname(this%VNUM),this%bccsound, &
156 this%fcsound, &
157 stat = err)
158 IF (err.NE.0) &
159 CALL this%Error("InitPhysics_euler", "Unable to allocate memory.")
160
164 this%DENSITY = 1 ! mass density !
165 this%pvarname(this%DENSITY) = "density"
166 this%cvarname(this%DENSITY) = "density"
167 this%XVELOCITY = 2 ! x-velocity !
168 this%XMOMENTUM = 2 ! x-momentum !
169 IF (this%VDIM.GE.2) THEN
170 this%YVELOCITY = 3 ! y-velocity !
171 this%YMOMENTUM = 3 ! y-momentum !
172 ELSE
173 this%YVELOCITY = 0 ! no y-velocity !
174 this%YMOMENTUM = 0 ! no y-momentum !
175 END IF
176 IF (this%VDIM.EQ.3) THEN
177 this%ZVELOCITY = 4 ! z-velocity !
178 this%ZMOMENTUM = 4 ! z-momentum !
179 ELSE
180 this%ZVELOCITY = 0 ! no z-velocity !
181 this%ZMOMENTUM = 0 ! no z-momentum !
182 END IF
183 this%PRESSURE = this%VNUM ! pressure !
184 this%ENERGY = this%VNUM ! total energy !
185 this%pvarname(this%PRESSURE) = "pressure"
186 this%cvarname(this%ENERGY) = "energy"
187
188 ! check which vector components are available and
189 ! set names shown in the data file
190 next_idx = 2
191 IF (btest(mesh%VECTOR_COMPONENTS,0)) THEN
192 this%pvarname(next_idx) = "xvelocity"
193 this%cvarname(next_idx) = "xmomentum"
194 next_idx = next_idx + 1
195 END IF
196 IF (btest(mesh%VECTOR_COMPONENTS,1)) THEN
197 this%pvarname(next_idx) = "yvelocity"
198 this%cvarname(next_idx) = "ymomentum"
199 next_idx = next_idx + 1
200 END IF
201 IF (btest(mesh%VECTOR_COMPONENTS,2)) THEN
202 this%pvarname(next_idx) = "zvelocity"
203 this%cvarname(next_idx) = "zmomentum"
204 END IF
205
206 ! not used in non-isotherml physics
207 this%csiso = 0.0
208
209 ! create new mesh arrays for sound speeds
210 this%bccsound = marray_base()
211 this%fcsound = marray_base(mesh%NFACES)
212
213 ! enable support for absorbing and farfield boundary conditions
214 this%supports_absorbing = .true.
215 this%supports_farfield = .true.
216
217 CALL this%SetOutput(mesh,config,io)
218
219 ! print some information
220 WRITE(info_str,'("spec. heat ratio: ",F4.2)') this%gamma
221 CALL this%Info(repeat(" ",12) // trim(info_str))
222 END SUBROUTINE initphysics
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 IF (ASSOCIATED(new_sv)) THEN
233 DEALLOCATE(new_sv)
234#ifdef DEBUG
235 CALL this%Warning("physics_euler::new_statevector","new statevector already associated")
236#endif
237 END IF
238 ALLOCATE(statevector_euler::new_sv)
239 new_sv = statevector_euler(this,flavour,num)
240 END SUBROUTINE new_statevector
241
243 PURE SUBROUTINE convert2primitive_all(this,cvar,pvar)
244 IMPLICIT NONE
245 !------------------------------------------------------------------------!
246 CLASS(physics_euler), INTENT(IN) :: this
247 CLASS(marray_compound), INTENT(INOUT) :: cvar,pvar
248 !------------------------------------------------------------------------!
249 SELECT TYPE(c => cvar)
250 TYPE IS (statevector_euler)
251 SELECT TYPE(p => pvar)
252 TYPE IS (statevector_euler)
253 IF (c%flavour.EQ.conservative.AND.p%flavour.EQ.primitive) THEN
254 ! conservative -> primitive
255 SELECT CASE(this%VDIM)
256 CASE(1)
257 CALL cons2prim(this%gamma,c%density%data1d(:),c%momentum%data2d(:,1), &
258 c%energy%data1d(:),p%density%data1d(:), &
259 p%velocity%data2d(:,1),p%pressure%data1d(:))
260 CASE(2)
261 CALL cons2prim(this%gamma,c%density%data1d(:),c%momentum%data2d(:,1), &
262 c%momentum%data2d(:,2),c%energy%data1d(:), &
263 p%density%data1d(:),p%velocity%data2d(:,1), &
264 p%velocity%data2d(:,2),p%pressure%data1d(:))
265 CASE(3)
266 CALL cons2prim(this%gamma,c%density%data1d(:),c%momentum%data2d(:,1), &
267 c%momentum%data2d(:,2),c%momentum%data2d(:,3), &
268 c%energy%data1d(:),p%density%data1d(:), &
269 p%velocity%data2d(:,1),p%velocity%data2d(:,2),&
270 p%velocity%data2d(:,3),p%pressure%data1d(:))
271 END SELECT
272 p%fargo_transformation_applied = c%fargo_transformation_applied
273 ELSE
274 ! do nothing
275 END IF
276 END SELECT
277 END SELECT
278 END SUBROUTINE convert2primitive_all
279
281 PURE SUBROUTINE convert2primitive_subset(this,i1,i2,j1,j2,k1,k2,cvar,pvar)
282 IMPLICIT NONE
283 !------------------------------------------------------------------------!
284 CLASS(physics_euler), INTENT(IN) :: this
285 INTEGER, INTENT(IN) :: i1,i2,j1,j2,k1,k2
286 CLASS(marray_compound), INTENT(INOUT) :: cvar,pvar
287 !------------------------------------------------------------------------!
288 SELECT TYPE(c => cvar)
289 TYPE IS (statevector_euler)
290 SELECT TYPE(p => pvar)
291 TYPE IS (statevector_euler)
292 IF (c%flavour.EQ.conservative.AND.p%flavour.EQ.primitive) THEN
293 SELECT CASE (c%density%RANK)
294 CASE(0) ! state vector contains cell center values
295 ! perform the transformation depending on dimensionality
296 SELECT CASE(this%VDIM)
297 CASE(1)
298 CALL cons2prim(this%gamma,c%density%data3d(i1:i2,j1:j2,k1:k2), &
299 c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
300 c%energy%data3d(i1:i2,j1:j2,k1:k2), &
301 p%density%data3d(i1:i2,j1:j2,k1:k2), &
302 p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
303 p%pressure%data3d(i1:i2,j1:j2,k1:k2))
304 CASE(2)
305 CALL cons2prim(this%gamma,c%density%data3d(i1:i2,j1:j2,k1:k2), &
306 c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
307 c%momentum%data4d(i1:i2,j1:j2,k1:k2,2), &
308 c%energy%data3d(i1:i2,j1:j2,k1:k2), &
309 p%density%data3d(i1:i2,j1:j2,k1:k2), &
310 p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
311 p%velocity%data4d(i1:i2,j1:j2,k1:k2,2), &
312 p%pressure%data3d(i1:i2,j1:j2,k1:k2))
313 CASE(3)
314 CALL cons2prim(this%gamma,c%density%data3d(i1:i2,j1:j2,k1:k2), &
315 c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
316 c%momentum%data4d(i1:i2,j1:j2,k1:k2,2), &
317 c%momentum%data4d(i1:i2,j1:j2,k1:k2,3), &
318 c%energy%data3d(i1:i2,j1:j2,k1:k2), &
319 p%density%data3d(i1:i2,j1:j2,k1:k2), &
320 p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
321 p%velocity%data4d(i1:i2,j1:j2,k1:k2,2),&
322 p%velocity%data4d(i1:i2,j1:j2,k1:k2,3), &
323 p%pressure%data3d(i1:i2,j1:j2,k1:k2))
324 END SELECT
325 CASE(1) ! state vector contains cell face / corner values
326 ! perform the transformation depending on dimensionality
327 SELECT CASE(this%VDIM)
328 CASE(1)
329 CALL cons2prim(this%gamma,c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
330 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
331 c%energy%data4d(i1:i2,j1:j2,k1:k2,:), &
332 p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
333 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
334 p%pressure%data4d(i1:i2,j1:j2,k1:k2,:))
335 CASE(2)
336 CALL cons2prim(this%gamma,c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
337 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
338 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,2), &
339 c%energy%data4d(i1:i2,j1:j2,k1:k2,:), &
340 p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
341 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
342 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,2), &
343 p%pressure%data4d(i1:i2,j1:j2,k1:k2,:))
344 CASE(3)
345 CALL cons2prim(this%gamma,c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
346 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
347 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,2), &
348 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,3), &
349 c%energy%data4d(i1:i2,j1:j2,k1:k2,:), &
350 p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
351 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
352 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,2),&
353 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,3), &
354 p%pressure%data4d(i1:i2,j1:j2,k1:k2,:))
355 END SELECT
356 CASE DEFAULT
357 ! do nothing
358 END SELECT
359 ELSE
360 ! do nothing
361 END IF
362 END SELECT
363 END SELECT
364 END SUBROUTINE convert2primitive_subset
365
367 PURE SUBROUTINE convert2conservative_all(this,pvar,cvar)
368 IMPLICIT NONE
369 !------------------------------------------------------------------------!
370 CLASS(physics_euler), INTENT(IN) :: this
371 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
372 !------------------------------------------------------------------------!
373 SELECT TYPE(p => pvar)
374 TYPE IS (statevector_euler)
375 SELECT TYPE(c => cvar)
376 TYPE IS (statevector_euler)
377 IF (p%flavour.EQ.primitive.AND.c%flavour.EQ.conservative) THEN
378 ! perform the transformation depending on dimensionality
379 SELECT CASE(this%VDIM)
380 CASE(1) ! 1D velocity / momentum
381 CALL prim2cons(this%gamma,p%density%data1d(:),p%velocity%data2d(:,1), &
382 p%pressure%data1d(:),c%density%data1d(:), &
383 c%momentum%data2d(:,1),c%energy%data1d(:))
384 CASE(2) ! 2D velocity / momentum
385 CALL prim2cons(this%gamma,p%density%data1d(:),p%velocity%data2d(:,1), &
386 p%velocity%data2d(:,2),p%pressure%data1d(:), &
387 c%density%data1d(:),c%momentum%data2d(:,1), &
388 c%momentum%data2d(:,2),c%energy%data1d(:))
389 CASE(3) ! 3D velocity / momentum
390 CALL prim2cons(this%gamma,p%density%data1d(:),p%velocity%data2d(:,1), &
391 p%velocity%data2d(:,2),p%velocity%data2d(:,3), &
392 p%pressure%data1d(:),c%density%data1d(:), &
393 c%momentum%data2d(:,1),c%momentum%data2d(:,2), &
394 c%momentum%data2d(:,3),c%energy%data1d(:))
395 END SELECT
396 c%fargo_transformation_applied = p%fargo_transformation_applied
397 ELSE
398 ! do nothing
399 END IF
400 END SELECT
401 END SELECT
402 END SUBROUTINE convert2conservative_all
403
405 PURE SUBROUTINE convert2conservative_subset(this,i1,i2,j1,j2,k1,k2,pvar,cvar)
406 IMPLICIT NONE
407 !------------------------------------------------------------------------!
408 CLASS(physics_euler), INTENT(IN) :: this
409 INTEGER, INTENT(IN) :: i1,i2,j1,j2,k1,k2
410 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
411 !------------------------------------------------------------------------!
412 SELECT TYPE(p => pvar)
413 TYPE IS (statevector_euler)
414 SELECT TYPE(c => cvar)
415 TYPE IS (statevector_euler)
416 IF (p%flavour.EQ.primitive.AND.c%flavour.EQ.conservative) THEN
417 SELECT CASE (p%density%RANK)
418 CASE(0) ! state vector contains cell center values
419 ! perform the transformation depending on dimensionality
420 SELECT CASE(this%VDIM)
421 CASE(1) ! 1D velocity / momentum
422 CALL prim2cons(this%gamma,p%density%data3d(i1:i2,j1:j2,k1:k2), &
423 p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
424 p%pressure%data3d(i1:i2,j1:j2,k1:k2), &
425 c%density%data3d(i1:i2,j1:j2,k1:k2), &
426 c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
427 c%energy%data3d(i1:i2,j1:j2,k1:k2))
428 CASE(2) ! 2D velocity / momentum
429 CALL prim2cons(this%gamma,p%density%data3d(i1:i2,j1:j2,k1:k2), &
430 p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
431 p%velocity%data4d(i1:i2,j1:j2,k1:k2,2), &
432 p%pressure%data3d(i1:i2,j1:j2,k1:k2), &
433 c%density%data3d(i1:i2,j1:j2,k1:k2), &
434 c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
435 c%momentum%data4d(i1:i2,j1:j2,k1:k2,2), &
436 c%energy%data3d(i1:i2,j1:j2,k1:k2))
437 CASE(3) ! 3D velocity / momentum
438 CALL prim2cons(this%gamma,p%density%data3d(i1:i2,j1:j2,k1:k2), &
439 p%velocity%data4d(i1:i2,j1:j2,k1:k2,1), &
440 p%velocity%data4d(i1:i2,j1:j2,k1:k2,2), &
441 p%velocity%data4d(i1:i2,j1:j2,k1:k2,3), &
442 p%pressure%data3d(i1:i2,j1:j2,k1:k2), &
443 c%density%data3d(i1:i2,j1:j2,k1:k2), &
444 c%momentum%data4d(i1:i2,j1:j2,k1:k2,1), &
445 c%momentum%data4d(i1:i2,j1:j2,k1:k2,2), &
446 c%momentum%data4d(i1:i2,j1:j2,k1:k2,3), &
447 c%energy%data3d(i1:i2,j1:j2,k1:k2))
448 END SELECT
449 CASE(1) ! state vector contains cell face / corner values
450 ! perform the transformation depending on dimensionality
451 SELECT CASE(this%VDIM)
452 CASE(1) ! 1D velocity / momentum
453 CALL prim2cons(this%gamma,p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
454 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
455 p%pressure%data4d(i1:i2,j1:j2,k1:k2,:), &
456 c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
457 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
458 c%energy%data4d(i1:i2,j1:j2,k1:k2,:))
459 CASE(2) ! 2D velocity / momentum
460 CALL prim2cons(this%gamma,p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
461 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
462 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,2), &
463 p%pressure%data4d(i1:i2,j1:j2,k1:k2,:), &
464 c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
465 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
466 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,2), &
467 c%energy%data4d(i1:i2,j1:j2,k1:k2,:))
468 CASE(3) ! 3D velocity / momentum
469 CALL prim2cons(this%gamma,p%density%data4d(i1:i2,j1:j2,k1:k2,:), &
470 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,1), &
471 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,2), &
472 p%velocity%data5d(i1:i2,j1:j2,k1:k2,:,3), &
473 p%pressure%data4d(i1:i2,j1:j2,k1:k2,:), &
474 c%density%data4d(i1:i2,j1:j2,k1:k2,:), &
475 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,1), &
476 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,2), &
477 c%momentum%data5d(i1:i2,j1:j2,k1:k2,:,3), &
478 c%energy%data4d(i1:i2,j1:j2,k1:k2,:))
479 END SELECT
480 CASE DEFAULT
481 ! do nothing
482 END SELECT
483 ELSE
484 ! do nothing
485 END IF
486 END SELECT
487 END SELECT
488 END SUBROUTINE convert2conservative_subset
489
491 !\todo NOT VERIFIED
492 PURE SUBROUTINE calcfluxesx(this,Mesh,nmin,nmax,prim,cons,xfluxes)
493 IMPLICIT NONE
494 !------------------------------------------------------------------------!
495 CLASS(physics_euler), INTENT(IN) :: this
496 CLASS(mesh_base), INTENT(IN) :: mesh
497 INTEGER, INTENT(IN) :: nmin,nmax
498 CLASS(marray_compound), INTENT(INOUT) :: prim,cons,xfluxes
499 !------------------------------------------------------------------------!
500 SELECT TYPE(p => prim)
501 TYPE IS(statevector_euler)
502 SELECT TYPE(c => cons)
503 TYPE IS(statevector_euler)
504 SELECT TYPE(f => xfluxes)
505 TYPE IS(statevector_euler)
506 SELECT CASE(this%VDIM)
507 CASE(1) ! 1D flux
508 CALL setflux(p%density%data2d(:,nmin:nmax), &
509 p%velocity%data3d(:,nmin:nmax,1), &
510 p%pressure%data2d(:,nmin:nmax), &
511 c%momentum%data3d(:,nmin:nmax,1), &
512 c%energy%data2d(:,nmin:nmax), &
513 f%density%data2d(:,nmin:nmax), &
514 f%momentum%data3d(:,nmin:nmax,1), &
515 f%energy%data2d(:,nmin:nmax))
516 CASE(2) ! 2D flux
517 CALL setflux(p%density%data2d(:,nmin:nmax), &
518 p%velocity%data3d(:,nmin:nmax,1), &
519 p%pressure%data2d(:,nmin:nmax), &
520 c%momentum%data3d(:,nmin:nmax,1), &
521 c%momentum%data3d(:,nmin:nmax,2), &
522 c%energy%data2d(:,nmin:nmax), &
523 f%density%data2d(:,nmin:nmax), &
524 f%momentum%data3d(:,nmin:nmax,1), &
525 f%momentum%data3d(:,nmin:nmax,2), &
526 f%energy%data2d(:,nmin:nmax))
527 CASE(3) ! 3D flux
528 CALL setflux(p%density%data2d(:,nmin:nmax), &
529 p%velocity%data3d(:,nmin:nmax,1), &
530 p%pressure%data2d(:,nmin:nmax), &
531 c%momentum%data3d(:,nmin:nmax,1), &
532 c%momentum%data3d(:,nmin:nmax,2), &
533 c%momentum%data3d(:,nmin:nmax,3), &
534 c%energy%data2d(:,nmin:nmax), &
535 f%density%data2d(:,nmin:nmax), &
536 f%momentum%data3d(:,nmin:nmax,1), &
537 f%momentum%data3d(:,nmin:nmax,2), &
538 f%momentum%data3d(:,nmin:nmax,3), &
539 f%energy%data2d(:,nmin:nmax))
540 END SELECT
541 END SELECT
542 END SELECT
543 END SELECT
544 END SUBROUTINE calcfluxesx
545
547 !\todo NOT VERIFIED
548 PURE SUBROUTINE calcfluxesy(this,Mesh,nmin,nmax,prim,cons,yfluxes)
549 IMPLICIT NONE
550 !------------------------------------------------------------------------!
551 CLASS(physics_euler), INTENT(IN) :: this
552 CLASS(mesh_base), INTENT(IN) :: mesh
553 INTEGER, INTENT(IN) :: nmin,nmax
554 CLASS(marray_compound), INTENT(INOUT) :: prim,cons,yfluxes
555 !------------------------------------------------------------------------!
556 SELECT TYPE(p => prim)
557 TYPE IS(statevector_euler)
558 SELECT TYPE(c => cons)
559 TYPE IS(statevector_euler)
560 SELECT TYPE(f => yfluxes)
561 TYPE IS(statevector_euler)
562 SELECT CASE(this%VDIM)
563 CASE(1) ! 1D flux
564 CALL setflux(p%density%data2d(:,nmin:nmax), &
565 p%velocity%data3d(:,nmin:nmax,1), &
566 p%pressure%data2d(:,nmin:nmax), &
567 c%momentum%data3d(:,nmin:nmax,1), &
568 c%energy%data2d(:,nmin:nmax), &
569 f%density%data2d(:,nmin:nmax), &
570 f%momentum%data3d(:,nmin:nmax,1), &
571 f%energy%data2d(:,nmin:nmax))
572 CASE(2) ! 2D flux
573 CALL setflux(p%density%data2d(:,nmin:nmax), &
574 p%velocity%data3d(:,nmin:nmax,2), &
575 p%pressure%data2d(:,nmin:nmax), &
576 c%momentum%data3d(:,nmin:nmax,2), &
577 c%momentum%data3d(:,nmin:nmax,1), &
578 c%energy%data2d(:,nmin:nmax), &
579 f%density%data2d(:,nmin:nmax), &
580 f%momentum%data3d(:,nmin:nmax,2), &
581 f%momentum%data3d(:,nmin:nmax,1), &
582 f%energy%data2d(:,nmin:nmax))
583 CASE(3) ! 3D flux
584 CALL setflux(p%density%data2d(:,nmin:nmax), &
585 p%velocity%data3d(:,nmin:nmax,2), &
586 p%pressure%data2d(:,nmin:nmax), &
587 c%momentum%data3d(:,nmin:nmax,2), &
588 c%momentum%data3d(:,nmin:nmax,1), &
589 c%momentum%data3d(:,nmin:nmax,3), &
590 c%energy%data2d(:,nmin:nmax), &
591 f%density%data2d(:,nmin:nmax), &
592 f%momentum%data3d(:,nmin:nmax,2), &
593 f%momentum%data3d(:,nmin:nmax,1), &
594 f%momentum%data3d(:,nmin:nmax,3), &
595 f%energy%data2d(:,nmin:nmax))
596 END SELECT
597 END SELECT
598 END SELECT
599 END SELECT
600 END SUBROUTINE calcfluxesy
601
603 !\todo NOT VERIFIED
604 PURE SUBROUTINE calcfluxesz(this,Mesh,nmin,nmax,prim,cons,zfluxes)
605 IMPLICIT NONE
606 !------------------------------------------------------------------------!
607 CLASS(physics_euler), INTENT(IN) :: this
608 CLASS(mesh_base), INTENT(IN) :: mesh
609 INTEGER, INTENT(IN) :: nmin,nmax
610 CLASS(marray_compound), INTENT(INOUT) :: prim,cons,zfluxes
611 !------------------------------------------------------------------------!
612 SELECT TYPE(p => prim)
613 TYPE IS(statevector_euler)
614 SELECT TYPE(c => cons)
615 TYPE IS(statevector_euler)
616 SELECT TYPE(f => zfluxes)
617 TYPE IS(statevector_euler)
618 SELECT CASE(this%VDIM)
619 CASE(1) ! 1D flux
620 CALL setflux(p%density%data2d(:,nmin:nmax), &
621 p%velocity%data3d(:,nmin:nmax,1), &
622 p%pressure%data2d(:,nmin:nmax), &
623 c%momentum%data3d(:,nmin:nmax,1), &
624 c%energy%data2d(:,nmin:nmax), &
625 f%density%data2d(:,nmin:nmax), &
626 f%momentum%data3d(:,nmin:nmax,1), &
627 f%energy%data2d(:,nmin:nmax))
628 CASE(2) ! 2D flux
629 CALL setflux(p%density%data2d(:,nmin:nmax), &
630 p%velocity%data3d(:,nmin:nmax,2), &
631 p%pressure%data2d(:,nmin:nmax), &
632 c%momentum%data3d(:,nmin:nmax,2), &
633 c%momentum%data3d(:,nmin:nmax,1), &
634 c%energy%data2d(:,nmin:nmax), &
635 f%density%data2d(:,nmin:nmax), &
636 f%momentum%data3d(:,nmin:nmax,2), &
637 f%momentum%data3d(:,nmin:nmax,1), &
638 f%energy%data2d(:,nmin:nmax))
639 CASE(3) ! 3D flux
640 CALL setflux(p%density%data2d(:,nmin:nmax), &
641 p%velocity%data3d(:,nmin:nmax,3), &
642 p%pressure%data2d(:,nmin:nmax), &
643 c%momentum%data3d(:,nmin:nmax,3), &
644 c%momentum%data3d(:,nmin:nmax,1), &
645 c%momentum%data3d(:,nmin:nmax,2), &
646 c%energy%data2d(:,nmin:nmax), &
647 f%density%data2d(:,nmin:nmax), &
648 f%momentum%data3d(:,nmin:nmax,3), &
649 f%momentum%data3d(:,nmin:nmax,1), &
650 f%momentum%data3d(:,nmin:nmax,2), &
651 f%energy%data2d(:,nmin:nmax))
652 END SELECT
653 END SELECT
654 END SELECT
655 END SELECT
656 END SUBROUTINE calcfluxesz
657
658! !> Reconstruction of the intermediate state for HLLC
659! !\todo NOT VERIFIED
660! PURE SUBROUTINE CalcIntermediateStateX(this,Mesh,prim,cons, &
661! amin,amax,cstar,astar)
662! IMPLICIT NONE
663! !------------------------------------------------------------------------!
664! CLASS(physics_euler), INTENT(IN) :: this
665! CLASS(mesh_base), INTENT(IN) :: Mesh
666! REAL, INTENT(IN), &
667! DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Mesh%NFACES,this%VNUM) &
668! :: prim,cons
669! REAL, INTENT(OUT), &
670! DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,this%VNUM) &
671! :: cstar
672! REAL, INTENT(IN), &
673! DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
674! :: amin,amax
675! REAL, INTENT(OUT), &
676! DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
677! :: astar
678! !------------------------------------------------------------------------!
679! INTEGER :: i,j,k
680! !------------------------------------------------------------------------!
681! DO k=Mesh%KMIN,Mesh%KMAX
682! DO j=Mesh%JMIN,Mesh%JMAX
683!!NEC$ IVDEP
684! DO i=Mesh%IMIN-1,Mesh%IMAX
685! CALL SetIntermediateState( &
686! prim(i,j,k,2,this%DENSITY),prim(i+1,j,k,1,this%DENSITY), &
687! prim(i,j,k,2,this%XVELOCITY),prim(i+1,j,k,1,this%XVELOCITY), &
688! prim(i,j,k,2,this%YVELOCITY),prim(i+1,j,k,1,this%YVELOCITY), &
689! prim(i,j,k,2,this%PRESSURE),prim(i+1,j,k,1,this%PRESSURE), &
690! cons(i,j,k,2,this%ENERGY),cons(i+1,j,k,1,this%ENERGY), &
691! amin(i,j,k),amax(i,j,k), &
692! cstar(i,j,k,this%DENSITY),cstar(i,j,k,this%XMOMENTUM), &
693! cstar(i,j,k,this%YMOMENTUM),cstar(i,j,k,this%ENERGY), &
694! astar(i,j,k))
695! END DO
696! END DO
697! END DO
698! END SUBROUTINE CalcIntermediateStateX
699!
700!
701! !> Reconstruction of the intermediate state for HLLC
702! !\todo NOT VERIFIED
703! PURE SUBROUTINE CalcIntermediateStateY(this,Mesh,prim,cons, &
704! bmin,bmax,cstar,bstar)
705! IMPLICIT NONE
706! !------------------------------------------------------------------------!
707! CLASS(physics_euler), INTENT(INOUT) :: this
708! CLASS(mesh_base), INTENT(IN) :: Mesh
709! REAL, INTENT(IN), &
710! DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Mesh%NFACES,this%VNUM) &
711! :: prim,cons
712! REAL, INTENT(OUT), &
713! DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,this%VNUM) &
714! :: cstar
715! REAL, INTENT(IN), &
716! DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
717! :: bmin,bmax
718! REAL, INTENT(OUT), &
719! DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
720! :: bstar
721! !------------------------------------------------------------------------!
722! INTEGER :: i,j,k
723! !------------------------------------------------------------------------!
724! DO k=Mesh%KGMIN,Mesh%KGMAX
725! DO j=Mesh%JMIN-1,Mesh%JMAX
726!!NEC$ IVDEP
727! DO i=Mesh%IGMIN,Mesh%IGMAX
728! CALL SetIntermediateState( &
729! prim(i,j,k,4,this%DENSITY),prim(i,j+1,k,3,this%DENSITY), &
730! prim(i,j,k,4,this%YVELOCITY),prim(i,j+1,k,3,this%YVELOCITY), &
731! prim(i,j,k,4,this%XVELOCITY),prim(i,j+1,k,3,this%XVELOCITY), &
732! prim(i,j,k,4,this%PRESSURE),prim(i,j+1,k,3,this%PRESSURE), &
733! cons(i,j,k,4,this%ENERGY),cons(i,j+1,k,3,this%ENERGY), &
734! bmin(i,j,k),bmax(i,j,k), &
735! cstar(i,j,k,this%DENSITY),cstar(i,j,k,this%YMOMENTUM), &
736! cstar(i,j,k,this%XMOMENTUM),cstar(i,j,k,this%ENERGY), &
737! bstar(i,j,k))
738! END DO
739! END DO
740! END DO
741! END SUBROUTINE CalcIntermediateStateY
742
743 PURE SUBROUTINE calculatecharsystemx(this,Mesh,i1,i2,pvar,lambda,xvar)
744 IMPLICIT NONE
745 !------------------------------------------------------------------------!
746 CLASS(physics_euler), INTENT(IN) :: this
747 CLASS(mesh_base), INTENT(IN) :: mesh
748 INTEGER, INTENT(IN) :: i1,i2
749 CLASS(marray_compound), INTENT(INOUT) :: pvar
750 REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM), &
751 INTENT(OUT) :: lambda,xvar
752 !------------------------------------------------------------------------!
753 INTEGER :: il,ir
754 !------------------------------------------------------------------------!
755 SELECT TYPE(p => pvar)
756 TYPE IS(statevector_euler)
757 il = min(i1,i2)
758 ir = max(i1,i2)
759 SELECT CASE(this%VDIM)
760 CASE(1) ! 1D
761 ! compute eigenvalues at i
762 CALL seteigenvalues1d(this%gamma, &
763 p%density%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
764 p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
765 p%pressure%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
766 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
767 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
768 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3))
769 ! compute characteristic variables using cell mean values of adjacent
770 ! cells to calculate derivatives and the isothermal speed of sound
771 ! at the intermediate cell face
772 CALL setcharvars1d(this%gamma, &
773 p%density%data3d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
774 p%density%data3d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
775 p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
776 p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
777 p%pressure%data3d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
778 p%pressure%data3d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
779 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
780 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
781 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
782 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
783 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3))
784 CASE(2) ! 2D
785 ! compute eigenvalues at i1
786 CALL seteigenvalues2d(this%gamma, &
787 p%density%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
788 p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
789 p%pressure%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
790 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
791 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
792 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
793 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4))
794 ! compute characteristic variables using cell mean values of adjacent
795 ! cells to calculate derivatives and the isothermal speed of sound
796 ! at the intermediate cell face
797 CALL setcharvars2d(this%gamma, &
798 p%density%data3d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
799 p%density%data3d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
800 p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
801 p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
802 p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
803 p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
804 p%pressure%data3d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
805 p%pressure%data3d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
806 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
807 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
808 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
809 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
810 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
811 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4))
812 CASE(3) ! 3D
813 ! compute eigenvalues at i1
814 CALL seteigenvalues3d(this%gamma, &
815 p%density%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
816 p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
817 p%pressure%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
818 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
819 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
820 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
821 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
822 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,5))
823 ! compute characteristic variables using cell mean values of adjacent
824 ! cells to calculate derivatives and the isothermal speed of sound
825 ! at the intermediate cell face
826 CALL setcharvars3d(this%gamma, &
827 p%density%data3d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
828 p%density%data3d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
829 p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
830 p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
831 p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
832 p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
833 p%velocity%data4d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
834 p%velocity%data4d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
835 p%pressure%data3d(il,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
836 p%pressure%data3d(ir,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
837 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
838 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,5), &
839 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
840 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
841 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
842 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
843 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,5))
844 END SELECT
845 END SELECT
846 END SUBROUTINE calculatecharsystemx
847
848 PURE SUBROUTINE calculatecharsystemy(this,Mesh,j1,j2,pvar,lambda,xvar)
849 IMPLICIT NONE
850 !------------------------------------------------------------------------!
851 CLASS(physics_euler), INTENT(IN) :: this
852 CLASS(mesh_base), INTENT(IN) :: mesh
853 INTEGER, INTENT(IN) :: j1,j2
854 CLASS(marray_compound), INTENT(INOUT) :: pvar
855 REAL, DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM), &
856 INTENT(OUT) :: lambda,xvar
857 !------------------------------------------------------------------------!
858 INTEGER :: jl,jr,vn,vt
859 !------------------------------------------------------------------------!
860 SELECT TYPE(p => pvar)
861 TYPE IS(statevector_euler)
862 jl = min(j1,j2)
863 jr = max(j1,j2)
864 SELECT CASE(this%VDIM)
865 CASE(1) ! 1D
866 ! compute eigenvalues at j
867 CALL seteigenvalues1d(this%gamma, &
868 p%density%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
869 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,2), &
870 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
871 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
872 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
873 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3))
874 ! compute characteristic variables using cell mean values of adjacent
875 ! cells to calculate derivatives and the isothermal speed of sound
876 ! at the intermediate cell face
877 CALL setcharvars1d(this%gamma, &
878 p%density%data3d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX), &
879 p%density%data3d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX), &
880 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,1), &
881 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,1), &
882 p%pressure%data3d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX), &
883 p%pressure%data3d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX), &
884 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
885 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
886 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
887 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
888 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3))
889 CASE(2) ! 2D
890 ! check which velocity component is along y-direction
891 SELECT CASE(mesh%VECTOR_COMPONENTS)
892 CASE(ior(vector_x,vector_y)) ! 2D velocities in x-y-plane
893 vt = 1
894 vn = 2
895 CASE(ior(vector_y,vector_z)) ! 2D velocities in y-z-plane
896 vt = 2
897 vn = 1
898 CASE DEFAULT
899 ! this should not happen
900 RETURN
901 END SELECT
902 ! compute eigenvalues at j
903 CALL seteigenvalues2d(this%gamma, &
904 p%density%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
905 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,vn), &
906 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
907 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
908 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
909 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
910 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4))
911 ! compute characteristic variables using cell mean values of adjacent
912 ! cells to calculate derivatives and the isothermal speed of sound
913 ! at the intermediate cell face
914 CALL setcharvars2d(this%gamma, &
915 p%density%data3d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX), &
916 p%density%data3d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX), &
917 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,vn), &
918 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,vn), &
919 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,vt), &
920 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,vt), &
921 p%pressure%data3d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX), &
922 p%pressure%data3d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX), &
923 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
924 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
925 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
926 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
927 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
928 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4))
929 CASE(3) ! 3D
930 ! compute eigenvalues at j
931 CALL seteigenvalues3d(this%gamma, &
932 p%density%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
933 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,2), &
934 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
935 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
936 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
937 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
938 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
939 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,5))
940 ! compute characteristic variables using cell mean values of adjacent
941 ! cells to calculate derivatives and the isothermal speed of sound
942 ! at the intermediate cell face
943 CALL setcharvars3d(this%gamma, &
944 p%density%data3d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX), &
945 p%density%data3d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX), &
946 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,2), &
947 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,2), &
948 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,1), &
949 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,1), &
950 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX,3), &
951 p%velocity%data4d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX,3), &
952 p%pressure%data3d(mesh%IMIN:mesh%IMAX,jl,mesh%KMIN:mesh%KMAX), &
953 p%pressure%data3d(mesh%IMIN:mesh%IMAX,jr,mesh%KMIN:mesh%KMAX), &
954 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
955 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,5), &
956 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
957 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
958 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
959 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
960 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,5))
961 END SELECT
962 END SELECT
963 END SUBROUTINE calculatecharsystemy
964
965
966 PURE SUBROUTINE calculatecharsystemz(this,Mesh,k1,k2,pvar,lambda,xvar)
967 IMPLICIT NONE
968 !------------------------------------------------------------------------!
969 CLASS(physics_euler), INTENT(IN) :: this
970 CLASS(mesh_base), INTENT(IN) :: mesh
971 INTEGER, INTENT(IN) :: k1,k2
972 CLASS(marray_compound), INTENT(INOUT) :: pvar
973 REAL, DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,this%VNUM), &
974 INTENT(OUT) :: lambda,xvar
975 !------------------------------------------------------------------------!
976 INTEGER :: kl,kr
977 !------------------------------------------------------------------------!
978 SELECT TYPE(p => pvar)
979 TYPE IS(statevector_euler)
980 kl = min(k1,k2)
981 kr = max(k1,k2)
982 SELECT CASE(this%VDIM)
983 CASE(1) ! 1D
984 ! compute eigenvalues at k
985 CALL seteigenvalues1d(this%gamma, &
986 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
987 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,3), &
988 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
989 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
990 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
991 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3))
992 ! compute characteristic variables using cell mean values of adjacent
993 ! cells to calculate derivatives and the isothermal speed of sound
994 ! at the intermediate cell face
995 CALL setcharvars1d(this%gamma, &
996 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl), &
997 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr), &
998 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,1), &
999 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,1), &
1000 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl), &
1001 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr), &
1002 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1003 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1004 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1005 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1006 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3))
1007 CASE(2) ! 2D
1008 ! compute eigenvalues at k
1009 CALL seteigenvalues2d(this%gamma, &
1010 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1011 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,3), &
1012 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1013 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1014 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1015 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1016 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4))
1017 ! compute characteristic variables using cell mean values of adjacent
1018 ! cells to calculate derivatives and the isothermal speed of sound
1019 ! at the intermediate cell face
1020 CALL setcharvars2d(this%gamma, &
1021 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl), &
1022 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr), &
1023 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,2), & ! 2nd component is vz
1024 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,2), &
1025 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,1), & ! 1st component: vx or vy
1026 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,1), &
1027 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl), &
1028 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr), &
1029 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1030 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1031 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1032 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1033 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1034 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4))
1035 CASE(3) ! 3D
1036 ! compute eigenvalues at k
1037 CALL seteigenvalues3d(this%gamma, &
1038 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1039 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,3), &
1040 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1041 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1042 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1043 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1044 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1045 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,5))
1046 ! compute characteristic variables using cell mean values of adjacent
1047 ! cells to calculate derivatives and the isothermal speed of sound
1048 ! at the intermediate cell face
1049 CALL setcharvars3d(this%gamma, &
1050 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl), &
1051 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr), &
1052 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,3), &
1053 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,3), &
1054 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,1), &
1055 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,1), &
1056 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl,2), &
1057 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr,2), &
1058 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kl), &
1059 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,kr), &
1060 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1061 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,5), &
1062 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1063 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1064 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1065 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1066 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,5))
1067 END SELECT
1068 END SELECT
1069 END SUBROUTINE calculatecharsystemz
1070
1072 PURE SUBROUTINE calculateboundarydatax(this,Mesh,i1,i2,xvar,pvar)
1073 IMPLICIT NONE
1074 !------------------------------------------------------------------------!
1075 CLASS(physics_euler), INTENT(IN) :: this
1076 CLASS(mesh_base), INTENT(IN) :: mesh
1077 INTEGER, INTENT(IN) :: i1,i2
1078 REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM), &
1079 INTENT(IN) :: xvar
1080 CLASS(marray_compound), INTENT(INOUT) :: pvar
1081 !------------------------------------------------------------------------!
1082 SELECT TYPE(p => pvar)
1083 TYPE IS(statevector_euler)
1084 SELECT CASE(this%VDIM)
1085 CASE(1) ! 1D
1086 CALL setboundarydata1d(i2-i1,this%gamma, &
1087 p%density%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1088 p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1089 p%pressure%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1090 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1091 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1092 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1093 p%density%data3d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1094 p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1095 p%pressure%data3d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
1096 CASE(2) ! 2D
1097 CALL setboundarydata2d(i2-i1,this%gamma, &
1098 p%density%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1099 p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1100 p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1101 p%pressure%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1102 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1103 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1104 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1105 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
1106 p%density%data3d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1107 p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1108 p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1109 p%pressure%data3d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
1110 CASE(3) ! 3D
1111 CALL setboundarydata3d(i2-i1,this%gamma, &
1112 p%density%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1113 p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1114 p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1115 p%velocity%data4d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1116 p%pressure%data3d(i1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1117 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1118 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1119 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1120 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
1121 xvar(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,5), &
1122 p%density%data3d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1123 p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1124 p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1125 p%velocity%data4d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1126 p%pressure%data3d(i2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
1127 END SELECT
1128 END SELECT
1129 END SUBROUTINE calculateboundarydatax
1130
1131 PURE SUBROUTINE calculateboundarydatay(this,Mesh,j1,j2,xvar,pvar)
1132 IMPLICIT NONE
1133 !------------------------------------------------------------------------!
1134 CLASS(physics_euler), INTENT(IN) :: this
1135 CLASS(mesh_base), INTENT(IN) :: mesh
1136 INTEGER, INTENT(IN) :: j1,j2
1137 REAL, DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM), &
1138 INTENT(IN) :: xvar
1139 CLASS(marray_compound), INTENT(INOUT) :: pvar
1140 !------------------------------------------------------------------------!
1141 INTEGER :: vt,vn
1142 !------------------------------------------------------------------------!
1143 SELECT TYPE(p => pvar)
1144 TYPE IS(statevector_euler)
1145 SELECT CASE(this%VDIM)
1146 CASE(1) ! 1D
1147 CALL setboundarydata1d(j2-j1,this%gamma, &
1148 p%density%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1149 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,1), &
1150 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1151 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1152 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1153 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1154 p%density%data3d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX), &
1155 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,1), &
1156 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX))
1157 CASE(2) ! 2D
1158 ! check which velocity component is along y-direction
1159 SELECT CASE(mesh%VECTOR_COMPONENTS)
1160 CASE(ior(vector_x,vector_y)) ! 2D velocities in x-y-plane
1161 vt = 1
1162 vn = 2
1163 CASE(ior(vector_y,vector_z)) ! 2D velocities in y-z-plane
1164 vt = 2
1165 vn = 1
1166 CASE DEFAULT
1167 ! this should not happen
1168 RETURN
1169 END SELECT
1170 CALL setboundarydata2d(j2-j1,this%gamma, &
1171 p%density%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1172 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,vn), &
1173 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,vt), &
1174 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1175 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1176 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1177 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1178 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
1179 p%density%data3d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX), &
1180 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,vn), &
1181 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,vt), &
1182 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX))
1183 CASE(3) ! 3D
1184 CALL setboundarydata3d(j2-j1,this%gamma, &
1185 p%density%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1186 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,2), &
1187 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,1), &
1188 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX,3), &
1189 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j1,mesh%KMIN:mesh%KMAX), &
1190 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1191 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1192 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1193 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
1194 xvar(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,5), &
1195 p%density%data3d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX), &
1196 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,2), &
1197 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,1), &
1198 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX,3), &
1199 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j2,mesh%KMIN:mesh%KMAX))
1200 END SELECT
1201 END SELECT
1202 END SUBROUTINE calculateboundarydatay
1203
1204 PURE SUBROUTINE calculateboundarydataz(this,Mesh,k1,k2,xvar,pvar)
1205 IMPLICIT NONE
1206 !------------------------------------------------------------------------!
1207 CLASS(physics_euler), INTENT(IN) :: this
1208 CLASS(mesh_base), INTENT(IN) :: mesh
1209 INTEGER, INTENT(IN) :: k1,k2
1210 REAL, DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,this%VNUM), &
1211 INTENT(IN) :: xvar
1212 CLASS(marray_compound), INTENT(INOUT) :: pvar
1213 !------------------------------------------------------------------------!
1214 SELECT TYPE(p => pvar)
1215 TYPE IS(statevector_euler)
1216 SELECT CASE(this%VDIM)
1217 CASE(1) ! 1D
1218 CALL setboundarydata1d(k2-k1,this%gamma, &
1219 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1220 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,1), &
1221 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1222 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1223 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1224 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1225 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2), &
1226 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,1), &
1227 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2))
1228 CASE(2) ! 2D
1229 CALL setboundarydata2d(k2-k1,this%gamma, &
1230 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1231 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,2), &
1232 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,1), &
1233 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1234 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1235 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1236 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1237 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1238 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2), &
1239 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,2), &
1240 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,1), &
1241 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2))
1242 CASE(3) ! 3D
1243 CALL setboundarydata3d(k2-k1,this%gamma, &
1244 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1245 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,3), &
1246 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,1), &
1247 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1,2), &
1248 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k1), &
1249 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1250 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1251 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1252 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1253 xvar(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,5), &
1254 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2), &
1255 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,3), &
1256 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,1), &
1257 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2,2), &
1258 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k2))
1259 END SELECT
1260 END SELECT
1261 END SUBROUTINE calculateboundarydataz
1262
1263
1265 !\todo NOT VERIFIED
1266 PURE SUBROUTINE calculateprim2riemannx(this,Mesh,i,pvar,lambda,Rinv)
1267 IMPLICIT NONE
1268 !------------------------------------------------------------------------!
1269 CLASS(physics_euler), INTENT(IN) :: this
1270 CLASS(mesh_base), INTENT(IN) :: mesh
1271 INTEGER, INTENT(IN) :: i
1272 CLASS(marray_compound), INTENT(IN) :: pvar
1273 REAL, INTENT(OUT), &
1274 DIMENSION(Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM) &
1275 :: lambda
1276 REAL, INTENT(OUT), &
1277 DIMENSION(Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM) &
1278 :: rinv
1279 !------------------------------------------------------------------------!
1280 SELECT TYPE(p => pvar)
1281 TYPE IS(statevector_euler)
1282 SELECT CASE(this%VDIM)
1283 CASE(1) ! 1D
1284 CALL seteigenvalues1d(this%gamma, &
1285 p%density%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1286 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1287 p%pressure%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX),&
1288 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1289 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1290 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3))
1291 ! compute Riemann invariants
1292 CALL prim2riemann1d(this%gamma, &
1293 p%density%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1294 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1295 p%pressure%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1296 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1297 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1298 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1299 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1300 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1301 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3))
1302 CASE(2) ! 2D
1303 CALL seteigenvalues2d(this%gamma, &
1304 p%density%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1305 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1306 p%pressure%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX),&
1307 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1308 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1309 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1310 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4))
1311 ! compute Riemann invariants
1312 CALL prim2riemann2d(this%gamma, &
1313 p%density%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1314 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1315 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1316 p%pressure%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1317 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1318 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1319 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1320 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
1321 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1322 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1323 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1324 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4))
1325 CASE(3) ! 3D
1326 ! compute eigenvalues at i
1327 CALL seteigenvalues3d(this%gamma, &
1328 p%density%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1329 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1330 p%pressure%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX),&
1331 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1332 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1333 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1334 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
1335 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,5))
1336 ! compute Riemann invariants
1337 CALL prim2riemann3d(this%gamma, &
1338 p%density%data3d( i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1339 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1340 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1341 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1342 p%pressure%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1343 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1344 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1345 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1346 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
1347 lambda(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,5), &
1348 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1349 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1350 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1351 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
1352 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,5))
1353 END SELECT
1354 END SELECT
1355 END SUBROUTINE calculateprim2riemannx
1356
1357
1359 !\todo NOT VERIFIED
1360 PURE SUBROUTINE calculateprim2riemanny(this,Mesh,j,pvar,lambda,Rinv)
1361 IMPLICIT NONE
1362 !------------------------------------------------------------------------!
1363 CLASS(physics_euler), INTENT(IN) :: this
1364 CLASS(mesh_base), INTENT(IN) :: mesh
1365 INTEGER, INTENT(IN) :: j
1366 CLASS(marray_compound), INTENT(IN) :: pvar
1367 REAL, INTENT(OUT), &
1368 DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM) &
1369 :: lambda
1370 REAL, INTENT(OUT), &
1371 DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM) &
1372 :: rinv
1373 !------------------------------------------------------------------------!
1374 SELECT TYPE(p => pvar)
1375 TYPE IS(statevector_euler)
1376 SELECT CASE(this%VDIM)
1377 CASE(1) ! 1D
1378 ! compute eigenvalues at j
1379 CALL seteigenvalues1d(this%gamma, &
1380 p%density%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1381 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,1), &
1382 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX),&
1383 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1384 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1385 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3))
1386 ! compute Riemann invariants
1387 CALL prim2riemann1d(this%gamma, &
1388 p%density%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1389 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,1), &
1390 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1391 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1392 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1393 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1394 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1395 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1396 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3))
1397 CASE(2) ! 2D
1398 ! compute eigenvalues at j
1399 CALL seteigenvalues2d(this%gamma, &
1400 p%density%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1401 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,2), &
1402 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX),&
1403 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1404 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1405 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1406 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4))
1407 ! compute Riemann invariants
1408 CALL prim2riemann2d(this%gamma, &
1409 p%density%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1410 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,2), &
1411 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,1), &
1412 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1413 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1414 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1415 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1416 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
1417 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1418 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1419 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1420 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4))
1421 CASE(3) ! 3D
1422 ! compute eigenvalues at j
1423 CALL seteigenvalues3d(this%gamma, &
1424 p%density%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1425 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,2), &
1426 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX),&
1427 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1428 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1429 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1430 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
1431 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,5))
1432 ! compute Riemann invariants
1433 CALL prim2riemann3d(this%gamma, &
1434 p%density%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1435 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,2), &
1436 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,3), &
1437 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,1), &
1438 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1439 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1440 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1441 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1442 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
1443 lambda(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,5), &
1444 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1445 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1446 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1447 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
1448 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,5))
1449 END SELECT
1450 END SELECT
1451 END SUBROUTINE calculateprim2riemanny
1452
1453
1455 !\todo NOT VERIFIED
1456 PURE SUBROUTINE calculateprim2riemannz(this,Mesh,k,pvar,lambda,Rinv)
1457 IMPLICIT NONE
1458 !------------------------------------------------------------------------!
1459 CLASS(physics_euler), INTENT(IN) :: this
1460 CLASS(mesh_base), INTENT(IN) :: mesh
1461 INTEGER, INTENT(IN) :: k
1462 CLASS(marray_compound), INTENT(IN) :: pvar
1463 REAL, INTENT(OUT), &
1464 DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,this%VNUM) &
1465 :: lambda
1466 REAL, INTENT(OUT), &
1467 DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,this%VNUM) &
1468 :: rinv
1469 !------------------------------------------------------------------------!
1470 SELECT TYPE(p => pvar)
1471 TYPE IS(statevector_euler)
1472 SELECT CASE(this%VDIM)
1473 CASE(1) ! 1D
1474 ! compute eigenvalues at k
1475 CALL seteigenvalues1d(this%gamma, &
1476 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1477 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,1), &
1478 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k),&
1479 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1480 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1481 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3))
1482 ! compute Riemann invariants
1483 CALL prim2riemann1d(this%gamma, &
1484 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1485 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,1), &
1486 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1487 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1488 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1489 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1490 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1491 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1492 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3))
1493 CASE(2) ! 2D
1494 ! compute eigenvalues at k
1495 CALL seteigenvalues2d(this%gamma, &
1496 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1497 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,2), &
1498 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k),&
1499 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1500 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1501 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1502 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4))
1503 ! compute Riemann invariants
1504 CALL prim2riemann2d(this%gamma, &
1505 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1506 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,2), &
1507 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,1), &
1508 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1509 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1510 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1511 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1512 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1513 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1514 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1515 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1516 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4))
1517
1518 CASE(3) ! 3D
1519 ! compute eigenvalues at k
1520 CALL seteigenvalues3d(this%gamma, &
1521 p%density%data3d( mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1522 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,3), &
1523 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k),&
1524 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1525 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1526 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1527 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1528 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,5))
1529
1530 ! compute Riemann invariants
1531 CALL prim2riemann3d(this%gamma, &
1532 p%density%data3d( mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1533 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,3), &
1534 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,1), &
1535 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,2), &
1536 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1537 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1538 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1539 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1540 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1541 lambda(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,5), &
1542 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1543 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1544 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1545 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1546 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,5))
1547
1548 END SELECT
1549 END SELECT
1550 END SUBROUTINE calculateprim2riemannz
1551
1552
1554 !\todo NOT VERIFIED
1555 PURE SUBROUTINE calculateriemann2primx(this,Mesh,i,Rinv,pvar)
1556 IMPLICIT NONE
1557 !------------------------------------------------------------------------!
1558 CLASS(physics_euler), INTENT(IN) :: this
1559 CLASS(mesh_base), INTENT(IN) :: mesh
1560 INTEGER, INTENT(IN) :: i
1561 REAL, INTENT(IN), &
1562 DIMENSION(Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM) &
1563 :: rinv
1564 CLASS(marray_compound), INTENT(INOUT) :: pvar
1565 !------------------------------------------------------------------------!
1566 SELECT TYPE(p => pvar)
1567 TYPE IS(statevector_euler)
1568 SELECT CASE(this%VDIM)
1569 CASE(1) ! 1D
1570 CALL riemann2prim1d(this%gamma, &
1571 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1572 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1573 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1574 p%density%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1575 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1576 p%pressure%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
1577 CASE(2) ! 2D
1578 CALL riemann2prim2d(this%gamma, &
1579 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1580 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1581 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1582 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
1583 p%density%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1584 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1585 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1586 p%pressure%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
1587 CASE(3) ! 3D
1588 CALL riemann2prim3d(this%gamma, &
1589 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1590 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1591 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1592 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4), &
1593 rinv(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,5), &
1594 p%density%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
1595 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1), &
1596 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2), &
1597 p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,3), &
1598 p%pressure%data3d(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
1599 END SELECT
1600 END SELECT
1601 END SUBROUTINE calculateriemann2primx
1602
1603
1605 !\todo NOT VERIFIED
1606 PURE SUBROUTINE calculateriemann2primy(this,Mesh,j,Rinv,pvar)
1607 IMPLICIT NONE
1608 !------------------------------------------------------------------------!
1609 CLASS(physics_euler), INTENT(IN) :: this
1610 CLASS(mesh_base), INTENT(IN) :: mesh
1611 INTEGER, INTENT(IN) :: j
1612 REAL, INTENT(IN), &
1613 DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%KMIN:Mesh%KMAX,this%VNUM) &
1614 :: rinv
1615 CLASS(marray_compound), INTENT(INOUT) :: pvar
1616 !------------------------------------------------------------------------!
1617 SELECT TYPE(p => pvar)
1618 TYPE IS(statevector_euler)
1619 SELECT CASE(this%VDIM)
1620 CASE(1) ! 1D
1621 CALL riemann2prim1d(this%gamma, &
1622 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1623 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1624 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1625 p%density%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1626 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,1), &
1627 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX))
1628 CASE(2) ! 2D
1629 CALL riemann2prim2d(this%gamma, &
1630 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1631 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1632 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1633 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
1634 p%density%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1635 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,2), &
1636 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,1), &
1637 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX))
1638 CASE(3) ! 3D
1639 CALL riemann2prim3d(this%gamma, &
1640 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,1), &
1641 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,2), &
1642 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,3), &
1643 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,4), &
1644 rinv(mesh%IMIN:mesh%IMAX,mesh%KMIN:mesh%KMAX,5), &
1645 p%density%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX), &
1646 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,2), &
1647 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,3), &
1648 p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,1), &
1649 p%pressure%data3d(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX))
1650 END SELECT
1651 END SELECT
1652 END SUBROUTINE calculateriemann2primy
1653
1654
1656 !\todo NOT VERIFIED
1657 PURE SUBROUTINE calculateriemann2primz(this,Mesh,k,Rinv,pvar)
1658 IMPLICIT NONE
1659 !------------------------------------------------------------------------!
1660 CLASS(physics_euler), INTENT(IN) :: this
1661 CLASS(mesh_base), INTENT(IN) :: mesh
1662 INTEGER, INTENT(IN) :: k
1663 REAL, INTENT(IN), &
1664 DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,this%VNUM) &
1665 :: rinv
1666 CLASS(marray_compound), INTENT(INOUT) :: pvar
1667 !------------------------------------------------------------------------!
1668 SELECT TYPE(p => pvar)
1669 TYPE IS(statevector_euler)
1670 SELECT CASE(this%VDIM)
1671 CASE(1) ! 1D
1672 CALL riemann2prim1d(this%gamma, &
1673 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1674 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1675 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1676 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1677 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,1), &
1678 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k))
1679 CASE(2) ! 2D
1680 CALL riemann2prim2d(this%gamma, &
1681 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1682 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1683 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1684 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1685 p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1686 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,2), &
1687 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,1), &
1688 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k))
1689 CASE(3) ! 3D
1690 CALL riemann2prim3d(this%gamma, &
1691 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1), &
1692 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2), &
1693 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,3), &
1694 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,4), &
1695 rinv(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,5), &
1696 p%density%data3d( mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k), &
1697 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,3), &
1698 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,1), &
1699 p%velocity%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,2), &
1700 p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k))
1701 END SELECT
1702 END SELECT
1703 END SUBROUTINE calculateriemann2primz
1704
1706 PURE SUBROUTINE geometricalsources(this,Mesh,pvar,cvar,sterm)
1707 USE geometry_generic_mod
1708 IMPLICIT NONE
1709 !------------------------------------------------------------------------!
1710 CLASS(physics_euler), INTENT(INOUT) :: this
1711 CLASS(mesh_base), INTENT(IN) :: mesh
1712 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar,sterm
1713 !------------------------------------------------------------------------!
1714 ! compute geometrical source only for non-cartesian mesh
1715 SELECT TYPE(mgeo => mesh%geometry)
1716 TYPE IS(geometry_cartesian)
1717 ! do nothing
1718 CLASS DEFAULT
1719 ! curvilinear geometry
1720 SELECT TYPE(p => pvar)
1721 TYPE IS(statevector_euler)
1722 SELECT TYPE(c => cvar)
1723 TYPE IS(statevector_euler)
1724 SELECT TYPE(s => sterm)
1725 TYPE IS(statevector_euler)
1726 ! no source terms
1727 s%density%data1d(:) = 0.0
1728 s%energy%data1d(:) = 0.0
1729 SELECT CASE(mesh%VECTOR_COMPONENTS)
1730 CASE(vector_x) ! 1D momentum in x-direction
1731 ! vy = vz = my = mz = 0
1732 s%momentum%data2d(:,1) = getgeometricalsourcex( &
1733 mesh%cxyx%data2d(:,2),mesh%cxzx%data2d(:,2), &
1734 mesh%cyxy%data2d(:,2),mesh%czxz%data2d(:,2), &
1735 p%velocity%data2d(:,1),0.0,0.0, &
1736 p%pressure%data1d(:), &
1737 0.0,0.0)
1738 CASE(vector_y) ! 1D momentum in y-direction
1739 ! vx = vz = mx = mz = 0
1740 s%momentum%data2d(:,1) = getgeometricalsourcey( &
1741 mesh%cxyx%data2d(:,2),mesh%cyxy%data2d(:,2), &
1742 mesh%cyzy%data2d(:,2),mesh%czyz%data2d(:,2), &
1743 0.0,p%velocity%data2d(:,1),0.0, &
1744 p%pressure%data1d(:), &
1745 0.0,0.0)
1746 CASE(vector_z) ! 1D momentum in z-direction
1747 ! vx = vy = mx = my = 0
1748 s%momentum%data2d(:,1) = getgeometricalsourcez( &
1749 mesh%cxzx%data2d(:,2),mesh%cyzy%data2d(:,2), &
1750 mesh%czxz%data2d(:,2),mesh%czyz%data2d(:,2), &
1751 0.0,0.0,p%velocity%data2d(:,1), &
1752 p%pressure%data1d(:), &
1753 0.0,0.0)
1754 CASE(ior(vector_x,vector_y)) ! 2D momentum in x-y-plane
1755 ! vz = mz = 0
1756 ! x-momentum
1757 s%momentum%data2d(:,1) = getgeometricalsourcex( &
1758 mesh%cxyx%data2d(:,2),mesh%cxzx%data2d(:,2), &
1759 mesh%cyxy%data2d(:,2),mesh%czxz%data2d(:,2), &
1760 p%velocity%data2d(:,1),p%velocity%data2d(:,2),0.0, &
1761 p%pressure%data1d(:), &
1762 c%momentum%data2d(:,2),0.0)
1763 ! y-momentum
1764 s%momentum%data2d(:,2) = getgeometricalsourcey( &
1765 mesh%cxyx%data2d(:,2),mesh%cyxy%data2d(:,2), &
1766 mesh%cyzy%data2d(:,2),mesh%czyz%data2d(:,2), &
1767 p%velocity%data2d(:,1),p%velocity%data2d(:,2),0.0, &
1768 p%pressure%data1d(:), &
1769 c%momentum%data2d(:,1),0.0)
1770 CASE(ior(vector_x,vector_z)) ! 2D momentum in x-z-plane
1771 ! vy = my = 0
1772 ! x-momentum
1773 s%momentum%data2d(:,1) = getgeometricalsourcex( &
1774 mesh%cxyx%data2d(:,2),mesh%cxzx%data2d(:,2), &
1775 mesh%cyxy%data2d(:,2),mesh%czxz%data2d(:,2), &
1776 p%velocity%data2d(:,1),0.0,p%velocity%data2d(:,2), &
1777 p%pressure%data1d(:), &
1778 0.0,c%momentum%data2d(:,2))
1779 ! z-momentum
1780 s%momentum%data2d(:,2) = getgeometricalsourcez( &
1781 mesh%cxzx%data2d(:,2),mesh%cyzy%data2d(:,2), &
1782 mesh%czxz%data2d(:,2),mesh%czyz%data2d(:,2), &
1783 p%velocity%data2d(:,1),0.0,p%velocity%data2d(:,2), &
1784 p%pressure%data1d(:), &
1785 c%momentum%data2d(:,1),0.0)
1786 CASE(ior(vector_y,vector_z)) ! 2D momentum in y-z-plane
1787 ! vx = mx = 0
1788 ! y-momentum
1789 s%momentum%data2d(:,1) = getgeometricalsourcey( &
1790 mesh%cxyx%data2d(:,2),mesh%cyxy%data2d(:,2), &
1791 mesh%cyzy%data2d(:,2),mesh%czyz%data2d(:,2), &
1792 0.0,p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
1793 p%pressure%data1d(:), &
1794 0.0,c%momentum%data2d(:,2))
1795 ! z-momentum
1796 s%momentum%data2d(:,2) = getgeometricalsourcez( &
1797 mesh%cxzx%data2d(:,2),mesh%cyzy%data2d(:,2), &
1798 mesh%czxz%data2d(:,2),mesh%czyz%data2d(:,2), &
1799 0.0,p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
1800 p%pressure%data1d(:), &
1801 0.0,c%momentum%data2d(:,1))
1802 CASE(ior(ior(vector_x,vector_y),vector_z)) ! 3D momentum
1803 ! x-momentum
1804 s%momentum%data2d(:,1) = getgeometricalsourcex( &
1805 mesh%cxyx%data2d(:,2),mesh%cxzx%data2d(:,2), &
1806 mesh%cyxy%data2d(:,2),mesh%czxz%data2d(:,2), &
1807 p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
1808 p%velocity%data2d(:,3),p%pressure%data1d(:), &
1809 c%momentum%data2d(:,2),c%momentum%data2d(:,3))
1810 ! y-momentum
1811 s%momentum%data2d(:,2) = getgeometricalsourcey( &
1812 mesh%cxyx%data2d(:,2),mesh%cyxy%data2d(:,2), &
1813 mesh%cyzy%data2d(:,2),mesh%czyz%data2d(:,2), &
1814 p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
1815 p%velocity%data2d(:,3),p%pressure%data1d(:), &
1816 c%momentum%data2d(:,1),c%momentum%data2d(:,3))
1817 ! z-momentum
1818 s%momentum%data2d(:,3) = getgeometricalsourcez( &
1819 mesh%cxzx%data2d(:,2),mesh%cyzy%data2d(:,2), &
1820 mesh%czxz%data2d(:,2),mesh%czyz%data2d(:,2), &
1821 p%velocity%data2d(:,1),p%velocity%data2d(:,2), &
1822 p%velocity%data2d(:,3),p%pressure%data1d(:), &
1823 c%momentum%data2d(:,1),c%momentum%data2d(:,2))
1824 CASE DEFAULT
1825 ! return NaN
1826 s%momentum%data1d(:) = 0.0 !NAN_DEFAULT_REAL
1827 END SELECT
1828 END SELECT
1829 END SELECT
1830 END SELECT
1831 ! reset ghost cell data
1832 IF (mesh%INUM.GT.1) THEN
1833 sterm%data4d(mesh%IGMIN:mesh%IMIN+mesh%IM1,:,:,:) = 0.0
1834 sterm%data4d(mesh%IMAX+mesh%IP1:mesh%IGMAX,:,:,:) = 0.0
1835 END IF
1836 IF (mesh%JNUM.GT.1) THEN
1837 sterm%data4d(:,mesh%JGMIN:mesh%JMIN+mesh%JM1,:,:) = 0.0
1838 sterm%data4d(:,mesh%JMAX+mesh%JP1:mesh%JGMAX,:,:) = 0.0
1839 END IF
1840 IF (mesh%KNUM.GT.1) THEN
1841 ! collapse the first 2 dimensions
1842 sterm%data3d(:,mesh%KGMIN:mesh%KMIN+mesh%KM1,:) = 0.0
1843 sterm%data3d(:,mesh%KMAX+mesh%KP1:mesh%KGMAX,:) = 0.0
1844 END IF
1845 END SELECT
1846 END SUBROUTINE geometricalsources
1847
1849 SUBROUTINE externalsources(this,accel,pvar,cvar,sterm)
1850 !------------------------------------------------------------------------!
1851 CLASS(physics_euler), INTENT(IN) :: this
1852 CLASS(marray_base), INTENT(IN) :: accel
1853 CLASS(marray_compound), INTENT(IN) :: pvar,cvar
1854 CLASS(marray_compound), INTENT(INOUT) :: sterm
1855 !------------------------------------------------------------------------!
1856 INTEGER :: m
1857 !------------------------------------------------------------------------!
1858 CALL this%physics_eulerisotherm%ExternalSources(accel,pvar,cvar,sterm)
1859 SELECT TYPE(p => pvar)
1860 TYPE IS(statevector_euler)
1861 SELECT TYPE(c => cvar)
1862 TYPE IS(statevector_euler)
1863 SELECT TYPE(s => sterm)
1864 TYPE IS(statevector_euler)
1865 SELECT CASE(this%VDIM)
1866 CASE(1)
1867 DO m=1,SIZE(s%energy%data1d)
1868 s%energy%data1d(m) = c%momentum%data2d(m,1) * accel%data2d(m,1)
1869 END DO
1870 CASE(2)
1871 DO m=1,SIZE(s%energy%data1d)
1872 s%energy%data1d(m) = c%momentum%data2d(m,1) * accel%data2d(m,1) &
1873 + c%momentum%data2d(m,2) * accel%data2d(m,2)
1874 END DO
1875 CASE(3)
1876 DO m=1,SIZE(s%energy%data1d)
1877 s%energy%data1d(m) = c%momentum%data2d(m,1) * accel%data2d(m,1) &
1878 + c%momentum%data2d(m,2) * accel%data2d(m,2) &
1879 + c%momentum%data2d(m,3) * accel%data2d(m,3)
1880 END DO
1881 END SELECT
1882 END SELECT
1883 END SELECT
1884 END SELECT
1885 END SUBROUTINE externalsources
1886
1887 SUBROUTINE viscositysources(this,Mesh,pvar,btxx,btxy,btxz,btyy,btyz,btzz,sterm)
1888 IMPLICIT NONE
1889 !------------------------------------------------------------------------!
1890 CLASS(physics_euler), INTENT(INOUT) :: this
1891 CLASS(mesh_base), INTENT(IN) :: Mesh
1892 CLASS(marray_compound), INTENT(INOUT) :: pvar,sterm
1893 REAL, INTENT(IN), &
1894 DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
1895 :: btxx,btxy,btxz,btyy,btyz,btzz
1896 !------------------------------------------------------------------------!
1897 CALL this%physics_eulerisotherm%ViscositySources(mesh,pvar,btxx,btxy,btxz,btyy,btyz,btzz,sterm)
1898 SELECT TYPE(p => pvar)
1899 CLASS IS(statevector_euler)
1900 SELECT TYPE(s => sterm)
1901 CLASS IS(statevector_euler)
1902 SELECT CASE(this%VDIM)
1903 ! CASE(1)
1904 ! this%tmp(:,:,:) = pvar(:,:,:,this%XVELOCITY)*btxx(:,:,:)
1905 !
1906 ! CALL Mesh%Divergence(this%tmp(:,:,:),sterm(:,:,:,this%ENERGY))
1907 CASE(2)
1908 !compute scalar product of v and tau (x-component)
1909 this%tmp(:,:,:) = p%velocity%data4d(:,:,:,1)*btxx(:,:,:) &
1910 + p%velocity%data4d(:,:,:,2)*btxy(:,:,:)
1911
1912 !compute scalar product of v and tau (y-component)
1913 this%tmp1(:,:,:) = p%velocity%data4d(:,:,:,1)*btxy(:,:,:) &
1914 + p%velocity%data4d(:,:,:,2)*btyy(:,:,:)
1915
1916 ! compute vector divergence of scalar product v and tau
1917 CALL mesh%Divergence(this%tmp(:,:,:),this%tmp1(:,:,:), &
1918 s%energy%data3d(:,:,:))
1919 CASE(3)
1920 !compute scalar product of v and tau (x-component)
1921 this%tmp(:,:,:) = p%velocity%data4d(:,:,:,1)*btxx(:,:,:) &
1922 + p%velocity%data4d(:,:,:,2)*btxy(:,:,:) &
1923 + p%velocity%data4d(:,:,:,3)*btxz(:,:,:)
1924
1925 !compute scalar product of v and tau (y-component)
1926 this%tmp1(:,:,:) = p%velocity%data4d(:,:,:,1)*btxy(:,:,:) &
1927 + p%velocity%data4d(:,:,:,2)*btyy(:,:,:) &
1928 + p%velocity%data4d(:,:,:,3)*btyz(:,:,:)
1929
1930 !compute scalar product of v and tau (z-component)
1931 this%tmp2(:,:,:) = p%velocity%data4d(:,:,:,1)*btxz(:,:,:) &
1932 + p%velocity%data4d(:,:,:,2)*btyz(:,:,:) &
1933 + p%velocity%data4d(:,:,:,3)*btzz(:,:,:)
1934 ! compute vector divergence of scalar product v and tau
1935 CALL mesh%Divergence(this%tmp(:,:,:),this%tmp1(:,:,:),this%tmp2(:,:,:), &
1936 s%energy%data3d(:,:,:))
1937 CASE DEFAULT
1938 ! return NaN
1939 s%data1d(:) = nan_default_real
1940 END SELECT
1941 END SELECT
1942 END SELECT
1943 END SUBROUTINE viscositysources
1944
1955 PURE SUBROUTINE addbackgroundvelocityx(this,Mesh,w,pvar,cvar)
1956 IMPLICIT NONE
1957 !------------------------------------------------------------------------!
1958 CLASS(physics_euler), INTENT(INOUT) :: this
1959 CLASS(mesh_base), INTENT(IN) :: mesh
1960 REAL, DIMENSION(Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
1961 INTENT(IN) :: w
1962 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
1963 !------------------------------------------------------------------------!
1964 INTEGER :: i,j,k
1965 !------------------------------------------------------------------------!
1966 SELECT TYPE(p => pvar)
1967 TYPE IS(statevector_euler)
1968 SELECT TYPE(c => cvar)
1969 TYPE IS(statevector_euler)
1970 IF (c%fargo_transformation_applied) THEN
1971 IF (.NOT.p%fargo_transformation_applied) RETURN ! should not happen
1972 DO k=mesh%KGMIN,mesh%KGMAX
1973 DO j=mesh%JGMIN,mesh%JGMAX
1974 DO i=mesh%IGMIN,mesh%IGMAX
1975 ! ATTENTION: don't change the order; on the RHS of the first
1976 ! assignment there must be the old momentum
1977 c%energy%data3d(i,j,k) = c%energy%data3d(i,j,k) &
1978 + w(j,k)*(c%momentum%data4d(i,j,k,1) &
1979 + 0.5*c%density%data3d(i,j,k)*w(j,k))
1980 p%velocity%data4d(i,j,k,1) = p%velocity%data4d(i,j,k,1) + w(j,k)
1981 c%momentum%data4d(i,j,k,1) = c%momentum%data4d(i,j,k,1) &
1982 + c%density%data3d(i,j,k)*w(j,k)
1983 END DO
1984 END DO
1985 END DO
1986 c%fargo_transformation_applied = .false.
1987 p%fargo_transformation_applied = .false.
1988 END IF
1989 END SELECT
1990 END SELECT
1991 END SUBROUTINE addbackgroundvelocityx
1992
2003 PURE SUBROUTINE addbackgroundvelocityy(this,Mesh,w,pvar,cvar)
2004 IMPLICIT NONE
2005 !------------------------------------------------------------------------!
2006 CLASS(physics_euler), INTENT(INOUT) :: this
2007 CLASS(mesh_base), INTENT(IN) :: mesh
2008 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%KGMIN:Mesh%KGMAX), &
2009 INTENT(IN) :: w
2010 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
2011 !------------------------------------------------------------------------!
2012 INTEGER :: i,j,k,v_idx
2013 !------------------------------------------------------------------------!
2014 SELECT TYPE(p => pvar)
2015 TYPE IS(statevector_euler)
2016 SELECT TYPE(c => cvar)
2017 TYPE IS(statevector_euler)
2018 IF (c%fargo_transformation_applied) THEN
2019 IF (.NOT.p%fargo_transformation_applied) RETURN ! should not happen
2020 ! check if x-component of velocity vector is available
2021 IF (btest(mesh%VECTOR_COMPONENTS,0)) THEN
2022 ! y-component is the 2nd component
2023 v_idx = 2
2024 ELSE
2025 ! no x-component -> y-component is the first component
2026 v_idx = 1
2027 END IF
2028 DO k=mesh%KGMIN,mesh%KGMAX
2029 DO j=mesh%JGMIN,mesh%JGMAX
2030 DO i=mesh%IGMIN,mesh%IGMAX
2031 ! ATTENTION: don't change the order; on the RHS of the first
2032 ! assignment there must be the old momentum
2033 c%energy%data3d(i,j,k) = c%energy%data3d(i,j,k) &
2034 + w(i,k)*(c%momentum%data4d(i,j,k,v_idx) &
2035 + 0.5*c%density%data3d(i,j,k)*w(i,k))
2036 p%velocity%data4d(i,j,k,v_idx) = p%velocity%data4d(i,j,k,v_idx) + w(i,k)
2037 c%momentum%data4d(i,j,k,v_idx) = c%momentum%data4d(i,j,k,v_idx) &
2038 + c%density%data3d(i,j,k)*w(i,k)
2039 END DO
2040 END DO
2041 END DO
2042 c%fargo_transformation_applied = .false.
2043 p%fargo_transformation_applied = .false.
2044 END IF
2045 END SELECT
2046 END SELECT
2047 END SUBROUTINE addbackgroundvelocityy
2048
2059 PURE SUBROUTINE addbackgroundvelocityz(this,Mesh,w,pvar,cvar)
2060 IMPLICIT NONE
2061 !------------------------------------------------------------------------!
2062 CLASS(physics_euler), INTENT(INOUT) :: this
2063 CLASS(mesh_base), INTENT(IN) :: mesh
2064 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX), &
2065 INTENT(IN) :: w
2066 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
2067 !------------------------------------------------------------------------!
2068 INTEGER :: i,j,k
2069 !------------------------------------------------------------------------!
2070 SELECT TYPE(p => pvar)
2071 TYPE IS(statevector_euler)
2072 SELECT TYPE(c => cvar)
2073 TYPE IS(statevector_euler)
2074 IF (c%fargo_transformation_applied) THEN
2075 IF (.NOT.p%fargo_transformation_applied) RETURN ! should not happen
2076 DO k=mesh%KGMIN,mesh%KGMAX
2077 DO j=mesh%JGMIN,mesh%JGMAX
2078 DO i=mesh%IGMIN,mesh%IGMAX
2079 ! ATTENTION: don't change the order; on the RHS of the first
2080 ! assignment there must be the old momentum
2081 c%energy%data3d(i,j,k) = c%energy%data3d(i,j,k) &
2082 + w(i,j)*(c%momentum%data4d(i,j,k,this%VDIM) &
2083 + 0.5*c%density%data3d(i,j,k)*w(i,j))
2084 p%velocity%data4d(i,j,k,this%VDIM) = p%velocity%data4d(i,j,k,this%VDIM) + w(i,j)
2085 c%momentum%data4d(i,j,k,this%VDIM) = c%momentum%data4d(i,j,k,this%VDIM) &
2086 + c%density%data3d(i,j,k)*w(i,j)
2087 END DO
2088 END DO
2089 END DO
2090 c%fargo_transformation_applied = .false.
2091 p%fargo_transformation_applied = .false.
2092 END IF
2093 END SELECT
2094 END SELECT
2095 END SUBROUTINE addbackgroundvelocityz
2096
2109 PURE SUBROUTINE subtractbackgroundvelocityx(this,Mesh,w,pvar,cvar)
2110 IMPLICIT NONE
2111 !------------------------------------------------------------------------!
2112 CLASS(physics_euler), INTENT(INOUT) :: this
2113 CLASS(mesh_base), INTENT(IN) :: mesh
2114 REAL,DIMENSION(Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
2115 INTENT(IN) :: w
2116 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
2117 !------------------------------------------------------------------------!
2118 INTEGER :: i,j,k
2119 !------------------------------------------------------------------------!
2120 SELECT TYPE(p => pvar)
2121 TYPE IS(statevector_euler)
2122 SELECT TYPE(c => cvar)
2123 TYPE IS(statevector_euler)
2124 IF (.NOT.c%fargo_transformation_applied) THEN
2125 IF (p%fargo_transformation_applied) RETURN ! should not happen
2126 DO k=mesh%KGMIN,mesh%KGMAX
2127 DO j=mesh%JGMIN,mesh%JGMAX
2128 DO i=mesh%IGMIN,mesh%IGMAX
2129 ! ATTENTION: don't change the order; on the RHS of the first
2130 ! assignment there must be the old momentum
2131 c%energy%data3d(i,j,k) = c%energy%data3d(i,j,k) &
2132 - w(j,k)*(c%momentum%data4d(i,j,k,1) &
2133 - 0.5*c%density%data3d(i,j,k)*w(j,k))
2134 p%velocity%data4d(i,j,k,1) = p%velocity%data4d(i,j,k,1) - w(j,k)
2135 c%momentum%data4d(i,j,k,1) = c%momentum%data4d(i,j,k,1) &
2136 - c%density%data3d(i,j,k)*w(j,k)
2137 END DO
2138 END DO
2139 END DO
2140 c%fargo_transformation_applied = .true.
2141 p%fargo_transformation_applied = .true.
2142 END IF
2143 END SELECT
2144 END SELECT
2145 END SUBROUTINE subtractbackgroundvelocityx
2146
2159 PURE SUBROUTINE subtractbackgroundvelocityy(this,Mesh,w,pvar,cvar)
2160 IMPLICIT NONE
2161 !------------------------------------------------------------------------!
2162 CLASS(physics_euler), INTENT(INOUT) :: this
2163 CLASS(mesh_base), INTENT(IN) :: mesh
2164 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%KGMIN:Mesh%KGMAX), &
2165 INTENT(IN) :: w
2166 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
2167 !------------------------------------------------------------------------!
2168 INTEGER :: i,j,k,v_idx
2169 !------------------------------------------------------------------------!
2170 SELECT TYPE(p => pvar)
2171 TYPE IS(statevector_euler)
2172 SELECT TYPE(c => cvar)
2173 TYPE IS(statevector_euler)
2174 IF (.NOT.c%fargo_transformation_applied) THEN
2175 IF (p%fargo_transformation_applied) RETURN ! should not happen
2176 ! check if x-component of velocity vector is available
2177 IF (btest(mesh%VECTOR_COMPONENTS,0)) THEN
2178 ! y-component is the 2nd component
2179 v_idx = 2
2180 ELSE
2181 ! no x-component -> y-component is the first component
2182 v_idx = 1
2183 END IF
2184 DO k=mesh%KGMIN,mesh%KGMAX
2185 DO j=mesh%JGMIN,mesh%JGMAX
2186 DO i=mesh%IGMIN,mesh%IGMAX
2187 ! ATTENTION: don't change the order; on the RHS of the first
2188 ! assignment there must be the old momentum
2189 c%energy%data3d(i,j,k) = c%energy%data3d(i,j,k) &
2190 - w(i,k)*(c%momentum%data4d(i,j,k,v_idx) &
2191 - 0.5*c%density%data3d(i,j,k)*w(i,k))
2192 p%velocity%data4d(i,j,k,v_idx) = p%velocity%data4d(i,j,k,v_idx) - w(i,k)
2193 c%momentum%data4d(i,j,k,v_idx) = c%momentum%data4d(i,j,k,v_idx) &
2194 - c%density%data3d(i,j,k)*w(i,k)
2195 END DO
2196 END DO
2197 END DO
2198 c%fargo_transformation_applied = .true.
2199 p%fargo_transformation_applied = .true.
2200 END IF
2201 END SELECT
2202 END SELECT
2203 END SUBROUTINE subtractbackgroundvelocityy
2204
2217 PURE SUBROUTINE subtractbackgroundvelocityz(this,Mesh,w,pvar,cvar)
2218 IMPLICIT NONE
2219 !------------------------------------------------------------------------!
2220 CLASS(physics_euler), INTENT(INOUT) :: this
2221 CLASS(mesh_base), INTENT(IN) :: mesh
2222 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX), &
2223 INTENT(IN) :: w
2224 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
2225 !------------------------------------------------------------------------!
2226 INTEGER :: i,j,k
2227 !------------------------------------------------------------------------!
2228 SELECT TYPE(p => pvar)
2229 TYPE IS(statevector_euler)
2230 SELECT TYPE(c => cvar)
2231 TYPE IS(statevector_euler)
2232 IF (.NOT.c%fargo_transformation_applied) THEN
2233 IF (p%fargo_transformation_applied) RETURN ! should not happen
2234 DO k=mesh%KGMIN,mesh%KGMAX
2235 DO j=mesh%JGMIN,mesh%JGMAX
2236 DO i=mesh%IGMIN,mesh%IGMAX
2237 ! ATTENTION: don't change the order; on the RHS of the first
2238 ! assignment there must be the old momentum
2239 c%energy%data3d(i,j,k) = c%energy%data3d(i,j,k) &
2240 - w(i,j)*(c%momentum%data4d(i,j,k,this%VDIM) &
2241 - 0.5*c%density%data3d(i,j,k)*w(i,j))
2242 p%velocity%data4d(i,j,k,this%VDIM) = p%velocity%data4d(i,j,k,this%VDIM) - w(i,j)
2243 c%momentum%data4d(i,j,k,this%VDIM) = c%momentum%data4d(i,j,k,this%VDIM) &
2244 - c%density%data3d(i,j,k)*w(i,j)
2245 END DO
2246 END DO
2247 END DO
2248 c%fargo_transformation_applied = .true.
2249 p%fargo_transformation_applied = .true.
2250 END IF
2251 END SELECT
2252 END SELECT
2253 END SUBROUTINE subtractbackgroundvelocityz
2254
2255
2267 PURE SUBROUTINE addfargosourcesx(this,Mesh,w,pvar,cvar,sterm)
2268 IMPLICIT NONE
2269 !------------------------------------------------------------------------!
2270 CLASS(physics_euler), INTENT(INOUT) :: this
2271 CLASS(mesh_base), INTENT(IN) :: mesh
2272 REAL, DIMENSION(Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), INTENT(IN) :: w
2273 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar,sterm
2274 !------------------------------------------------------------------------!
2275 INTEGER :: i,j,k
2276 !------------------------------------------------------------------------!
2277 SELECT TYPE(c => cvar)
2278 CLASS IS(statevector_euler)
2279 SELECT TYPE(p => pvar)
2280 TYPE IS(statevector_euler)
2281 SELECT TYPE(s => sterm)
2282 CLASS IS(statevector_euler)
2283 ! ATTENTION: fargo sources are added to the given sterm
2284 SELECT CASE(mesh%VECTOR_COMPONENTS)
2285 CASE(ior(vector_x,vector_y))
2286 ! 2D (x,y) momentum vector
2287 DO k=mesh%KMIN,mesh%KMAX
2288 DO j=mesh%JMIN,mesh%JMAX
2289 DO i=mesh%IMIN,mesh%IMAX
2290 this%tmp(i,j,k) = c%momentum%data4d(i,j,k,2) * 0.5 * (w(j+1,k)-w(j-1,k)) &
2291 / mesh%dly%data3d(i,j,k)
2292 ! x-momentum source
2293 s%momentum%data4d(i,j,k,1) = s%momentum%data4d(i,j,k,1) - this%tmp(i,j,k)
2294 ! add geometrical terms
2295 this%tmp(i,j,k) = this%tmp(i,j,k) + getgeometricalsourcex( &
2296 mesh%cxyx%data4d(i,j,k,2),mesh%cxzx%data4d(i,j,k,2), &
2297 mesh%cyxy%data4d(i,j,k,2),mesh%czxz%data4d(i,j,k,2), &
2298 p%velocity%data4d(i,j,k,1),p%velocity%data4d(i,j,k,2),0.0, &
2299 0.0,c%momentum%data4d(i,j,k,2),0.0)
2300 ! energy source
2301 s%energy%data3d(i,j,k) = s%energy%data3d(i,j,k) &
2302 - p%velocity%data4d(i,j,k,1)*this%tmp(i,j,k)
2303 END DO
2304 END DO
2305 END DO
2306 CASE(ior(vector_x,vector_z))
2307 ! 2D (x,z) momentum vector
2308 DO k=mesh%KMIN,mesh%KMAX
2309 DO j=mesh%JMIN,mesh%JMAX
2310 DO i=mesh%IMIN,mesh%IMAX
2311 this%tmp(i,j,k) = c%momentum%data4d(i,j,k,2) * 0.5 * (w(j,k+1)-w(j,k-1)) &
2312 / mesh%dlz%data3d(i,j,k)
2313 ! x-momentum source
2314 s%momentum%data4d(i,j,k,1) = s%momentum%data4d(i,j,k,1) - this%tmp(i,j,k)
2315 ! add geometrical terms
2316 this%tmp(i,j,k) = this%tmp(i,j,k) + getgeometricalsourcex( &
2317 mesh%cxyx%data4d(i,j,k,2),mesh%cxzx%data4d(i,j,k,2), &
2318 mesh%cyxy%data4d(i,j,k,2),mesh%czxz%data4d(i,j,k,2), &
2319 p%velocity%data4d(i,j,k,1),0.0,p%velocity%data4d(i,j,k,2), &
2320 0.0,0.0,c%momentum%data4d(i,j,k,2))
2321 ! energy source
2322 s%energy%data3d(i,j,k) = s%energy%data3d(i,j,k) &
2323 - p%velocity%data4d(i,j,k,1)*this%tmp(i,j,k)
2324 END DO
2325 END DO
2326 END DO
2327 CASE(ior(ior(vector_x,vector_y),vector_z))
2328 ! 3D momentum vector
2329 DO k=mesh%KMIN,mesh%KMAX
2330 DO j=mesh%JMIN,mesh%JMAX
2331 DO i=mesh%IMIN,mesh%IMAX
2332 this%tmp(i,j,k) = c%momentum%data4d(i,j,k,2) * 0.5 * (w(j+1,k)-w(j-1,k)) / mesh%dly%data3d(i,j,k) &
2333 + c%momentum%data4d(i,j,k,3) * 0.5 * (w(j,k+1)-w(j,k-1)) / mesh%dlz%data3d(i,j,k)
2334 ! x-momentum source
2335 s%momentum%data4d(i,j,k,1) = s%momentum%data4d(i,j,k,1) - this%tmp(i,j,k)
2336 ! add geometrical terms
2337 this%tmp(i,j,k) = this%tmp(i,j,k) + getgeometricalsourcex( &
2338 mesh%cxyx%data4d(i,j,k,2),mesh%cxzx%data4d(i,j,k,2), &
2339 mesh%cyxy%data4d(i,j,k,2),mesh%czxz%data4d(i,j,k,2), &
2340 p%velocity%data4d(i,j,k,1),p%velocity%data4d(i,j,k,2), &
2341 p%velocity%data4d(i,j,k,3), &
2342 0.0,c%momentum%data4d(i,j,k,2),c%momentum%data4d(i,j,k,3))
2343 ! energy source
2344 s%energy%data3d(i,j,k) = s%energy%data3d(i,j,k) &
2345 - p%velocity%data4d(i,j,k,1)*this%tmp(i,j,k)
2346 END DO
2347 END DO
2348 END DO
2349 CASE DEFAULT
2350 ! do nothing in any other case
2351 END SELECT
2352 END SELECT
2353 END SELECT
2354 END SELECT
2355 END SUBROUTINE addfargosourcesx
2356
2368 PURE SUBROUTINE addfargosourcesy(this,Mesh,w,pvar,cvar,sterm)
2369 IMPLICIT NONE
2370 !------------------------------------------------------------------------!
2371 CLASS(physics_euler), INTENT(INOUT) :: this
2372 CLASS(mesh_base), INTENT(IN) :: mesh
2373 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%KGMIN:Mesh%KGMAX), INTENT(IN) :: w
2374 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar,sterm
2375 !------------------------------------------------------------------------!
2376 INTEGER :: i,j,k
2377 !------------------------------------------------------------------------!
2378 SELECT TYPE(c => cvar)
2379 CLASS IS(statevector_euler)
2380 SELECT TYPE(p => pvar)
2381 TYPE IS(statevector_euler)
2382 SELECT TYPE(s => sterm)
2383 CLASS IS(statevector_euler)
2384 SELECT CASE(mesh%VECTOR_COMPONENTS)
2385 CASE(ior(vector_x,vector_y))
2386 ! 2D (x,y) momentum vector
2387 DO k=mesh%KMIN,mesh%KMAX
2388 DO j=mesh%JMIN,mesh%JMAX
2389 DO i=mesh%IMIN,mesh%IMAX
2390 this%tmp(i,j,k) = c%momentum%data4d(i,j,k,1) * 0.5 * (w(i+1,k)-w(i-1,k)) &
2391 / mesh%dlx%data3d(i,j,k)
2392 ! y-momentum source
2393 s%momentum%data4d(i,j,k,2) = s%momentum%data4d(i,j,k,2) - this%tmp(i,j,k)
2394 ! add geometrical terms
2395 this%tmp(i,j,k) = this%tmp(i,j,k) + getgeometricalsourcey( &
2396 mesh%cxyx%data4d(i,j,k,2),mesh%cyxy%data4d(i,j,k,2), &
2397 mesh%cyzy%data4d(i,j,k,2),mesh%czyz%data4d(i,j,k,2), &
2398 p%velocity%data4d(i,j,k,1),p%velocity%data4d(i,j,k,2),0.0, &
2399 0.0,c%momentum%data4d(i,j,k,1),0.0)
2400 ! energy source
2401 s%energy%data3d(i,j,k) = s%energy%data3d(i,j,k) &
2402 - p%velocity%data4d(i,j,k,2)*this%tmp(i,j,k)
2403 END DO
2404 END DO
2405 END DO
2406 CASE(ior(vector_y,vector_z))
2407 ! 2D (y,z) momentum vector
2408 DO k=mesh%KMIN,mesh%KMAX
2409 DO j=mesh%JMIN,mesh%JMAX
2410 DO i=mesh%IMIN,mesh%IMAX
2411 this%tmp(i,j,k) = c%momentum%data4d(i,j,k,2) * 0.5 * (w(i,k+1)-w(i,k-1)) &
2412 / mesh%dlz%data3d(i,j,k)
2413 ! y-momentum source
2414 s%momentum%data4d(i,j,k,1) = s%momentum%data4d(i,j,k,1) - this%tmp(i,j,k)
2415 ! add geometrical terms
2416 this%tmp(i,j,k) = this%tmp(i,j,k) + getgeometricalsourcey( &
2417 mesh%cxyx%data4d(i,j,k,2),mesh%cyxy%data4d(i,j,k,2), &
2418 mesh%cyzy%data4d(i,j,k,2),mesh%czyz%data4d(i,j,k,2), &
2419 0.0,p%velocity%data4d(i,j,k,1),p%velocity%data4d(i,j,k,2), &
2420 0.0,0.0,c%momentum%data4d(i,j,k,2))
2421 ! energy source
2422 s%energy%data3d(i,j,k) = s%energy%data3d(i,j,k) &
2423 - p%velocity%data4d(i,j,k,1)*this%tmp(i,j,k)
2424 END DO
2425 END DO
2426 END DO
2427 CASE(ior(ior(vector_x,vector_y),vector_z))
2428 ! 3D momentum vector
2429 DO k=mesh%KMIN,mesh%KMAX
2430 DO j=mesh%JMIN,mesh%JMAX
2431 DO i=mesh%IMIN,mesh%IMAX
2432 this%tmp(i,j,k) = c%momentum%data4d(i,j,k,1) * 0.5 * (w(i+1,k)-w(i-1,k)) / mesh%dlx%data3d(i,j,k) &
2433 + c%momentum%data4d(i,j,k,3) * 0.5 * (w(j,k+1)-w(j,k-1)) / mesh%dlz%data3d(i,j,k)
2434 ! y-momentum source
2435 s%momentum%data4d(i,j,k,2) = s%momentum%data4d(i,j,k,2) - this%tmp(i,j,k)
2436 ! add geometrical terms
2437 this%tmp(i,j,k) = this%tmp(i,j,k) + getgeometricalsourcey( &
2438 mesh%cxyx%data4d(i,j,k,2),mesh%cyxy%data4d(i,j,k,2), &
2439 mesh%cyzy%data4d(i,j,k,2),mesh%czyz%data4d(i,j,k,2), &
2440 p%velocity%data4d(i,j,k,1),p%velocity%data4d(i,j,k,2), &
2441 p%velocity%data4d(i,j,k,3),0.0, &
2442 c%momentum%data4d(i,j,k,1),c%momentum%data4d(i,j,k,3))
2443 ! energy source
2444 s%energy%data3d(i,j,k) = s%energy%data3d(i,j,k) &
2445 - p%velocity%data4d(i,j,k,2)*this%tmp(i,j,k)
2446 END DO
2447 END DO
2448 END DO
2449 CASE DEFAULT
2450 ! do nothing in any other case
2451 END SELECT
2452 END SELECT
2453 END SELECT
2454 END SELECT
2455 END SUBROUTINE addfargosourcesy
2456
2468 PURE SUBROUTINE addfargosourcesz(this,Mesh,w,pvar,cvar,sterm)
2469 IMPLICIT NONE
2470 !------------------------------------------------------------------------!
2471 CLASS(physics_euler), INTENT(INOUT) :: this
2472 CLASS(mesh_base), INTENT(IN) :: mesh
2473 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX), INTENT(IN) :: w
2474 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar,sterm
2475 !------------------------------------------------------------------------!
2476 INTEGER :: i,j,k
2477 !------------------------------------------------------------------------!
2478 SELECT TYPE(c => cvar)
2479 CLASS IS(statevector_euler)
2480 SELECT TYPE(p => pvar)
2481 TYPE IS(statevector_euler)
2482 SELECT TYPE(s => sterm)
2483 CLASS IS(statevector_euler)
2484 SELECT CASE(mesh%VECTOR_COMPONENTS)
2485 CASE(ior(vector_x,vector_z))
2486 ! 2D (x,z) momentum vector
2487 DO k=mesh%KMIN,mesh%KMAX
2488 DO j=mesh%JMIN,mesh%JMAX
2489 DO i=mesh%IMIN,mesh%IMAX
2490 this%tmp(i,j,k) = c%momentum%data4d(i,j,k,1) * 0.5 * (w(i+1,j)-w(i-1,j)) / mesh%dlx%data3d(i,j,k)
2491 ! z-momentum source
2492 s%momentum%data4d(i,j,k,2) = s%momentum%data4d(i,j,k,2) - this%tmp(i,j,k)
2493 ! add geometrical terms
2494 this%tmp(i,j,k) = this%tmp(i,j,k) + getgeometricalsourcez( &
2495 mesh%cxzx%data4d(i,j,k,2),mesh%cyzy%data4d(i,j,k,2), &
2496 mesh%czxz%data4d(i,j,k,2),mesh%czyz%data4d(i,j,k,2), &
2497 p%velocity%data4d(i,j,k,1),0.0,p%velocity%data4d(i,j,k,2), &
2498 0.0,c%momentum%data4d(i,j,k,1),0.0)
2499 ! energy source
2500 s%energy%data3d(i,j,k) = s%energy%data3d(i,j,k) &
2501 - p%velocity%data4d(i,j,k,2)*this%tmp(i,j,k)
2502 END DO
2503 END DO
2504 END DO
2505 CASE(ior(vector_y,vector_z))
2506 ! 2D (y,z) momentum vector
2507 DO k=mesh%KMIN,mesh%KMAX
2508 DO j=mesh%JMIN,mesh%JMAX
2509 DO i=mesh%IMIN,mesh%IMAX
2510 this%tmp(i,j,k) = c%momentum%data4d(i,j,k,1) * 0.5 * (w(i,j+1)-w(i,j-1)) / mesh%dly%data3d(i,j,k)
2511 ! z-momentum source
2512 s%momentum%data4d(i,j,k,2) = s%momentum%data4d(i,j,k,2) - this%tmp(i,j,k)
2513 ! add geometrical terms
2514 this%tmp(i,j,k) = this%tmp(i,j,k) + getgeometricalsourcez( &
2515 mesh%cxzx%data4d(i,j,k,2),mesh%cyzy%data4d(i,j,k,2), &
2516 mesh%czxz%data4d(i,j,k,2),mesh%czyz%data4d(i,j,k,2), &
2517 0.0,p%velocity%data4d(i,j,k,1),p%velocity%data4d(i,j,k,2), &
2518 0.0,0.0,c%momentum%data4d(i,j,k,1))
2519 ! energy source
2520 s%energy%data3d(i,j,k) = s%energy%data3d(i,j,k) &
2521 - p%velocity%data4d(i,j,k,2)*this%tmp(i,j,k)
2522 END DO
2523 END DO
2524 END DO
2525 CASE(ior(ior(vector_x,vector_y),vector_z))
2526 ! 3D momentum vector
2527 DO k=mesh%KMIN,mesh%KMAX
2528 DO j=mesh%JMIN,mesh%JMAX
2529 DO i=mesh%IMIN,mesh%IMAX
2530 this%tmp(i,j,k) = c%momentum%data4d(i,j,k,1) * 0.5 * (w(i+1,j)-w(i-1,j)) / mesh%dlx%data3d(i,j,k) &
2531 + c%momentum%data4d(i,j,k,2) * 0.5 * (w(j,j+1)-w(j,j-1)) / mesh%dly%data3d(i,j,k)
2532 ! z-momentum source
2533 s%momentum%data4d(i,j,k,3) = s%momentum%data4d(i,j,k,3) - this%tmp(i,j,k)
2534 ! add geometrical terms
2535 this%tmp(i,j,k) = this%tmp(i,j,k) + getgeometricalsourcez( &
2536 mesh%cxzx%data4d(i,j,k,2),mesh%cyzy%data4d(i,j,k,2), &
2537 mesh%czxz%data4d(i,j,k,2),mesh%czyz%data4d(i,j,k,2), &
2538 p%velocity%data4d(i,j,k,1),p%velocity%data4d(i,j,k,2), &
2539 p%velocity%data4d(i,j,k,3), &
2540 0.0,c%momentum%data4d(i,j,k,1),c%momentum%data4d(i,j,k,2))
2541 ! energy source
2542 s%energy%data3d(i,j,k) = s%energy%data3d(i,j,k) &
2543 - p%velocity%data4d(i,j,k,3)*this%tmp(i,j,k)
2544 END DO
2545 END DO
2546 END DO
2547 CASE DEFAULT
2548 ! do nothing in any other case
2549 END SELECT
2550 END SELECT
2551 END SELECT
2552 END SELECT
2553 END SUBROUTINE addfargosourcesz
2554
2555 PURE SUBROUTINE updatesoundspeed(this,pvar)
2556 IMPLICIT NONE
2557 !------------------------------------------------------------------------!
2558 CLASS(physics_euler), INTENT(INOUT) :: this
2559 CLASS(statevector_euler),INTENT(INOUT) :: pvar
2560 !------------------------------------------------------------------------!
2561 INTEGER :: i
2562 !------------------------------------------------------------------------!
2563 IF (pvar%density%RANK.EQ.0) THEN
2564 DO concurrent(i=1:SIZE(this%bccsound%data1d))
2565 this%bccsound%data1d(i) = getsoundspeed(this%gamma, &
2566 pvar%density%data1d(i),pvar%pressure%data1d(i))
2567 END DO
2568 ELSE
2569 DO concurrent(i=1:SIZE(this%fcsound%data1d))
2570 this%fcsound%data1d(i) = getsoundspeed(this%gamma, &
2571 pvar%density%data1d(i),pvar%pressure%data1d(i))
2572 END DO
2573 END IF
2574 END SUBROUTINE updatesoundspeed
2575
2577 SUBROUTINE finalize(this)
2578 IMPLICIT NONE
2579 !------------------------------------------------------------------------!
2580 CLASS(physics_euler), INTENT(INOUT) :: this
2581 !------------------------------------------------------------------------!
2582 CALL this%physics_eulerisotherm%Finalize()
2583 END SUBROUTINE finalize
2584
2585!----------------------------------------------------------------------------!
2587
2593 FUNCTION createstatevector(Physics,flavour,num) RESULT(new_sv)
2594 IMPLICIT NONE
2595 !-------------------------------------------------------------------!
2596 CLASS(physics_euler), INTENT(IN) :: physics
2597 INTEGER, OPTIONAL, INTENT(IN) :: flavour,num
2598 TYPE(statevector_euler) :: new_sv
2599 !-------------------------------------------------------------------!
2600 LOGICAL :: success = .false.
2601 INTEGER :: err
2602 !-------------------------------------------------------------------!
2603 ! call inherited function
2604 new_sv = statevector_eulerisotherm(physics,flavour,num)
2605 ! add entries specific for euler physics
2606 ALLOCATE(new_sv%pressure,stat=err)
2607 IF (err.EQ.0) THEN
2608 IF (PRESENT(num)) THEN
2609 new_sv%pressure = marray_base(num) ! num scalars
2610 ELSE
2611 new_sv%pressure = marray_base() ! one scalar
2612 END IF
2613 success = new_sv%AppendMArray(new_sv%pressure)
2614 SELECT CASE(new_sv%flavour)
2615 CASE(primitive)
2616 new_sv%energy => null()
2617 CASE(conservative)
2618 new_sv%energy => new_sv%pressure
2619 new_sv%pressure => null()
2620 CASE DEFAULT
2621#ifdef DEBUG
2622 print *,"ERROR in physics_euler::CreateStateVector: unknown flavour"
2623#endif
2624 END SELECT
2625 END IF
2626 IF (.NOT.success) &
2627 CALL physics%Error("physics_euler::CreateStateVector", &
2628 "state vector initialization failed")
2629 END FUNCTION createstatevector
2630
2632 SUBROUTINE assignmarray_0(this,ma)
2633 IMPLICIT NONE
2634 !------------------------------------------------------------------------!
2635 CLASS(statevector_euler),INTENT(INOUT) :: this
2636 CLASS(marray_base),INTENT(IN) :: ma
2637 !------------------------------------------------------------------------!
2638#if DEBUG > 2
2639 print *,"DEBUG INFO in physics_euler::AssignMArray_0: assigning 2 state vectors"
2640#endif
2641 CALL this%statevector_eulerisotherm%AssignMArray_0(ma)
2642 IF (SIZE(this%data1d).LE.0) RETURN ! empty compound
2643 SELECT TYPE(src => ma)
2644 CLASS IS(statevector_euler)
2645#if DEBUG > 2
2646 print *,"DEBUG INFO in physics_euler::AssignMArray_0: restoring pressure/energy pointers"
2647#endif
2648 SELECT CASE(this%flavour)
2649 CASE(primitive)
2650 ! pressure is the third item
2651 this%pressure => this%GetItem(this%NextItem(this%NextItem(this%FirstItem())))
2652 this%energy => null()
2653 CASE(conservative)
2654 ! energy is the third item
2655 this%energy => this%GetItem(this%NextItem(this%NextItem(this%FirstItem())))
2656 this%pressure => null()
2657 CASE DEFAULT
2658 ! error, this should not happen
2659 END SELECT
2660 CLASS DEFAULT
2661 ! do nothing: ma may not be of type euler, i.e. during initialization (see CreateStatevector)
2662 END SELECT
2663 END SUBROUTINE assignmarray_0
2664
2666 SUBROUTINE finalize_statevector(this)
2667 IMPLICIT NONE
2668 !------------------------------------------------------------------------!
2669 TYPE(statevector_euler),INTENT(INOUT) :: this
2670 !------------------------------------------------------------------------!
2671#if DEBUG > 2
2672 print *,"DEBUG INFO in physics_euler::Finalize: cleanup statevector"
2673#endif
2674 NULLIFY(this%pressure,this%energy)
2675 END SUBROUTINE finalize_statevector
2676
2677
2678!----------------------------------------------------------------------------!
2680
2681
2682 ELEMENTAL FUNCTION getsoundspeed(gamma,density,pressure) RESULT(cs)
2683 IMPLICIT NONE
2684 !------------------------------------------------------------------------!
2685 REAL, INTENT(IN) :: gamma,density,pressure
2686 REAL :: cs
2687 !------------------------------------------------------------------------!
2688 cs = sqrt(max(2.0*tiny(cs),gamma*pressure/density))
2689 END FUNCTION getsoundspeed
2690
2693 ELEMENTAL SUBROUTINE setroeaverages(gamma,rhoL,rhoR,ul,uR,vL,vR,pL,pR,eL,eR,u,cs)
2694 IMPLICIT NONE
2695 !------------------------------------------------------------------------!
2696 REAL, INTENT(IN) :: gamma,rhol,rhor,ul,ur,vl,vr,pl,pr,el,er
2697 REAL, INTENT(OUT) :: u,cs
2698 !------------------------------------------------------------------------!
2699 REAL :: sqrtrhol,sqrtrhor,invsqrtrho,v,hl,hr,h
2700 !------------------------------------------------------------------------!
2701 sqrtrhol = sqrt(rhol)
2702 sqrtrhor = sqrt(rhor)
2703 ! density
2704 invsqrtrho = 1./ (sqrtrhol + sqrtrhor)
2705 ! velocities
2706 u = (sqrtrhol*ul + sqrtrhor*ur) * invsqrtrho
2707 v = (sqrtrhol*vl + sqrtrhor*vr) * invsqrtrho
2708 ! enthalpy
2709 hl = (el + pl) / rhol
2710 hr = (er + pr) / rhor
2711 h = (sqrtrhol * hl + sqrtrhor * hr) * invsqrtrho
2712 ! sound speed
2713 cs = sqrt((gamma-1.)*(h-0.5*(u**2+v**2)))
2714 END SUBROUTINE setroeaverages
2715
2717 ELEMENTAL SUBROUTINE setwavespeeds(cs,v,minwav,maxwav)
2718 IMPLICIT NONE
2719 !------------------------------------------------------------------------!
2720 REAL, INTENT(IN) :: cs,v
2721 REAL, INTENT(OUT) :: minwav,maxwav
2722 !------------------------------------------------------------------------!
2723 ! minimal and maximal wave speeds
2724 minwav = min(0.,v-cs)
2725 maxwav = max(0.,v+cs)
2726 END SUBROUTINE setwavespeeds
2727
2729 ELEMENTAL SUBROUTINE seteigenvalues1d(gamma,rho,v,P,l1,l2,l3)
2730 IMPLICIT NONE
2731 !------------------------------------------------------------------------!
2732 REAL, INTENT(IN) :: gamma,rho,v,p
2733 REAL, INTENT(OUT) :: l1,l2,l3
2734 !------------------------------------------------------------------------!
2735 REAL :: cs
2736 !------------------------------------------------------------------------!
2737 ! adiabatic sound speed
2738 cs = getsoundspeed(gamma,rho,p)
2739 ! call subroutine for isothermal case with the adiabatic sound speed
2740 l1 = v - cs
2741 l2 = v
2742 l3 = v + cs
2743 END SUBROUTINE seteigenvalues1d
2744
2746 ELEMENTAL SUBROUTINE seteigenvalues2d(gamma,rho,v,P,l1,l2,l3,l4)
2747 IMPLICIT NONE
2748 !------------------------------------------------------------------------!
2749 REAL, INTENT(IN) :: gamma,rho,v,p
2750 REAL, INTENT(OUT) :: l1,l2,l3,l4
2751 !------------------------------------------------------------------------!
2752 CALL seteigenvalues1d(gamma,rho,v,p,l1,l2,l4)
2753 l3 = v
2754 END SUBROUTINE seteigenvalues2d
2755
2757 ELEMENTAL SUBROUTINE seteigenvalues3d(gamma,rho,v,P,l1,l2,l3,l4,l5)
2758 IMPLICIT NONE
2759 !------------------------------------------------------------------------!
2760 REAL, INTENT(IN) :: gamma,rho,v,p
2761 REAL, INTENT(OUT) :: l1,l2,l3,l4,l5
2762 !------------------------------------------------------------------------!
2763 CALL seteigenvalues2d(gamma,rho,v,p,l1,l2,l3,l5)
2764 l4 = v
2765 END SUBROUTINE seteigenvalues3d
2766
2767 ! \todo NOT VERIFIED
2768 !! only for HLLC fluxes
2769 ELEMENTAL SUBROUTINE setintermediatestate(rhoL,rhoR,uL,uR,vl,vR,&
2770 pL,pR,eL,eR,amin,amax,rho,mu,mv,e,a)
2771 IMPLICIT NONE
2772 !------------------------------------------------------------------------!
2773 REAL, INTENT(IN) :: rhol,rhor,ul,ur,vl,vr,pl,pr,el,er,amin,amax
2774 REAL, INTENT(OUT) :: rho,mu,mv,e,a
2775 !------------------------------------------------------------------------!
2776 REAL :: ql,qr
2777 !------------------------------------------------------------------------!
2778 ql = rhol * (ul-amin)
2779 qr = rhor * (ur-amax)
2780 ! wave speed in intermediate region
2781 a = (pr - pl + qr*ur - ql*ul) / (qr - ql)
2782 ! set left or right state vector depending on wave speed "a"
2783 IF (a.GT.0.0) THEN
2784 ! left state
2785 rho = ql / (a - amin)
2786 mu = rho * a
2787 mv = rho * vl
2788 e = rho * (el/rhol + (a - ul) * (a - pl/ql))
2789 ELSE
2790 ! right state
2791 rho = qr / (a - amax)
2792 mu = rho * a
2793 mv = rho * vr
2794 e = rho * (er/rhor + (a - ur) * (a - pr/qr))
2795 END IF
2796 END SUBROUTINE setintermediatestate
2797
2799 ELEMENTAL SUBROUTINE setcharvars1d(gamma,rho1,rho2,u1,u2,P1,P2,l1,l3, &
2800 xvar1,xvar2,xvar3)
2801 IMPLICIT NONE
2802 !------------------------------------------------------------------------!
2803 REAL, INTENT(IN) :: gamma,rho1,rho2,u1,u2,p1,p2,l1,l3
2804 REAL, INTENT(OUT) :: xvar1,xvar2,xvar3
2805 !------------------------------------------------------------------------!
2806 REAL :: gammadu_cs,dlnp
2807 !------------------------------------------------------------------------!
2808 dlnp = log(p2/p1) ! = LOG(P2)-LOG(P1)
2809 ! gamma/cs = 2*gamma / (v+cs - (v-cs))
2810 gammadu_cs = 2*gamma*(u2-u1) / (l3-l1)
2811 ! characteristic variables
2812 xvar1 = dlnp - gammadu_cs
2813 xvar2 = dlnp - gamma * log(rho2/rho1)
2814 xvar3 = dlnp + gammadu_cs
2815 END SUBROUTINE setcharvars1d
2816
2818 ELEMENTAL SUBROUTINE setcharvars2d(gamma,rho1,rho2,u1,u2,v1,v2,P1,P2,l1,l4, &
2819 xvar1,xvar2,xvar3,xvar4)
2820 IMPLICIT NONE
2821 !------------------------------------------------------------------------!
2822 REAL, INTENT(IN) :: gamma,rho1,rho2,u1,u2,v1,v2,p1,p2,l1,l4
2823 REAL, INTENT(OUT) :: xvar1,xvar2,xvar3,xvar4
2824 !------------------------------------------------------------------------!
2825 CALL setcharvars1d(gamma,rho1,rho2,u1,u2,p1,p2,l1,l4,xvar1,xvar2,xvar4)
2826 xvar3 = (v2-v1)
2827 END SUBROUTINE setcharvars2d
2828
2830 ELEMENTAL SUBROUTINE setcharvars3d(gamma,rho1,rho2,u1,u2,v1,v2,w1,w2,P1,P2, &
2831 l1,l5,xvar1,xvar2,xvar3,xvar4,xvar5)
2832 IMPLICIT NONE
2833 !------------------------------------------------------------------------!
2834 REAL, INTENT(IN) :: gamma,rho1,rho2,u1,u2,v1,v2,w1,w2,p1,p2,l1,l5
2835 REAL, INTENT(OUT) :: xvar1,xvar2,xvar3,xvar4,xvar5
2836 !------------------------------------------------------------------------!
2837 CALL setcharvars2d(gamma,rho1,rho2,u1,u2,v1,v2,p1,p2,l1,l5,xvar1,xvar2,xvar3,xvar5)
2838 xvar4 = (w2-w1)
2839 END SUBROUTINE setcharvars3d
2840
2843 ELEMENTAL SUBROUTINE setboundarydata1d(delta,gamma,rho1,u1,P1,xvar1, &
2844 xvar2,xvar3,rho2,u2,P2)
2845 IMPLICIT NONE
2846 !------------------------------------------------------------------------!
2847 INTEGER, INTENT(IN) :: delta
2848 REAL, INTENT(IN) :: gamma,rho1,u1,p1,xvar1,xvar2,xvar3
2849 REAL, INTENT(OUT) :: rho2,u2,p2
2850 !------------------------------------------------------------------------!
2851 REAL :: dlnp,cs
2852 !------------------------------------------------------------------------!
2853 dlnp = 0.5*(xvar3+xvar1)
2854 ! extrapolate boundary values using characteristic variables
2855 rho2 = rho1 * exp(delta/gamma*(dlnp-xvar2))
2856 p2 = p1 * exp(delta*dlnp)
2857 ! compute sound speed with arithmetic mean of density and pressure
2858 ! (the factor 0.5 cancels out, because cs depends on the quotient P/rho)
2859 cs = getsoundspeed(gamma,rho1+rho2,p1+p2)
2860 u2 = u1 + 0.5*delta*cs/gamma * (xvar3-xvar1)
2861 END SUBROUTINE setboundarydata1d
2862
2865 ELEMENTAL SUBROUTINE setboundarydata2d(delta,gamma,rho1,u1,v1,P1,xvar1, &
2866 xvar2,xvar3,xvar4,rho2,u2,v2,P2)
2867 IMPLICIT NONE
2868 !------------------------------------------------------------------------!
2869 INTEGER, INTENT(IN) :: delta
2870 REAL, INTENT(IN) :: gamma,rho1,u1,v1,p1,xvar1,xvar2,xvar3,xvar4
2871 REAL, INTENT(OUT) :: rho2,u2,v2,p2
2872 !------------------------------------------------------------------------!
2873 CALL setboundarydata1d(delta,gamma,rho1,u1,p1,xvar1,xvar2,xvar4,rho2,u2,p2)
2874 v2 = v1 + delta*xvar3
2875 END SUBROUTINE setboundarydata2d
2876
2879 ELEMENTAL SUBROUTINE setboundarydata3d(delta,gamma,rho1,u1,v1,w1,P1,xvar1, &
2880 xvar2,xvar3,xvar4,xvar5,rho2,u2,v2,w2,P2)
2881 IMPLICIT NONE
2882 !------------------------------------------------------------------------!
2883 INTEGER, INTENT(IN) :: delta
2884 REAL, INTENT(IN) :: gamma,rho1,u1,v1,w1,p1,xvar1,xvar2,xvar3,xvar4,xvar5
2885 REAL, INTENT(OUT) :: rho2,u2,v2,w2,p2
2886 !------------------------------------------------------------------------!
2887 CALL setboundarydata2d(delta,gamma,rho1,u1,v1,p1,xvar1,xvar2,xvar3,xvar5,rho2,u2,v2,p2)
2888 w2 = w1 + delta*xvar4
2889 END SUBROUTINE setboundarydata3d
2890
2891 ! \todo NOT VERIFIED
2892 !! only for farfield boundary conditions
2893 ELEMENTAL SUBROUTINE prim2riemann1d(gamma,rho,vx,p,l1,l2,l3,Rminus,Rs,Rplus)
2894 IMPLICIT NONE
2895 !------------------------------------------------------------------------!
2896 REAL, INTENT(IN) :: gamma,rho,vx,p,l1,l2,l3
2897 REAL, INTENT(OUT) :: rminus,rs,rplus
2898 !------------------------------------------------------------------------!
2899 REAL :: cs
2900 !------------------------------------------------------------------------!
2901 cs = l3-l2 ! l2 = v, l3 = v+cs
2902 ! compute 1st Riemann invariant (R+)
2903 rplus = vx + 2./(gamma-1.0) * cs
2904 ! compute 2st Riemann invariant (R-)
2905 rminus = vx - 2./(gamma-1.0) * cs
2906 ! compute entropy
2907 rs = p/rho**gamma
2908 END SUBROUTINE prim2riemann1d
2909
2910
2911 ! \todo NOT VERIFIED
2912 !! only for farfield boundary conditions
2913 ELEMENTAL SUBROUTINE prim2riemann2d(gamma,rho,vx,vy,p,&
2914 l1,l2,l3,l4,Rminus,Rs,Rvt,Rplus)
2915 IMPLICIT NONE
2916 !------------------------------------------------------------------------!
2917 REAL, INTENT(IN) :: gamma,rho,vx,vy,p,l1,l2,l3,l4
2918 REAL, INTENT(OUT) :: rminus,rs,rvt,rplus
2919 !------------------------------------------------------------------------!
2920 REAL :: cs
2921 !------------------------------------------------------------------------!
2922 CALL prim2riemann1d(gamma,rho,vx,p,l1,l2,l4,rminus,rs,rplus)
2923 ! tangential velocities
2924 rvt = vy
2925 END SUBROUTINE prim2riemann2d
2926
2927
2928 ! \todo NOT VERIFIED
2929 !! only for farfield boundary conditions
2930 ELEMENTAL SUBROUTINE prim2riemann3d(gamma,rho,vx,vy,vz,p,&
2931 l1,l2,l3,l4,l5,Rminus,Rs,Rvt,Rwt,Rplus)
2932 IMPLICIT NONE
2933 !------------------------------------------------------------------------!
2934 REAL, INTENT(IN) :: gamma,rho,vx,vy,vz,p,l1,l2,l3,l4,l5
2935 REAL, INTENT(OUT) :: rminus,rs,rvt,rwt,rplus
2936 !------------------------------------------------------------------------!
2937 REAL :: cs
2938 !------------------------------------------------------------------------!
2939 CALL prim2riemann2d(gamma,rho,vx,vy,p,l1,l2,l3,l5,rminus,rs,rvt,rplus)
2940 rwt = vz
2941 END SUBROUTINE prim2riemann3d
2942
2943 ! \todo NOT VERIFIED
2944 !! only for farfield boundary conditions
2945 ELEMENTAL SUBROUTINE riemann2prim1d(gamma,Rminus,Rs,Rplus,rho,vx,p)
2946 IMPLICIT NONE
2947 !------------------------------------------------------------------------!
2948 REAL, INTENT(IN) :: gamma,rminus,rs,rplus
2949 REAL, INTENT(OUT) :: rho,vx,p
2950 !------------------------------------------------------------------------!
2951 REAL :: cs2gam
2952 !------------------------------------------------------------------------!
2953 ! normal velocity
2954 vx = 0.5*(rplus+rminus)
2955 ! cs**2 / gamma
2956 cs2gam = (0.25*(gamma-1.0)*(rplus-rminus))**2 / gamma
2957 ! density
2958 rho = (cs2gam/rs)**(1./(gamma-1.0))
2959 ! pressure
2960 p = cs2gam * rho
2961 END SUBROUTINE riemann2prim1d
2962
2963 ! \todo NOT VERIFIED
2964 !! only for farfield boundary conditions
2965 ELEMENTAL SUBROUTINE riemann2prim2d(gamma,Rminus,Rs,Rvt,Rplus,&
2966 rho,vx,vy,p)
2967 IMPLICIT NONE
2968 !------------------------------------------------------------------------!
2969 REAL, INTENT(IN) :: gamma,rminus,rs,rvt,rplus
2970 REAL, INTENT(OUT) :: rho,vx,vy,p
2971 !------------------------------------------------------------------------!
2972 REAL :: cs2gam
2973 !------------------------------------------------------------------------!
2974 ! tangential velocity
2975 vy = rvt
2976
2977 CALL riemann2prim1d(gamma,rminus,rs,rplus,rho,vx,p)
2978 END SUBROUTINE riemann2prim2d
2979
2980
2981 ! \todo NOT VERIFIED
2982 !! only for farfield boundary conditions
2983 ELEMENTAL SUBROUTINE riemann2prim3d(gamma,Rminus,Rs,Rvt,Rwt,Rplus,&
2984 rho,vx,vy,vz,p)
2985 IMPLICIT NONE
2986 !------------------------------------------------------------------------!
2987 REAL, INTENT(IN) :: gamma,rminus,rs,rvt,rwt,rplus
2988 REAL, INTENT(OUT) :: rho,vx,vy,vz,p
2989 !------------------------------------------------------------------------!
2990 REAL :: cs2gam
2991 !------------------------------------------------------------------------!
2992 ! tangential velocity
2993 vz = rwt
2994
2995 CALL riemann2prim2d(gamma,rminus,rs,rvt,rplus,rho,vx,vy,p)
2996 END SUBROUTINE riemann2prim3d
2997
2999 ELEMENTAL SUBROUTINE setflux1d(rho,u,P,mu,E,f1,f2,f3)
3000 IMPLICIT NONE
3001 !------------------------------------------------------------------------!
3002 REAL, INTENT(IN) :: rho,u,p,mu,e
3003 REAL, INTENT(OUT) :: f1,f2,f3
3004 !------------------------------------------------------------------------!
3005 f1 = rho*u
3006 f2 = mu*u + p
3007 f3 = (e+p)*u
3008 END SUBROUTINE setflux1d
3009
3011 ELEMENTAL SUBROUTINE setflux2d(rho,u,P,mu,mv,E,f1,f2,f3,f4)
3012 IMPLICIT NONE
3013 !------------------------------------------------------------------------!
3014 REAL, INTENT(IN) :: rho,u,p,mu,mv,e
3015 REAL, INTENT(OUT) :: f1,f2,f3,f4
3016 !------------------------------------------------------------------------!
3017 CALL setflux1d(rho,u,p,mu,e,f1,f2,f4)
3018 f3 = mv*u
3019 END SUBROUTINE setflux2d
3020
3022 ELEMENTAL SUBROUTINE setflux3d(rho,u,P,mu,mv,mw,E,f1,f2,f3,f4,f5)
3023 IMPLICIT NONE
3024 !------------------------------------------------------------------------!
3025 REAL, INTENT(IN) :: rho,u,p,mu,mv,mw,e
3026 REAL, INTENT(OUT) :: f1,f2,f3,f4,f5
3027 !------------------------------------------------------------------------!
3028 CALL setflux2d(rho,u,p,mu,mv,e,f1,f2,f3,f5)
3029 f4 = mw*u
3030 END SUBROUTINE setflux3d
3031
3033 ELEMENTAL SUBROUTINE cons2prim1d(gamma,rho_in,mu,E,rho_out,u,P)
3034 IMPLICIT NONE
3035 !------------------------------------------------------------------------!
3036 REAL, INTENT(IN) :: gamma,rho_in,mu,e
3037 REAL, INTENT(OUT) :: rho_out,u,p
3038 !------------------------------------------------------------------------!
3039 REAL :: inv_rho
3040 !------------------------------------------------------------------------!
3041 inv_rho = 1./rho_in
3042 rho_out = rho_in
3043 u = mu * inv_rho
3044 p = (gamma-1.)*(e - 0.5 * inv_rho * mu*mu)
3045 END SUBROUTINE cons2prim1d
3046
3048 ELEMENTAL SUBROUTINE cons2prim2d(gamma,rho_in,mu,mv,E,rho_out,u,v,P)
3049 IMPLICIT NONE
3050 !------------------------------------------------------------------------!
3051 REAL, INTENT(IN) :: gamma,rho_in,mu,mv,e
3052 REAL, INTENT(OUT) :: rho_out,u,v,p
3053 !------------------------------------------------------------------------!
3054 REAL :: inv_rho
3055 !------------------------------------------------------------------------!
3056 inv_rho = 1./rho_in
3057 rho_out = rho_in
3058 u = mu * inv_rho
3059 v = mv * inv_rho
3060 p = (gamma-1.)*(e - 0.5 * inv_rho * (mu*mu + mv*mv))
3061 END SUBROUTINE cons2prim2d
3062
3064 ELEMENTAL SUBROUTINE cons2prim3d(gamma,rho_in,mu,mv,mw,E,rho_out,u,v,w,P)
3065 IMPLICIT NONE
3066 !------------------------------------------------------------------------!
3067 REAL, INTENT(IN) :: gamma,rho_in,mu,mv,mw,e
3068 REAL, INTENT(OUT) :: rho_out,u,v,w,p
3069 !------------------------------------------------------------------------!
3070 REAL :: inv_rho
3071 !------------------------------------------------------------------------!
3072 inv_rho = 1./rho_in
3073 rho_out = rho_in
3074 u = mu * inv_rho
3075 v = mv * inv_rho
3076 w = mw * inv_rho
3077 p = (gamma-1.)*(e - 0.5 * inv_rho * (mu*mu + mv*mv + mw*mw))
3078 END SUBROUTINE cons2prim3d
3079
3081 ELEMENTAL SUBROUTINE prim2cons1d(gamma,rho_in,u,P,rho_out,mu,E)
3082 IMPLICIT NONE
3083 !------------------------------------------------------------------------!
3084 REAL, INTENT(IN) :: gamma,rho_in,u,p
3085 REAL, INTENT(OUT) :: rho_out,mu,e
3086 !------------------------------------------------------------------------!
3087 rho_out = rho_in
3088 mu = rho_in * u
3089 e = p/(gamma-1.) + 0.5 * rho_in * u*u
3090 END SUBROUTINE prim2cons1d
3091
3093 ELEMENTAL SUBROUTINE prim2cons2d(gamma,rho_in,u,v,P,rho_out,mu,mv,E)
3094 IMPLICIT NONE
3095 !------------------------------------------------------------------------!
3096 REAL, INTENT(IN) :: gamma,rho_in,u,v,p
3097 REAL, INTENT(OUT) :: rho_out,mu,mv,e
3098 !------------------------------------------------------------------------!
3099 rho_out = rho_in
3100 mu = rho_in * u
3101 mv = rho_in * v
3102 e = p/(gamma-1.) + 0.5 * rho_in * (u*u + v*v)
3103 END SUBROUTINE prim2cons2d
3104
3106 ELEMENTAL SUBROUTINE prim2cons3d(gamma,rho_in,u,v,w,P,rho_out,mu,mv,mw,E)
3107 IMPLICIT NONE
3108 !------------------------------------------------------------------------!
3109 REAL, INTENT(IN) :: gamma,rho_in,u,v,w,p
3110 REAL, INTENT(OUT) :: rho_out,mu,mv,mw,e
3111 !------------------------------------------------------------------------!
3112 rho_out = rho_in
3113 mu = rho_in * u
3114 mv = rho_in * v
3115 mw = rho_in * w
3116 e = p/(gamma-1.) + 0.5 * rho_in * (u*u + v*v + w*w)
3117 END SUBROUTINE prim2cons3d
3118
3128 ELEMENTAL FUNCTION getgeometricalsourcex(cxyx,cxzx,cyxy,czxz,vx,vy,vz,P,my,mz)
3129 IMPLICIT NONE
3130 !------------------------------------------------------------------------!
3131 REAL, INTENT(IN) :: cxyx,cxzx,cyxy,czxz,vx,vy,vz,p,my,mz
3132 REAL :: getgeometricalsourcex
3133 !------------------------------------------------------------------------!
3134 getgeometricalsourcex = my * (cyxy*vy - cxyx*vx) + mz * (czxz*vz - cxzx*vx) &
3135 + (cyxy + czxz) * p
3136 END FUNCTION getgeometricalsourcex
3137
3139 ELEMENTAL FUNCTION getgeometricalsourcey(cxyx,cyxy,cyzy,czyz,vx,vy,vz,P,mx,mz)
3140 IMPLICIT NONE
3141 !------------------------------------------------------------------------!
3142 REAL, INTENT(IN) :: cxyx,cyxy,cyzy,czyz,vx,vy,vz,p,mx,mz
3143 REAL :: getgeometricalsourcey
3144 !------------------------------------------------------------------------!
3145 getgeometricalsourcey = mz * (czyz*vz - cyzy*vy) + mx * (cxyx*vx - cyxy*vy) &
3146 + (cxyx + czyz) * p
3147 END FUNCTION getgeometricalsourcey
3148
3150 ELEMENTAL FUNCTION getgeometricalsourcez(cxzx,cyzy,czxz,czyz,vx,vy,vz,P,mx,my)
3151 IMPLICIT NONE
3152 !------------------------------------------------------------------------!
3153 REAL, INTENT(IN) :: cxzx,cyzy,czxz,czyz,vx,vy,vz,p,mx,my
3154 REAL :: getgeometricalsourcez
3155 !------------------------------------------------------------------------!
3156 getgeometricalsourcez = mx * (cxzx*vx - czxz*vz) + my * (cyzy*vy - czyz*vz) &
3157 + (cxzx + cyzy) * p
3158 END FUNCTION getgeometricalsourcez
3159
3160END MODULE physics_euler_mod
physics_base_mod::new_statevector
Definition: physics_base.f90:195
physics_euler_mod::cons2prim
Definition: physics_euler.f90:119
physics_euler_mod::prim2cons
Definition: physics_euler.f90:122
physics_euler_mod::setflux
Definition: physics_euler.f90:116
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
geometry_generic_mod
constructor for geometry class
Definition: geometry_generic.f90:35
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
logging_base_mod::nan_default_real
real, save nan_default_real
NaN real constant.
Definition: logging_base.f90:85
logging_base_mod::finalize
subroutine finalize(this)
Destructor of logging_base class.
Definition: logging_base.f90:160
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
marray_base_mod::assignmarray_0
subroutine assignmarray_0(this, ma)
assigns one mesh array to another mesh array
Definition: marray_base.f90:346
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
mesh_base_mod::vector_y
integer, parameter vector_y
Definition: mesh_base.f90:109
mesh_base_mod::vector_z
integer, parameter vector_z
Definition: mesh_base.f90:109
mesh_base_mod::vector_x
integer, parameter vector_x
flags to check which vector components are enabled
Definition: mesh_base.f90:109
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
physics_base_mod::primitive
@, public primitive
Definition: physics_base.f90:67
physics_base_mod::euler
integer, parameter, public euler
Definition: physics_base.f90:498
physics_base_mod::conservative
@, public conservative
Definition: physics_base.f90:67
physics_euler_mod
physics module for 1D,2D and 3D non-isothermal Euler equations
Definition: physics_euler.f90:41
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:3000
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:3049
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:3023
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:1073
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:244
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:2832
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:3034
physics_euler_mod::calcfluxesz
pure subroutine calcfluxesz(this, Mesh, nmin, nmax, prim, cons, zfluxes)
Calculate Fluxes in z-direction.
Definition: physics_euler.f90:605
physics_euler_mod::calculatecharsystemy
pure subroutine calculatecharsystemy(this, Mesh, j1, j2, pvar, lambda, xvar)
Definition: physics_euler.f90:849
physics_euler_mod::updatesoundspeed
pure subroutine updatesoundspeed(this, pvar)
Definition: physics_euler.f90:2556
physics_euler_mod::addfargosourcesx
pure subroutine addfargosourcesx(this, Mesh, w, pvar, cvar, sterm)
sources terms for fargo advection along x-direction
Definition: physics_euler.f90:2268
physics_euler_mod::setwavespeeds
elemental subroutine setwavespeeds(cs, v, minwav, maxwav)
set minimal and maximal wave speeds
Definition: physics_euler.f90:2718
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:3107
physics_euler_mod::subtractbackgroundvelocityz
pure subroutine subtractbackgroundvelocityz(this, Mesh, w, pvar, cvar)
Substracts a background velocity field for fargo routines.
Definition: physics_euler.f90:2218
physics_euler_mod::calculateprim2riemanny
pure subroutine calculateprim2riemanny(this, Mesh, j, pvar, lambda, Rinv)
Conversion from primitive to riemann invariants for farfield boundaries.
Definition: physics_euler.f90:1361
physics_euler_mod::addfargosourcesz
pure subroutine addfargosourcesz(this, Mesh, w, pvar, cvar, sterm)
sources terms for fargo advection along z-direction
Definition: physics_euler.f90:2469
physics_euler_mod::calculateboundarydatay
pure subroutine calculateboundarydatay(this, Mesh, j1, j2, xvar, pvar)
Definition: physics_euler.f90:1132
physics_euler_mod::geometricalsources
pure subroutine geometricalsources(this, Mesh, pvar, cvar, sterm)
Calculates geometrical sources.
Definition: physics_euler.f90:1707
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:3094
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:3012
physics_euler_mod::addbackgroundvelocityz
pure subroutine addbackgroundvelocityz(this, Mesh, w, pvar, cvar)
Adds a background velocity field for fargo routines.
Definition: physics_euler.f90:2060
physics_euler_mod::problem_name
character(len=32), parameter problem_name
Definition: physics_euler.f90:52
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:3129
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:3140
physics_euler_mod::addbackgroundvelocityy
pure subroutine addbackgroundvelocityy(this, Mesh, w, pvar, cvar)
Adds a background velocity field for fargo routines.
Definition: physics_euler.f90:2004
physics_euler_mod::riemann2prim1d
elemental subroutine riemann2prim1d(gamma, Rminus, Rs, Rplus, rho, vx, p)
Definition: physics_euler.f90:2946
physics_euler_mod::getsoundspeed
elemental real function getsoundspeed(gamma, density, pressure)
Definition: physics_euler.f90:2683
physics_euler_mod::calcfluxesy
pure subroutine calcfluxesy(this, Mesh, nmin, nmax, prim, cons, yfluxes)
Calculate Fluxes in y-direction.
Definition: physics_euler.f90:549
physics_euler_mod::calculateriemann2primy
pure subroutine calculateriemann2primy(this, Mesh, j, Rinv, pvar)
Convert Riemann invariants to primitives for farfield boundaries.
Definition: physics_euler.f90:1607
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:2801
physics_euler_mod::riemann2prim3d
elemental subroutine riemann2prim3d(gamma, Rminus, Rs, Rvt, Rwt, Rplus, rho, vx, vy, vz, p)
Definition: physics_euler.f90:2985
physics_euler_mod::finalize_statevector
subroutine finalize_statevector(this)
actual destructor of the statevector_eulerisotherm type
Definition: physics_euler.f90:2667
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:2758
physics_euler_mod::calculateriemann2primx
pure subroutine calculateriemann2primx(this, Mesh, i, Rinv, pvar)
Convert Riemann invariants to primitives for farfield boundaries.
Definition: physics_euler.f90:1556
physics_euler_mod::subtractbackgroundvelocityy
pure subroutine subtractbackgroundvelocityy(this, Mesh, w, pvar, cvar)
Substracts a background velocity field for fargo routines.
Definition: physics_euler.f90:2160
physics_euler_mod::initphysics
subroutine initphysics(this, Mesh, config, IO)
constructor of physics_euler class
Definition: physics_euler.f90:136
physics_euler_mod::calculateprim2riemannx
pure subroutine calculateprim2riemannx(this, Mesh, i, pvar, lambda, Rinv)
Conversion from primitive to riemann invariants for farfield boundaries.
Definition: physics_euler.f90:1267
physics_euler_mod::prim2riemann3d
elemental subroutine prim2riemann3d(gamma, rho, vx, vy, vz, p, l1, l2, l3, l4, l5, Rminus, Rs, Rvt, Rwt, Rplus)
Definition: physics_euler.f90:2932
physics_euler_mod::calcfluxesx
pure subroutine calcfluxesx(this, Mesh, nmin, nmax, prim, cons, xfluxes)
Calculate Fluxes in x-direction.
Definition: physics_euler.f90:493
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:2820
physics_euler_mod::calculateriemann2primz
pure subroutine calculateriemann2primz(this, Mesh, k, Rinv, pvar)
Convert Riemann invariants to primitives for farfield boundaries.
Definition: physics_euler.f90:1658
physics_euler_mod::calculateprim2riemannz
pure subroutine calculateprim2riemannz(this, Mesh, k, pvar, lambda, Rinv)
Conversion from primitive to riemann invariants for farfield boundaries.
Definition: physics_euler.f90:1457
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:368
physics_euler_mod::viscositysources
subroutine viscositysources(this, Mesh, pvar, btxx, btxy, btxz, btyy, btyz, btzz, sterm)
Definition: physics_euler.f90:1888
physics_euler_mod::riemann2prim2d
elemental subroutine riemann2prim2d(gamma, Rminus, Rs, Rvt, Rplus, rho, vx, vy, p)
Definition: physics_euler.f90:2967
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:3082
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:2881
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:3065
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:3151
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:282
physics_euler_mod::prim2riemann2d
elemental subroutine prim2riemann2d(gamma, rho, vx, vy, p, l1, l2, l3, l4, Rminus, Rs, Rvt, Rplus)
Definition: physics_euler.f90:2915
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:2730
physics_euler_mod::addfargosourcesy
pure subroutine addfargosourcesy(this, Mesh, w, pvar, cvar, sterm)
sources terms for fargo advection along y-direction
Definition: physics_euler.f90:2369
physics_euler_mod::externalsources
subroutine externalsources(this, accel, pvar, cvar, sterm)
compute momentum and energy sources given an external force
Definition: physics_euler.f90:1850
physics_euler_mod::addbackgroundvelocityx
pure subroutine addbackgroundvelocityx(this, Mesh, w, pvar, cvar)
Adds a background velocity field for fargo routines.
Definition: physics_euler.f90:1956
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:406
physics_euler_mod::subtractbackgroundvelocityx
pure subroutine subtractbackgroundvelocityx(this, Mesh, w, pvar, cvar)
Substracts a background velocity field for fargo routines.
Definition: physics_euler.f90:2110
physics_euler_mod::calculateboundarydataz
pure subroutine calculateboundarydataz(this, Mesh, k1, k2, xvar, pvar)
Definition: physics_euler.f90:1205
physics_euler_mod::createstatevector
type(statevector_euler) function createstatevector(Physics, flavour, num)
Constructor of statevector_euler.
Definition: physics_euler.f90:2594
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:2747
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:2771
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:2845
physics_euler_mod::setroeaverages
elemental subroutine setroeaverages(gamma, rhoL, rhoR, ul, uR, vL, vR, pL, pR, eL, eR, u, cs)
Definition: physics_euler.f90:2694
physics_euler_mod::prim2riemann1d
elemental subroutine prim2riemann1d(gamma, rho, vx, p, l1, l2, l3, Rminus, Rs, Rplus)
Definition: physics_euler.f90:2894
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:2867
physics_euler_mod::calculatecharsystemx
pure subroutine calculatecharsystemx(this, Mesh, i1, i2, pvar, lambda, xvar)
Definition: physics_euler.f90:744
physics_euler_mod::calculatecharsystemz
pure subroutine calculatecharsystemz(this, Mesh, k1, k2, pvar, lambda, xvar)
Definition: physics_euler.f90:967
physics_eulerisotherm_mod
physics module for 1D,2D and 3D isothermal Euler equations
Definition: physics_eulerisotherm.f90:44
common_dict::dict_typ
Definition: common_dict.f90:123
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:69
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
physics_euler_mod::physics_euler
Definition: physics_euler.f90:54
physics_euler_mod::statevector_euler
Definition: physics_euler.f90:105
physics_eulerisotherm_mod::physics_eulerisotherm
Definition: physics_eulerisotherm.f90:56
physics_eulerisotherm_mod::statevector_eulerisotherm
Definition: physics_eulerisotherm.f90:122
Generated by   doxygen 1.9.4