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