timedisc_base.f90
Go to the documentation of this file.
1!r#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: timedisc_base.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!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
11!# #
12!# This program is free software; you can redistribute it and/or modify #
13!# it under the terms of the GNU General Public License as published by #
14!# the Free Software Foundation; either version 2 of the License, or (at #
15!# your option) any later version. #
16!# #
17!# This program is distributed in the hope that it will be useful, but #
18!# WITHOUT ANY WARRANTY; without even the implied warranty of #
19!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
20!# NON INFRINGEMENT. See the GNU General Public License for more #
21!# details. #
22!# #
23!# You should have received a copy of the GNU General Public License #
24!# along with this program; if not, write to the Free Software #
25!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
26!# #
27!#############################################################################
49!----------------------------------------------------------------------------!
56!----------------------------------------------------------------------------!
57MODULE timedisc_base_mod
58 USE logging_base_mod
59 USE boundary_generic_mod
60 USE mesh_base_mod
61 USE marray_compound_mod
62 USE marray_base_mod
63 USE physics_base_mod
64 USE physics_generic_mod
65 USE sources_base_mod
66 USE sources_generic_mod
67 USE sources_gravity_mod
68 USE fluxes_base_mod
69 USE reconstruction_base_mod
70 USE common_dict
71#ifdef PARALLEL
72#ifdef HAVE_MPI_MOD
73 USE mpi
74#endif
75#endif
76 IMPLICIT NONE
77#ifdef PARALLEL
78#ifdef HAVE_MPIF_H
79 include 'mpif.h'
80#endif
81#endif
82!----------------------------------------------------------------------------!
83PRIVATE
84 TYPE, ABSTRACT, EXTENDS (logging_base) :: timedisc_base
86 CLASS(boundary_generic), ALLOCATABLE :: boundary
87 CLASS(marray_compound), POINTER &
88 :: pvar => null(),cvar => null(), &
89 ptmp => null(),ctmp => null(), &
90 cold => null(), &
91 src => null(),geo_src => null(), &
92 rhs => null(), &
93 cerr => null(),cerr_max => null(), &
94 solution => null()
95 INTEGER :: order
96 REAL :: cfl
97 REAL :: dt
98 REAL :: dtold
99 REAL :: dtmin
100 INTEGER :: dtcause,dtmincause
101 REAL :: dtmax
102 REAL,POINTER :: dtmean, dtstddev
103 INTEGER :: dtaccept
104 REAL,POINTER :: time
105 REAL :: stoptime
106 REAL :: dtlimit
107 REAL :: tol_rel
108 REAL :: maxerrold
109 LOGICAL :: break
110 LOGICAL :: always_update_bccsound
112 INTEGER :: maxiter
113 INTEGER :: n_adj
114 INTEGER :: m
115 INTEGER :: degree
116 INTEGER :: rhstype
117 INTEGER :: checkdatabm
118 REAL :: pmin,rhomin,tmin
120 TYPE(marray_base) :: dq, flux, dql
122 REAL :: err_n, h_n
123 REAL, DIMENSION(:), POINTER :: tol_abs
124 REAL :: beta
125
127 REAL, DIMENSION(:,:,:,:), POINTER :: phi,oldphi_s,&
128 newphi_s
129 REAL, DIMENSION(:), POINTER :: gamma
130 INTEGER :: pc
132 REAL, DIMENSION(:,:,:,:), POINTER :: xfluxdydz=>null(),yfluxdzdx=>null(),zfluxdxdy=>null()
133 REAL, DIMENSION(:,:,:,:), POINTER :: amax=>null()
134 REAL, DIMENSION(:,:), POINTER :: bflux=>null()
135 LOGICAL :: write_error
136 INTEGER, DIMENSION(:,:), POINTER :: shift=>null()
137 REAL, DIMENSION(:,:), POINTER :: buf=>null()
138 REAL, DIMENSION(:,:), POINTER :: w=>null()
139 REAL, DIMENSION(:,:), POINTER :: wfactor=>null()
140 REAL, DIMENSION(:,:), POINTER :: delxy =>null()
141
142 CONTAINS
143
144 PROCEDURE :: inittimedisc
145 PROCEDURE :: setoutput
146 PROCEDURE :: integrationstep
147 PROCEDURE :: computerhs
148 PROCEDURE :: adjusttimestep
149 PROCEDURE :: calctimestep
150 PROCEDURE :: acceptsolution
151 PROCEDURE :: rejectsolution
152 PROCEDURE :: checkdata
153 PROCEDURE :: computeerror
154 PROCEDURE :: getcentrifugalvelocity
155 PROCEDURE :: getorder
156 PROCEDURE :: getcfl
157 PROCEDURE :: fargoadvectionx
158 PROCEDURE :: fargoadvectiony
159 PROCEDURE :: fargoadvectionz
160 PROCEDURE :: calcbackgroundvelocity
161 procedure(solveode), DEFERRED :: solveode
162 PROCEDURE :: finalize_base
163 procedure(finalize), DEFERRED :: finalize
164 END TYPE timedisc_base
165!----------------------------------------------------------------------------!
166 abstract INTERFACE
167 SUBROUTINE solveode(this,Mesh,Physics,Sources,Fluxes,time,dt,err)
168 IMPORT timedisc_base, mesh_base, physics_base, fluxes_base, sources_base, sources_list
169 IMPLICIT NONE
170 CLASS(timedisc_base), INTENT(INOUT) :: this
171 CLASS(mesh_base), INTENT(IN) :: Mesh
172 CLASS(physics_base), INTENT(INOUT) :: Physics
173 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
174 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
175 REAL, INTENT(IN) :: time
176 REAL, INTENT(INOUT) :: dt, err
177 END SUBROUTINE
178 SUBROUTINE finalize(this)
179 IMPORT timedisc_base
180 IMPLICIT NONE
181 !------------------------------------------------------------------------!
182 CLASS(timedisc_base) :: this
183 END SUBROUTINE
184
185 END INTERFACE
186
187 INTEGER, PARAMETER :: modified_euler = 1
188 INTEGER, PARAMETER :: rk_fehlberg = 2
189 INTEGER, PARAMETER :: cash_karp = 3
190 INTEGER, PARAMETER :: dormand_prince = 4
191 INTEGER, PARAMETER :: ssprk = 5
192 !--------------------------------------------------------------------------!
193 INTEGER, PARAMETER :: dtcause_cfl = 0 ! smallest ts due to cfl cond. !
194 INTEGER, PARAMETER :: dtcause_erradj = -1 ! smallest ts due to err adj. !
195 INTEGER, PARAMETER :: dtcause_smallerr = -2 ! smallest ts due to err
196 INTEGER, PARAMETER :: dtcause_fileio = -4
197 ! CheckData_modeuler Bit Mask constants
198 INTEGER, PARAMETER :: check_nothing = int(b'0')
199 INTEGER, PARAMETER :: check_all = not(check_nothing)
200 INTEGER, PARAMETER :: check_csound = int(b'1')
201 INTEGER, PARAMETER :: check_pmin = int(b'10')
202 INTEGER, PARAMETER :: check_rhomin = int(b'100')
203 INTEGER, PARAMETER :: check_invalid = int(b'1000')
204 INTEGER, PARAMETER :: check_tmin = int(b'10000')
205 !--------------------------------------------------------------------------!
206 CHARACTER(LEN=40), PARAMETER, DIMENSION(3) :: fargo_method = (/ &
207 "dynamic background velocity ", &
208 "user supplied fixed background velocity ", &
209 "shearing box " /)
210 PUBLIC :: &
211 ! types
212 timedisc_base, &
213 ! constants
214 modified_euler, rk_fehlberg, cash_karp, dormand_prince, &
215 ssprk, &
216 dtcause_cfl,dtcause_erradj,dtcause_smallerr, dtcause_fileio, &
217 check_all, check_nothing, check_csound, check_rhomin, check_pmin, &
218 check_invalid, check_tmin
219 !--------------------------------------------------------------------------!
220
221CONTAINS
222
223
224 SUBROUTINE inittimedisc(this,Mesh,Physics,config,IO,ttype,tname)
225 USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
226 USE physics_euler_mod, ONLY : physics_euler
227 USE geometry_cylindrical_mod, ONLY: geometry_cylindrical
228 USE geometry_logcylindrical_mod, ONLY: geometry_logcylindrical
229 USE geometry_spherical_mod, ONLY: geometry_spherical
230 USE geometry_cartesian_mod, ONLY: geometry_cartesian
231 IMPLICIT NONE
232 !------------------------------------------------------------------------!
233 CLASS(timedisc_base), INTENT(INOUT) :: this
234 CLASS(mesh_base), INTENT(INOUT) :: Mesh
235 CLASS(physics_base), INTENT(IN) :: Physics
236 TYPE(dict_typ), POINTER :: config,IO
237 INTEGER, INTENT(IN) :: ttype
238 CHARACTER(LEN=32), INTENT(IN) :: tname
239 !------------------------------------------------------------------------!
240 INTEGER :: err, d, i,j,k
241 CHARACTER(LEN=32) :: order_str,cfl_str,stoptime_str,dtmax_str,beta_str
242 CHARACTER(LEN=32) :: info_str,shear_direction
243 INTEGER :: method
244 !------------------------------------------------------------------------!
245 CALL this%InitLogging(ttype,tname)
246
247 IF (.NOT.physics%Initialized().OR..NOT.mesh%Initialized()) &
248 CALL this%Error("InitTimedisc","physics and/or mesh module uninitialized")
249
250 ! allocate memory for data structures needed in all timedisc modules
251 ALLOCATE( &
252 this%xfluxdydz(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
253 this%yfluxdzdx(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
254 this%zfluxdxdy(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
255 this%amax(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
256 this%tol_abs(physics%VNUM), &
257 this%dtmean, this%dtstddev, &
258 this%time, &
259 stat = err)
260 IF (err.NE.0) THEN
261 CALL this%Error("InitTimedisc", "Unable to allocate memory.")
262 END IF
263
264 ! initialize state vectors
265 CALL physics%new_statevector(this%pvar,primitive)
266 CALL physics%new_statevector(this%ptmp,primitive)
267 CALL physics%new_statevector(this%cvar,conservative)
268 CALL physics%new_statevector(this%ctmp,conservative)
269 CALL physics%new_statevector(this%cold,conservative)
270 CALL physics%new_statevector(this%geo_src,conservative)
271 CALL physics%new_statevector(this%src,conservative)
272 CALL physics%new_statevector(this%rhs,conservative)
273
274 ! initialize all variables
275 this%pvar%data1d(:) = 0.
276 this%ptmp%data1d(:) = 0.
277 this%ctmp%data1d(:) = 0.
278 this%cvar%data1d(:) = 0.
279 this%cold%data1d(:) = 0.
280 this%geo_src%data1d(:) = 0.
281 this%src%data1d(:) = 0.
282 this%rhs%data1d(:) = 0.
283 this%xfluxdydz = 0.
284 this%yfluxdzdx = 0.
285 this%zfluxdxdy = 0.
286 this%amax = 0.
287 this%tol_abs = 0.
288 this%dtmean = 0.
289 this%dtstddev = 0.
290 this%break = .false.
291 this%n_adj = 0
292 this%maxerrold = 0.
293 this%dtaccept = 0
294
295 CALL getattr(config, "starttime", this%time, 0.0)
296 CALL getattr(config, "method", method)
297 CALL getattr(config, "stoptime", this%stoptime)
298 this%dt = this%stoptime
299 this%dtold = this%dt
300 this%dtmin = this%dt
301 this%dtmincause = -99
302
303 ! set default values
304 ! CFL number
305 CALL getattr(config, "cfl", this%cfl, 0.4)
306
307 ! time step minimum
308 CALL getattr(config, "dtlimit", this%dtlimit, epsilon(this%dtlimit)*this%stoptime)
309
310 ! time step maximum in units of [CFL timestep] (Used in Dumka)
311 CALL getattr(config, "dtmax", this%dtmax, 5.0)
312
313 ! maximum iterations, maxiter <= 0 means no iteration limit
314 CALL getattr(config, "maxiter", this%maxiter, 0)
315
316 ! relative tolerance for adaptive step size control
317 CALL getattr(config, "tol_rel", this%tol_rel, 0.01)
318
319 ! absolute tolerance for adaptive step size control
320 this%tol_abs(:) = 0.001
321 CALL getattr(config, "tol_abs", this%tol_abs, this%tol_abs)
322
323 ! rhs type. 0 = default, 1 = conserve angular momentum
324 CALL getattr(config, "rhstype", this%rhstype, 0)
325! IF (this%rhstype.NE.0.AND.Physics%VDIM.NE.2) THEN
326! CALL this%Error("InitTimedisc", "Alternative rhstype only works in 2D.")
327! END IF
328
329 ! time step friction parameter for PI-Controller
330 CALL getattr(config, "beta", this%beta, 0.08)
331
332 ! enable/disable update of bccsound at each substep
333 ! Usually bccsound is only needed at the beginning of each integration step
334 ! to determine the next time step. Hence it is not necessary to update
335 ! bccsound at substeps when using higher order integration schemes. Since
336 ! non-isothermal physics involve evaluation of the SQRT function to compute
337 ! the speed of sound it might speed up the simulation a bit if one disables
338 ! the update at substeps.
339 ! Diabling the update is not possible if bccsound is needed somewhere else,
340 ! e.g. some source terms depend on bccsound.
341 CALL getattr(config, "always_update_bccsound", d, 1)
342 IF (d.EQ.0) THEN
343 this%always_update_bccsound = .false.
344 ELSE
345 this%always_update_bccsound = .true.
346 END IF
347
348 ! check data bit mask
349 ! \todo{expected argument list...}
350 CALL getattr(config, "checkdata", this%checkdatabm, check_all)
351
352 ! pressure minimum to check if CHECK_PMIN is active
353 CALL getattr(config, "pmin", this%pmin, tiny(this%pmin))
354
355 ! temperature minimum to check if CHECK_TMIN is active
356 CALL getattr(config, "tmin", this%tmin, tiny(this%tmin))
357
358 ! density minimum to check if CHECK_RHOMIN is active
359 CALL getattr(config, "rhomin", this%rhomin, tiny(this%rhomin))
360
361 CALL getattr(config, "output/"//trim(physics%pvarname(1)), d, 1)
362
363 CALL this%SetOutput(mesh,physics,config,io)
364
365 ! prepare data arrays for fargo transport
366 ! check physics first; fargo is currently supported for physics_euler(isotherm) and derived types
367 SELECT TYPE(phys => physics)
368 CLASS IS(physics_eulerisotherm)
369 ! check for dimensionality of velocity vector
370 IF (mesh%fargo%GetType().GT.0) THEN
371 IF (physics%VDIM.LT.2) &
372 CALL this%Error("timedisc_base::InitTimedisc","FARGO transport with 1D velocity vector not supported!")
373 END IF
374 ! check for fargo transport direction and allocate some arrays
375 ! if direction is not in {1,2,3} skip this (fargo disabled)
376 ! see mesh_base::InitMesh for general initialization of fargo transport
377 SELECT CASE(mesh%fargo%GetDirection())
378 CASE(1) ! fargo transport along x-direction
379 ALLOCATE( &
380 this%w(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
381 this%delxy(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
382 this%shift(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
383#ifdef PARALLEL
384 this%buf(physics%VNUM+physics%PNUM,1:mesh%MININUM), & !!! NOT CHANGED, because not clear were used !
385#endif
386 stat = err)
387 CASE(2) ! fargo transport along y-direction
388 ALLOCATE( &
389 this%w(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
390 this%delxy(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
391 this%shift(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
392#ifdef PARALLEL
393 this%buf(physics%VNUM+physics%PNUM,1:mesh%MINJNUM), & !!! NOT CHANGED, because not clear were used !
394#endif
395 stat = err)
396 CASE(3) ! fargo transport along z-direction
397 ALLOCATE( &
398 this%w(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX), &
399 this%delxy(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX), &
400 this%shift(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX), &
401#ifdef PARALLEL
402 this%buf(physics%VNUM+physics%PNUM,1:mesh%MINJNUM), & !!! NOT CHANGED, because not clear where used !
403#endif
404 stat = err)
405 CASE DEFAULT
406 err = 0
407 END SELECT
408 ! additional array used for computation of background velocity
409 IF (err.EQ.0.AND.mesh%fargo%GetType().EQ.1) ALLOCATE(this%wfactor,mold=this%w,stat=err)
410 IF (err.NE.0) &
411 CALL this%Error("timedisc_base::InitTimedisc", "Unable to allocate memory for fargo advection.")
412
413 ! initialize background velocity field w
414 SELECT CASE(mesh%fargo%GetType())
415 CASE(0) ! fargo disabled -> do nothing
416 CASE(1,2) ! fargo transport with fixed user supplied or dynamical velocity
417 CALL this%Warning("timedisc_base::InitTimedisc","FARGO transport is experimental, use with care!")
418 ! fargo advection type 1: w is computed in each time step (see CalcBackgroundVelocity below)
419 ! fargo advection type 2: w is provided by the user, e.g. in InitData
420 this%w(:,:) = 0.0
421 this%shift(:,:) = 0.0
422 this%dq = marray_base()
423 this%dq%data1d(:) = 0.0
424 this%dql = marray_base()
425 this%dql%data1d(:) = 0.0
426 this%flux = marray_base()
427 this%flux%data1d(:) = 0.0
428 IF (ASSOCIATED(this%wfactor)) THEN
429 ! special geometrical factor for computation of dynamic background velocity;
430 ! this factor is the ratio of the shortest width (in fargo transport direction)
431 ! of a cell over all faces (perpendicular to the fargo transport direction) and
432 ! the extend of the cell (in fargo transport direction) through the middle of the
433 ! cell (see CalcBackgroundVelocity)
434 SELECT CASE(mesh%fargo%GetDirection())
435 CASE(1) ! fargo transport along x-direction
436 DO k=mesh%KGMIN,mesh%KGMAX
437 DO j=mesh%JGMIN,mesh%JGMAX
438 this%wfactor(j,k) = min(minval(mesh%hx%faces(mesh%IMIN,j,k,south:north)),&
439 minval(mesh%hx%faces(mesh%IMIN,j,k,bottom:top))) &
440 / mesh%hx%bcenter(mesh%IMIN,j,k)
441 END DO
442 END DO
443 CASE(2) ! fargo transport along y-direction
444 DO k=mesh%KGMIN,mesh%KGMAX
445 DO i=mesh%IGMIN,mesh%IGMAX
446 this%wfactor(i,k) = min(minval(mesh%hy%faces(i,mesh%JMIN,k,west:east)),&
447 minval(mesh%hy%faces(i,mesh%JMIN,k,bottom:top))) &
448 / mesh%hy%bcenter(i,mesh%JMIN,k)
449 END DO
450 END DO
451 CASE(3) ! fargo transport along z-direction
452 DO j=mesh%JGMIN,mesh%JGMAX
453 DO i=mesh%IGMIN,mesh%IGMAX
454 this%wfactor(i,j) = min(minval(mesh%hz%faces(i,j,mesh%KMIN,west:east)),&
455 minval(mesh%hz%faces(i,j,mesh%KMIN,south:north))) &
456 / mesh%hz%bcenter(i,j,mesh%KMIN)
457 END DO
458 END DO
459 END SELECT
460 END IF
461 CASE(3) ! fixed background velocity in shearing box
462 IF(mesh%shear_dir.EQ.2) THEN
463 this%w(:,:) = -mesh%Q*mesh%omega*mesh%bcenter(:,mesh%JMIN,:,1) !-Q*Omega*x
464 this%shift(:,:) = 0.0
465 shear_direction = "west<->east"
466 ELSE IF(mesh%shear_dir.EQ.1) THEN
467 this%w(:,:) = mesh%Q*mesh%omega*mesh%bcenter(mesh%IMIN,:,:,2) !Q*Omega*y
468 this%shift(:,:) = 0.0
469 shear_direction = "south<->north"
470 END IF
471 this%dq = marray_base()
472 this%dq%data1d(:) = 0.0
473 this%dql = marray_base()
474 this%dql%data1d(:) = 0.0
475 this%flux = marray_base()
476 this%flux%data1d(:) = 0.0
477 CASE DEFAULT
478 CALL this%Error("timedisc_base::InitTimedisc","unknown fargo advection scheme")
479 END SELECT
480 CLASS DEFAULT
481 CALL this%Error("timedisc_base::InitTimedisc","selected physics doesn't support fargo transport")
482 END SELECT
483
484 ! print some information
485 WRITE (order_str, '(I0)') this%GetOrder()
486 WRITE (cfl_str, '(F4.2)') this%GetCFL()
487 WRITE (stoptime_str, '(ES12.4)') this%stoptime
488 WRITE (dtmax_str, '(ES12.4)') this%dtmax
489 WRITE (beta_str, '(ES12.4)') this%beta
490 CALL this%Info(" TIMEDISC-> ODE solver: " //trim(this%GetName()))
491 CALL this%Info(" order: " //trim(order_str))
492 CALL this%Info(" CFL number: " //trim(cfl_str))
493 CALL this%Info(" dtmax: " //trim(dtmax_str))
494 CALL this%Info(" stoptime: " //trim(stoptime_str))
495 CALL this%Info(" beta: " //trim(beta_str))
496 ! adaptive step size control
497 IF (this%tol_rel.LT.1.0.AND.this%order.GT.1) THEN
498 ! create state vector to store the error
499 CALL physics%new_statevector(this%cerr,conservative)
500 this%cerr%data1d(:) = 0.
501 ! create selection for the internal region
502 WRITE (info_str,'(ES9.1)') this%tol_rel*100
503 CALL this%Info(" step size control: enabled")
504 CALL this%Info(" rel. precision: "//trim(info_str)//" %")
505 ELSE
506 WRITE (info_str,'(A)') "disabled"
507 END IF
508 CALL this%Info(" fargo: " //trim(mesh%fargo%GetName()))
509 IF(mesh%fargo%GetType().GT.0) &
510 CALL this%Info(" transport dir.: " //trim(mesh%fargo%GetDirectionName()))
511 SELECT CASE(this%rhstype)
512 CASE(0)
513 ! special rhs disabled, print nothing
514 CASE(1)
515 IF (.NOT.ASSOCIATED(this%w)) THEN
516 SELECT CASE(mesh%Geometry%GetAzimuthIndex())
517 CASE(2)
518 ALLOCATE(this%w(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX),stat = err)
519 CASE(3)
520 ALLOCATE(this%w(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX),stat = err)
521 END SELECT
522 IF (err.NE.0) THEN
523 CALL this%Error("InitTimedisc", "Unable to allocate memory for special RHS time stepping.")
524 END IF
525 this%w(:,:) = 0.0
526 END IF
527 CALL this%Info(" special rhs: enabled")
528 CASE DEFAULT
529 CALL this%Error("InitTimedisc","unknown rhstype")
530 END SELECT
531 END SUBROUTINE inittimedisc
532
533
534 SUBROUTINE setoutput(this,Mesh,Physics,config,IO)
535 IMPLICIT NONE
536 !-----------------------------------------------------------------------!
537 CLASS(timedisc_base), INTENT(INOUT) :: this
538 CLASS(mesh_base), INTENT(IN) :: Mesh
539 CLASS(physics_base), INTENT(IN) :: Physics
540 TYPE(dict_typ), POINTER :: config,IO
541 !------------------------------------------------------------------------!
542 INTEGER :: valwrite,i
543 CHARACTER(LEN=128) :: key,pvar_key
544 LOGICAL :: writeSolution
545 !------------------------------------------------------------------------!
546 valwrite = 0
547 CALL getattr(config, "output/error", valwrite, 0)
548 IF((valwrite.EQ.1).AND.this%tol_rel.LE.1.0) THEN
549 this%write_error = .true.
550 CALL physics%new_statevector(this%cerr_max,conservative)
551 this%cerr_max%data1d(:) = 0.
552 ELSE
553 this%write_error = .false.
554 END IF
555
556 valwrite = 0
557 CALL getattr(config, "output/solution", valwrite, 0)
558 IF(valwrite.EQ.1) THEN
559 writesolution = .true.
560 CALL physics%new_statevector(this%solution,primitive)
561 ELSE
562 NULLIFY(this%solution)
563 writesolution = .false.
564 END IF
565
566 CALL getattr(config, "output/time", valwrite, 1)
567 IF(valwrite.EQ.1) &
568 CALL setattr(io, "time", ref(this%time))
569
570 DO i=1, physics%VNUM
571 !prim
572 key = trim(physics%pvarname(i))
573 pvar_key = key
574 valwrite = 0
575 CALL getattr(config, "output/" // trim(key), valwrite, 1)
576 ! second argument is important if pvarname is used twice
577 IF (valwrite.EQ.1) THEN
578 CALL setattr(io, trim(key), &
579 this%pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
580 END IF
581
582 IF(writesolution) THEN
583 CALL setattr(io, trim(key)//"_ref", &
584 this%solution%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
585 END IF
586
587 !cons
588 key = trim(physics%cvarname(i))
589 IF(key.EQ.pvar_key) key = trim(key) // "_cvar"
590 valwrite = 0
591 CALL getattr(config, "output/" // trim(key), valwrite, 0)
592 ! second argument is important if pvarname is used twice
593 IF (valwrite.EQ.1) THEN
594 CALL setattr(io, trim(key), &
595 this%cvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
596 END IF
597
598 key = trim(physics%cvarname(i))
599 IF(this%write_error) THEN
600 CALL setattr(io, trim(key) // "_error", &
601 this%cerr_max%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
602 END IF
603
604 ! write geometrical sources
605 CALL getattr(config, "output/" // "geometrical_sources", valwrite, 0)
606 IF (valwrite.EQ.1) THEN
607 CALL setattr(io, trim(key)//"_geo_src", &
608 this%geo_src%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
609 END IF
610
611 ! write external sources
612 CALL getattr(config, "output/" // "external_sources", valwrite, 0)
613 IF (valwrite.EQ.1) THEN
614 CALL setattr(io, trim(key)//"_src", &
615 this%src%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
616 END IF
617
618 ! write right hand side
619 CALL getattr(config, "output/" // "rhs", valwrite, 0)
620 IF (valwrite.EQ.1) THEN
621 CALL setattr(io, trim(key)//"_rhs", &
622 this%rhs%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
623 END IF
624
625 ! write fluxes
626 ! ATTENTION: this are the numerical fluxes devided by dy or dx respectively
627 CALL getattr(config, "output/" // "fluxes", valwrite, 0)
628 IF (valwrite.EQ.1) THEN
629 CALL setattr(io, trim(key)//"_xfluxdydz", &
630 this%xfluxdydz(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
631 CALL setattr(io, trim(key)//"_yfluxdzdx", &
632 this%yfluxdzdx(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
633 CALL setattr(io, trim(key)//"_zfluxdxdy", &
634 this%zfluxdxdy(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
635 END IF
636
637 END DO
638 ! write bfluxes
639 CALL getattr(config, "output/" // "bflux", valwrite, 0)
640 IF(valwrite.EQ.1) THEN
641 ALLOCATE(this%bflux(physics%VNUM,6))
642 CALL setattr(io, "bflux", this%bflux)
643 ELSE
644 NULLIFY(this%bflux)
645 END IF
646 END SUBROUTINE setoutput
647
648
649 SUBROUTINE integrationstep(this,Mesh,Physics,Sources,Fluxes,iter,config,IO)
650 !------------------------------------------------------------------------!
651 CLASS(timedisc_base), INTENT(INOUT) :: this
652 CLASS(mesh_base), INTENT(INOUT) :: Mesh
653 CLASS(physics_base), INTENT(INOUT) :: Physics
654 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
655 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
656 INTEGER :: iter
657 TYPE(dict_typ), POINTER :: config,IO
658 !------------------------------------------------------------------------!
659 REAL :: err,dtold,dt,time
660 !------------------------------------------------------------------------!
661 INTENT(INOUT) :: iter
662 !------------------------------------------------------------------------!
663 ! determine the background velocity if fargo advection type 1 is enabled
664 IF (mesh%fargo%GetType().EQ.1) &
665 CALL this%CalcBackgroundVelocity(mesh,physics,this%pvar,this%cvar)
666 ! transform to comoving frame if fargo is enabled
667 SELECT CASE(mesh%fargo%GetDirection())
668 CASE(1)
669 CALL physics%SubtractBackgroundVelocityX(mesh,this%w,this%pvar,this%cvar)
670 CASE(2)
671 CALL physics%SubtractBackgroundVelocityY(mesh,this%w,this%pvar,this%cvar)
672 CASE(3)
673 CALL physics%SubtractBackgroundVelocityZ(mesh,this%w,this%pvar,this%cvar)
674 END SELECT
675
676 time = this%time
677 dt = this%dt
678 ! save dtmin if it is not determined by the fileio, i.e. time
679 ! between two successive write operations
680 IF (dt.LT.this%dtmin .AND. this%dtcause .NE. dtcause_fileio) THEN
681 this%dtmin = dt
682 this%dtmincause = this%dtcause
683 END IF
684
685 ! loop until new solution for time+dt is accepted
686 timestep: DO WHILE (.NOT.this%break.AND.time+dt.LE.this%time+this%dt)
687 dtold = dt
688
689 CALL this%SolveODE(mesh,physics,sources,fluxes,time,dt,err)
690
691 ! check truncation error and restart if necessary,
692 ! but first check data integrity of err
693 IF (err.NE.err) THEN
694 ! err = NaN
695 CALL this%Warning("timedisc_base::IntegrationStep","max error returned by SolveODE is NaN")
696 this%break = .true.
697 ELSE IF (err.GT.huge(err)) THEN
698 ! err = Inf
699 CALL this%Warning("timedisc_base::IntegrationStep","max error returned by SolveODE is Inf")
700 this%break = .true.
701 ELSE IF (err.GE.1.0) THEN
702 ! err >= 1.0
703 CALL this%RejectSolution(mesh,physics,sources,fluxes,time,dt)
704 ELSE
705 ! err < 1.0
706 CALL this%AcceptSolution(mesh,physics,sources,fluxes,time,dtold,iter)
707 END IF
708
709 ! abort integration if lower time step limit reached
710 IF (dt.LT.this%dtlimit) this%break = .true.
711
712 END DO timestep
713
714 IF (.NOT.this%break) THEN
715 ! Save true advanced time (for fargo linear advection)
716 this%dt = time - this%time
717
718 this%time = time
719 this%dtold = dt
720
721 ! perform the fargo advection step if enabled
722 IF (mesh%fargo%GetType().GT.0) THEN
723 ! ATTENTION: Boundary conditions are applied to data with subtracted background
724 ! velocity here which may yield errornous results for the real boundaries
725 ! but not for the periodic boundaries which is essential for the
726 ! fargo advection step.
727 ! Since we call ComputeRHS after the fargo step, boundary conditions
728 ! are applied again and this time with added background velocity,
729 ! at least if fargo method is set to 1 or 2 (see ComputeRHS).
730 CALL this%boundary%CenterBoundary(mesh,physics,this%time,this%pvar,this%cvar)
731 SELECT CASE(mesh%fargo%GetDirection())
732 CASE(1)
733 CALL this%FargoAdvectionX(fluxes,mesh,physics,sources)
734 CASE(2)
735 CALL this%FargoAdvectionY(fluxes,mesh,physics,sources)
736 CASE(3)
737 CALL this%FargoAdvectionZ(fluxes,mesh,physics,sources)
738 END SELECT
739 ! convert conservative to primitive variables
740 CALL physics%Convert2Primitive(this%cvar,this%pvar)
741 ! Calculate RHS after the Advection Step
742 CALL this%ComputeRHS(mesh,physics,sources,fluxes,this%time,0.,this%pvar,this%cvar,&
743 this%checkdatabm,this%rhs)
744 ELSE
745 CALL this%ComputeRHS(mesh,physics,sources,fluxes,this%time,dtold,this%pvar,this%cvar,&
746 this%checkdatabm,this%rhs)
747 END IF
748 this%cold = this%cvar
749 END IF
750 END SUBROUTINE integrationstep
751
752
759 FUNCTION adjusttimestep(this,maxerr,dtold) RESULT(dtnew)
760 IMPLICIT NONE
761 !------------------------------------------------------------------------!
762 CLASS(timedisc_base) :: this
763 REAL :: maxerr,dtold,dtnew
764 !------------------------------------------------------------------------!
765 REAL :: dttmp
766 !------------------------------------------------------------------------!
767 INTENT(IN) :: maxerr,dtold
768 INTENT(INOUT) :: this
769 !------------------------------------------------------------------------!
770 ! adaptive step size control disabled
771 IF (this%tol_rel.GE.1.0) THEN
772 dtnew = huge(dtnew)
773 RETURN
774 END IF
775
776 ! time step estimate with maxerr determined from the comparison of numerical
777 ! results of two explicit schemes with different order where GetOrder(this)
778 ! returns the order of the higher order scheme (P-Controller)
779
780 IF(maxerr.GT.0.) THEN
781 dttmp = 0.9*dtold*exp(-log(maxerr)/getorder(this))
782 ELSE
783 ! ths limit for maxerr->0 is dttmp -> oo. But this cut off to
784 ! 4*dtold later in this function.
785 dttmp = 4.*dtold
786 END IF
787
788 ! this controls the increase/decrease of time steps based on the
789 ! old value of maxerr from the previous time step (PI-Controller)
790 IF (this%maxerrold.GT.0.) THEN
791 dtnew = dttmp * exp(-this%beta*log(this%maxerrold/maxerr))
792 ELSE
793 dtnew = dttmp
794 END IF
795
796 ! adjust the time step based on the error estimate
797 IF (maxerr.LT.1.0) THEN
798 ! increase time step
799 dtnew = min(dtnew,4.*dtold) ! enlarge at most by a factor of 4
800 ! it is possible: maxerr < 1 => 0.9*dt < dtnew => maybe dtnew < dt
801 IF (dtnew .LT. dtold) this%dtcause = dtcause_smallerr
802 ELSE
803 ! If maxerrold >> maxerr dtnew can become larger than dtold even
804 ! if the timestep is rejected. If this is the case fall back to
805 ! the simple P-Controller.
806 IF (dtnew.GT.dtold) dtnew = dttmp
807 ! decrease time step
808 dtnew = max(dtnew,0.25*dtold) ! reduce at most by a factor of 1/4
809 END IF
810
811 ! store maxerr for next time step
812 this%maxerrold = maxerr
813
814 END FUNCTION adjusttimestep
815
816
818 REAL function calctimestep(this,mesh,physics,sources,fluxes,time,dtcause) result(dt)
819 USE physics_euler_mod, ONLY : physics_euler, statevector_euler
820 IMPLICIT NONE
821 !------------------------------------------------------------------------!
822 CLASS(timedisc_base), INTENT(INOUT) :: this
823 CLASS(mesh_base), INTENT(IN) :: mesh
824 CLASS(physics_base), INTENT(INOUT) :: physics
825 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: sources
826 CLASS(fluxes_base), INTENT(INOUT) :: fluxes
827 REAL, INTENT(IN) :: time
828 INTEGER, INTENT(INOUT) :: dtcause
829 !------------------------------------------------------------------------!
830 REAL :: invdt
831 REAL :: dt_cfl, dt_src
832 REAL :: invdt_x,invdt_y,invdt_z
833 !------------------------------------------------------------------------!
834 ! CFL condition:
835 ! maximal wave speeds in each direction
836 IF (.NOT.this%always_update_bccsound) THEN
837 SELECT TYPE(phys => physics)
838 CLASS IS(physics_euler)
839 SELECT TYPE(pvar => this%pvar)
840 CLASS IS(statevector_euler)
841 CALL phys%UpdateSoundSpeed(pvar)
842 END SELECT
843 END SELECT
844 END IF
845
846 SELECT TYPE(phys => physics)
847 TYPE IS(physics_eulerisotherm)
848 SELECT TYPE(pvar => this%pvar)
849 CLASS IS(statevector_eulerisotherm)
850 CALL phys%CalculateWaveSpeeds(mesh,pvar,fluxes%minwav,fluxes%maxwav)
851 END SELECT
852 TYPE IS(physics_euler)
853 SELECT TYPE(pvar => this%pvar)
854 CLASS IS(statevector_euler)
855 CALL phys%CalculateWaveSpeeds(mesh,pvar,fluxes%minwav,fluxes%maxwav)
856 END SELECT
857 END SELECT
858
859 ! compute maximum of inverse time for CFL condition
860 IF ((mesh%JNUM.EQ.1).AND.(mesh%KNUM.EQ.1)) THEN
861 ! 1D, only x-direction
862 invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlx%data1d(:))
863 ELSE IF ((mesh%INUM.EQ.1).AND.(mesh%KNUM.EQ.1)) THEN
864 ! 1D, only y-direction
865 invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dly%data1d(:))
866 ELSE IF ((mesh%INUM.EQ.1).AND.(mesh%JNUM.EQ.1)) THEN
867 ! 1D, only z-direction
868 invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlz%data1d(:))
869 ELSE IF ((mesh%INUM.GT.1).AND.(mesh%JNUM.GT.1).AND.(mesh%KNUM.EQ.1)) THEN
870 ! 2D, x-y-plane
871 invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlx%data1d(:) &
872 + max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) / mesh%dly%data1d(:))
873 ELSE IF ((mesh%INUM.GT.1).AND.(mesh%KNUM.GT.1).AND.(mesh%JNUM.EQ.1)) THEN
874 ! 2D, x-z-plane
875 invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlx%data1d(:) &
876 + max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) / mesh%dlz%data1d(:))
877 ELSE IF ((mesh%JNUM.GT.1).AND.(mesh%KNUM.GT.1).AND.(mesh%INUM.EQ.1)) THEN
878 ! 2D, y-z-plane
879 invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dly%data1d(:) &
880 + max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) / mesh%dlz%data1d(:))
881 ELSE
882 ! full 3D
883 !TODO: Achtung: Hier wurde fuer eine bessere Symmetrie fuer jede Richtung ein eigenes invdt
884 ! berechnet. Dies koennte jedoch einen Verlust an Stabilitaet bewirken. Hier muesste mal eine
885 ! Stabilitaetsanalyse gemacht werden
886 invdt_x = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlx%data1d(:))
887 invdt_y = maxval(max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) / mesh%dly%data1d(:))
888 invdt_z = maxval(max(fluxes%maxwav%data2d(:,3),-fluxes%minwav%data2d(:,3)) / mesh%dlz%data1d(:))
889 invdt = max(invdt_y,invdt_z,invdt_x)
890 ! invdt = MAXVAL(MAX(Fluxes%maxwav(:,:,:,3),-Fluxes%minwav(:,:,:,3)) / Mesh%dlz(:,:,:) &
891 ! + MAX(Fluxes%maxwav(:,:,:,2),-Fluxes%minwav(:,:,:,2)) / Mesh%dly(:,:,:) &
892 ! + MAX(Fluxes%maxwav(:,:,:,1),-Fluxes%minwav(:,:,:,1)) / Mesh%dlx(:,:,:))
893 END IF
894
895 ! largest time step due to CFL condition
896 dt_cfl = this%cfl / invdt
897 ! time step limitation due to cfl -> dtcause = 0
898 dtcause = 0
899
900 ! check for sources
901 IF (ALLOCATED(sources)) THEN
902 ! initialize this to be sure dt_src > 0
903 dt_src = dt_cfl
904 CALL sources%CalcTimestep(mesh,physics,fluxes,this%pvar,this%cvar,time,dt_src,dtcause)
905 dt = min(dt_cfl,dt_src)
906 ELSE
907 dt = dt_cfl
908 END IF
909 END FUNCTION calctimestep
910
911 SUBROUTINE acceptsolution(this,Mesh,Physics,Sources,Fluxes,time,dt,iter)
912 IMPLICIT NONE
913 !------------------------------------------------------------------------!
914 CLASS(timedisc_base), INTENT(INOUT) :: this
915 CLASS(mesh_base), INTENT(IN) :: Mesh
916 CLASS(physics_base), INTENT(INOUT) :: Physics
917 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
918 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
919 REAL :: time,dt
920 INTEGER :: iter
921 !------------------------------------------------------------------------!
922 REAL :: dtmeanold
923 REAL, DIMENSION(:), POINTER :: bflux,bfold
924 !------------------------------------------------------------------------!
925 INTENT(INOUT) :: time,dt
926 !------------------------------------------------------------------------!
927 time=time+dt
928 this%dtaccept = this%dtaccept + 1
929 dtmeanold = this%dtmean
930 this%dtmean = this%dtmean + (dt - this%dtmean)/this%dtaccept
931 this%dtstddev = this%dtstddev + (dt - dtmeanold)*(dt-this%dtmean)
932 IF (time.LT.this%time+this%dt) THEN
933 CALL this%ComputeRHS(mesh,physics,sources,fluxes,time,dt,this%pvar,this%cvar,&
934 this%checkdatabm,this%rhs)
935 this%cold = this%cvar
936 END IF
937 IF (mesh%INUM.GT.1) THEN
938 ! collapse the 4D arrays to improve vector performance;
939 ! maybe we can switch to the old style assignments if
940 ! automatic collapsing of NEC nfort compiler works
941 bflux(1:SIZE(fluxes%bxflux)) => fluxes%bxflux
942 bfold(1:SIZE(fluxes%bxfold)) => fluxes%bxfold
943 bfold(:) = bflux(:)
944! Fluxes%bxfold(:,:,:,:) = Fluxes%bxflux(:,:,:,:)
945 END IF
946 IF (mesh%JNUM.GT.1) THEN
947 bflux(1:SIZE(fluxes%byflux)) => fluxes%byflux
948 bfold(1:SIZE(fluxes%byfold)) => fluxes%byfold
949 bfold(:) = bflux(:)
950! Fluxes%byfold(:,:,:,:) = Fluxes%byflux(:,:,:,:)
951 END IF
952 IF (mesh%KNUM.GT.1) THEN
953 bflux(1:SIZE(fluxes%bzflux)) => fluxes%bzflux
954 bfold(1:SIZE(fluxes%bzfold)) => fluxes%bzfold
955 bfold(:) = bflux(:)
956! Fluxes%bzfold(:,:,:,:) = Fluxes%bzflux(:,:,:,:)
957 END IF
958 iter = iter + 1
959 END SUBROUTINE acceptsolution
960
961 SUBROUTINE rejectsolution(this,Mesh,Physics,Sources,Fluxes,time,dt)
962 IMPLICIT NONE
963 !------------------------------------------------------------------------!
964 CLASS(timedisc_base), INTENT(INOUT) :: this
965 CLASS(mesh_base), INTENT(IN) :: Mesh
966 CLASS(physics_base), INTENT(INOUT) :: Physics
967 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
968 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
969 REAL :: time,dt
970 !------------------------------------------------------------------------!
971 INTEGER :: i,j,k,l,m
972 !------------------------------------------------------------------------!
973 INTENT(IN) :: time
974 INTENT(INOUT) :: dt
975 !------------------------------------------------------------------------!
976 this%cvar = this%cold
977 ! This data has already been checked before in AcceptSolution
978 CALL this%ComputeRHS(mesh,physics,sources,fluxes,time,dt,this%pvar, &
979 this%cvar,check_nothing,this%rhs)
980!NEC$ SHORTLOOP
981 DO m=1,physics%VNUM
982!NEC$ UNROLL(2)
983 DO l=1,2
984 DO k=mesh%KGMIN,mesh%KGMAX
985 DO j=mesh%JGMIN,mesh%JGMAX
986 fluxes%bxflux(j,k,l,m) = fluxes%bxfold(j,k,l,m)
987 END DO
988 END DO
989 DO i=mesh%IGMIN,mesh%IGMAX
990 DO k=mesh%KGMIN,mesh%KGMAX
991 fluxes%byflux(k,i,l,m) = fluxes%byfold(k,i,l,m)
992 END DO
993 END DO
994 DO j=mesh%JGMIN,mesh%JGMAX
995 DO i=mesh%IGMIN,mesh%IGMAX
996 fluxes%bzflux(i,j,l,m) = fluxes%bzfold(i,j,l,m)
997 END DO
998 END DO
999 END DO
1000 END DO
1001 ! count adjustments for information
1002 this%n_adj = this%n_adj + 1
1003 this%dtcause = dtcause_erradj
1004 ! only save dtmin if the reason is not the fileio (automatically satisfied)
1005 IF (dt.LT.this%dtmin) THEN
1006 this%dtmin = dt
1007 this%dtmincause = this%dtcause
1008 END IF
1009 END SUBROUTINE rejectsolution
1010
1011
1012 FUNCTION computeerror(this,Mesh,Physics,cvar_high,cvar_low) RESULT(maxerr)
1013 IMPLICIT NONE
1014 !------------------------------------------------------------------------!
1015 CLASS(timedisc_base), INTENT(INOUT) :: this
1016 CLASS(mesh_base), INTENT(IN) :: mesh
1017 CLASS(physics_base), INTENT(IN) :: physics
1018 CLASS(marray_compound),INTENT(INOUT):: cvar_high,cvar_low
1019 REAL :: maxerr
1020 !------------------------------------------------------------------------!
1021 INTEGER :: i,l
1022#ifdef PARALLEL
1023 INTEGER :: ierror
1024#endif
1025 REAL :: rel_err(physics%vnum)
1026 !------------------------------------------------------------------------!
1027!NEC$ SHORTLOOP
1028 DO l=1,physics%VNUM
1029 DO concurrent(i=1:SIZE(this%cerr%data2d,dim=1))
1030 ! compute the local error (including ghost zones)
1031 this%cerr%data2d(i,l) = abs(cvar_high%data2d(i,l)-cvar_low%data2d(i,l)) &
1032 / (this%tol_rel*abs(cvar_high%data2d(i,l)) + this%tol_abs(l))
1033 END DO
1034 ! determine the global maximum on the whole grid except for ghost zones
1035 rel_err(l) = maxval(this%cerr%data2d(:,l),mask=mesh%without_ghost_zones%mask1d(:))
1036 ! store the maximum between two output time steps
1037 IF (this%write_error) THEN
1038 DO concurrent(i=1:SIZE(this%cerr_max%data2d,dim=1))
1039 this%cerr_max%data2d(i,l) = max(this%cerr_max%data2d(i,l),this%cerr%data2d(i,l))
1040 END DO
1041 END IF
1042 END DO
1043
1044 ! compute the maximum of all variables
1045 maxerr = maxval(rel_err(:))
1046
1047#ifdef PARALLEL
1048 ! find the global maximum of all processes
1049 CALL mpi_allreduce(mpi_in_place,maxerr,1,default_mpi_real,mpi_max,&
1050 mesh%comm_cart,ierror)
1051#endif
1052
1053 END FUNCTION computeerror
1054
1055
1056
1073 SUBROUTINE computerhs(this,Mesh,Physics,Sources,Fluxes,time,dt,pvar,cvar,checkdatabm,rhs)
1074 USE physics_euler_mod, ONLY : physics_euler,statevector_euler
1075 USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
1076 USE geometry_spherical_mod, ONLY: geometry_spherical
1077 USE geometry_cylindrical_mod, ONLY: geometry_cylindrical
1078 IMPLICIT NONE
1079 !------------------------------------------------------------------------!
1080 CLASS(timedisc_base), INTENT(INOUT) :: this
1081 CLASS(mesh_base), INTENT(IN) :: Mesh
1082 CLASS(physics_base), INTENT(INOUT) :: Physics
1083 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
1084 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
1085 REAL, INTENT(IN) :: time, dt
1086 CLASS(marray_compound),INTENT(INOUT):: pvar,cvar,rhs
1087 INTEGER, INTENT(IN) :: checkdatabm
1088 !------------------------------------------------------------------------!
1089 INTEGER :: i,j,k,l
1090 LOGICAL :: have_potential
1091 REAL :: t
1092 REAL :: phi,wp
1093 REAL, DIMENSION(:,:,:,:), POINTER :: pot
1094 CLASS(sources_base), POINTER :: sp => null()
1095 CLASS(sources_gravity), POINTER :: grav => null()
1096 !------------------------------------------------------------------------!
1097 t = time
1098 SELECT CASE(mesh%fargo%GetType())
1099 CASE(1,2)
1100 ! transform to real velocity, i.e. v_residual + w_background,
1101 ! before setting the boundary conditions
1102 SELECT CASE(mesh%fargo%GetDirection())
1103 CASE(1)
1104 CALL physics%AddBackgroundVelocityX(mesh,this%w,pvar,cvar)
1105 CASE(2)
1106 CALL physics%AddBackgroundVelocityY(mesh,this%w,pvar,cvar)
1107 CASE(3)
1108 CALL physics%AddBackgroundVelocityZ(mesh,this%w,pvar,cvar)
1109 END SELECT
1110 CASE(3)
1111 ! boundary conditions are set for residual velocity in shearing box simulations
1112 ! usually it's not necessary to subtract the background velocity here, but
1113 ! in case someone adds it before, we subtract it here; the subroutine checks,
1114 ! if the velocities have already been transformed
1115 SELECT CASE(mesh%fargo%GetDirection())
1116 CASE(1)
1117 CALL physics%SubtractBackgroundVelocityX(mesh,this%w,pvar,cvar)
1118 CASE(2)
1119 CALL physics%SubtractBackgroundVelocityY(mesh,this%w,pvar,cvar)
1120 CASE(3)
1121 CALL physics%SubtractBackgroundVelocityZ(mesh,this%w,pvar,cvar)
1122 END SELECT
1123 ! ATTENTION: the time must be the initial time of the whole time step
1124 ! not the time of a substep if fargo type 3 (shearing box)
1125 ! is enabled
1126 t = this%time
1127 CASE DEFAULT
1128 ! fargo disabled (do nothing)
1129 END SELECT
1130
1131 ! set boundary values and convert conservative to primitive variables
1132 ! ATTENTION: if fargo type 3 is enabled, background velocity is subtracted
1133 ! afterwards, otherwise not
1134 CALL this%boundary%CenterBoundary(mesh,physics,t,pvar,cvar)
1135
1136 IF(iand(checkdatabm,check_tmin).NE.check_nothing.AND.&
1137 this%tmin.GT.1.e-10) THEN
1138 SELECT TYPE(p => pvar)
1139 CLASS IS(statevector_euler)
1140 SELECT TYPE(c => cvar)
1141 CLASS IS(statevector_euler)
1142 ! If temperature is below TMIN limit pressure to density*RG/MU*TMIN.
1143 p%pressure%data1d(:) = max(p%pressure%data1d(:), &
1144 p%density%data1d(:)*physics%Constants%RG/physics%MU*this%TMIN)
1145 CALL physics%Convert2Conservative(p,c)
1146 END SELECT
1147 END SELECT
1148 END IF
1149
1150 IF(iand(checkdatabm,check_pmin).NE.check_nothing.AND.&
1151 physics%PRESSURE.GT.0) THEN
1152 ! Check if the pressure is below pmin. If it is, increase the pressure
1153 ! to reach pmin
1154 SELECT TYPE(p => pvar)
1155 CLASS IS(statevector_euler)
1156 SELECT TYPE(c => cvar)
1157 CLASS IS(statevector_euler)
1158 ! If temperature is below PMIN limit pressure to PMIN
1159 p%pressure%data1d(:) &
1160 = max(p%pressure%data1d(:),this%pmin)
1161 CALL physics%Convert2Conservative(p,c)
1162 END SELECT
1163 END SELECT
1164 END IF
1165
1166 ! update the speed of sound (non-isotherml physics only)
1167 IF (this%always_update_bccsound) THEN
1168 SELECT TYPE(phys => physics)
1169 CLASS IS(physics_euler)
1170 SELECT TYPE(pvar => this%pvar)
1171 CLASS IS(statevector_euler)
1172 CALL phys%UpdateSoundSpeed(pvar)
1173 END SELECT
1174 END SELECT
1175 END IF
1176
1177 ! check for illegal data
1178 IF(checkdatabm.NE.check_nothing) &
1179 CALL this%CheckData(mesh,physics,fluxes,pvar,cvar,checkdatabm)
1180
1181 ! get geometrical sources
1182 CALL physics%GeometricalSources(mesh,pvar,cvar,this%geo_src)
1183
1184 ! get source terms due to external forces if present
1185 IF (ALLOCATED(sources)) &
1186 CALL sources%ExternalSources(mesh,physics,fluxes,sources, &
1187 t,dt,pvar,cvar,this%src)
1188
1189 ! if fargo advection is enabled additional source terms occur;
1190 ! furthermore computation of numerical fluxes should always be
1191 ! carried out with residual velocity
1192 SELECT CASE(mesh%fargo%GetType())
1193 CASE(1,2)
1194 ! add fargo source terms to geometrical source terms and
1195 ! subtract background velocity
1196 ! zero geometry sources for cartesian simulations first
1197 IF (mesh%geometry%GetType().EQ.cartesian) this%geo_src%data1d(:) = 0.0
1198 SELECT CASE(mesh%fargo%GetDirection())
1199 CASE(1)
1200 CALL physics%AddFargoSourcesX(mesh,this%w,pvar,cvar,this%geo_src)
1201 CALL physics%SubtractBackgroundVelocityX(mesh,this%w,pvar,cvar)
1202 CASE(2)
1203 CALL physics%AddFargoSourcesY(mesh,this%w,pvar,cvar,this%geo_src)
1204 CALL physics%SubtractBackgroundVelocityY(mesh,this%w,pvar,cvar)
1205 CASE(3)
1206 CALL physics%AddFargoSourcesZ(mesh,this%w,pvar,cvar,this%geo_src)
1207 CALL physics%SubtractBackgroundVelocityZ(mesh,this%w,pvar,cvar)
1208 END SELECT
1209 CASE(3)
1210 ! background velocity field has already been subtracted (do nothing);
1211 ! fargo specific source terms are handled in the shearing box source
1212 ! term module
1213 CASE DEFAULT
1214 ! fargo disabled (do nothing)
1215 END SELECT
1216
1217 ! get the numerical fluxes
1218 CALL fluxes%CalculateFluxes(mesh,physics,pvar,cvar,this%xfluxdydz, &
1219 this%yfluxdzdx,this%zfluxdxdy)
1220
1221 ! compute the right hand side for boundary flux computation;
1222 ! this is probably wrong for the special rhs (rhstype=1, see above)
1223!NEC$ SHORTLOOP
1224 rhsbdflux: DO l=1,physics%VNUM+physics%PNUM
1225!NEC$ IVDEP
1226 DO k=mesh%KMIN,mesh%KMAX
1227!NEC$ IVDEP
1228 DO j=mesh%JMIN,mesh%JMAX
1229 ! western and eastern boundary fluxes
1230 rhs%data4d(mesh%IMIN-mesh%Ip1,j,k,l) = mesh%dy*mesh%dz * this%xfluxdydz(mesh%IMIN-mesh%Ip1,j,k,l)
1231 rhs%data4d(mesh%IMAX+mesh%Ip1,j,k,l) = -mesh%dy*mesh%dz * this%xfluxdydz(mesh%IMAX,j,k,l)
1232 END DO
1233 END DO
1234!NEC$ IVDEP
1235 DO k=mesh%KMIN,mesh%KMAX
1236!NEC$ IVDEP
1237 DO i=mesh%IMIN,mesh%IMAX
1238 ! southern and northern boundary fluxes
1239 rhs%data4d(i,mesh%JMIN-mesh%Jp1,k,l) = mesh%dz*mesh%dx * this%yfluxdzdx(i,mesh%JMIN-mesh%Jp1,k,l)
1240 rhs%data4d(i,mesh%JMAX+mesh%Jp1,k,l) = -mesh%dz*mesh%dx * this%yfluxdzdx(i,mesh%JMAX,k,l)
1241 END DO
1242 END DO
1243!NEC$ IVDEP
1244 DO j=mesh%JMIN,mesh%JMAX
1245!NEC$ IVDEP
1246 DO i=mesh%IMIN,mesh%IMAX
1247 ! bottomer and upper boundary fluxes
1248 rhs%data4d(i,j,mesh%KMIN-mesh%Kp1,l) = mesh%dx*mesh%dy * this%zfluxdxdy(i,j,mesh%KMIN-mesh%Kp1,l)
1249 rhs%data4d(i,j,mesh%KMAX+mesh%Kp1,l) = -mesh%dx*mesh%dy * this%zfluxdxdy(i,j,mesh%KMAX,l)
1250 END DO
1251 END DO
1252 END DO rhsbdflux
1253
1254 SELECT CASE(this%rhstype)
1255 CASE(0)
1256!NEC$ SHORTLOOP
1257 DO l=1,physics%VNUM+physics%PNUM
1258 DO k=mesh%KMIN,mesh%KMAX
1259 DO j=mesh%JMIN,mesh%JMAX
1260!NEC$ IVDEP
1261 DO i=mesh%IMIN,mesh%IMAX
1262 ! update right hand side of ODE
1263 rhs%data4d(i,j,k,l) = &
1264 mesh%dydzdV%data3d(i,j,k)*(this%xfluxdydz(i,j,k,l) - this%xfluxdydz(i-mesh%Ip1,j,k,l)) &
1265 + mesh%dzdxdV%data3d(i,j,k)*(this%yfluxdzdx(i,j,k,l) - this%yfluxdzdx(i,j-mesh%Jp1,k,l)) &
1266 + mesh%dxdydV%data3d(i,j,k)*(this%zfluxdxdy(i,j,k,l) - this%zfluxdxdy(i,j,k-mesh%Kp1,l)) &
1267 - this%geo_src%data4d(i,j,k,l) - this%src%data4d(i,j,k,l)
1268 END DO
1269 END DO
1270 END DO
1271 END DO
1272
1273 CASE(1)
1277 IF (ALLOCATED(sources)) &
1278 sp => sources%GetSourcesPointer(gravity)
1279
1280 NULLIFY(pot)
1281 have_potential = .false.
1282 IF(ASSOCIATED(grav)) THEN
1283 IF(.NOT.grav%addtoenergy) THEN
1284 pot => mesh%RemapBounds(grav%pot%data4d(:,:,:,:))
1285 have_potential=.true.
1286 END IF
1287 END IF
1288
1289 phi = 0.
1290
1291 SELECT CASE(mesh%geometry%GetAzimuthIndex())
1292 CASE(2) ! Cylindrical
1293 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX
1294 DO j=mesh%JMIN-mesh%JP1,mesh%JMAX
1295!NEC$ IVDEP
1296 DO i=mesh%IMIN-mesh%IP1,mesh%IMAX
1297 wp = 0.5*(this%w(i,k)+this%w(i+1,k)) + mesh%hy%faces(i+1,j,k,1)*mesh%OMEGA
1298 IF(physics%PRESSURE.GT.0) THEN
1299 this%xfluxdydz(i,j,k,physics%ENERGY) &
1300 = this%xfluxdydz(i,j,k,physics%ENERGY) &
1301 + wp * (0.5 * wp * this%xfluxdydz(i,j,k,physics%DENSITY) &
1302 + this%xfluxdydz(i,j,k,physics%YMOMENTUM))
1303 IF (have_potential) THEN
1304 this%xfluxdydz(i,j,k,physics%ENERGY) = this%xfluxdydz(i,j,k,physics%ENERGY) &
1305 + pot(i,j,k,2)*this%xfluxdydz(i,j,k,physics%DENSITY)
1306 END IF
1307 END IF
1308
1309 this%xfluxdydz(i,j,k,physics%YMOMENTUM) &
1310 = (this%xfluxdydz(i,j,k,physics%YMOMENTUM) &
1311 + wp * this%xfluxdydz(i,j,k,physics%DENSITY)) * mesh%hy%faces(i+mesh%IP1,j,k,1)
1312
1313 wp = this%w(i,k) + mesh%hy%bcenter(i,j,k)*mesh%OMEGA
1314
1315 IF(physics%PRESSURE.GT.0) THEN
1316 this%yfluxdzdx(i,j,k,physics%ENERGY) &
1317 = this%yfluxdzdx(i,j,k,physics%ENERGY) &
1318 + wp * ( 0.5 * wp * this%yfluxdzdx(i,j,k,physics%DENSITY) &
1319 + this%yfluxdzdx(i,j,k,physics%YMOMENTUM))
1320 IF (have_potential) THEN
1321 this%yfluxdzdx(i,j,k,physics%ENERGY) = this%yfluxdzdx(i,j,k,physics%ENERGY) &
1322 + pot(i,j,k,3)*this%yfluxdzdx(i,j,k,physics%DENSITY)
1323 END IF
1324 END IF
1325
1326 this%yfluxdzdx(i,j,k,physics%YMOMENTUM) &
1327 = this%yfluxdzdx(i,j,k,physics%YMOMENTUM) &
1328 + wp * this%yfluxdzdx(i,j,k,physics%DENSITY)
1329
1330 wp = this%w(i,k) + mesh%hy%bcenter(i,j,k)*mesh%OMEGA
1331
1332 IF(physics%PRESSURE.GT.0) THEN
1333 this%zfluxdxdy(i,j,k,physics%ENERGY) &
1334 = this%zfluxdxdy(i,j,k,physics%ENERGY) &
1335 + wp * ( 0.5 * wp * this%zfluxdxdy(i,j,k,physics%DENSITY) &
1336 + this%zfluxdxdy(i,j,k,physics%YMOMENTUM))
1337 IF (have_potential) THEN
1338 this%zfluxdxdy(i,j,k,physics%ENERGY) = this%zfluxdxdy(i,j,k,physics%ENERGY) &
1339 + pot(i,j,k,4)*this%zfluxdxdy(i,j,k,physics%DENSITY)
1340 END IF
1341 END IF
1342
1343 this%zfluxdxdy(i,j,k,physics%YMOMENTUM) &
1344 = (this%zfluxdxdy(i,j,k,physics%YMOMENTUM) &
1345 + wp * this%zfluxdxdy(i,j,k,physics%DENSITY)) * mesh%hy%faces(i,j,k+mesh%KP1,1)
1346
1347 END DO
1348 END DO
1349 END DO
1350
1351 ! set isothermal sound speeds
1352 SELECT TYPE (phys => physics)
1353 CLASS IS (physics_eulerisotherm)
1354 DO k=mesh%KMIN,mesh%KMAX
1355 DO j=mesh%JMIN,mesh%JMAX
1356!NEC$ IVDEP
1357 DO i=mesh%IMIN,mesh%IMAX
1358 wp = pvar%data4d(i,j,k,physics%YVELOCITY) + this%w(i,k) + mesh%hy%bcenter(i,j,k)*mesh%OMEGA
1359
1360 this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1361 = pvar%data4d(i,j,k,physics%DENSITY) * wp * (wp*mesh%cyxy%bcenter(i,j,k) - &
1362 pvar%data4d(i,j,k,physics%XVELOCITY) * mesh%cxyx%bcenter(i,j,k))
1363
1364 !cyxy and cyzy are eleminated due to the modified divergence operator
1365 this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1366 = cvar%data4d(i,j,k,physics%XMOMENTUM) * ( pvar%data4d(i,j,k,physics%XVELOCITY) * mesh%cxyx%center(i,j,k))
1367
1368 IF(physics%ZMOMENTUM.GT.0) THEN
1369 this%geo_src%data4d(i,j,k,physics%XMOMENTUM) = this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1370 + cvar%data4d(i,j,k,physics%ZVELOCITY) * (pvar%data4d(i,j,k,physics%ZVELOCITY)*mesh%czxz%bcenter(i,j,k) - &
1371 pvar%data4d(i,j,k,physics%XVELOCITY) * mesh%cxzx%bcenter(i,j,k))
1372 this%geo_src%data4d(i,j,k,physics%YMOMENTUM) = this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1373 + cvar%data4d(i,j,k,physics%ZMOMENTUM) * ( pvar%data4d(i,j,k,physics%ZVELOCITY) * mesh%czyz%center(i,j,k))
1374 this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) = cvar%data4d(i,j,k,physics%XMOMENTUM) &
1375 * ( pvar%data4d(i,j,k,physics%XVELOCITY) * mesh%cxzx%center(i,j,k) &
1376 - pvar%data4d(i,j,k,physics%ZVELOCITY) * mesh%czxz%center(i,j,k) ) &
1377 + pvar%data4d(i,j,k,physics%DENSITY)*wp* ( wp*mesh%cyzy%center(i,j,k) &
1378 - pvar%data4d(i,j,k,physics%ZVELOCITY) * mesh%czyz%center(i,j,k) )
1379 END IF
1380
1381 IF(physics%PRESSURE.GT.0) THEN
1382 this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1383 = this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1384 +pvar%data4d(i,j,k,physics%PRESSURE) &
1385 *( mesh%cyxy%center(i,j,k) + mesh%czxz%center(i,j,k) )
1386 this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1387 = this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1388 + pvar%data4d(i,j,k,physics%PRESSURE) &
1389 *( mesh%cxyx%center(i,j,k) + mesh%czyz%center(i,j,k) )
1390 IF(physics%ZMOMENTUM.GT.0) &
1391 this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) &
1392 = this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) &
1393 + pvar%data4d(i,j,k,physics%PRESSURE) &
1394 *( mesh%cxzx%center(i,j,k) + mesh%cyzy%center(i,j,k) )
1395
1396 this%geo_src%data4d(i,j,k,physics%ENERGY) = 0.
1397 ELSE
1398
1399 this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1400 = this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1401 +pvar%data4d(i,j,k,physics%DENSITY)*phys%bccsound%data3d(i,j,k)**2 &
1402 * ( mesh%cyxy%center(i,j,k) + mesh%czxz%center(i,j,k) )
1403 this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1404 = this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1405 + pvar%data4d(i,j,k,physics%DENSITY)*phys%bccsound%data3d(i,j,k)**2 &
1406 *( mesh%cxyx%center(i,j,k) + mesh%czyz%center(i,j,k) )
1407 IF(physics%ZMOMENTUM.GT.0) &
1408 this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) &
1409 = this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) &
1410 + pvar%data4d(i,j,k,physics%DENSITY)*phys%bccsound%data3d(i,j,k)**2 &
1411 *( mesh%cxzx%center(i,j,k) + mesh%cyzy%center(i,j,k) )
1412
1413 END IF
1414 END DO
1415 END DO
1416 END DO
1417 CLASS DEFAULT
1418 CALL this%Error("timedisc_base::ComputeRHS","physics not supported with rhstype=1")
1419 END SELECT
1420
1421
1422!NEC$ NOVECTOR
1423 DO l=1,physics%VNUM+physics%PNUM
1424!NEC$ OUTERLOOP_UNROLL(8)
1425 DO k=mesh%KMIN,mesh%KMAX
1426 DO j=mesh%JMIN,mesh%JMAX
1427!NEC$ IVDEP
1428 DO i=mesh%IMIN,mesh%IMAX
1429 ! update right hand side of ODE
1430 IF(l.EQ.physics%YMOMENTUM) THEN
1431 rhs%data4d(i,j,k,l) = mesh%dydzdV%data3d(i,j,k)*(this%xfluxdydz(i,j,k,l) - this%xfluxdydz(i-mesh%IP1,j,k,l)) &
1432 / mesh%hy%center(i,j,k) &
1433 + mesh%dzdxdV%data3d(i,j,k)*(this%yfluxdzdx(i,j,k,l) - this%yfluxdzdx(i,j-mesh%JP1,k,l)) &
1434 + mesh%dxdydV%data3d(i,j,k)*(this%zfluxdxdy(i,j,k,l) - this%zfluxdxdy(i,j,k-mesh%KP1,l)) &
1435 / mesh%hy%center(i,j,k)
1436 ELSE
1437 rhs%data4d(i,j,k,l) = mesh%dydzdV%data3d(i,j,k)*(this%xfluxdydz(i,j,k,l) - this%xfluxdydz(i-mesh%IP1,j,k,l)) &
1438 + mesh%dzdxdV%data3d(i,j,k)*(this%yfluxdzdx(i,j,k,l) - this%yfluxdzdx(i,j-mesh%JP1,k,l)) &
1439 + mesh%dxdydV%data3d(i,j,k)*(this%zfluxdxdy(i,j,k,l) - this%zfluxdxdy(i,j,k-mesh%KP1,l))
1440 END IF
1441 END DO
1442 END DO
1443 END DO
1444 END DO
1445
1446 DO k=mesh%KMIN,mesh%KMAX
1447 DO j=mesh%JMIN,mesh%JMAX
1448!NEC$ IVDEP
1449 DO i=mesh%IMIN,mesh%IMAX
1450 wp = this%w(i,k) + mesh%radius%bcenter(i,j,k) * mesh%OMEGA
1451 rhs%data4d(i,j,k,physics%YMOMENTUM) = rhs%data4d(i,j,k,physics%YMOMENTUM) &
1452 - wp * rhs%data4d(i,j,k,physics%DENSITY)
1453 IF(physics%PRESSURE.GT.0) THEN
1454 rhs%data4d(i,j,k,physics%ENERGY) = rhs%data4d(i,j,k,physics%ENERGY) &
1455 - wp*( rhs%data4d(i,j,k,physics%YMOMENTUM) &
1456 + 0.5 * wp* rhs%data4d(i,j,k,physics%DENSITY))
1457 IF (have_potential) THEN
1458 rhs%data4d(i,j,k,physics%ENERGY) = &
1459 rhs%data4d(i,j,k,physics%ENERGY) - pot(i,j,k,1) * rhs%data4d(i,j,k,physics%DENSITY)
1460 END IF
1461 END IF
1462 END DO
1463 END DO
1464 END DO
1465
1466 CASE(3) !spherical, tancylindrical
1467 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX
1468 DO j=mesh%JMIN-mesh%JP1,mesh%JMAX
1469!NEC$ IVDEP
1470 DO i=mesh%IMIN-mesh%IP1,mesh%IMAX
1471 wp = 0.5*(this%w(i,j)+this%w(i+1,j)) + mesh%hz%faces(i+1,j,k,1)*mesh%OMEGA
1472 IF(physics%PRESSURE.GT.0) THEN
1473 this%xfluxdydz(i,j,k,physics%ENERGY) &
1474 = this%xfluxdydz(i,j,k,physics%ENERGY) &
1475 + wp * (0.5 * wp * this%xfluxdydz(i,j,k,physics%DENSITY) &
1476 + this%xfluxdydz(i,j,k,physics%ZMOMENTUM))
1477 IF (have_potential) THEN
1478 this%xfluxdydz(i,j,k,physics%ENERGY) = this%xfluxdydz(i,j,k,physics%ENERGY) &
1479 + pot(i,j,k,2)*this%xfluxdydz(i,j,k,physics%DENSITY)
1480 END IF
1481 END IF
1482
1483 this%xfluxdydz(i,j,k,physics%ZMOMENTUM) &
1484 = (this%xfluxdydz(i,j,k,physics%ZMOMENTUM) &
1485 + wp * this%xfluxdydz(i,j,k,physics%DENSITY)) * mesh%hz%faces(i+1,j,k,1)
1486
1487 wp = this%w(i,j) + mesh%hz%bcenter(i,j,k)*mesh%OMEGA
1488
1489 IF(physics%PRESSURE.GT.0) THEN
1490 this%yfluxdzdx(i,j,k,physics%ENERGY) &
1491 = this%yfluxdzdx(i,j,k,physics%ENERGY) &
1492 + wp * ( 0.5 * wp * this%yfluxdzdx(i,j,k,physics%DENSITY) &
1493 + this%yfluxdzdx(i,j,k,physics%ZMOMENTUM))
1494 IF (have_potential) THEN
1495 this%yfluxdzdx(i,j,k,physics%ENERGY) = this%yfluxdzdx(i,j,k,physics%ENERGY) &
1496 + pot(i,j,k,3)*this%yfluxdzdx(i,j,k,physics%DENSITY)
1497 END IF
1498 END IF
1499
1500 this%yfluxdzdx(i,j,k,physics%ZMOMENTUM) &
1501 = (this%yfluxdzdx(i,j,k,physics%ZMOMENTUM) &
1502 + wp * this%yfluxdzdx(i,j,k,physics%DENSITY)) * mesh%hz%faces(i,j+1,k,1)
1503
1504 wp = this%w(i,j) + mesh%hz%bcenter(i,j,k)*mesh%OMEGA
1505
1506 IF(physics%PRESSURE.GT.0) THEN
1507 this%zfluxdxdy(i,j,k,physics%ENERGY) &
1508 = this%zfluxdxdy(i,j,k,physics%ENERGY) &
1509 + wp * ( 0.5 * wp * this%zfluxdxdy(i,j,k,physics%DENSITY) &
1510 + this%zfluxdxdy(i,j,k,physics%ZMOMENTUM))
1511 IF (have_potential) THEN
1512 this%zfluxdxdy(i,j,k,physics%ENERGY) = this%zfluxdxdy(i,j,k,physics%ENERGY) &
1513 + pot(i,j,k,4)*this%zfluxdxdy(i,j,k,physics%DENSITY)
1514 END IF
1515 END IF
1516
1517 this%zfluxdxdy(i,j,k,physics%ZMOMENTUM) &
1518 = this%zfluxdxdy(i,j,k,physics%ZMOMENTUM) &
1519 + wp * this%zfluxdxdy(i,j,k,physics%DENSITY)
1520
1521 END DO
1522 END DO
1523 END DO
1524
1525 ! set isothermal sound speeds
1526 SELECT TYPE (phys => physics)
1527 CLASS IS (physics_eulerisotherm)
1528 DO k=mesh%KMIN,mesh%KMAX
1529 DO j=mesh%JMIN,mesh%JMAX
1530!NEC$ IVDEP
1531 DO i=mesh%IMIN,mesh%IMAX
1532 wp = pvar%data4d(i,j,k,physics%ZVELOCITY) + this%w(i,j) + mesh%hz%bcenter(i,j,k)*mesh%OMEGA
1533
1534 this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1535 = cvar%data4d(i,j,k,physics%YMOMENTUM) * (pvar%data4d(i,j,k,physics%YVELOCITY)*mesh%cyxy%bcenter(i,j,k) - &
1536 pvar%data4d(i,j,k,physics%XVELOCITY) * mesh%cxyx%bcenter(i,j,k)) &
1537 + pvar%data4d(i,j,k,physics%DENSITY)*wp * (wp*mesh%czxz%bcenter(i,j,k) - &
1538 pvar%data4d(i,j,k,physics%XVELOCITY) * mesh%cxzx%bcenter(i,j,k))
1539
1540 this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1541 = cvar%data4d(i,j,k,physics%XMOMENTUM) &
1542 * ( pvar%data4d(i,j,k,physics%XVELOCITY) * mesh%cxyx%center(i,j,k) &
1543 - pvar%data4d(i,j,k,physics%YVELOCITY) * mesh%cyxy%center(i,j,k) ) &
1544 + pvar%data4d(i,j,k,physics%DENSITY) * wp &
1545 * ( wp * mesh%czyz%center(i,j,k) - pvar%data4d(i,j,k,physics%YVELOCITY) * mesh%cyzy%center(i,j,k) )
1546
1547 ! czxz and czyz are eleminated due to the modified divergence operator
1548 IF(physics%ZMOMENTUM.GT.0) &
1549 this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) &
1550 = cvar%data4d(i,j,k,physics%XMOMENTUM) * ( pvar%data4d(i,j,k,physics%XVELOCITY) * mesh%cxzx%center(i,j,k) ) &
1551 + cvar%data4d(i,j,k,physics%YMOMENTUM) * ( pvar%data4d(i,j,k,physics%YVELOCITY) * mesh%cyzy%center(i,j,k) )
1552
1553 IF(physics%PRESSURE.GT.0) THEN
1554 this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1555 = this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1556 +pvar%data4d(i,j,k,physics%PRESSURE) &
1557 *( mesh%cyxy%center(i,j,k) + mesh%czxz%center(i,j,k) )
1558 this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1559 = this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1560 + pvar%data4d(i,j,k,physics%PRESSURE) &
1561 *( mesh%cxyx%center(i,j,k) + mesh%czyz%center(i,j,k) )
1562 IF(physics%ZMOMENTUM.GT.0) &
1563 this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) &
1564 = this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) &
1565 + pvar%data4d(i,j,k,physics%PRESSURE) &
1566 *( mesh%cxzx%center(i,j,k) + mesh%cyzy%center(i,j,k) )
1567
1568 this%geo_src%data4d(i,j,k,physics%ENERGY) = 0.
1569 ELSE
1570
1571 this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1572 = this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1573 +pvar%data4d(i,j,k,physics%DENSITY)*phys%bccsound%data3d(i,j,k)**2 &
1574 * ( mesh%cyxy%center(i,j,k) + mesh%czxz%center(i,j,k) )
1575 this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1576 = this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1577 + pvar%data4d(i,j,k,physics%DENSITY)*phys%bccsound%data3d(i,j,k)**2 &
1578 *( mesh%cxyx%center(i,j,k) + mesh%czyz%center(i,j,k) )
1579 IF(physics%ZMOMENTUM.GT.0) &
1580 this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) &
1581 = this%geo_src%data4d(i,j,k,physics%ZMOMENTUM) &
1582 + pvar%data4d(i,j,k,physics%DENSITY)*phys%bccsound%data3d(i,j,k)**2 &
1583 *( mesh%cxzx%center(i,j,k) + mesh%cyzy%center(i,j,k) )
1584
1585 END IF
1586 END DO
1587 END DO
1588 END DO
1589
1590 CLASS DEFAULT
1591 ! abort
1592 CALL this%Error("timedisc_base::ComputeRHS","physics not supported with rhstype=1")
1593 END SELECT
1594
1595
1596!NEC$ NOVECTOR
1597 DO l=1,physics%VNUM+physics%PNUM
1598!NEC$ OUTERLOOP_UNROLL(8)
1599 DO k=mesh%KMIN,mesh%KMAX
1600 DO j=mesh%JMIN,mesh%JMAX
1601!NEC$ IVDEP
1602 DO i=mesh%IMIN,mesh%IMAX
1603 ! update right hand side of ODE
1604 IF(l.EQ.physics%ZMOMENTUM) THEN
1605 rhs%data4d(i,j,k,l) = mesh%dydzdV%data3d(i,j,k)*(this%xfluxdydz(i,j,k,l) - this%xfluxdydz(i-mesh%IP1,j,k,l)) &
1606 / mesh%hz%center(i,j,k) &
1607 + mesh%dzdxdV%data3d(i,j,k)*(this%yfluxdzdx(i,j,k,l) - this%yfluxdzdx(i,j-mesh%JP1,k,l)) &
1608 / mesh%hz%center(i,j,k) &
1609 + mesh%dxdydV%data3d(i,j,k)*(this%zfluxdxdy(i,j,k,l) - this%zfluxdxdy(i,j,k-mesh%KP1,l))
1610 ELSE
1611 rhs%data4d(i,j,k,l) = mesh%dydzdV%data3d(i,j,k)*(this%xfluxdydz(i,j,k,l) - this%xfluxdydz(i-mesh%IP1,j,k,l)) &
1612 + mesh%dzdxdV%data3d(i,j,k)*(this%yfluxdzdx(i,j,k,l) - this%yfluxdzdx(i,j-mesh%JP1,k,l)) &
1613 + mesh%dxdydV%data3d(i,j,k)*(this%zfluxdxdy(i,j,k,l) - this%zfluxdxdy(i,j,k-mesh%KP1,l))
1614 END IF
1615 END DO
1616 END DO
1617 END DO
1618 END DO
1619
1620 DO k=mesh%KMIN,mesh%KMAX
1621 DO j=mesh%JMIN,mesh%JMAX
1622!NEC$ IVDEP
1623 DO i=mesh%IMIN,mesh%IMAX
1624 wp = this%w(i,j) + mesh%radius%bcenter(i,j,k) * mesh%OMEGA
1625 rhs%data4d(i,j,k,physics%ZMOMENTUM) = rhs%data4d(i,j,k,physics%ZMOMENTUM) &
1626 - wp * rhs%data4d(i,j,k,physics%DENSITY)
1627 IF(physics%PRESSURE.GT.0) THEN
1628 rhs%data4d(i,j,k,physics%ENERGY) = rhs%data4d(i,j,k,physics%ENERGY) &
1629 - wp*( rhs%data4d(i,j,k,physics%ZMOMENTUM) &
1630 + 0.5 * wp* rhs%data4d(i,j,k,physics%DENSITY))
1631 IF (have_potential) THEN
1632 rhs%data4d(i,j,k,physics%ENERGY) = &
1633 rhs%data4d(i,j,k,physics%ENERGY) - pot(i,j,k,1) * rhs%data4d(i,j,k,physics%DENSITY)
1634 END IF
1635 END IF
1636 END DO
1637 END DO
1638 END DO
1639
1640 CASE DEFAULT
1641 CALL this%Error("timedisc_base::ComputeRHS","rhstype 1 is not supported for this geometry")
1642 END SELECT
1643
1644
1645
1646!NEC$ NOVECTOR
1647 DO l=1,physics%VNUM+physics%PNUM
1648!NEC$ OUTERLOOP_UNROLL(8)
1649 DO k=mesh%KMIN,mesh%KMAX
1650 DO j=mesh%JMIN,mesh%JMAX
1651!NEC$ IVDEP
1652 DO i=mesh%IMIN,mesh%IMAX
1653 ! update right hand side of ODE
1654 rhs%data4d(i,j,k,l) = rhs%data4d(i,j,k,l) - this%geo_src%data4d(i,j,k,l) - this%src%data4d(i,j,k,l)
1655 END DO
1656 END DO
1657 END DO
1658 END DO
1659 CASE DEFAULT
1660 CALL this%Error("timedisc_base::ComputeRHS","only rhstype 0 or 1 supported")
1661 END SELECT
1662 END SUBROUTINE computerhs
1663
1666 SUBROUTINE checkdata(this,Mesh,Physics,Fluxes,pvar,cvar,checkdatabm)
1667 USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
1668 USE physics_euler_mod, ONLY : physics_euler
1669 IMPLICIT NONE
1670 !------------------------------------------------------------------------!
1671 CLASS(timedisc_base), INTENT(INOUT) :: this
1672 CLASS(mesh_base), INTENT(IN) :: Mesh
1673 CLASS(physics_base), INTENT(INOUT) :: Physics
1674 CLASS(fluxes_base), INTENT(IN) :: Fluxes
1675 CLASS(marray_compound),INTENT(INOUT):: pvar,cvar
1676 INTEGER, INTENT(IN) :: checkdatabm
1677 !------------------------------------------------------------------------!
1678#ifdef PARALLEL
1679 INTEGER :: err
1680#endif
1681 REAL :: val
1682 !------------------------------------------------------------------------!
1683 IF(iand(checkdatabm,check_csound).NE.check_nothing) THEN
1684 ! check for speed of sound
1685 SELECT TYPE(phys => physics)
1686 CLASS IS(physics_eulerisotherm)
1687 val = minval(phys%bccsound%data1d(:))
1688 IF(val.LE.0.) THEN
1689 ! warn now and stop after file output
1690 CALL this%Warning("CheckData","Illegal speed of sound value less than 0.")
1691 this%break = .true.
1692 END IF
1693 CLASS DEFAULT
1694 CALL this%Warning("CheckData","check speed of sound selected, but bccsound not defined")
1695 END SELECT
1696 END IF
1697 IF((iand(checkdatabm,check_pmin).NE.check_nothing)) THEN
1698 ! check for non-isothermal physics with pressure defined
1699 SELECT TYPE(p => pvar)
1700 CLASS IS(statevector_euler)
1701 val = minval(p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX, &
1702 mesh%KMIN:mesh%KMAX))
1703 IF(val.LT.this%pmin) THEN
1704 ! warn now and stop after file output
1705 CALL this%Warning("CheckData","Pressure below allowed pmin value.")
1706 this%break = .true.
1707 END IF
1708 END SELECT
1709 END IF
1710
1711 IF(iand(checkdatabm,check_rhomin).NE.check_nothing) THEN
1712 ! check for physics with density defined
1713 SELECT TYPE(p => pvar)
1714 CLASS IS(statevector_eulerisotherm)
1715 val = minval(p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX, &
1716 mesh%KMIN:mesh%KMAX))
1717 IF(val.LT.this%rhomin) THEN
1718 ! warn now and stop after file output
1719 CALL this%Warning("CheckData","Density below allowed rhomin value.")
1720 this%break = .true.
1721 END IF
1722 END SELECT
1723 END IF
1724
1725 IF(iand(checkdatabm,check_invalid).NE.check_nothing) THEN
1726 IF(any((cvar%data1d.NE.cvar%data1d).OR.(pvar%data1d.NE.pvar%data1d))) THEN
1727 ! warn now and stop after file output
1728 CALL this%Warning("CheckData","Found NaN in pvar or cvar.")
1729 this%break = .true.
1730 END IF
1731 IF(any((cvar%data1d.GT.huge(cvar%data1d)).OR.pvar%data1d.GT.huge(pvar%data1d))) THEN
1732 ! warn now and stop after file output
1733 CALL this%Warning("CheckData","Found Infinity in pvar or cvar.")
1734 this%break = .true.
1735 END IF
1736 END IF
1737
1738#ifdef PARALLEL
1739 CALL mpi_allreduce(mpi_in_place, this%break, 1, mpi_logical, mpi_lor, &
1740 mesh%comm_cart, err)
1741#endif
1742 END SUBROUTINE checkdata
1743
1744 PURE FUNCTION getorder(this) RESULT(odr)
1745 IMPLICIT NONE
1746 !------------------------------------------------------------------------!
1747 CLASS(timedisc_base), INTENT(IN) :: this
1748 INTEGER :: odr
1749 !------------------------------------------------------------------------!
1750 odr = this%order
1751 END FUNCTION getorder
1752
1753 PURE FUNCTION getcfl(this) RESULT(cfl)
1754 IMPLICIT NONE
1755 !------------------------------------------------------------------------!
1756 CLASS(timedisc_base), INTENT(IN) :: this
1757 REAL :: cfl
1758 !------------------------------------------------------------------------!
1759 cfl = this%CFL
1760 END FUNCTION getcfl
1761
1762
1781 SUBROUTINE fargoadvectionx(this,Fluxes,Mesh,Physics,Sources)
1782 IMPLICIT NONE
1783 !------------------------------------------------------------------------!
1784 CLASS(timedisc_base), INTENT(INOUT) :: this
1785 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
1786 CLASS(mesh_base), INTENT(IN) :: Mesh
1787 CLASS(physics_base), INTENT(INOUT) :: Physics
1788 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
1789 !------------------------------------------------------------------------!
1790 INTEGER :: i,j,k,l
1791#ifdef PARALLEL
1792 CHARACTER(LEN=80) :: str
1793 INTEGER :: status(MPI_STATUS_SIZE)
1794 INTEGER :: ierror
1795 REAL :: mpi_buf(2*Mesh%GNUM)
1796 REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX) :: buf
1797#endif
1798 !------------------------------------------------------------------------!
1799 ! determine step size of integer shift and length of remaining transport step
1800 ! first compute the whole step
1801 this%delxy(:,:) = this%w(:,:) * this%dt / mesh%dlx%data3d(mesh%IMIN,:,:)
1802
1803#ifdef PARALLEL
1804 ! make sure all MPI processes use the same step if domain is decomposed
1805 ! along the x-direction (can be different due to round-off errors)
1806 IF (mesh%dims(1).GT.1) THEN
1807 CALL mpi_allreduce(mpi_in_place,this%delxy,(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1), &
1808 default_mpi_real,mpi_min,mesh%Icomm,ierror)
1809 END IF
1810#endif
1811
1812 ! then subdivide into integer shift and remaining linear advection step
1813!NEC$ COLLAPSE
1814 DO k = mesh%KGMIN,mesh%KGMAX
1815!NEC$ IVDEP
1816 DO j = mesh%JGMIN,mesh%JGMAX
1817 this%shift(j,k) = nint(this%delxy(j,k))
1818 this%delxy(j,k) = this%delxy(j,k)-dble(this%shift(j,k))
1819 END DO
1820 END DO
1821
1822!NEC$ SHORTLOOP
1823 DO l=1,physics%VNUM+physics%PNUM
1824 this%dq%data3d(mesh%IGMIN,:,:) = this%cvar%data4d(mesh%IGMIN+1,:,:,l) - this%cvar%data4d(mesh%IGMIN,:,:,l)
1825!NEC$ IVDEP
1826 DO i=mesh%IGMIN+1,mesh%IGMAX-1
1827 this%dq%data3d(i,:,:) = this%cvar%data4d(i+1,:,:,l)-this%cvar%data4d(i,:,:,l)
1828 ! apply minmod limiter
1829 WHERE (sign(1.0,this%dq%data3d(i-1,:,:))*sign(1.0,this%dq%data3d(i,:,:)).GT.0.)
1830 this%dql%data3d(i,:,:) = sign(min(abs(this%dq%data3d(i-1,:,:)),abs(this%dq%data3d(i,:,:))),this%dq%data3d(i-1,:,:))
1831 ELSEWHERE
1832 this%dql%data3d(i,:,:) = 0.
1833 END WHERE
1834 END DO
1835!NEC$ IVDEP
1836 DO i=mesh%IMIN-1,mesh%IMAX
1837 WHERE(this%delxy(:,:).GT.0.)
1838 this%flux%data3d(i,:,:) = this%cvar%data4d(i,:,:,l) + .5 * this%dql%data3d(i,:,:) * (1. - this%delxy(:,:))
1839 ELSEWHERE
1840 this%flux%data3d(i,:,:) = this%cvar%data4d(i+1,:,:,l) - .5*this%dql%data3d(i+1,:,:)*(1. + this%delxy(:,:))
1841 END WHERE
1842 this%cvar%data4d(i,:,:,l) = this%cvar%data4d(i,:,:,l) - this%delxy(:,:)*(this%flux%data3d(i,:,:) - this%flux%data3d(i-1,:,:))
1843 END DO
1844 END DO
1845
1846!#ifdef PARALLEL
1847! ! We only need to do something, if we (also) are dealing with domain decomposition in
1848! ! the second (phi) direction
1849! IF(PAR_DIMS.GT.1) THEN
1850! DO i=GMIN1,GMAX1
1851! IF(this%shift(i,k).GT.0) THEN
1852! DO l=1,Physics%VNUM+Physics%PNUM
1853! IF(Mesh%SN_shear) THEN
1854! this%buf(k,1:this%shift(i)) = this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k)
1855! ELSE
1856! this%buf(k,1:this%shift(i)) = this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k)
1857! END IF
1858! END DO
1859! CALL MPI_Sendrecv_replace(&
1860! this%buf,&
1861! this%shift(i)*Physics%VNUM, &
1862! DEFAULT_MPI_REAL, &
1863! Mesh%neighbor(EAST), i+Mesh%GNUM, &
1864! Mesh%neighbor(WEST), i+Mesh%GNUM, &
1865! Mesh%comm_cart, status, ierror)
1866! DO k=1,Physics%VNUM
1867! IF(Mesh%SN_shear) THEN
1868! this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k) = this%buf(k,1:this%shift(i))
1869! ELSE
1870! this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k) = this%buf(k,1:this%shift(i))
1871! END IF
1872! END DO
1873! ELSE IF(this%shift(i).LT.0) THEN
1874! DO k=1,Physics%VNUM
1875! IF(Mesh%SN_shear) THEN
1876! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k)
1877! ELSE
1878! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k)
1879! END IF
1880! END DO
1881! CALL MPI_Sendrecv_replace(&
1882! this%buf,&
1883! -this%shift(i)*Physics%VNUM, &
1884! DEFAULT_MPI_REAL, &
1885! Mesh%neighbor(EAST), i+Mesh%GNUM, &
1886! Mesh%neighbor(WEST), i+Mesh%GNUM, &
1887! Mesh%comm_cart, status, ierror)
1888! DO k=1,Physics%VNUM
1889! IF (Mesh%SN_shear) THEN
1890! this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k) = this%buf(k,1:-this%shift(i))
1891! ELSE
1892! this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k) = this%buf(k,1:-this%shift(i))
1893! END IF
1894! END DO
1895! END IF
1896! END DO
1897! END DO
1898! END IF
1899!#endif
1900
1901 ! Integer shift along x- or y-direction
1902!NEC$ SHORTLOOP
1903 DO l=1,physics%VNUM+physics%PNUM
1904!NEC$ IVDEP
1905 DO k=mesh%KGMIN,mesh%KGMAX
1906!NEC$ IVDEP
1907 DO j=mesh%JGMIN,mesh%JGMAX
1908 this%cvar%data4d(mesh%IMIN:mesh%IMAX,j,k,l) &
1909 = cshift(this%cvar%data4d(mesh%IMIN:mesh%IMAX,j,k,l),-this%shift(j,k),1)
1910 END DO
1911 END DO
1912 END DO
1913
1914 END SUBROUTINE fargoadvectionx
1915
1916
1935 SUBROUTINE fargoadvectiony(this,Fluxes,Mesh,Physics,Sources)
1936 IMPLICIT NONE
1937 !------------------------------------------------------------------------!
1938 CLASS(timedisc_base), INTENT(INOUT) :: this
1939 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
1940 CLASS(mesh_base), INTENT(IN) :: Mesh
1941 CLASS(physics_base), INTENT(INOUT) :: Physics
1942 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
1943 !------------------------------------------------------------------------!
1944 INTEGER :: i,j,k,l
1945#ifdef PARALLEL
1946 CHARACTER(LEN=80) :: str
1947 INTEGER :: status(MPI_STATUS_SIZE)
1948 INTEGER :: ierror
1949 REAL :: mpi_buf(2*Mesh%GNUM)
1950 REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX) :: buf
1951#endif
1952 !------------------------------------------------------------------------!
1953 ! determine step size of integer shift and length of remaining transport step
1954 ! first compute the whole step
1955 this%delxy(:,:) = this%w(:,:) * this%dt / mesh%dly%data3d(:,mesh%JMIN,:)
1956
1957#ifdef PARALLEL
1958 ! make sure all MPI processes use the same step if domain is decomposed
1959 ! along the y-direction (can be different due to round-off errors)
1960 IF (mesh%dims(2).GT.1) THEN
1961 CALL mpi_allreduce(mpi_in_place,this%delxy,(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1), &
1962 default_mpi_real,mpi_min,mesh%Jcomm,ierror)
1963 END IF
1964#endif
1965
1966 ! then subdivide into integer shift and remaining linear advection step
1967!NEC$ COLLAPSE
1968 DO k = mesh%KGMIN,mesh%KGMAX
1969!NEC$ IVDEP
1970 DO i = mesh%IGMIN,mesh%IGMAX
1971 this%shift(i,k) = nint(this%delxy(i,k))
1972 this%delxy(i,k) = this%delxy(i,k)-dble(this%shift(i,k))
1973 END DO
1974 END DO
1975
1976!NEC$ SHORTLOOP
1977 DO l=1,physics%VNUM+physics%PNUM
1978 this%dq%data3d(:,mesh%JGMIN,:) = this%cvar%data4d(:,mesh%JGMIN+1,:,l) - this%cvar%data4d(:,mesh%JGMIN,:,l)
1979!NEC$ IVDEP
1980 DO j=mesh%JGMIN+1,mesh%JGMAX-1
1981 this%dq%data3d(:,j,:) = this%cvar%data4d(:,j+1,:,l)-this%cvar%data4d(:,j,:,l)
1982 ! apply minmod limiter
1983 WHERE (sign(1.0,this%dq%data3d(:,j-1,:))*sign(1.0,this%dq%data3d(:,j,:)).GT.0.)
1984 this%dql%data3d(:,j,:) = sign(min(abs(this%dq%data3d(:,j-1,:)),abs(this%dq%data3d(:,j,:))),this%dq%data3d(:,j-1,:))
1985 ELSEWHERE
1986 this%dql%data3d(:,j,:) = 0.
1987 END WHERE
1988 END DO
1989!NEC$ IVDEP
1990 DO j=mesh%JMIN-1,mesh%JMAX
1991 WHERE(this%delxy(:,:).GT.0.)
1992 this%flux%data3d(:,j,:) = this%cvar%data4d(:,j,:,l) + .5 * this%dql%data3d(:,j,:) * (1. - this%delxy(:,:))
1993 ELSEWHERE
1994 this%flux%data3d(:,j,:) = this%cvar%data4d(:,j+1,:,l) - .5*this%dql%data3d(:,j+1,:)*(1. + this%delxy(:,:))
1995 END WHERE
1996 this%cvar%data4d(:,j,:,l) = this%cvar%data4d(:,j,:,l) - this%delxy(:,:)*(this%flux%data3d(:,j,:) - this%flux%data3d(:,j-1,:))
1997 END DO
1998 END DO
1999
2000
2001!#ifdef PARALLEL
2002! ! We only need to do something, if we (also) are dealing with domain decomposition in
2003! ! the second (phi) direction
2004! IF(PAR_DIMS.GT.1) THEN
2005! DO i=GMIN1,GMAX1
2006! IF(this%shift(i,k).GT.0) THEN
2007! DO l=1,Physics%VNUM+Physics%PNUM
2008! IF(Mesh%SN_shear) THEN
2009! this%buf(k,1:this%shift(i)) = this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k)
2010! ELSE
2011! this%buf(k,1:this%shift(i)) = this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k)
2012! END IF
2013! END DO
2014! CALL MPI_Sendrecv_replace(&
2015! this%buf,&
2016! this%shift(i)*Physics%VNUM, &
2017! DEFAULT_MPI_REAL, &
2018! Mesh%neighbor(EAST), i+Mesh%GNUM, &
2019! Mesh%neighbor(WEST), i+Mesh%GNUM, &
2020! Mesh%comm_cart, status, ierror)
2021! DO k=1,Physics%VNUM
2022! IF(Mesh%SN_shear) THEN
2023! this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k) = this%buf(k,1:this%shift(i))
2024! ELSE
2025! this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k) = this%buf(k,1:this%shift(i))
2026! END IF
2027! END DO
2028! ELSE IF(this%shift(i).LT.0) THEN
2029! DO k=1,Physics%VNUM
2030! IF(Mesh%SN_shear) THEN
2031! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k)
2032! ELSE
2033! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k)
2034! END IF
2035! END DO
2036! CALL MPI_Sendrecv_replace(&
2037! this%buf,&
2038! -this%shift(i)*Physics%VNUM, &
2039! DEFAULT_MPI_REAL, &
2040! Mesh%neighbor(EAST), i+Mesh%GNUM, &
2041! Mesh%neighbor(WEST), i+Mesh%GNUM, &
2042! Mesh%comm_cart, status, ierror)
2043! DO k=1,Physics%VNUM
2044! IF (Mesh%SN_shear) THEN
2045! this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k) = this%buf(k,1:-this%shift(i))
2046! ELSE
2047! this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k) = this%buf(k,1:-this%shift(i))
2048! END IF
2049! END DO
2050! END IF
2051! END DO
2052! END DO
2053! END IF
2054!#endif
2055
2056 ! Integer shift along x- or y-direction
2057!NEC$ SHORTLOOP
2058 DO l=1,physics%VNUM+physics%PNUM
2059!NEC$ IVDEP
2060 DO k=mesh%KGMIN,mesh%KGMAX
2061!NEC$ IVDEP
2062 DO i=mesh%IGMIN,mesh%IGMAX
2063 this%cvar%data4d(i,mesh%JMIN:mesh%JMAX,k,l) &
2064 = cshift(this%cvar%data4d(i,mesh%JMIN:mesh%JMAX,k,l),-this%shift(i,k),1)
2065 END DO
2066 END DO
2067 END DO
2068
2069 END SUBROUTINE fargoadvectiony
2070
2089 SUBROUTINE fargoadvectionz(this,Fluxes,Mesh,Physics,Sources)
2090 IMPLICIT NONE
2091 !------------------------------------------------------------------------!
2092 CLASS(timedisc_base), INTENT(INOUT) :: this
2093 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
2094 CLASS(mesh_base), INTENT(IN) :: Mesh
2095 CLASS(physics_base), INTENT(INOUT) :: Physics
2096 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
2097 !------------------------------------------------------------------------!
2098 INTEGER :: i,j,k,l
2099#ifdef PARALLEL
2100 CHARACTER(LEN=80) :: str
2101 INTEGER :: status(MPI_STATUS_SIZE)
2102 INTEGER :: ierror
2103 REAL :: mpi_buf(2*Mesh%GNUM)
2104 REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX) :: buf
2105#endif
2106 !------------------------------------------------------------------------!
2107 ! determine step size of integer shift and length of remaining transport step
2108 ! first compute the whole step
2109 this%delxy(:,:) = this%w(:,:) * this%dt / mesh%dlz%data3d(:,:,mesh%KMIN)
2110
2111#ifdef PARALLEL
2112 ! make sure all MPI processes use the same step if domain is decomposed
2113 ! along the z-direction (can be different due to round-off errors)
2114 IF (mesh%dims(3).GT.1) THEN
2115 CALL mpi_allreduce(mpi_in_place,this%delxy,(mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1), &
2116 default_mpi_real,mpi_min,mesh%Kcomm,ierror)
2117 END IF
2118#endif
2119
2120 ! then subdivide into integer shift and remaining linear advection step
2121!NEC$ COLLAPSE
2122 DO j = mesh%JGMIN,mesh%JGMAX
2123!NEC$ IVDEP
2124 DO i = mesh%IGMIN,mesh%IGMAX
2125 this%shift(i,j) = nint(this%delxy(i,j))
2126 this%delxy(i,j) = this%delxy(i,j)-dble(this%shift(i,j))
2127 END DO
2128 END DO
2129! PRINT *,MINVAL(this%w),MAXVAL(this%w),MINVAL(this%shift),MAXVAL(this%shift),MINVAL(this%delxy),MAXVAL(this%delxy)
2130
2131!NEC$ SHORTLOOP
2132 DO l=1,physics%VNUM+physics%PNUM
2133 this%dq%data3d(:,:,mesh%KGMIN) = this%cvar%data4d(:,:,mesh%KGMIN+1,l) - this%cvar%data4d(:,:,mesh%KGMIN,l)
2134!NEC$ IVDEP
2135 DO k=mesh%KGMIN+1,mesh%KGMAX-1
2136 this%dq%data3d(:,:,k) = this%cvar%data4d(:,:,k+1,l)-this%cvar%data4d(:,:,k,l)
2137 ! apply minmod limiter
2138 WHERE (this%dq%data3d(:,:,k-1)*this%dq%data3d(:,:,k).LE.0.)
2139 ! different left and right sign -> no gradient
2140 this%dql%data3d(:,:,k) = 0.
2141 ELSEWHERE
2142 ! both positive/negative -> take the minium of the absolute value
2143 this%dql%data3d(:,:,k) = sign(min(abs(this%dq%data3d(:,:,k-1)),abs(this%dq%data3d(:,:,k))),this%dq%data3d(:,:,k))
2144 END WHERE
2145! WHERE (SIGN(1.0,this%dq%data3d(:,:,k-1))*SIGN(1.0,this%dq%data3d(:,:,k)).GT.0.)
2146! this%dql%data3d(:,:,k) = SIGN(MIN(ABS(this%dq%data3d(:,:,k-1)),ABS(this%dq%data3d(:,:,k))),this%dq%data3d(:,:,k))
2147! ELSEWHERE
2148! this%dql%data3d(:,:,k) = 0.
2149! END WHERE
2150 END DO
2151! PRINT '(2(ES14.5,3(I4)))',MINVAL(this%dql%data3d),MINLOC(this%dql%data3d),MAXVAL(this%dql%data3d),MAXLOC(this%dql%data3d)
2152!NEC$ IVDEP
2153 DO k=mesh%KMIN-1,mesh%KMAX
2154! IF (k.EQ.Mesh%KMAX) THEN
2155! this%flux%data3d(:,:,k) = this%flux%data3d(:,:,Mesh%KMIN-1)
2156! ELSE
2157 WHERE(this%delxy(:,:).GT.0.)
2158 this%flux%data3d(:,:,k) = this%cvar%data4d(:,:,k,l) + .5 * this%dql%data3d(:,:,k) * (1. - this%delxy(:,:))
2159 ELSEWHERE
2160 this%flux%data3d(:,:,k) = this%cvar%data4d(:,:,k+1,l) - .5*this%dql%data3d(:,:,k+1)*(1. + this%delxy(:,:))
2161 END WHERE
2162! END IF
2163 IF (k.LT.mesh%KMIN) CONTINUE
2164 this%cvar%data4d(:,:,k,l) = this%cvar%data4d(:,:,k,l) - this%delxy(:,:)*(this%flux%data3d(:,:,k) - this%flux%data3d(:,:,k-1))
2165! IF (l.GT.1) STOP
2166! PRINT *,k,this%flux%data3d(Mesh%IMIN,Mesh%JMIN,k),this%dq%data3d(Mesh%IMIN,Mesh%JMIN,k),this%dql%data3d(Mesh%IMIN,Mesh%JMIN,k)
2167 END DO
2168 END DO
2169!#ifdef PARALLEL
2170! ! We only need to do something, if we (also) are dealing with domain decomposition in
2171! ! the second (phi) direction
2172! IF(PAR_DIMS.GT.1) THEN
2173! DO i=GMIN1,GMAX1
2174! IF(this%shift(i,k).GT.0) THEN
2175! DO l=1,Physics%VNUM+Physics%PNUM
2176! IF(Mesh%SN_shear) THEN
2177! this%buf(k,1:this%shift(i)) = this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k)
2178! ELSE
2179! this%buf(k,1:this%shift(i)) = this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k)
2180! END IF
2181! END DO
2182! CALL MPI_Sendrecv_replace(&
2183! this%buf,&
2184! this%shift(i)*Physics%VNUM, &
2185! DEFAULT_MPI_REAL, &
2186! Mesh%neighbor(EAST), i+Mesh%GNUM, &
2187! Mesh%neighbor(WEST), i+Mesh%GNUM, &
2188! Mesh%comm_cart, status, ierror)
2189! DO k=1,Physics%VNUM
2190! IF(Mesh%SN_shear) THEN
2191! this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k) = this%buf(k,1:this%shift(i))
2192! ELSE
2193! this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k) = this%buf(k,1:this%shift(i))
2194! END IF
2195! END DO
2196! ELSE IF(this%shift(i).LT.0) THEN
2197! DO k=1,Physics%VNUM
2198! IF(Mesh%SN_shear) THEN
2199! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k)
2200! ELSE
2201! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k)
2202! END IF
2203! END DO
2204! CALL MPI_Sendrecv_replace(&
2205! this%buf,&
2206! -this%shift(i)*Physics%VNUM, &
2207! DEFAULT_MPI_REAL, &
2208! Mesh%neighbor(EAST), i+Mesh%GNUM, &
2209! Mesh%neighbor(WEST), i+Mesh%GNUM, &
2210! Mesh%comm_cart, status, ierror)
2211! DO k=1,Physics%VNUM
2212! IF (Mesh%SN_shear) THEN
2213! this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k) = this%buf(k,1:-this%shift(i))
2214! ELSE
2215! this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k) = this%buf(k,1:-this%shift(i))
2216! END IF
2217! END DO
2218! END IF
2219! END DO
2220! END DO
2221! END IF
2222!#endif
2223
2224 ! Integer shift along x- or y-direction
2225!NEC$ SHORTLOOP
2226 DO l=1,physics%VNUM+physics%PNUM
2227!NEC$ IVDEP
2228 DO j=mesh%JGMIN,mesh%JGMAX
2229!NEC$ IVDEP
2230 DO i=mesh%IGMIN,mesh%IGMAX
2231 this%cvar%data4d(i,j,mesh%KMIN:mesh%KMAX,l) &
2232 = cshift(this%cvar%data4d(i,j,mesh%KMIN:mesh%KMAX,l),-this%shift(i,j),1)
2233 END DO
2234 END DO
2235 END DO
2236
2237 END SUBROUTINE fargoadvectionz
2238
2239
2243 SUBROUTINE calcbackgroundvelocity(this,Mesh,Physics,pvar,cvar)
2244 IMPLICIT NONE
2245 !------------------------------------------------------------------------!
2246 CLASS(timedisc_base), INTENT(INOUT) :: this
2247 CLASS(mesh_base), INTENT(IN) :: Mesh
2248 CLASS(physics_base), INTENT(INOUT) :: Physics
2249 CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
2250 !------------------------------------------------------------------------!
2251 REAL :: wi
2252 INTEGER :: i,j,k,v_idx
2253#ifdef PARALLEL
2254 INTEGER :: ierror
2255#endif
2256 !------------------------------------------------------------------------!
2257 ! currently only eulerisotherm and derived physics are supported
2258 ! if fargo is enabled (see InitTimedisc)
2259 SELECT TYPE(p => pvar)
2260 CLASS IS(statevector_eulerisotherm)
2261 v_idx = mesh%fargo%GetDirection() ! index of fargo transport velocity corresponds with transport direction
2262 SELECT CASE(mesh%fargo%GetDirection())
2263 CASE(1)
2264 ! transform back to real, i.e. not co-moving, quantities
2265 CALL physics%AddBackgroundVelocityX(mesh,this%w,pvar,cvar)
2266 DO k=mesh%KGMIN,mesh%KGMAX
2267 DO j=mesh%JGMIN,mesh%JGMAX
2268 ! some up all xvelocities along the x-direction
2269 wi = sum(p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,k,v_idx)*this%wfactor(j,k))
2270#ifdef PARALLEL
2271 ! extend the sum over all partitions
2272 IF(mesh%dims(1).GT.1) THEN
2273 CALL mpi_allreduce(mpi_in_place, wi, 1, default_mpi_real, mpi_sum, &
2274 mesh%Icomm, ierror)
2275 END IF
2276#endif
2277 ! set new background velocity to the arithmetic mean of the
2278 ! xvelocity field along the x-direction
2279 this%w(j,k) = wi / mesh%INUM
2280 END DO
2281 END DO
2282 CASE(2)
2283 ! transform back to real, i.e. not co-moving, quantities
2284 CALL physics%AddBackgroundVelocityY(mesh,this%w,pvar,cvar)
2285 ! if x-velocity is suppressed, i.e. zero bit not set -> first velocity component is y-velocity
2286 IF (.NOT.btest(mesh%VECTOR_COMPONENTS,0)) v_idx = 1
2287 DO k=mesh%KGMIN,mesh%KGMAX
2288 DO i=mesh%IGMIN,mesh%IGMAX
2289 ! some up all yvelocities along the y-direction
2290 wi = sum(p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,k,v_idx)*this%wfactor(i,k))
2291#ifdef PARALLEL
2292 ! extend the sum over all partitions
2293 IF(mesh%dims(2).GT.1) THEN
2294 CALL mpi_allreduce(mpi_in_place, wi, 1, default_mpi_real, mpi_sum, &
2295 mesh%Jcomm, ierror)
2296 END IF
2297#endif
2298 this%w(i,k) = wi / mesh%JNUM
2299 END DO
2300 END DO
2301 CASE(3)
2302 ! transform back to real, i.e. not co-moving, quantities
2303 CALL physics%AddBackgroundVelocityZ(mesh,this%w,pvar,cvar)
2304 ! last velocity component should be the z-velocity
2305 v_idx = physics%VDIM ! could be < 3
2306 DO j=mesh%JGMIN,mesh%JGMAX
2307 DO i=mesh%IGMIN,mesh%IGMAX
2308 ! some up all zvelocities along the z-direction
2309 wi = sum(p%velocity%data4d(i,j,mesh%KMIN:mesh%KMAX,v_idx)*this%wfactor(i,j))
2310#ifdef PARALLEL
2311 ! extend the sum over all partitions
2312 IF(mesh%dims(3).GT.1) THEN
2313 CALL mpi_allreduce(mpi_in_place, wi, 1, default_mpi_real, mpi_sum, &
2314 mesh%Kcomm, ierror)
2315 END IF
2316#endif
2317 ! set new background velocity to the arithmetic mean of the
2318 ! yvelocity field along the y-direction
2319 this%w(i,j) = wi / mesh%KNUM
2320 END DO
2321 END DO
2322 END SELECT
2323 CLASS DEFAULT
2324 CALL this%Error("timedisc_base::CalcBackgroundVelocity","physics currently not supported with fargo transport")
2325 END SELECT
2326 END SUBROUTINE calcbackgroundvelocity
2327
2328
2355 FUNCTION getcentrifugalvelocity(this,Mesh,Physics,Fluxes,Sources,&
2356 dir_omega_,accel_,centrot) RESULT(velocity)
2357 USE physics_euler_mod, ONLY: physics_euler
2358 USE physics_eulerisotherm_mod, ONLY: physics_eulerisotherm
2359 IMPLICIT NONE
2360 !------------------------------------------------------------------------!
2361 CLASS(timedisc_base), INTENT(INOUT) :: this
2362 CLASS(mesh_base), INTENT(IN) :: mesh
2363 CLASS(physics_base), INTENT(INOUT) :: physics
2364 CLASS(fluxes_base), INTENT(INOUT) :: fluxes
2365 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: sources
2366 REAL, OPTIONAL, INTENT(IN) :: dir_omega_(3)
2367 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VDIM), &
2368 OPTIONAL, INTENT(IN) :: accel_
2369 REAL, OPTIONAL, INTENT(IN) :: centrot(3)
2370 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VDIM) &
2371 :: velocity
2372 !------------------------------------------------------------------------!
2373 CLASS(marray_base), ALLOCATABLE &
2374 :: accel,dist_axis_projected,ephi_projected,eomega,tmpvec, &
2375 posvec,tmp
2376 REAL :: omega2,dir_omega(3)
2377 INTEGER :: k,l,err
2378 !------------------------------------------------------------------------!
2379 ALLOCATE(accel,dist_axis_projected,ephi_projected,posvec,eomega,tmpvec,tmp, &
2380 stat=err)
2381 IF (err.NE.0) CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
2382 "Unable to allocate memory.")
2383 ! initialize marrays
2384 tmp = marray_base()
2385 tmpvec = marray_base(3)
2386 eomega = marray_base(3)
2387 posvec = marray_base(3)
2388 dist_axis_projected = marray_base(physics%VDIM)
2389 ephi_projected = marray_base(physics%VDIM)
2390 accel = marray_base(physics%VDIM)
2391
2392 ! do some sanity checks on the input data
2393 IF (PRESENT(dir_omega_)) THEN
2394 dir_omega(:) = dir_omega_(:)
2395 omega2 = sum(dir_omega(:)*dir_omega(:))
2396 ! omega must not be the zero vector
2397 IF (.NOT. (omega2.GT.0.0)) &
2398 CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
2399 "omega must not be the zero vector")
2400 ! normalize dir_omega
2401 dir_omega(:) = dir_omega(:) / sqrt(omega2)
2402 ELSE
2403 ! default is rotation in a plane perpendicular to ê_z
2404 dir_omega(1:3) = (/ 0.0, 0.0, 1.0 /)
2405 END IF
2406
2407 ! compute (local) curvilinear vector components of the unit
2408 ! vectors pointing in the direction of the rotational axis
2409 tmpvec%data2d(:,1) = dir_omega(1)
2410 tmpvec%data2d(:,2) = dir_omega(2)
2411 tmpvec%data2d(:,3) = dir_omega(3)
2412 CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,tmpvec%data4d,eomega%data4d)
2413
2414 ! 1. Get the curvilinear components of all vectors pointing
2415 ! from the center of rotation into each cell using cartesian
2416 ! coordinates.
2417 ! The cartesian coordinates are the components of the position
2418 ! vector with respect to the cartesian standard orthonormal
2419 ! basis {ê_x, ê_y, ê_z}.
2420 IF (PRESENT(centrot)) THEN
2421 IF ((mesh%rotsym.GT.0.0).AND.any(centrot(:).NE.0.0)) &
2422 CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
2423 "if rotational symmetry is enabled center of rotation must be the origin")
2424 ! Translate the position vector to the center of rotation.
2425 DO k=1,3
2426 tmpvec%data4d(:,:,:,k) = mesh%bccart(:,:,:,k) - centrot(k)
2427 END DO
2428 ! compute curvilinear components of translated position vectors
2429 CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,tmpvec%data4d,posvec%data4d)
2430 ELSE
2431 ! default center of rotation is the origin of the mesh
2432 ! -> use position vectors with respect to the origin
2433 posvec%data4d = mesh%posvec%bcenter
2434 END IF
2435
2436 ! 2. Compute vectors pointing from axis of rotation into each cell
2437 ! and the azimuthal unit vector ephi
2438 tmp%data1d(:) = sum(posvec%data2d(:,:)*eomega%data2d(:,:),dim=2)
2439 posvec%data2d(:,1) = posvec%data2d(:,1) - tmp%data1d(:)*eomega%data2d(:,1)
2440 posvec%data2d(:,2) = posvec%data2d(:,2) - tmp%data1d(:)*eomega%data2d(:,2)
2441 posvec%data2d(:,3) = posvec%data2d(:,3) - tmp%data1d(:)*eomega%data2d(:,3)
2442
2443 ! distance to axis
2444 tmp%data1d(:) = sqrt(sum(posvec%data2d(:,:)*posvec%data2d(:,:),dim=2))
2445 ! use tmpvec for temporary storage of ephi
2446 CALL cross_product(eomega%data2d(:,1),eomega%data2d(:,2),eomega%data2d(:,3), &
2447 posvec%data2d(:,1)/tmp%data1d(:),posvec%data2d(:,2)/tmp%data1d(:), &
2448 posvec%data2d(:,3)/tmp%data1d(:),tmpvec%data2d(:,1),tmpvec%data2d(:,2),tmpvec%data2d(:,3))
2449
2450 ! project it onto the mesh
2451 SELECT TYPE(phys => physics)
2452 CLASS IS(physics_eulerisotherm)
2453 l = 1
2454 DO k=1,3
2455 ! check if vector component is used
2456 IF (btest(mesh%VECTOR_COMPONENTS,k-1)) THEN
2457 dist_axis_projected%data2d(:,l) = posvec%data2d(:,k)
2458 ephi_projected%data2d(:,l) = tmpvec%data2d(:,k)
2459 l = l + 1
2460 END IF
2461 END DO
2462 CLASS DEFAULT
2463 CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
2464 "physics not supported")
2465 END SELECT
2466
2467 ! 3. determine the acceleration
2468 IF(PRESENT(accel_)) THEN
2469 accel%data4d(:,:,:,:) = accel_(:,:,:,:)
2470 ELSE
2471 ! Works only for centrot = (0,0,0)
2472 IF(PRESENT(centrot)) THEN
2473 IF(any(centrot(:).NE.0.0)) &
2474 CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
2475 "You are not allowed to define centrot without accel.")
2476 END IF
2477 ! This may not work for physics with unusual conservative variables.
2478 ! We assume: conservative momentum = density * velocity
2479 ! sanity check if someone derives new physics from the euler modules
2480 ! which does not have this property
2481 SELECT TYPE(phys => physics)
2482 TYPE IS(physics_eulerisotherm)
2483 ! supported -> do nothing
2484 TYPE IS(physics_euler)
2485 ! supported -> do nothing
2486 CLASS DEFAULT
2487 CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
2488 "physics not supported")
2489 END SELECT
2490 SELECT TYPE(phys => physics)
2491 CLASS IS(physics_eulerisotherm)
2492 ! evaluate the right hand side for given (preliminary) initial data
2493 CALL this%ComputeRHS(mesh,phys,sources,fluxes,this%time,0.0, &
2494 this%pvar,this%cvar,this%checkdatabm,this%rhs)
2495 SELECT TYPE(pvar => this%pvar)
2496 CLASS IS(statevector_eulerisotherm)
2497 SELECT TYPE(rhs => this%rhs)
2498 CLASS IS(statevector_eulerisotherm)
2499 DO k=1,phys%VDIM
2500 accel%data2d(:,k) = -rhs%momentum%data2d(:,k) / pvar%density%data1d(:)
2501 END DO
2502 END SELECT
2503 END SELECT
2504 END SELECT
2505 END IF
2506
2507 ! 4. compute absolute value of azimuthal velocity
2508 tmp%data1d(:) = sqrt(max(0.0,-sum(accel%data2d(:,:)*dist_axis_projected%data2d(:,:),dim=2)))
2509
2510 ! 5. Compute components of the azimuthal velocity vector
2511 DO k=1,physics%VDIM
2512 velocity(:,:,:,k) = tmp%data3d(:,:,:) * ephi_projected%data4d(:,:,:,k)
2513 END DO
2514
2515 DEALLOCATE(accel,dist_axis_projected,ephi_projected,posvec,eomega,tmpvec,tmp)
2516
2517 CONTAINS
2518
2519 ELEMENTAL SUBROUTINE cross_product(ax,ay,az,bx,by,bz,aXbx,aXby,aXbz)
2520 IMPLICIT NONE
2521 !------------------------------------------------------------------------!
2522 REAL, INTENT(IN) :: ax,ay,az,bx,by,bz
2523 REAL, INTENT(OUT) :: axbx,axby,axbz
2524 !------------------------------------------------------------------------!
2525 axbx = ay*bz - az*by
2526 axby = az*bx - ax*bz
2527 axbz = ax*by - ay*bx
2528 END SUBROUTINE cross_product
2529
2530 END FUNCTION getcentrifugalvelocity
2531
2532 SUBROUTINE finalize_base(this)
2533 IMPLICIT NONE
2534 !------------------------------------------------------------------------!
2535 CLASS(timedisc_base) :: this
2536 !------------------------------------------------------------------------!
2537 IF (.NOT.this%Initialized()) &
2538 CALL this%Error("CloseTimedisc","not initialized")
2539 ! call boundary destructor
2540 CALL this%Boundary%Finalize()
2541
2542 DEALLOCATE( &
2543 this%pvar,this%cvar,this%ptmp,this%ctmp,this%cold, &
2544 this%geo_src,this%src,this%rhs, &
2545 this%xfluxdydz,this%yfluxdzdx,this%zfluxdxdy,this%amax,this%tol_abs,&
2546 this%dtmean,this%dtstddev,this%time)
2547
2548 IF (ASSOCIATED(this%cerr)) DEALLOCATE(this%cerr)
2549 IF (ASSOCIATED(this%cerr_max)) DEALLOCATE(this%cerr_max)
2550 IF (ASSOCIATED(this%w)) DEALLOCATE(this%w)
2551 IF (ASSOCIATED(this%wfactor)) DEALLOCATE(this%wfactor)
2552 IF (ASSOCIATED(this%delxy))DEALLOCATE(this%delxy)
2553 IF (ASSOCIATED(this%shift))DEALLOCATE(this%shift)
2554#ifdef PARALLEL
2555 IF(ASSOCIATED(this%buf)) DEALLOCATE(this%buf)
2556#endif
2557 IF(ASSOCIATED(this%bflux)) DEALLOCATE(this%bflux)
2558 IF(ASSOCIATED(this%solution)) DEALLOCATE(this%solution)
2559 END SUBROUTINE finalize_base
2560
2561
2562END MODULE timedisc_base_mod
pc
elemental real function pc(tau, rho_s0, cs_inf, gamma)
Definition: agndisk.f90:515
boundary_generic_mod
Generic boundary module.
Definition: boundary_generic.f90:44
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
fluxes_base_mod
base module for numerical flux functions
Definition: fluxes_base.f90:39
geometry_cartesian_mod
defines properties of a 3D cartesian mesh
Definition: geometry_cartesian.f90:36
geometry_cylindrical_mod
defines properties of a 3D cylindrical mesh
Definition: geometry_cylindrical.f90:36
geometry_logcylindrical_mod
defines properties of a 3D logcylindrical mesh
Definition: geometry_logcylindrical.f90:33
geometry_spherical_mod
defines properties of a 3D spherical mesh
Definition: geometry_spherical.f90:42
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
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_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::south
integer, parameter south
named constant for southern boundary
Definition: mesh_base.f90:101
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::conservative
@, public conservative
Definition: physics_base.f90:67
physics_euler_mod
physics module for 1D,2D and 3D non-isothermal Euler equations
Definition: physics_euler.f90:41
physics_eulerisotherm_mod
physics module for 1D,2D and 3D isothermal Euler equations
Definition: physics_eulerisotherm.f90:44
physics_generic_mod
constructor for physics class
Definition: physics_generic.f90:35
reconstruction_base_mod
base module for reconstruction process
Definition: reconstruction_base.f90:35
sources_base_mod
generic source terms module providing functionaly common to all source terms
Definition: sources_base.f90:42
sources_generic_mod
module to manage list of source terms
Definition: sources_generic.f90:37
sources_gravity_mod
generic gravity terms module providing functionaly common to all gravity terms
Definition: sources_gravity.f90:41
timedisc_base_mod
Definition: timedisc_base.f90:57
timedisc_base_mod::dtcause_smallerr
integer, parameter, public dtcause_smallerr
Definition: timedisc_base.f90:195
timedisc_base_mod::check_tmin
integer, parameter, public check_tmin
Definition: timedisc_base.f90:204
timedisc_base_mod::finalize_base
subroutine finalize_base(this)
Definition: timedisc_base.f90:2533
timedisc_base_mod::check_csound
integer, parameter, public check_csound
Definition: timedisc_base.f90:200
timedisc_base_mod::acceptsolution
subroutine acceptsolution(this, Mesh, Physics, Sources, Fluxes, time, dt, iter)
Definition: timedisc_base.f90:912
timedisc_base_mod::check_invalid
integer, parameter, public check_invalid
Definition: timedisc_base.f90:203
timedisc_base_mod::cash_karp
integer, parameter, public cash_karp
Definition: timedisc_base.f90:189
timedisc_base_mod::getcfl
pure real function getcfl(this)
Definition: timedisc_base.f90:1754
timedisc_base_mod::fargo_method
character(len=40), dimension(3), parameter fargo_method
Definition: timedisc_base.f90:206
timedisc_base_mod::computerhs
subroutine computerhs(this, Mesh, Physics, Sources, Fluxes, time, dt, pvar, cvar, checkdatabm, rhs)
compute the RHS of the spatially discretized PDE
Definition: timedisc_base.f90:1074
timedisc_base_mod::check_nothing
integer, parameter, public check_nothing
Definition: timedisc_base.f90:198
timedisc_base_mod::calctimestep
real function calctimestep(this, Mesh, Physics, Sources, Fluxes, time, dtcause)
Determines the CFL time step and time step limits due to source terms.
Definition: timedisc_base.f90:819
timedisc_base_mod::setoutput
subroutine setoutput(this, Mesh, Physics, config, IO)
Definition: timedisc_base.f90:535
timedisc_base_mod::check_rhomin
integer, parameter, public check_rhomin
Definition: timedisc_base.f90:202
timedisc_base_mod::dormand_prince
integer, parameter, public dormand_prince
Definition: timedisc_base.f90:190
timedisc_base_mod::checkdata
subroutine checkdata(this, Mesh, Physics, Fluxes, pvar, cvar, checkdatabm)
compute the RHS of the spatially discretized PDE
Definition: timedisc_base.f90:1667
timedisc_base_mod::fargoadvectiony
subroutine fargoadvectiony(this, Fluxes, Mesh, Physics, Sources)
Calculates the linear transport step in Fargo Scheme along y-axis .
Definition: timedisc_base.f90:1936
timedisc_base_mod::inittimedisc
subroutine inittimedisc(this, Mesh, Physics, config, IO, ttype, tname)
Definition: timedisc_base.f90:225
timedisc_base_mod::dtcause_erradj
integer, parameter, public dtcause_erradj
Definition: timedisc_base.f90:194
timedisc_base_mod::rejectsolution
subroutine rejectsolution(this, Mesh, Physics, Sources, Fluxes, time, dt)
Definition: timedisc_base.f90:962
timedisc_base_mod::check_all
integer, parameter, public check_all
Definition: timedisc_base.f90:199
timedisc_base_mod::integrationstep
subroutine integrationstep(this, Mesh, Physics, Sources, Fluxes, iter, config, IO)
Definition: timedisc_base.f90:650
timedisc_base_mod::modified_euler
integer, parameter, public modified_euler
Definition: timedisc_base.f90:187
timedisc_base_mod::dtcause_fileio
integer, parameter, public dtcause_fileio
smallest ts due to fileio
Definition: timedisc_base.f90:196
timedisc_base_mod::rk_fehlberg
integer, parameter, public rk_fehlberg
Definition: timedisc_base.f90:188
timedisc_base_mod::getorder
pure integer function getorder(this)
Definition: timedisc_base.f90:1745
timedisc_base_mod::computeerror
real function computeerror(this, Mesh, Physics, cvar_high, cvar_low)
Definition: timedisc_base.f90:1013
timedisc_base_mod::dtcause_cfl
integer, parameter, public dtcause_cfl
Definition: timedisc_base.f90:193
timedisc_base_mod::getcentrifugalvelocity
real function, dimension(mesh%igmin:mesh%igmax, mesh%jgmin:mesh%jgmax, mesh%kgmin:mesh%kgmax, physics%vdim) getcentrifugalvelocity(this, Mesh, Physics, Fluxes, Sources, dir_omega_, accel_, centrot)
Compute velocity leading to a centrifugal acceleration with respect to some given axis of rotation wh...
Definition: timedisc_base.f90:2357
timedisc_base_mod::check_pmin
integer, parameter, public check_pmin
Definition: timedisc_base.f90:201
timedisc_base_mod::fargoadvectionx
subroutine fargoadvectionx(this, Fluxes, Mesh, Physics, Sources)
Calculates the linear transport step in Fargo Scheme along x-axis .
Definition: timedisc_base.f90:1782
timedisc_base_mod::adjusttimestep
real function adjusttimestep(this, maxerr, dtold)
adjust the time step
Definition: timedisc_base.f90:760
timedisc_base_mod::fargoadvectionz
subroutine fargoadvectionz(this, Fluxes, Mesh, Physics, Sources)
Calculates the linear transport step in Fargo Scheme along z-axis .
Definition: timedisc_base.f90:2090
timedisc_base_mod::calcbackgroundvelocity
subroutine calcbackgroundvelocity(this, Mesh, Physics, pvar, cvar)
Calculates new background velocity for fargo advection.
Definition: timedisc_base.f90:2244
timedisc_base_mod::ssprk
integer, parameter, public ssprk
Definition: timedisc_base.f90:191
boundary_generic_mod::boundary_generic
Definition: boundary_generic.f90:67
common_dict::dict_typ
Definition: common_dict.f90:123
fluxes_base_mod::fluxes_base
Definition: fluxes_base.f90:61
geometry_cartesian_mod::geometry_cartesian
Definition: geometry_cartesian.f90:41
geometry_cylindrical_mod::geometry_cylindrical
Definition: geometry_cylindrical.f90:41
geometry_logcylindrical_mod::geometry_logcylindrical
Definition: geometry_logcylindrical.f90:39
geometry_spherical_mod::geometry_spherical
Definition: geometry_spherical.f90:49
logging_base_mod::logging_base
common data structure
Definition: logging_base.f90:56
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_euler_mod::physics_euler
Definition: physics_euler.f90:54
physics_euler_mod::statevector_euler
Definition: physics_euler.f90:105
physics_eulerisotherm_mod::physics_eulerisotherm
Definition: physics_eulerisotherm.f90:56
physics_eulerisotherm_mod::statevector_eulerisotherm
Definition: physics_eulerisotherm.f90:122
sources_base_mod::sources_base
Definition: sources_base.f90:57
sources_generic_mod::sources_list
container class to manage the list of source terms
Definition: sources_generic.f90:58
timedisc_base_mod::timedisc_base
Definition: timedisc_base.f90:84
cross_product
elemental subroutine cross_product(ax, ay, az, bx, by, bz, aXbx, aXby, aXbz)
Definition: timedisc_base.f90:2520
Generated by   doxygen 1.9.4