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!----------------------------------------------------------------------------!
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
155 PROCEDURE :: getorder
156 PROCEDURE :: getcfl
157 PROCEDURE :: fargoadvectionx
158 PROCEDURE :: fargoadvectiony
159 PROCEDURE :: fargoadvectionz
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)
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
213 ! constants
215 ssprk, &
219 !--------------------------------------------------------------------------!
220
221CONTAINS
222
223
224 SUBROUTINE inittimedisc(this,Mesh,Physics,config,IO,ttype,tname)
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)
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)
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)
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)
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)
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)
2497 SELECT TYPE(rhs => this%rhs)
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
elemental real function pc(tau, rho_s0, cs_inf, gamma)
Definition: agndisk.f90:515
Generic boundary module.
Dictionary for generic data types.
Definition: common_dict.f90:61
type(logging_base), save this
base module for numerical flux functions
Definition: fluxes_base.f90:39
defines properties of a 3D cartesian mesh
defines properties of a 3D cylindrical mesh
defines properties of a 3D logcylindrical mesh
defines properties of a 3D spherical mesh
Basic fosite module.
subroutine finalize(this)
Destructor of logging_base class.
base class for mesh arrays
Definition: marray_base.f90:36
derived class for compound of mesh arrays
basic mesh module
Definition: mesh_base.f90:72
integer, parameter east
named constant for eastern boundary
Definition: mesh_base.f90:101
integer, parameter bottom
named constant for bottom boundary
Definition: mesh_base.f90:101
integer, parameter south
named constant for southern boundary
Definition: mesh_base.f90:101
integer, parameter top
named constant for top boundary
Definition: mesh_base.f90:101
integer, parameter north
named constant for northern boundary
Definition: mesh_base.f90:101
integer, parameter west
named constant for western boundary
Definition: mesh_base.f90:101
Basic physics module.
@, public primitive
@, public conservative
physics module for 1D,2D and 3D non-isothermal Euler equations
physics module for 1D,2D and 3D isothermal Euler equations
constructor for physics class
base module for reconstruction process
generic source terms module providing functionaly common to all source terms
module to manage list of source terms
generic gravity terms module providing functionaly common to all gravity terms
integer, parameter, public dtcause_smallerr
integer, parameter, public check_tmin
subroutine finalize_base(this)
integer, parameter, public check_csound
subroutine acceptsolution(this, Mesh, Physics, Sources, Fluxes, time, dt, iter)
integer, parameter, public check_invalid
integer, parameter, public cash_karp
pure real function getcfl(this)
character(len=40), dimension(3), parameter fargo_method
subroutine computerhs(this, Mesh, Physics, Sources, Fluxes, time, dt, pvar, cvar, checkdatabm, rhs)
compute the RHS of the spatially discretized PDE
integer, parameter, public check_nothing
real function calctimestep(this, Mesh, Physics, Sources, Fluxes, time, dtcause)
Determines the CFL time step and time step limits due to source terms.
subroutine setoutput(this, Mesh, Physics, config, IO)
integer, parameter, public check_rhomin
integer, parameter, public dormand_prince
subroutine checkdata(this, Mesh, Physics, Fluxes, pvar, cvar, checkdatabm)
compute the RHS of the spatially discretized PDE
subroutine fargoadvectiony(this, Fluxes, Mesh, Physics, Sources)
Calculates the linear transport step in Fargo Scheme along y-axis .
subroutine inittimedisc(this, Mesh, Physics, config, IO, ttype, tname)
integer, parameter, public dtcause_erradj
subroutine rejectsolution(this, Mesh, Physics, Sources, Fluxes, time, dt)
integer, parameter, public check_all
subroutine integrationstep(this, Mesh, Physics, Sources, Fluxes, iter, config, IO)
integer, parameter, public modified_euler
integer, parameter, public dtcause_fileio
smallest ts due to fileio
integer, parameter, public rk_fehlberg
pure integer function getorder(this)
real function computeerror(this, Mesh, Physics, cvar_high, cvar_low)
integer, parameter, public dtcause_cfl
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...
integer, parameter, public check_pmin
subroutine fargoadvectionx(this, Fluxes, Mesh, Physics, Sources)
Calculates the linear transport step in Fargo Scheme along x-axis .
real function adjusttimestep(this, maxerr, dtold)
adjust the time step
subroutine fargoadvectionz(this, Fluxes, Mesh, Physics, Sources)
Calculates the linear transport step in Fargo Scheme along z-axis .
subroutine calcbackgroundvelocity(this, Mesh, Physics, pvar, cvar)
Calculates new background velocity for fargo advection.
integer, parameter, public ssprk
common data structure
basic mesh array class
Definition: marray_base.f90:69
mesh data structure
Definition: mesh_base.f90:122
container class to manage the list of source terms
elemental subroutine cross_product(ax, ay, az, bx, by, bz, aXbx, aXby, aXbz)