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-2018 #
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 !#############################################################################
44 
45 !----------------------------------------------------------------------------!
52 !----------------------------------------------------------------------------!
53 MODULE timedisc_base_mod
54  USE logging_base_mod
55  USE boundary_generic_mod
56  USE mesh_base_mod
57  USE marray_compound_mod
58  USE marray_base_mod
59  USE physics_base_mod
60  USE physics_generic_mod
61  USE sources_base_mod
62  USE sources_gravity_mod
63  USE fluxes_base_mod
64  USE reconstruction_base_mod
65  USE common_dict
66 #ifdef PARALLEL
67 #ifdef HAVE_MPI_MOD
68  USE mpi
69 #endif
70 #endif
71  IMPLICIT NONE
72 #ifdef PARALLEL
73 #ifdef HAVE_MPIF_H
74  include 'mpif.h'
75 #endif
76 #endif
77 !----------------------------------------------------------------------------!
78 PRIVATE
79  TYPE, ABSTRACT, EXTENDS (logging_base) :: timedisc_base
81  CLASS(boundary_generic), ALLOCATABLE :: boundary
82  CLASS(marray_compound), POINTER &
83  :: pvar,cvar,ptmp,ctmp,cold, & !< prim/cons state vectors
84  src,geo_src, & !< source terms
85  rhs, & !< ODE right hand side
86  cerr,cerr_max, & !< error control & output
87  solution
88  INTEGER :: order
89  REAL :: cfl
90  REAL :: dt
91  REAL :: dtold
92  REAL :: dtmin
93  INTEGER :: dtcause,dtmincause
94  REAL :: dtmax
95  REAL,POINTER :: dtmean, dtstddev
96  INTEGER :: dtaccept
97  REAL,POINTER :: time
98  REAL :: stoptime
99  REAL :: dtlimit
100  REAL :: tol_rel
101  REAL :: maxerrold
102  LOGICAL :: break
103  LOGICAL :: always_update_bccsound
105  INTEGER :: maxiter
106  INTEGER :: n_adj
107  INTEGER :: m
108  INTEGER :: degree
109  INTEGER :: rhstype
110  INTEGER :: checkdatabm
111  REAL :: pmin,rhomin,tmin
113  TYPE(marray_base) :: dq, flux, dql
115  REAL :: err_n, h_n
116  REAL, DIMENSION(:), POINTER :: tol_abs
117  REAL :: beta
118 
120  REAL, DIMENSION(:,:,:,:), POINTER :: phi,oldphi_s,&
121  newphi_s
122  REAL, DIMENSION(:), POINTER :: gamma
123  INTEGER :: pc
125  REAL, DIMENSION(:,:,:,:), POINTER :: xfluxdydz,yfluxdzdx,zfluxdxdy
126  REAL, DIMENSION(:,:,:,:), POINTER :: amax
127  REAL, DIMENSION(:,:), POINTER :: bflux
128  LOGICAL :: write_error
129  INTEGER, DIMENSION(:,:), POINTER :: shift=>null()
130  REAL, DIMENSION(:,:), POINTER :: buf=>null()
131  REAL, DIMENSION(:,:), POINTER :: w=>null()
132  REAL, DIMENSION(:,:), POINTER :: delxy =>null()
133 
134  CONTAINS
135 
136  PROCEDURE :: inittimedisc
137  PROCEDURE :: setoutput
138  PROCEDURE :: integrationstep
139  PROCEDURE :: computerhs
140  PROCEDURE :: adjusttimestep
141  PROCEDURE :: calctimestep
142  PROCEDURE :: acceptsolution
143  PROCEDURE :: rejectsolution
144  PROCEDURE :: checkdata
145  PROCEDURE :: computeerror
146  PROCEDURE :: getcentrifugalvelocity
147  PROCEDURE :: getorder
148  PROCEDURE :: getcfl
149  PROCEDURE :: fargoadvectionx
150  PROCEDURE :: fargoadvectiony
151  PROCEDURE :: calcbackgroundvelocity
152  procedure(solveode), DEFERRED :: solveode
153  PROCEDURE :: finalize_base
154  procedure(finalize), DEFERRED :: finalize
155  END TYPE timedisc_base
156 !----------------------------------------------------------------------------!
157  abstract INTERFACE
158  SUBROUTINE solveode(this,Mesh,Physics,Sources,Fluxes,time,dt,err)
159  IMPORT timedisc_base, mesh_base, physics_base, fluxes_base, sources_base
160  IMPLICIT NONE
161  CLASS(timedisc_base), INTENT(INOUT) :: this
162  CLASS(mesh_base), INTENT(IN) :: Mesh
163  CLASS(physics_base), INTENT(INOUT) :: Physics
164  CLASS(sources_base), POINTER :: Sources
165  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
166  REAL, INTENT(IN) :: time
167  REAL, INTENT(INOUT) :: dt, err
168  END SUBROUTINE
169  SUBROUTINE finalize(this)
170  IMPORT timedisc_base
171  IMPLICIT NONE
172  !------------------------------------------------------------------------!
173  CLASS(timedisc_base) :: this
174  END SUBROUTINE
175 
176  END INTERFACE
177 
178  INTEGER, PARAMETER :: modified_euler = 1
179  INTEGER, PARAMETER :: rk_fehlberg = 2
180  INTEGER, PARAMETER :: cash_karp = 3
181  INTEGER, PARAMETER :: dormand_prince = 4
182  INTEGER, PARAMETER :: ssprk = 5
183  !--------------------------------------------------------------------------!
184  INTEGER, PARAMETER :: dtcause_cfl = 0 ! smallest ts due to cfl cond. !
185  INTEGER, PARAMETER :: dtcause_erradj = -1 ! smallest ts due to err adj. !
186  INTEGER, PARAMETER :: dtcause_smallerr = -2 ! smallest ts due to err
187  INTEGER, PARAMETER :: dtcause_fileio = -4
188  ! CheckData_modeuler Bit Mask constants
189  INTEGER, PARAMETER :: check_nothing = int(b'0')
190  INTEGER, PARAMETER :: check_all = not(check_nothing)
191  INTEGER, PARAMETER :: check_csound = int(b'1')
192  INTEGER, PARAMETER :: check_pmin = int(b'10')
193  INTEGER, PARAMETER :: check_rhomin = int(b'100')
194  INTEGER, PARAMETER :: check_invalid = int(b'1000')
195  INTEGER, PARAMETER :: check_tmin = int(b'10000')
196  !--------------------------------------------------------------------------!
197  CHARACTER(LEN=40), PARAMETER, DIMENSION(3) :: fargo_method = (/ &
198  "dynamic background velocity ", &
199  "user supplied fixed background velocity ", &
200  "shearing box " /)
201  PUBLIC :: &
202  ! types
203  timedisc_base, &
204  ! constants
205  modified_euler, rk_fehlberg, cash_karp, dormand_prince, &
206  ssprk, &
207  dtcause_cfl,dtcause_erradj,dtcause_smallerr, dtcause_fileio, &
208  check_all, check_nothing, check_csound, check_rhomin, check_pmin, &
209  check_invalid, check_tmin
210  !--------------------------------------------------------------------------!
211 
212 CONTAINS
213 
214 
215  SUBROUTINE inittimedisc(this,Mesh,Physics,config,IO,ttype,tname)
216  USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
217  USE physics_euler_mod, ONLY : physics_euler
218  USE geometry_cylindrical_mod, ONLY: geometry_cylindrical
219  USE geometry_logcylindrical_mod, ONLY: geometry_logcylindrical
220  USE geometry_cartesian_mod, ONLY: geometry_cartesian
221  IMPLICIT NONE
222  !------------------------------------------------------------------------!
223  CLASS(timedisc_base), INTENT(INOUT) :: this
224  CLASS(mesh_base), INTENT(INOUT) :: Mesh
225  CLASS(physics_base), INTENT(IN) :: Physics
226  TYPE(Dict_TYP), POINTER :: config,IO
227  INTEGER, INTENT(IN) :: ttype
228  CHARACTER(LEN=32), INTENT(IN) :: tname
229  !------------------------------------------------------------------------!
230  INTEGER :: err, d
231  CHARACTER(LEN=32) :: order_str,cfl_str,stoptime_str,dtmax_str,beta_str
232  CHARACTER(LEN=32) :: info_str,shear_direction
233  INTEGER :: method
234  !------------------------------------------------------------------------!
235  CALL this%InitLogging(ttype,tname)
236 
237  IF (.NOT.physics%Initialized().OR..NOT.mesh%Initialized()) &
238  CALL this%Error("InitTimedisc","physics and/or mesh module uninitialized")
239 
240  ! allocate memory for data structures needed in all timedisc modules
241  ALLOCATE( &
242  this%xfluxdydz(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
243  this%yfluxdzdx(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
244  this%zfluxdxdy(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,physics%VNUM), &
245  this%amax(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
246  this%tol_abs(physics%VNUM), &
247  this%dtmean, this%dtstddev, &
248  this%time, &
249  stat = err)
250  IF (err.NE.0) THEN
251  CALL this%Error("InitTimedisc", "Unable to allocate memory.")
252  END IF
253 
254  ! initialize state vectors
255  CALL physics%new_statevector(this%pvar,primitive)
256  CALL physics%new_statevector(this%ptmp,primitive)
257  CALL physics%new_statevector(this%cvar,conservative)
258  CALL physics%new_statevector(this%ctmp,conservative)
259  CALL physics%new_statevector(this%cold,conservative)
260  CALL physics%new_statevector(this%geo_src,conservative)
261  CALL physics%new_statevector(this%src,conservative)
262  CALL physics%new_statevector(this%rhs,conservative)
263  NULLIFY(this%cerr,this%cerr_max)
264 
265  ! initialize all variables
266  this%pvar%data1d(:) = 0.
267  this%ptmp%data1d(:) = 0.
268  this%ctmp%data1d(:) = 0.
269  this%cvar%data1d(:) = 0.
270  this%cold%data1d(:) = 0.
271  this%geo_src%data1d(:) = 0.
272  this%src%data1d(:) = 0.
273  this%rhs%data1d(:) = 0.
274  this%xfluxdydz = 0.
275  this%yfluxdzdx = 0.
276  this%zfluxdxdy = 0.
277  this%amax = 0.
278  this%tol_abs = 0.
279  this%dtmean = 0.
280  this%dtstddev = 0.
281  this%break = .false.
282  this%n_adj = 0
283  this%maxerrold = 0.
284  this%dtaccept = 0
285 
286  CALL getattr(config, "starttime", this%time, 0.0)
287  CALL getattr(config, "method", method)
288  CALL getattr(config, "stoptime", this%stoptime)
289  this%dt = this%stoptime
290  this%dtold = this%dt
291  this%dtmin = this%dt
292 
293  ! set default values
294  ! CFL number
295  CALL getattr(config, "cfl", this%cfl, 0.4)
296 
297  ! time step minimum
298  CALL getattr(config, "dtlimit", this%dtlimit, epsilon(this%dtlimit)*this%stoptime)
299 
300  ! time step maximum in units of [CFL timestep] (Used in Dumka)
301  CALL getattr(config, "dtmax", this%dtmax, 5.0)
302 
303  ! maximum iterations, maxiter <= 0 means no iteration limit
304  CALL getattr(config, "maxiter", this%maxiter, 0)
305 
306  ! relative tolerance for adaptive step size control
307  CALL getattr(config, "tol_rel", this%tol_rel, 0.01)
308 
309  ! absolute tolerance for adaptive step size control
310  this%tol_abs(:) = 0.001
311  CALL getattr(config, "tol_abs", this%tol_abs, this%tol_abs)
312 
313  ! rhs type. 0 = default, 1 = conserve angular momentum
314  CALL getattr(config, "rhstype", this%rhstype, 0)
315  IF (this%rhstype.NE.0.AND.physics%VDIM.NE.2) THEN
316  CALL this%Error("InitTimedisc", "Alternative rhstype only works in 2D.")
317  END IF
318 
319  ! time step friction parameter for PI-Controller
320  CALL getattr(config, "beta", this%beta, 0.08)
321 
322  ! enable/disable update of bccsound at each substep
323  ! Usually bccsound is only needed at the beginning of each integration step
324  ! to determine the next time step. Hence it is not necessary to update
325  ! bccsound at substeps when using higher order integration schemes. Since
326  ! non-isothermal physics involve evaluation of the SQRT function to compute
327  ! the speed of sound it might speed up the simulation a bit if one disables
328  ! the update at substeps.
329  ! Diabling the update is not possible if bccsound is needed somewhere else,
330  ! e.g. some source terms depend on bccsound.
331  CALL getattr(config, "always_update_bccsound", d, 1)
332  IF (d.EQ.0) THEN
333  this%always_update_bccsound = .false.
334  ELSE
335  this%always_update_bccsound = .true.
336  END IF
337 
338  ! check data bit mask
339  ! \todo{expected argument list...}
340  CALL getattr(config, "checkdata", this%checkdatabm, check_all)
341 
342  ! pressure minimum to check if CHECK_PMIN is active
343  CALL getattr(config, "pmin", this%pmin, tiny(this%pmin))
344 
345  ! temperature minimum to check if CHECK_TMIN is active
346  CALL getattr(config, "tmin", this%tmin, tiny(this%tmin))
347 
348  ! density minimum to check if CHECK_RHOMIN is active
349  CALL getattr(config, "rhomin", this%rhomin, tiny(this%rhomin))
350 
351  CALL getattr(config, "output/"//trim(physics%pvarname(1)), d, 1)
352 
353  CALL this%SetOutput(mesh,physics,config,io)
354 
355  ! check if fargo can be used
356  IF(mesh%use_fargo.EQ.1) THEN
357  ! check physics
358  SELECT TYPE(phys => physics)
359  CLASS IS(physics_eulerisotherm)
360  ! check geometry
361  SELECT TYPE(geo=>mesh%Geometry)
362  TYPE IS(geometry_cylindrical)
363  IF (mesh%JNUM.LT.2) CALL this%Error("InitTimedisc", &
364  "fargo advection needs more the one cell in y-direction")
365  ! allocate data arrays used for fargo
366  ALLOCATE( &
367  this%w(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
368  this%delxy(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
369  this%shift(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
370 #ifdef PARALLEL
371  this%buf(physics%VNUM+physics%PNUM,1:mesh%MINJNUM), & !!! NOT CHANGED, because not clear were used !
372 #endif
373  stat = err)
374  IF (err.NE.0) THEN
375  CALL this%Error("InitTimedisc", "Unable to allocate memory for fargo advection.")
376  END IF
377 
378  TYPE IS(geometry_logcylindrical)
379  IF (mesh%JNUM.LT.2) CALL this%Error("InitTimedisc", &
380  "fargo advection needs more the one cell in y-direction")
381  ! allocate data arrays used for fargo
382  ALLOCATE( &
383  this%w(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
384  this%delxy(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
385  this%shift(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
386 #ifdef PARALLEL
387  this%buf(physics%VNUM+physics%PNUM,1:mesh%MINJNUM), & !!! NOT CHANGED, because not clear were used !
388 #endif
389  stat = err)
390  IF (err.NE.0) THEN
391  CALL this%Error("InitTimedisc", "Unable to allocate memory for fargo advection.")
392  END IF
393  TYPE IS(geometry_cartesian) ! in cartesian fargo shift can be chosen in either x- or y-direction
394  IF(mesh%shear_dir.EQ.2) THEN
395  ALLOCATE( &
396  this%w(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
397  this%delxy(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
398  this%shift(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX), &
399 #ifdef PARALLEL
400  this%buf(physics%VNUM+physics%PNUM,1:mesh%MINJNUM), & !!! SEE ABOVE
401 #endif
402  stat = err)
403  IF (err.NE.0) THEN
404  CALL this%Error("InitTimedisc", "Unable to allocate memory for fargo advection.")
405  END IF
406  ELSE IF(mesh%shear_dir.EQ.1) THEN
407  ALLOCATE( &
408  this%w(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
409  this%delxy(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
410  this%shift(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
411 #ifdef PARALLEL
412  this%buf(physics%VNUM+physics%PNUM,1:mesh%MININUM), & !!! SEE ABOVE
413 #endif
414  stat = err)
415  IF (err.NE.0) THEN
416  CALL this%Error("InitTimedisc", "Unable to allocate memory for fargo advection.")
417  END IF
418  ELSE
419  CALL this%Error("InitTimedisc", "Unable to allocate memory for fargo advection.")
420  END IF
421  CLASS DEFAULT
422  ! geometry not supported -> disable fargo
423  mesh%FARGO = 0
424  CALL this%Warning("InitTimedisc", &
425  "fargo has been disabled, because the geometry is not supported.")
426  END SELECT
427  CLASS DEFAULT
428  ! geometry not supported -> disable fargo
429  mesh%FARGO = 0
430  CALL this%Warning("InitTimedisc","fargo has been disabled, because the physics is not supported.")
431  END SELECT
432  ! initialize background velocity field w
433  SELECT CASE(mesh%FARGO)
434  CASE(1,2) ! set to 0;
435  ! fargo advection type 1: w is computed in each time step (see FargoCalcVelocity)
436  ! fargo advection type 2: w is provided by the user, e.g. in InitData
437  this%w(:,:) = 0.0
438  this%shift(:,:) = 0.0
439  this%dq = marray_base()
440  this%dq%data1d(:) = 0.0
441  this%dql = marray_base()
442  this%dql%data1d(:) = 0.0
443  this%flux = marray_base()
444  this%flux%data1d(:) = 0.0
445  CASE(3) ! fixed background velocity in shearing box
446  IF(mesh%shear_dir.EQ.2) THEN
447  this%w(:,:) = -mesh%Q*mesh%omega*mesh%bcenter(:,mesh%JMIN,:,1) !-Q*Omega*x
448  this%shift(:,:) = 0.0
449  shear_direction = "west<->east"
450  ELSE IF(mesh%shear_dir.EQ.1) THEN
451  this%w(:,:) = mesh%Q*mesh%omega*mesh%bcenter(mesh%IMIN,:,:,2) !Q*Omega*y
452  this%shift(:,:) = 0.0
453  shear_direction = "south<->north"
454  END IF
455  this%dq = marray_base()
456  this%dq%data1d(:) = 0.0
457  this%dql = marray_base()
458  this%dql%data1d(:) = 0.0
459  this%flux = marray_base()
460  this%flux%data1d(:) = 0.0
461  END SELECT
462  ELSE IF (mesh%use_fargo.EQ.0) THEN
463  ! do nothing
464  ELSE
465  CALL this%Error("InitTimedisc","unknown fargo advection scheme")
466  END IF
467 
468  ! print some information
469  WRITE (order_str, '(I0)') this%GetOrder()
470  WRITE (cfl_str, '(F4.2)') this%GetCFL()
471  WRITE (stoptime_str, '(ES10.4)') this%stoptime
472  WRITE (dtmax_str, '(ES10.4)') this%dtmax
473  WRITE (beta_str, '(ES10.4)') this%beta
474  CALL this%Info(" TIMEDISC-> ODE solver: " //trim(this%GetName()))
475  CALL this%Info(" order: " //trim(order_str))
476  CALL this%Info(" CFL number: " //trim(cfl_str))
477  CALL this%Info(" dtmax: " //trim(dtmax_str))
478  CALL this%Info(" stoptime: " //trim(stoptime_str))
479  CALL this%Info(" beta: " //trim(beta_str))
480  ! adaptive step size control
481  IF (this%tol_rel.LT.1.0.AND.this%order.GT.1) THEN
482  ! create state vector to store the error
483  CALL physics%new_statevector(this%cerr,conservative)
484  this%cerr%data1d(:) = 0.
485  ! create selection for the internal region
486  WRITE (info_str,'(ES7.1)') this%tol_rel*100
487  CALL this%Info(" step size control: enabled")
488  CALL this%Info(" rel. precision: "//trim(info_str)//" %")
489  ELSE
490  WRITE (info_str,'(A)') "disabled"
491  END IF
492  IF (mesh%FARGO.NE.0) &
493  CALL this%Info(" fargo: " //trim(fargo_method(mesh%FARGO)))
494  IF(mesh%FARGO.EQ.3) &
495  CALL this%Info(" shear-direction: " //trim(shear_direction))
496  SELECT CASE(this%rhstype)
497  CASE(0)
498  ! special rhs disabled, print nothing
499  CASE(1)
500  IF (.NOT.ASSOCIATED(this%w)) THEN
501  ALLOCATE(this%w(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX),stat = err)
502  IF (err.NE.0) THEN
503  CALL this%Error("InitTimedisc", "Unable to allocate memory for special RHS time stepping.")
504  END IF
505  this%w(:,:) = 0.0
506  END IF
507  CALL this%Info(" special rhs: enabled")
508  CASE DEFAULT
509  CALL this%Error("InitTimedisc","unknown rhstype")
510  END SELECT
511  END SUBROUTINE inittimedisc
512 
513 
514  SUBROUTINE setoutput(this,Mesh,Physics,config,IO)
515  IMPLICIT NONE
516  !-----------------------------------------------------------------------!
517  CLASS(timedisc_base), INTENT(INOUT) :: this
518  CLASS(mesh_base), INTENT(IN) :: Mesh
519  CLASS(physics_base), INTENT(IN) :: Physics
520  TYPE(Dict_TYP), POINTER :: config,IO
521  !------------------------------------------------------------------------!
522  INTEGER :: valwrite,i
523  CHARACTER(LEN=128) :: key,pvar_key
524  LOGICAL :: writeSolution
525  !------------------------------------------------------------------------!
526  valwrite = 0
527  CALL getattr(config, "output/error", valwrite, 0)
528  IF((valwrite.EQ.1).AND.this%tol_rel.LE.1.0) THEN
529  this%write_error = .true.
530  CALL physics%new_statevector(this%cerr_max,conservative)
531  this%cerr_max%data1d(:) = 0.
532  ELSE
533  this%write_error = .false.
534  END IF
535 
536  valwrite = 0
537  CALL getattr(config, "output/solution", valwrite, 0)
538  IF(valwrite.EQ.1) THEN
539  writesolution = .true.
540  CALL physics%new_statevector(this%solution,primitive)
541  ELSE
542  NULLIFY(this%solution)
543  writesolution = .false.
544  END IF
545 
546  CALL getattr(config, "output/time", valwrite, 1)
547  IF(valwrite.EQ.1) &
548  CALL setattr(io, "time", ref(this%time))
549 
550  DO i=1, physics%VNUM
551  !prim
552  key = trim(physics%pvarname(i))
553  pvar_key = key
554  valwrite = 0
555  CALL getattr(config, "output/" // trim(key), valwrite, 1)
556  ! second argument is important if pvarname is used twice
557  IF (valwrite.EQ.1) THEN
558  CALL setattr(io, trim(key), &
559  this%pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
560  END IF
561 
562  IF(writesolution) THEN
563  CALL setattr(io, trim(key)//"_solution", &
564  this%solution%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
565  END IF
566 
567  !cons
568  key = trim(physics%cvarname(i))
569  IF(key.EQ.pvar_key) key = trim(key) // "_cvar"
570  valwrite = 0
571  CALL getattr(config, "output/" // trim(key), valwrite, 0)
572  ! second argument is important if pvarname is used twice
573  IF (valwrite.EQ.1) THEN
574  CALL setattr(io, trim(key), &
575  this%cvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
576  END IF
577 
578  key = trim(physics%cvarname(i))
579  IF(this%write_error) THEN
580  CALL setattr(io, trim(key) // "_error", &
581  this%cerr_max%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
582  END IF
583 
584  ! write geometrical sources
585  CALL getattr(config, "output/" // "geometrical_sources", valwrite, 0)
586  IF (valwrite.EQ.1) THEN
587  CALL setattr(io, trim(key)//"_geo_src", &
588  this%geo_src%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
589  END IF
590 
591  ! write external sources
592  CALL getattr(config, "output/" // "external_sources", valwrite, 0)
593  IF (valwrite.EQ.1) THEN
594  CALL setattr(io, trim(key)//"_src", &
595  this%src%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
596  END IF
597 
598  ! write right hand side
599  CALL getattr(config, "output/" // "rhs", valwrite, 0)
600  IF (valwrite.EQ.1) THEN
601  CALL setattr(io, trim(key)//"_rhs", &
602  this%rhs%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
603  END IF
604 
605  ! write fluxes
606  ! ATTENTION: this are the numerical fluxes devided by dy or dx respectively
607  CALL getattr(config, "output/" // "fluxes", valwrite, 0)
608  IF (valwrite.EQ.1) THEN
609  CALL setattr(io, trim(key)//"_xfluxdydz", &
610  this%xfluxdydz(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
611  CALL setattr(io, trim(key)//"_yfluxdzdx", &
612  this%yfluxdzdx(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
613  CALL setattr(io, trim(key)//"_zfluxdxdy", &
614  this%zfluxdxdy(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,i))
615  END IF
616 
617  END DO
618  ! write bfluxes
619  CALL getattr(config, "output/" // "bflux", valwrite, 0)
620  IF(valwrite.EQ.1) THEN
621  ALLOCATE(this%bflux(physics%VNUM,6))
622  CALL setattr(io, "bflux", this%bflux)
623  ELSE
624  NULLIFY(this%bflux)
625  END IF
626  END SUBROUTINE setoutput
627 
628 
629  SUBROUTINE integrationstep(this,Mesh,Physics,Sources,Fluxes,iter,config,IO)
630  !------------------------------------------------------------------------!
631  CLASS(timedisc_base), INTENT(INOUT) :: this
632  CLASS(mesh_base), INTENT(INOUT) :: Mesh
633  CLASS(physics_base), INTENT(INOUT) :: Physics
634  CLASS(sources_base), POINTER :: Sources
635  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
636  INTEGER :: iter
637  TYPE(Dict_TYP), POINTER :: config,IO
638  !------------------------------------------------------------------------!
639  REAL :: err,dtold,dt,time
640  !------------------------------------------------------------------------!
641  INTENT(INOUT) :: iter
642  !------------------------------------------------------------------------!
643  ! transform to selenoidal velocities if fargo is enabled
644  IF (mesh%FARGO.GT.0) THEN
645  IF (mesh%shear_dir.EQ.1.AND.mesh%FARGO.EQ.3) THEN
646  CALL physics%SubtractBackgroundVelocityX(mesh,this%w,this%pvar,this%cvar)
647  ELSE
648  CALL physics%SubtractBackgroundVelocityY(mesh,this%w,this%pvar,this%cvar)
649  END IF
650  END IF
651 
652 
653  time = this%time
654  dt = this%dt
655  IF (dt.LT.this%dtmin .AND. this%dtcause .NE. dtcause_fileio) THEN
656  ! only save dtmin if the reasion is not the fileio
657  this%dtmin = dt
658  this%dtmincause = this%dtcause
659  END IF
660 !NEC$ NOVECTOR
661  timestep: DO WHILE (time+dt.LE.this%time+this%dt)
662  dtold = dt
663 
664  CALL this%SolveODE(mesh,physics,sources,fluxes,time,dt,err)
665  ! check truncation error and restart if necessary
666  IF (err.LT.1.0) THEN
667  CALL this%AcceptSolution(mesh,physics,sources,fluxes,time,dtold,iter)
668  ELSE
669  CALL this%RejectSolution(mesh,physics,sources,fluxes,time,dt)
670  END IF
671 
672  IF (dt.LT.this%dtlimit) &
673  this%break = .true.
674  ! Break if dt.LT.this%dtlimit or CheckData failed
675  IF(this%break) THEN
676  ! Do not attempt to fargo shift anymore
677  mesh%FARGO = 0
678  EXIT timestep
679  END IF
680  END DO timestep
681 
682  ! Save true advanced time (for fargo linear advection)
683  this%dt = time - this%time
684 
685  this%time = time
686  this%dtold = dt
687 
688  ! perform the fargo advection step if enabled
689  IF (mesh%FARGO.GT.0) THEN
690  IF (mesh%shear_dir.EQ.1) THEN
691  CALL this%FargoAdvectionX(fluxes,mesh,physics,sources)
692  ELSE
693  CALL this%FargoAdvectionY(fluxes,mesh,physics,sources)
694  END IF
695  END IF
696  END SUBROUTINE integrationstep
697 
698 
705  FUNCTION adjusttimestep(this,maxerr,dtold) RESULT(dtnew)
706  IMPLICIT NONE
707  !------------------------------------------------------------------------!
708  CLASS(timedisc_base) :: this
709  REAL :: maxerr,dtold,dtnew
710  !------------------------------------------------------------------------!
711  REAL :: dttmp
712  !------------------------------------------------------------------------!
713  INTENT(IN) :: maxerr,dtold
714  INTENT(INOUT) :: this
715  !------------------------------------------------------------------------!
716  ! adaptive step size control disabled
717  IF (this%tol_rel.GE.1.0) THEN
718  dtnew = huge(dtnew)
719  RETURN
720  END IF
721 
722  ! time step estimate with maxerr determined from the comparison of numerical
723  ! results of two explicit schemes with different order where GetOrder(this)
724  ! returns the order of the higher order scheme (P-Controller)
725 
726  IF(maxerr.GT.0.) THEN
727  dttmp = 0.9*dtold*exp(-log(maxerr)/getorder(this))
728  ELSE
729  ! ths limit for maxerr->0 is dttmp -> oo. But this cut off to
730  ! 4*dtold later in this function.
731  dttmp = 4.*dtold
732  END IF
733 
734  ! this controls the increase/decrease of time steps based on the
735  ! old value of maxerr from the previous time step (PI-Controller)
736  IF (this%maxerrold.GT.0.) THEN
737  dtnew = dttmp * exp(-this%beta*log(this%maxerrold/maxerr))
738  ELSE
739  dtnew = dttmp
740  END IF
741 
742  ! adjust the time step based on the error estimate
743  IF (maxerr.LT.1.0) THEN
744  ! increase time step
745  dtnew = min(dtnew,4.*dtold) ! enlarge at most by a factor of 4
746  ! it is possible: maxerr < 1 => 0.9*dt < dtnew => maybe dtnew < dt
747  IF (dtnew .LT. dtold) this%dtcause = dtcause_smallerr
748  ELSE
749  ! If maxerrold >> maxerr dtnew can become larger than dtold even
750  ! if the timestep is rejected. If this is the case fall back to
751  ! the simple P-Controller.
752  IF (dtnew.GT.dtold) dtnew = dttmp
753  ! decrease time step
754  dtnew = max(dtnew,0.25*dtold) ! reduce at most by a factor of 1/4
755  END IF
756 
757  ! store maxerr for next time step
758  this%maxerrold = maxerr
759 
760  END FUNCTION adjusttimestep
761 
762 
764  REAL FUNCTION calctimestep(this,Mesh,Physics,Sources,Fluxes,time,dtcause) RESULT(dt)
765  USE physics_euler_mod, ONLY : physics_euler, statevector_euler
766  IMPLICIT NONE
767  !------------------------------------------------------------------------!
768  CLASS(timedisc_base), INTENT(INOUT) :: this
769  CLASS(mesh_base), INTENT(IN) :: mesh
770  CLASS(physics_base), INTENT(INOUT) :: physics
771  CLASS(sources_base), POINTER :: sources
772  CLASS(fluxes_base), INTENT(INOUT) :: fluxes
773  REAL, INTENT(IN) :: time
774  INTEGER, INTENT(INOUT) :: dtcause
775  !------------------------------------------------------------------------!
776  REAL :: invdt
777  REAL :: dt_cfl, dt_src
778  REAL :: invdt_x,invdt_y,invdt_z
779  !------------------------------------------------------------------------!
780  ! CFL condition:
781  ! maximal wave speeds in each direction
782  IF (.NOT.this%always_update_bccsound) THEN
783  SELECT TYPE(phys => physics)
784  CLASS IS(physics_euler)
785  SELECT TYPE(pvar => this%pvar)
786  CLASS IS(statevector_euler)
787  CALL phys%UpdateSoundSpeed(pvar)
788  END SELECT
789  END SELECT
790  END IF
791 
792  SELECT TYPE(phys => physics)
793  TYPE IS(physics_eulerisotherm)
794  SELECT TYPE(pvar => this%pvar)
795  CLASS IS(statevector_eulerisotherm)
796  CALL phys%CalculateWaveSpeeds(mesh,pvar,fluxes%minwav,fluxes%maxwav)
797  END SELECT
798  TYPE IS(physics_euler)
799  SELECT TYPE(pvar => this%pvar)
800  CLASS IS(statevector_euler)
801  CALL phys%CalculateWaveSpeeds(mesh,pvar,fluxes%minwav,fluxes%maxwav)
802  END SELECT
803  END SELECT
804 
805  ! compute maximum of inverse time for CFL condition
806  IF ((mesh%JNUM.EQ.1).AND.(mesh%KNUM.EQ.1)) THEN
807  ! 1D, only x-direction
808  invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlx%data1d(:))
809  ELSE IF ((mesh%INUM.EQ.1).AND.(mesh%KNUM.EQ.1)) THEN
810  ! 1D, only y-direction
811  invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dly%data1d(:))
812  ELSE IF ((mesh%INUM.EQ.1).AND.(mesh%JNUM.EQ.1)) THEN
813  ! 1D, only z-direction
814  invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlz%data1d(:))
815  ELSE IF ((mesh%INUM.GT.1).AND.(mesh%JNUM.GT.1).AND.(mesh%KNUM.EQ.1)) THEN
816  ! 2D, x-y-plane
817  invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlx%data1d(:) &
818  + max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) / mesh%dly%data1d(:))
819  ELSE IF ((mesh%INUM.GT.1).AND.(mesh%KNUM.GT.1).AND.(mesh%JNUM.EQ.1)) THEN
820  ! 2D, x-z-plane
821  invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlx%data1d(:) &
822  + max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) / mesh%dlz%data1d(:))
823  ELSE IF ((mesh%JNUM.GT.1).AND.(mesh%KNUM.GT.1).AND.(mesh%INUM.EQ.1)) THEN
824  ! 2D, y-z-plane
825  invdt = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dly%data1d(:) &
826  + max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) / mesh%dlz%data1d(:))
827  ELSE
828  ! full 3D
829  !TODO: Achtung: Hier wurde fuer eine bessere Symmetrie fuer jede Richtung ein eigenes invdt
830  ! berechnet. Dies koennte jedoch einen Verlust an Stabilitaet bewirken. Hier muesste mal eine
831  ! Stabilitaetsanalyse gemacht werden
832  invdt_x = maxval(max(fluxes%maxwav%data2d(:,1),-fluxes%minwav%data2d(:,1)) / mesh%dlx%data1d(:))
833  invdt_y = maxval(max(fluxes%maxwav%data2d(:,2),-fluxes%minwav%data2d(:,2)) / mesh%dly%data1d(:))
834  invdt_z = maxval(max(fluxes%maxwav%data2d(:,3),-fluxes%minwav%data2d(:,3)) / mesh%dlz%data1d(:))
835  invdt = max(invdt_y,invdt_z,invdt_x)
836  ! invdt = MAXVAL(MAX(Fluxes%maxwav(:,:,:,3),-Fluxes%minwav(:,:,:,3)) / Mesh%dlz(:,:,:) &
837  ! + MAX(Fluxes%maxwav(:,:,:,2),-Fluxes%minwav(:,:,:,2)) / Mesh%dly(:,:,:) &
838  ! + MAX(Fluxes%maxwav(:,:,:,1),-Fluxes%minwav(:,:,:,1)) / Mesh%dlx(:,:,:))
839  END IF
840 
841  ! largest time step due to CFL condition
842  dt_cfl = this%cfl / invdt
843  ! time step limitation due to cfl -> dtcause = 0
844  dtcause = 0
845 
846  ! check for sources
847  IF (ASSOCIATED(sources)) THEN
848  ! initialize this to be sure dt_src > 0
849  dt_src = dt_cfl
850  CALL sources%CalcTimestep(mesh,physics,fluxes,this%pvar,this%cvar,time,dt_src,dtcause)
851  dt = min(dt_cfl,dt_src)
852  ELSE
853  dt = dt_cfl
854  END IF
855  END FUNCTION calctimestep
856 
857  SUBROUTINE acceptsolution(this,Mesh,Physics,Sources,Fluxes,time,dt,iter)
858  IMPLICIT NONE
859  !------------------------------------------------------------------------!
860  CLASS(timedisc_base), INTENT(INOUT) :: this
861  CLASS(mesh_base), INTENT(IN) :: Mesh
862  CLASS(physics_base), INTENT(INOUT) :: Physics
863  CLASS(sources_base), POINTER :: Sources
864  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
865  REAL :: time,dt
866  INTEGER :: iter
867  !------------------------------------------------------------------------!
868  REAL :: dtmeanold
869  REAL, DIMENSION(:), POINTER :: bflux,bfold
870  !------------------------------------------------------------------------!
871  INTENT(INOUT) :: time,dt
872  !------------------------------------------------------------------------!
873  time=time+dt
874  this%dtaccept = this%dtaccept + 1
875  dtmeanold = this%dtmean
876  this%dtmean = this%dtmean + (dt - this%dtmean)/this%dtaccept
877  this%dtstddev = this%dtstddev + (dt - dtmeanold)*(dt-this%dtmean)
878  CALL this%ComputeRHS(mesh,physics,sources,fluxes,time,dt,this%pvar,&
879  this%cvar,this%checkdatabm,this%rhs)
880  this%cold%data1d(:) = this%cvar%data1d(:)
881  IF (mesh%INUM.GT.1) THEN
882  ! collapse the 4D arrays to improve vector performance;
883  ! maybe we can switch to the old style assignments if
884  ! automatic collapsing of NEC nfort compiler works
885  bflux(1:SIZE(fluxes%bxflux)) => fluxes%bxflux
886  bfold(1:SIZE(fluxes%bxfold)) => fluxes%bxfold
887  bfold(:) = bflux(:)
888 ! Fluxes%bxfold(:,:,:,:) = Fluxes%bxflux(:,:,:,:)
889  END IF
890  IF (mesh%JNUM.GT.1) THEN
891  bflux(1:SIZE(fluxes%byflux)) => fluxes%byflux
892  bfold(1:SIZE(fluxes%byfold)) => fluxes%byfold
893  bfold(:) = bflux(:)
894 ! Fluxes%byfold(:,:,:,:) = Fluxes%byflux(:,:,:,:)
895  END IF
896  IF (mesh%KNUM.GT.1) THEN
897  bflux(1:SIZE(fluxes%bzflux)) => fluxes%bzflux
898  bfold(1:SIZE(fluxes%bzfold)) => fluxes%bzfold
899  bfold(:) = bflux(:)
900 ! Fluxes%bzfold(:,:,:,:) = Fluxes%bzflux(:,:,:,:)
901  END IF
902  iter = iter + 1
903  END SUBROUTINE acceptsolution
904 
905  SUBROUTINE rejectsolution(this,Mesh,Physics,Sources,Fluxes,time,dt)
906  IMPLICIT NONE
907  !------------------------------------------------------------------------!
908  CLASS(timedisc_base), INTENT(INOUT) :: this
909  CLASS(mesh_base), INTENT(IN) :: Mesh
910  CLASS(physics_base), INTENT(INOUT) :: Physics
911  CLASS(sources_base), POINTER :: Sources
912  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
913  REAL :: time,dt
914  !------------------------------------------------------------------------!
915  INTENT(IN) :: time
916  INTENT(INOUT) :: dt
917  !------------------------------------------------------------------------!
918  this%cvar%data1d(:) = this%cold%data1d(:)
919  ! This data has already been checked at AcceptSolution
920  CALL this%ComputeRHS(mesh,physics,sources,fluxes,time,dt,this%pvar, &
921  this%cvar,check_nothing,this%rhs)
922  fluxes%bxflux(:,:,:,:) = fluxes%bxfold(:,:,:,:)
923  fluxes%byflux(:,:,:,:) = fluxes%byfold(:,:,:,:)
924  fluxes%bzflux(:,:,:,:) = fluxes%bzfold(:,:,:,:)
925  ! count adjustments for information
926  this%n_adj = this%n_adj + 1
927  this%dtcause = dtcause_erradj
928  ! only save dtmin if the reasion is not the fileio (automatically satisfied)
929  IF (dt.LT.this%dtmin) THEN
930  this%dtmin = dt
931  this%dtmincause = this%dtcause
932  END IF
933  END SUBROUTINE rejectsolution
934 
935 
936  FUNCTION computeerror(this,Mesh,Physics,cvar_high,cvar_low) RESULT(maxerr)
937  IMPLICIT NONE
938  !------------------------------------------------------------------------!
939  CLASS(timedisc_base), INTENT(INOUT) :: this
940  CLASS(mesh_base), INTENT(IN) :: mesh
941  CLASS(physics_base), INTENT(IN) :: physics
942  CLASS(marray_compound),INTENT(INOUT):: cvar_high,cvar_low
943  REAL :: maxerr
944  !------------------------------------------------------------------------!
945  INTEGER :: l
946 #ifdef PARALLEL
947  INTEGER :: ierror
948 #endif
949  REAL :: rel_err(physics%vnum)
950  !------------------------------------------------------------------------!
951 !NEC$ SHORTLOOP
952  DO l=1,physics%VNUM
953  ! compute the local error (including ghost zones)
954  this%cerr%data2d(:,l) = abs(cvar_high%data2d(:,l)-cvar_low%data2d(:,l)) &
955  / (this%tol_rel*abs(cvar_high%data2d(:,l)) + this%tol_abs(l))
956  ! determine the global maximum on the whole grid except for ghost zones
957  rel_err(l) = maxval(this%cerr%data2d(:,l),mask=mesh%without_ghost_zones%mask1d(:))
958  ! store the maximum between two output time steps
959  IF (this%write_error) &
960  this%cerr_max%data2d(:,l) = max(this%cerr_max%data2d(:,l),this%cerr%data2d(:,l))
961  END DO
962 
963  ! compute the maximum of all variables
964  maxerr = maxval(rel_err(:))
965 
966 #ifdef PARALLEL
967  ! find the global maximum of all processes
968  CALL mpi_allreduce(mpi_in_place,maxerr,1,default_mpi_real,mpi_max,&
969  mesh%comm_cart,ierror)
970 #endif
971 
972  END FUNCTION computeerror
973 
974 
975 
992  SUBROUTINE computerhs(this,Mesh,Physics,Sources,Fluxes,time,dt,pvar,cvar,checkdatabm,rhs)
993  USE physics_euler_mod, ONLY : physics_euler,statevector_euler
994  USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
995  IMPLICIT NONE
996  !------------------------------------------------------------------------!
997  CLASS(timedisc_base), INTENT(INOUT) :: this
998  CLASS(mesh_base), INTENT(IN) :: Mesh
999  CLASS(physics_base), INTENT(INOUT) :: Physics
1000  CLASS(sources_base), POINTER :: Sources
1001  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
1002  REAL, INTENT(IN) :: time, dt
1003  CLASS(marray_compound),INTENT(INOUT):: pvar,cvar,rhs
1004  INTEGER, INTENT(IN) :: checkdatabm
1005  !------------------------------------------------------------------------!
1006  INTEGER :: i,j,k,l
1007  LOGICAL :: have_potential
1008  REAL :: t
1009  REAL :: phi,wp
1010  REAL, DIMENSION(:,:,:,:), POINTER :: pot
1011  CLASS(sources_base), POINTER :: sp
1012  CLASS(sources_gravity), POINTER :: grav
1013  !------------------------------------------------------------------------!
1014  t = time
1015  SELECT CASE(mesh%FARGO)
1016  CASE(1,2)
1017  ! transform to real velocity, i.e. v_residual + w_background,
1018  ! before setting the boundary conditions
1019  CALL physics%AddBackgroundVelocityY(mesh,this%w,pvar,cvar)
1020  CASE(3)
1021  ! boundary conditions are set for residual velocity in shearing box simulations
1022  ! usually it's not necessary to subtract the background velocity here, but
1023  ! in case someone adds it before, we subtract it here; the subroutine checks,
1024  ! if the velocities have already been transformed
1025  IF (mesh%shear_dir.EQ.1) THEN
1026  CALL physics%SubtractBackgroundVelocityX(mesh,this%w,pvar,cvar)
1027  ELSE IF (mesh%shear_dir.EQ.2) THEN
1028  CALL physics%SubtractBackgroundVelocityY(mesh,this%w,pvar,cvar)
1029  END IF
1030  ! ATTENTION: the time must be the initial time of the whole time step
1031  ! not the time of a substep
1032  t = this%time
1033  CASE DEFAULT
1034  ! fargo disabled (do nothing)
1035  END SELECT
1036 
1037  ! set boundary values and convert conservative to primitive variables
1038  CALL this%boundary%CenterBoundary(mesh,physics,t,pvar,cvar)
1039 
1040  IF(iand(checkdatabm,check_tmin).NE.check_nothing.AND.&
1041  this%tmin.GT.1.e-10) THEN
1042  SELECT TYPE(p => pvar)
1043  CLASS IS(statevector_euler)
1044  SELECT TYPE(c => cvar)
1045  CLASS IS(statevector_euler)
1046  ! If temperature is below TMIN limit pressure to density*RG/MU*TMIN.
1047  p%pressure%data1d(:) = max(p%pressure%data1d(:), &
1048  p%density%data1d(:)*physics%Constants%RG/physics%MU*this%TMIN)
1049  CALL physics%Convert2Conservative(p,c)
1050  END SELECT
1051  END SELECT
1052  END IF
1053 
1054  IF(iand(checkdatabm,check_pmin).NE.check_nothing.AND.&
1055  physics%PRESSURE.GT.0) THEN
1056  ! Check if the pressure is below pmin. If it is, increase the pressure
1057  ! to reach pmin
1058  SELECT TYPE(p => pvar)
1059  CLASS IS(statevector_euler)
1060  SELECT TYPE(c => cvar)
1061  CLASS IS(statevector_euler)
1062  ! If temperature is below PMIN limit pressure to PMIN
1063  p%pressure%data1d(:) &
1064  = max(p%pressure%data1d(:),this%pmin)
1065  CALL physics%Convert2Conservative(p,c)
1066  END SELECT
1067  END SELECT
1068  END IF
1069 
1070  ! update the speed of sound (non-isotherml physics only)
1071  IF (this%always_update_bccsound) THEN
1072  SELECT TYPE(phys => physics)
1073  CLASS IS(physics_euler)
1074  SELECT TYPE(pvar => this%pvar)
1075  CLASS IS(statevector_euler)
1076  CALL phys%UpdateSoundSpeed(pvar)
1077  END SELECT
1078  END SELECT
1079  END IF
1080 
1081  ! check for illegal data
1082  IF(checkdatabm.NE.check_nothing) &
1083  CALL this%CheckData(mesh,physics,fluxes,pvar,cvar,checkdatabm)
1084 
1085  ! get geometrical sources
1086  CALL physics%GeometricalSources(mesh,pvar,cvar,this%geo_src)
1087 
1088  ! get source terms due to external forces if present
1089  IF (ASSOCIATED(sources)) &
1090  CALL sources%ExternalSources(mesh,fluxes,physics, &
1091  t,dt,pvar,cvar,this%src)
1092 
1093  ! if fargo advection is enabled additional source terms occur;
1094  ! furthermore computation of numerical fluxes should always be
1095  ! carried out with residual velocity
1096  SELECT CASE(mesh%FARGO)
1097  CASE(1,2)
1098  ! add fargo source terms to geometrical source terms
1099  CALL physics%AddFargoSources(mesh,this%w,pvar,cvar,this%geo_src)
1100  ! subtract background velocity
1101  CALL physics%SubtractBackgroundVelocityY(mesh,this%w,pvar,cvar)
1102  CASE(3)
1103  ! background velocity field has already been subtracted (do nothing);
1104  ! fargo specific source terms are handled in the shearing box source
1105  ! term module
1106  CASE DEFAULT
1107  ! fargo disabled (do nothing)
1108  END SELECT
1109 
1110  ! get the numerical fluxes
1111  CALL fluxes%CalculateFluxes(mesh,physics,pvar,cvar,this%xfluxdydz, &
1112  this%yfluxdzdx,this%zfluxdxdy)
1113 
1114  ! compute the right hand side for boundary flux computation;
1115  ! this is probably wrong for the special rhs (rhstype=1, see above)
1116 !NEC$ SHORTLOOP
1117  DO l=1,physics%VNUM+physics%PNUM
1118 !NEC$ IVDEP
1119  DO k=mesh%KMIN,mesh%KMAX
1120 !NEC$ IVDEP
1121  DO j=mesh%JMIN,mesh%JMAX
1122  ! western and eastern boundary fluxes
1123  rhs%data4d(mesh%IMIN-mesh%Ip1,j,k,l) = mesh%dy*mesh%dz * this%xfluxdydz(mesh%IMIN-mesh%Ip1,j,k,l)
1124  rhs%data4d(mesh%IMAX+mesh%Ip1,j,k,l) = -mesh%dy*mesh%dz * this%xfluxdydz(mesh%IMAX,j,k,l)
1125  END DO
1126  END DO
1127 !NEC$ IVDEP
1128  DO k=mesh%KMIN,mesh%KMAX
1129 !NEC$ IVDEP
1130  DO i=mesh%IMIN,mesh%IMAX
1131  ! southern and northern boundary fluxes
1132  rhs%data4d(i,mesh%JMIN-mesh%Jp1,k,l) = mesh%dz*mesh%dx * this%yfluxdzdx(i,mesh%JMIN-mesh%Jp1,k,l)
1133  rhs%data4d(i,mesh%JMAX+mesh%Jp1,k,l) = -mesh%dz*mesh%dx * this%yfluxdzdx(i,mesh%JMAX,k,l)
1134  END DO
1135  END DO
1136 !NEC$ IVDEP
1137  DO j=mesh%JMIN,mesh%JMAX
1138 !NEC$ IVDEP
1139  DO i=mesh%IMIN,mesh%IMAX
1140  ! bottomer and upper boundary fluxes
1141  rhs%data4d(i,j,mesh%KMIN-mesh%Kp1,l) = mesh%dx*mesh%dy * this%zfluxdxdy(i,j,mesh%KMIN-mesh%Kp1,l)
1142  rhs%data4d(i,j,mesh%KMAX+mesh%Kp1,l) = -mesh%dx*mesh%dy * this%zfluxdxdy(i,j,mesh%KMAX,l)
1143  END DO
1144  END DO
1145  END DO
1146 
1147  SELECT CASE(this%rhstype)
1148  CASE(0)
1149 !NEC$ SHORTLOOP
1150  DO l=1,physics%VNUM+physics%PNUM
1151  DO k=mesh%KMIN,mesh%KMAX
1152  DO j=mesh%JMIN,mesh%JMAX
1153 !NEC$ IVDEP
1154  DO i=mesh%IMIN,mesh%IMAX
1155  ! update right hand side of ODE
1156  rhs%data4d(i,j,k,l) = &
1157  mesh%dydzdV%data3d(i,j,k)*(this%xfluxdydz(i,j,k,l) - this%xfluxdydz(i-mesh%Ip1,j,k,l)) &
1158  + mesh%dzdxdV%data3d(i,j,k)*(this%yfluxdzdx(i,j,k,l) - this%yfluxdzdx(i,j-mesh%Jp1,k,l)) &
1159  + mesh%dxdydV%data3d(i,j,k)*(this%zfluxdxdy(i,j,k,l) - this%zfluxdxdy(i,j,k-mesh%Kp1,l)) &
1160  - this%geo_src%data4d(i,j,k,l) - this%src%data4d(i,j,k,l)
1161  END DO
1162  END DO
1163  END DO
1164  END DO
1165  CASE(1)
1169  sp => sources
1170  DO
1171  IF (ASSOCIATED(sp).EQV..false.) RETURN
1172  SELECT TYPE(sp)
1173  CLASS IS(sources_gravity)
1174  grav => sp
1175  EXIT
1176  END SELECT
1177  sp => sp%next
1178  END DO
1179 
1180  NULLIFY(pot)
1181  have_potential = .false.
1182  IF(ASSOCIATED(grav)) THEN
1183  IF(.NOT.grav%addtoenergy) THEN
1184  pot => mesh%RemapBounds(grav%pot%data4d(:,:,:,:))
1185  have_potential=.true.
1186  END IF
1187  END IF
1188 
1189  phi = 0.
1190 
1191  DO k=mesh%KMIN-mesh%KP1,mesh%KMAX
1192  DO j=mesh%JMIN-mesh%JP1,mesh%JMAX
1193 !NEC$ IVDEP
1194  DO i=mesh%IMIN-mesh%IP1,mesh%IMAX
1195  wp = 0.5*(this%w(i,k)+this%w(i+1,k)) + mesh%hy%faces(i+1,j,k,1)*mesh%OMEGA
1196  IF(physics%PRESSURE.GT.0) THEN
1197  this%xfluxdydz(i,j,k,physics%ENERGY) &
1198  = this%xfluxdydz(i,j,k,physics%ENERGY) &
1199  + wp * (0.5 * wp * this%xfluxdydz(i,j,k,physics%DENSITY) &
1200  + this%xfluxdydz(i,j,k,physics%YMOMENTUM))
1201  IF (have_potential) THEN
1202  this%xfluxdydz(i,j,k,physics%ENERGY) = this%xfluxdydz(i,j,k,physics%ENERGY) + pot(i,j,k,2)*this%xfluxdydz(i,j,k,physics%DENSITY)
1203  END IF
1204  END IF
1205 
1206  this%xfluxdydz(i,j,k,physics%YMOMENTUM) &
1207  = (this%xfluxdydz(i,j,k,physics%YMOMENTUM) &
1208  + wp * this%xfluxdydz(i,j,k,physics%DENSITY)) * mesh%hy%faces(i+1,j,k,1)
1209 
1210  wp = this%w(i,k) + mesh%radius%bcenter(i,j,k)*mesh%OMEGA
1211  IF(physics%PRESSURE.GT.0) THEN
1212  this%yfluxdzdx(i,j,k,physics%ENERGY) &
1213  = this%yfluxdzdx(i,j,k,physics%ENERGY) &
1214  + wp * ( 0.5 * wp * this%yfluxdzdx(i,j,k,physics%DENSITY) &
1215  + this%yfluxdzdx(i,j,k,physics%YMOMENTUM))
1216  IF (have_potential) THEN
1217  this%yfluxdzdx(i,j,k,physics%ENERGY) = this%yfluxdzdx(i,j,k,physics%ENERGY) + pot(i,j,k,3)*this%yfluxdzdx(i,j,k,physics%DENSITY)
1218  END IF
1219  END IF
1220 
1221  this%yfluxdzdx(i,j,k,physics%YMOMENTUM) &
1222  = this%yfluxdzdx(i,j,k,physics%YMOMENTUM) &
1223  + wp * this%yfluxdzdx(i,j,k,physics%DENSITY)
1224  END DO
1225  END DO
1226  END DO
1227 
1228  ! set isothermal sound speeds
1229  SELECT TYPE (phys => physics)
1230  CLASS IS (physics_eulerisotherm)
1231  DO k=mesh%KMIN,mesh%KMAX
1232  DO j=mesh%JMIN,mesh%JMAX
1233 !NEC$ IVDEP
1234  DO i=mesh%IMIN,mesh%IMAX
1235  wp = pvar%data4d(i,j,k,physics%YVELOCITY) + this%w(i,k) + mesh%radius%bcenter(i,j,k)*mesh%OMEGA
1236  this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1237  = pvar%data4d(i,j,k,physics%DENSITY) * wp**2 * mesh%cyxy%bcenter(i,j,k) &
1238  + pvar%data4d(i,j,k,physics%DENSITY) * pvar%data4d(i,j,k,physics%XVELOCITY) * wp * mesh%cxyx%bcenter(i,j,k)
1239 
1240  this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1241  = cvar%data4d(i,j,k,physics%XMOMENTUM) &
1242  * ( pvar%data4d(i,j,k,physics%XVELOCITY) * mesh%cxyx%center(i,j,k) )
1243 
1244  IF(physics%PRESSURE.GT.0) THEN
1245  this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1246  = this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1247  +pvar%data4d(i,j,k,physics%PRESSURE) &
1248  *( mesh%cyxy%center(i,j,k) + mesh%czxz%center(i,j,k) )
1249  this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1250  = this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1251  + pvar%data4d(i,j,k,physics%PRESSURE) &
1252  *( mesh%cxyx%center(i,j,k) + mesh%czyz%center(i,j,k) )
1253  this%geo_src%data4d(i,j,k,physics%ENERGY) = 0.
1254  ELSE
1255 
1256  this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1257  = this%geo_src%data4d(i,j,k,physics%XMOMENTUM) &
1258  +pvar%data4d(i,j,k,physics%DENSITY)*phys%bccsound%data3d(i,j,k)**2 &
1259  * ( mesh%cyxy%center(i,j,k) + mesh%czxz%center(i,j,k) )
1260  this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1261  = this%geo_src%data4d(i,j,k,physics%YMOMENTUM) &
1262  + pvar%data4d(i,j,k,physics%DENSITY)*phys%bccsound%data3d(i,j,k)**2 &
1263  *( mesh%cxyx%center(i,j,k) + mesh%czyz%center(i,j,k) )
1264  END IF
1265  END DO
1266  END DO
1267  END DO
1268  CLASS DEFAULT
1269  ! abort
1270  END SELECT
1271 
1272 
1273 !NEC$ NOVECTOR
1274  DO l=1,physics%VNUM+physics%PNUM
1275 !NEC$ OUTERLOOP_UNROLL(8)
1276  DO k=mesh%KMIN,mesh%KMAX
1277  DO j=mesh%JMIN,mesh%JMAX
1278 !NEC$ IVDEP
1279  DO i=mesh%IMIN,mesh%IMAX
1280  ! update right hand side of ODE
1281  IF(l.EQ.physics%YMOMENTUM) THEN
1282  rhs%data4d(i,j,k,l) = mesh%dydzdV%data3d(i,j,k)*(this%xfluxdydz(i,j,k,l) - this%xfluxdydz(i-1,j,k,l)) &
1283  / mesh%hy%center(i,j,k)
1284  ELSE
1285  rhs%data4d(i,j,k,l) = mesh%dydzdV%data3d(i,j,k)*(this%xfluxdydz(i,j,k,l) - this%xfluxdydz(i-1,j,k,l))
1286  END IF
1287  ! time step update
1288  rhs%data4d(i,j,k,l) = rhs%data4d(i,j,k,l) &
1289  + mesh%dzdxdV%data3d(i,j,k)*(this%yfluxdzdx(i,j,k,l) - this%yfluxdzdx(i,j-1,k,l))
1290  END DO
1291  END DO
1292  END DO
1293  END DO
1294 
1295  DO k=mesh%KMIN,mesh%KMAX
1296  DO j=mesh%JMIN,mesh%JMAX
1297 !NEC$ IVDEP
1298  DO i=mesh%IMIN,mesh%IMAX
1299  wp = this%w(i,k) + mesh%radius%bcenter(i,j,k) * mesh%OMEGA
1300  rhs%data4d(i,j,k,physics%YMOMENTUM) = rhs%data4d(i,j,k,physics%YMOMENTUM) &
1301  - wp * rhs%data4d(i,j,k,physics%DENSITY)
1302  IF(physics%PRESSURE.GT.0) THEN
1303  rhs%data4d(i,j,k,physics%ENERGY) = rhs%data4d(i,j,k,physics%ENERGY) &
1304  - wp*( rhs%data4d(i,j,k,physics%YMOMENTUM) &
1305  + 0.5 * wp* rhs%data4d(i,j,k,physics%DENSITY))
1306  IF (have_potential) THEN
1307  rhs%data4d(i,j,k,physics%ENERGY) = &
1308  rhs%data4d(i,j,k,physics%ENERGY) - pot(i,j,k,1) * rhs%data4d(i,j,k,physics%DENSITY)
1309  END IF
1310  END IF
1311  END DO
1312  END DO
1313  END DO
1314 
1315 !NEC$ NOVECTOR
1316  DO l=1,physics%VNUM+physics%PNUM
1317 !NEC$ OUTERLOOP_UNROLL(8)
1318  DO k=mesh%KMIN,mesh%KMAX
1319  DO j=mesh%JMIN,mesh%JMAX
1320 !NEC$ IVDEP
1321  DO i=mesh%IMIN,mesh%IMAX
1322  ! update right hand side of ODE
1323  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)
1324  END DO
1325  END DO
1326  END DO
1327  END DO
1328  END SELECT
1329  END SUBROUTINE computerhs
1330 
1333  SUBROUTINE checkdata(this,Mesh,Physics,Fluxes,pvar,cvar,checkdatabm)
1334  USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
1335  USE physics_euler_mod, ONLY : physics_euler
1336  IMPLICIT NONE
1337  !------------------------------------------------------------------------!
1338  CLASS(timedisc_base), INTENT(INOUT) :: this
1339  CLASS(mesh_base), INTENT(IN) :: Mesh
1340  CLASS(physics_base), INTENT(INOUT) :: Physics
1341  CLASS(fluxes_base), INTENT(IN) :: Fluxes
1342  CLASS(marray_compound),INTENT(INOUT):: pvar,cvar
1343  INTEGER, INTENT(IN) :: checkdatabm
1344  !------------------------------------------------------------------------!
1345 #ifdef PARALLEL
1346  INTEGER :: err
1347 #endif
1348  REAL :: val
1349  !------------------------------------------------------------------------!
1350  IF(iand(checkdatabm,check_csound).NE.check_nothing) THEN
1351  ! check for speed of sound
1352  SELECT TYPE(phys => physics)
1353  CLASS IS(physics_eulerisotherm)
1354  val = minval(phys%bccsound%data1d(:))
1355  IF(val.LE.0.) THEN
1356  ! warn now and stop after file output
1357  CALL this%Warning("CheckData","Illegal speed of sound value less than 0.")
1358  this%break = .true.
1359  END IF
1360  CLASS DEFAULT
1361  CALL this%Warning("CheckData","check speed of sound selected, but bccsound not defined")
1362  END SELECT
1363  END IF
1364  IF((iand(checkdatabm,check_pmin).NE.check_nothing)) THEN
1365  ! check for non-isothermal physics with pressure defined
1366  SELECT TYPE(p => pvar)
1367  CLASS IS(statevector_euler)
1368  val = minval(p%pressure%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX, &
1369  mesh%KMIN:mesh%KMAX))
1370  IF(val.LT.this%pmin) THEN
1371  ! warn now and stop after file output
1372  CALL this%Warning("CheckData","Pressure below allowed pmin value.")
1373  this%break = .true.
1374  END IF
1375  END SELECT
1376  END IF
1377 
1378  IF(iand(checkdatabm,check_rhomin).NE.check_nothing) THEN
1379  ! check for physics with density defined
1380  SELECT TYPE(p => pvar)
1381  CLASS IS(statevector_eulerisotherm)
1382  val = minval(p%density%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX, &
1383  mesh%KMIN:mesh%KMAX))
1384  IF(val.LT.this%rhomin) THEN
1385  ! warn now and stop after file output
1386  CALL this%Warning("CheckData","Density below allowed rhomin value.")
1387  this%break = .true.
1388  END IF
1389  END SELECT
1390  END IF
1391 
1392  IF(iand(checkdatabm,check_invalid).NE.check_nothing) THEN
1393  IF(any((cvar%data1d.NE.cvar%data1d).OR.(pvar%data1d.NE.pvar%data1d))) THEN
1394  ! warn now and stop after file output
1395  CALL this%Warning("CheckData","Found NaN in pvar or cvar.")
1396  this%break = .true.
1397  END IF
1398  IF(any((cvar%data1d.GT.huge(cvar%data1d)).OR.pvar%data1d.GT.huge(pvar%data1d))) THEN
1399  ! warn now and stop after file output
1400  CALL this%Warning("CheckData","Found Infinity in pvar or cvar.")
1401  this%break = .true.
1402  END IF
1403  END IF
1404 
1405 #ifdef PARALLEL
1406  CALL mpi_allreduce(mpi_in_place, this%break, 1, mpi_logical, mpi_lor, &
1407  mesh%comm_cart, err)
1408 #endif
1409  END SUBROUTINE checkdata
1410 
1411  PURE FUNCTION getorder(this) RESULT(odr)
1412  IMPLICIT NONE
1413  !------------------------------------------------------------------------!
1414  CLASS(timedisc_base), INTENT(IN) :: this
1415  INTEGER :: odr
1416  !------------------------------------------------------------------------!
1417  odr = this%order
1418  END FUNCTION getorder
1419 
1420  PURE FUNCTION getcfl(this) RESULT(cfl)
1421  IMPLICIT NONE
1422  !------------------------------------------------------------------------!
1423  CLASS(timedisc_base), INTENT(IN) :: this
1424  REAL :: cfl
1425  !------------------------------------------------------------------------!
1426  cfl = this%CFL
1427  END FUNCTION getcfl
1428 
1429 
1448  SUBROUTINE fargoadvectionx(this,Fluxes,Mesh,Physics,Sources)
1449  IMPLICIT NONE
1450  !------------------------------------------------------------------------!
1451  CLASS(timedisc_base), INTENT(INOUT) :: this
1452  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
1453  CLASS(mesh_base), INTENT(IN) :: Mesh
1454  CLASS(physics_base), INTENT(INOUT) :: Physics
1455  CLASS(sources_base), POINTER :: Sources
1456  !------------------------------------------------------------------------!
1457  INTEGER :: i,j,k,l
1458 #ifdef PARALLEL
1459  CHARACTER(LEN=80) :: str
1460  INTEGER :: status(MPI_STATUS_SIZE)
1461  INTEGER :: ierror
1462  REAL :: mpi_buf(2*Mesh%GNUM)
1463  REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX) :: buf
1464 #endif
1465  !------------------------------------------------------------------------!
1466  ! determine step size of integer shift and length of remaining transport step
1467  ! first compute the whole step
1468  this%delxy(:,:) = this%w(:,:) * this%dt / mesh%dlx%data3d(mesh%IMIN,:,:)
1469 
1470 #ifdef PARALLEL
1471  ! make sure all MPI processes use the same step if domain is decomposed
1472  ! along the x- or y-direction (can be different due to round-off errors)
1473  IF (mesh%dims(1).GT.1) THEN
1474  CALL mpi_allreduce(mpi_in_place,this%delxy,(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1), &
1475  default_mpi_real,mpi_min,mesh%Icomm,ierror)
1476  END IF
1477 #endif
1478 
1479  ! then subdivide into integer shift and remaining linear advection step
1480 !NEC$ COLLAPSE
1481  DO k = mesh%KGMIN,mesh%KGMAX
1482 !NEC$ IVDEP
1483  DO j = mesh%JGMIN,mesh%JGMAX
1484  this%shift(j,k) = nint(this%delxy(j,k))
1485  this%delxy(j,k) = this%delxy(j,k)-dble(this%shift(j,k))
1486  END DO
1487  END DO
1488 
1489 !NEC$ SHORTLOOP
1490  DO l=1,physics%VNUM+physics%PNUM
1491  this%dq%data3d(mesh%IGMIN,:,:) = this%cvar%data4d(mesh%IGMIN+1,:,:,l) - this%cvar%data4d(mesh%IGMIN,:,:,l)
1492 !NEC$ IVDEP
1493  DO i=mesh%IGMIN+1,mesh%IGMAX-1
1494  this%dq%data3d(i,:,:) = this%cvar%data4d(i+1,:,:,l)-this%cvar%data4d(i,:,:,l)
1495  ! apply minmod limiter
1496  WHERE (sign(1.0,this%dq%data3d(i-1,:,:))*sign(1.0,this%dq%data3d(i,:,:)).GT.0)
1497  this%dql%data3d(i,:,:) = sign(min(abs(this%dq%data3d(i-1,:,:)),abs(this%dq%data3d(i,:,:))),this%dq%data3d(i-1,:,:))
1498  ELSEWHERE
1499  this%dql%data3d(i,:,:) = 0.
1500  END WHERE
1501  END DO
1502 !NEC$ IVDEP
1503  DO i=mesh%IMIN-1,mesh%IMAX
1504  WHERE(this%delxy(:,:).GT.0.)
1505  this%flux%data3d(i,:,:) = this%cvar%data4d(i,:,:,l) + .5 * this%dql%data3d(i,:,:) * (1. - this%delxy(:,:))
1506  ELSEWHERE
1507  this%flux%data3d(i,:,:) = this%cvar%data4d(i+1,:,:,l) - .5*this%dql%data3d(i+1,:,:)*(1. + this%delxy(:,:))
1508  END WHERE
1509  this%cvar%data4d(i,:,:,l) = this%cvar%data4d(i,:,:,l) - this%delxy(:,:)*(this%flux%data3d(i,:,:) - this%flux%data3d(i-1,:,:))
1510  END DO
1511  END DO
1512 
1513 !#ifdef PARALLEL
1514 ! ! We only need to do something, if we (also) are dealing with domain decomposition in
1515 ! ! the second (phi) direction
1516 ! IF(PAR_DIMS.GT.1) THEN
1517 ! DO i=GMIN1,GMAX1
1518 ! IF(this%shift(i,k).GT.0) THEN
1519 ! DO l=1,Physics%VNUM+Physics%PNUM
1520 ! IF(Mesh%SN_shear) THEN
1521 ! this%buf(k,1:this%shift(i)) = this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k)
1522 ! ELSE
1523 ! this%buf(k,1:this%shift(i)) = this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k)
1524 ! END IF
1525 ! END DO
1526 ! CALL MPI_Sendrecv_replace(&
1527 ! this%buf,&
1528 ! this%shift(i)*Physics%VNUM, &
1529 ! DEFAULT_MPI_REAL, &
1530 ! Mesh%neighbor(EAST), i+Mesh%GNUM, &
1531 ! Mesh%neighbor(WEST), i+Mesh%GNUM, &
1532 ! Mesh%comm_cart, status, ierror)
1533 ! DO k=1,Physics%VNUM
1534 ! IF(Mesh%SN_shear) THEN
1535 ! this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k) = this%buf(k,1:this%shift(i))
1536 ! ELSE
1537 ! this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k) = this%buf(k,1:this%shift(i))
1538 ! END IF
1539 ! END DO
1540 ! ELSE IF(this%shift(i).LT.0) THEN
1541 ! DO k=1,Physics%VNUM
1542 ! IF(Mesh%SN_shear) THEN
1543 ! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k)
1544 ! ELSE
1545 ! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k)
1546 ! END IF
1547 ! END DO
1548 ! CALL MPI_Sendrecv_replace(&
1549 ! this%buf,&
1550 ! -this%shift(i)*Physics%VNUM, &
1551 ! DEFAULT_MPI_REAL, &
1552 ! Mesh%neighbor(EAST), i+Mesh%GNUM, &
1553 ! Mesh%neighbor(WEST), i+Mesh%GNUM, &
1554 ! Mesh%comm_cart, status, ierror)
1555 ! DO k=1,Physics%VNUM
1556 ! IF (Mesh%SN_shear) THEN
1557 ! this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k) = this%buf(k,1:-this%shift(i))
1558 ! ELSE
1559 ! this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k) = this%buf(k,1:-this%shift(i))
1560 ! END IF
1561 ! END DO
1562 ! END IF
1563 ! END DO
1564 ! END DO
1565 ! END IF
1566 !#endif
1567 
1568  ! Integer shift along x- or y-direction
1569 !NEC$ SHORTLOOP
1570  DO l=1,physics%VNUM+physics%PNUM
1571 !NEC$ IVDEP
1572  DO k=mesh%KGMIN,mesh%KGMAX
1573 !NEC$ IVDEP
1574  DO j=mesh%JGMIN,mesh%JGMAX
1575  this%cvar%data4d(mesh%IMIN:mesh%IMAX,j,k,l) &
1576  = cshift(this%cvar%data4d(mesh%IMIN:mesh%IMAX,j,k,l),-this%shift(j,k),1)
1577  END DO
1578  END DO
1579  END DO
1580 
1581  ! convert conservative to primitive variables
1582  CALL physics%Convert2Primitive(this%cvar,this%pvar)
1583  ! Calculate RHS after the Advection Step
1584  CALL this%ComputeRHS(mesh,physics,sources,fluxes,this%time,0.,this%pvar,this%cvar,&
1585  this%checkdatabm,this%rhs)
1586 
1587  this%cold%data1d(:) = this%cvar%data1d(:)
1588 
1589 
1590  END SUBROUTINE fargoadvectionx
1591 
1592 
1611  SUBROUTINE fargoadvectiony(this,Fluxes,Mesh,Physics,Sources)
1612  IMPLICIT NONE
1613  !------------------------------------------------------------------------!
1614  CLASS(timedisc_base), INTENT(INOUT) :: this
1615  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
1616  CLASS(mesh_base), INTENT(IN) :: Mesh
1617  CLASS(physics_base), INTENT(INOUT) :: Physics
1618  CLASS(sources_base), POINTER :: Sources
1619  !------------------------------------------------------------------------!
1620  INTEGER :: i,j,k,l
1621 #ifdef PARALLEL
1622  CHARACTER(LEN=80) :: str
1623  INTEGER :: status(MPI_STATUS_SIZE)
1624  INTEGER :: ierror
1625  REAL :: mpi_buf(2*Mesh%GNUM)
1626  REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX) :: buf
1627 #endif
1628  !------------------------------------------------------------------------!
1629  ! determine step size of integer shift and length of remaining transport step
1630  ! first compute the whole step
1631  this%delxy(:,:) = this%w(:,:) * this%dt / mesh%dlx%data3d(:,mesh%JMIN,:)
1632 
1633 #ifdef PARALLEL
1634  ! make sure all MPI processes use the same step if domain is decomposed
1635  ! along the y-direction (can be different due to round-off errors)
1636  IF (mesh%dims(2).GT.1) THEN
1637  CALL mpi_allreduce(mpi_in_place,this%delxy,(mesh%IGMAX-mesh%IGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1), &
1638  default_mpi_real,mpi_min,mesh%Jcomm,ierror)
1639  END IF
1640 #endif
1641 
1642  ! then subdivide into integer shift and remaining linear advection step
1643 !NEC$ COLLAPSE
1644  DO k = mesh%KGMIN,mesh%KGMAX
1645 !NEC$ IVDEP
1646  DO i = mesh%IGMIN,mesh%IGMAX
1647  this%shift(i,k) = nint(this%delxy(i,k))
1648  this%delxy(i,k) = this%delxy(i,k)-dble(this%shift(i,k))
1649  END DO
1650  END DO
1651 
1652 !NEC$ SHORTLOOP
1653  DO l=1,physics%VNUM+physics%PNUM
1654  this%dq%data3d(:,mesh%JGMIN,:) = this%cvar%data4d(:,mesh%JGMIN+1,:,l) - this%cvar%data4d(:,mesh%JGMIN,:,l)
1655 !NEC$ IVDEP
1656  DO j=mesh%JGMIN+1,mesh%JGMAX-1
1657  this%dq%data3d(:,j,:) = this%cvar%data4d(:,j+1,:,l)-this%cvar%data4d(:,j,:,l)
1658  ! apply minmod limiter
1659  WHERE (sign(1.0,this%dq%data3d(:,j-1,:))*sign(1.0,this%dq%data3d(:,j,:)).GT.0)
1660  this%dql%data3d(:,j,:) = sign(min(abs(this%dq%data3d(:,j-1,:)),abs(this%dq%data3d(:,j,:))),this%dq%data3d(:,j-1,:))
1661  ELSEWHERE
1662  this%dql%data3d(:,j,:) = 0.
1663  END WHERE
1664  END DO
1665 !NEC$ IVDEP
1666  DO j=mesh%JMIN-1,mesh%JMAX
1667  WHERE(this%delxy(:,:).GT.0.)
1668  this%flux%data3d(:,j,:) = this%cvar%data4d(:,j,:,l) + .5 * this%dql%data3d(:,j,:) * (1. - this%delxy(:,:))
1669  ELSEWHERE
1670  this%flux%data3d(:,j,:) = this%cvar%data4d(:,j+1,:,l) - .5*this%dql%data3d(:,j+1,:)*(1. + this%delxy(:,:))
1671  END WHERE
1672  this%cvar%data4d(:,j,:,l) = this%cvar%data4d(:,j,:,l) - this%delxy(:,:)*(this%flux%data3d(:,j,:) - this%flux%data3d(:,j-1,:))
1673  END DO
1674  END DO
1675 
1676 
1677 !#ifdef PARALLEL
1678 ! ! We only need to do something, if we (also) are dealing with domain decomposition in
1679 ! ! the second (phi) direction
1680 ! IF(PAR_DIMS.GT.1) THEN
1681 ! DO i=GMIN1,GMAX1
1682 ! IF(this%shift(i,k).GT.0) THEN
1683 ! DO l=1,Physics%VNUM+Physics%PNUM
1684 ! IF(Mesh%SN_shear) THEN
1685 ! this%buf(k,1:this%shift(i)) = this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k)
1686 ! ELSE
1687 ! this%buf(k,1:this%shift(i)) = this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k)
1688 ! END IF
1689 ! END DO
1690 ! CALL MPI_Sendrecv_replace(&
1691 ! this%buf,&
1692 ! this%shift(i)*Physics%VNUM, &
1693 ! DEFAULT_MPI_REAL, &
1694 ! Mesh%neighbor(EAST), i+Mesh%GNUM, &
1695 ! Mesh%neighbor(WEST), i+Mesh%GNUM, &
1696 ! Mesh%comm_cart, status, ierror)
1697 ! DO k=1,Physics%VNUM
1698 ! IF(Mesh%SN_shear) THEN
1699 ! this%cvar%data4d(MAX2-this%shift(i)+1:MAX2,i,k) = this%buf(k,1:this%shift(i))
1700 ! ELSE
1701 ! this%cvar%data4d(i,MAX2-this%shift(i)+1:MAX2,k) = this%buf(k,1:this%shift(i))
1702 ! END IF
1703 ! END DO
1704 ! ELSE IF(this%shift(i).LT.0) THEN
1705 ! DO k=1,Physics%VNUM
1706 ! IF(Mesh%SN_shear) THEN
1707 ! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k)
1708 ! ELSE
1709 ! this%buf(k,1:-this%shift(i)) = this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k)
1710 ! END IF
1711 ! END DO
1712 ! CALL MPI_Sendrecv_replace(&
1713 ! this%buf,&
1714 ! -this%shift(i)*Physics%VNUM, &
1715 ! DEFAULT_MPI_REAL, &
1716 ! Mesh%neighbor(EAST), i+Mesh%GNUM, &
1717 ! Mesh%neighbor(WEST), i+Mesh%GNUM, &
1718 ! Mesh%comm_cart, status, ierror)
1719 ! DO k=1,Physics%VNUM
1720 ! IF (Mesh%SN_shear) THEN
1721 ! this%cvar%data4d(MIN2:MIN2-this%shift(i)-1,i,k) = this%buf(k,1:-this%shift(i))
1722 ! ELSE
1723 ! this%cvar%data4d(i,MIN2:MIN2-this%shift(i)-1,k) = this%buf(k,1:-this%shift(i))
1724 ! END IF
1725 ! END DO
1726 ! END IF
1727 ! END DO
1728 ! END DO
1729 ! END IF
1730 !#endif
1731 
1732  ! Integer shift along x- or y-direction
1733 !NEC$ SHORTLOOP
1734  DO l=1,physics%VNUM+physics%PNUM
1735 !NEC$ IVDEP
1736  DO k=mesh%KGMIN,mesh%KGMAX
1737 !NEC$ IVDEP
1738  DO i=mesh%IGMIN,mesh%IGMAX
1739  this%cvar%data4d(i,mesh%JMIN:mesh%JMAX,k,l) &
1740  = cshift(this%cvar%data4d(i,mesh%JMIN:mesh%JMAX,k,l),-this%shift(i,k),1)
1741  END DO
1742  END DO
1743  END DO
1744 
1745  ! convert conservative to primitive variables
1746  CALL physics%Convert2Primitive(this%cvar,this%pvar)
1747  ! Calculate RHS after the Advection Step
1748  CALL this%ComputeRHS(mesh,physics,sources,fluxes,this%time,0.,this%pvar,this%cvar,&
1749  this%checkdatabm,this%rhs)
1750 
1751  this%cold%data1d(:) = this%cvar%data1d(:)
1752 
1753 
1754  END SUBROUTINE fargoadvectiony
1755 
1756 
1760  SUBROUTINE calcbackgroundvelocity(this,Mesh,Physics,pvar,cvar,w)
1761  IMPLICIT NONE
1762  !------------------------------------------------------------------------!
1763  CLASS(timedisc_base), INTENT(INOUT) :: this
1764  CLASS(mesh_base), INTENT(IN) :: Mesh
1765  CLASS(physics_base), INTENT(INOUT) :: Physics
1766  CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
1767  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%KGMIN:Mesh%KGMAX), &
1768  INTENT(OUT) :: w
1769  !------------------------------------------------------------------------!
1770  REAL :: wi
1771  INTEGER :: i,k
1772 #ifdef PARALLEL
1773  INTEGER :: ierror
1774 #endif
1775  !------------------------------------------------------------------------!
1776  ! make sure we are using true, i.e. non-selenoidal, quantities
1777  CALL physics%AddBackgroundVelocityY(mesh,w,pvar,cvar)
1778  SELECT TYPE(p => pvar)
1779  CLASS IS(statevector_eulerisotherm)
1780  DO k=mesh%KGMIN,mesh%KGMAX
1781  DO i=mesh%IGMIN,mesh%IGMAX
1782  ! some up all yvelocities along the y-direction
1783  wi = sum(p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,k,2))
1784 #ifdef PARALLEL
1785  ! extend the sum over all partitions
1786  IF(mesh%dims(2).GT.1) THEN
1787  CALL mpi_allreduce(mpi_in_place, wi, 1, default_mpi_real, mpi_sum, &
1788  mesh%Jcomm, ierror)
1789  END IF
1790 #endif
1791  ! set new background velocity to the arithmetic mean of the
1792  ! yvelocity field along the y-direction
1793  w(i,k) = wi / mesh%JNUM
1794  END DO
1795  END DO
1796  END SELECT
1797  END SUBROUTINE calcbackgroundvelocity
1798 
1799 
1826  FUNCTION getcentrifugalvelocity(this,Mesh,Physics,Fluxes,Sources,&
1827  dir_omega_,accel_,centrot) RESULT(velocity)
1828  USE physics_euler_mod, ONLY: physics_euler
1829  USE physics_eulerisotherm_mod, ONLY: physics_eulerisotherm
1830  IMPLICIT NONE
1831  !------------------------------------------------------------------------!
1832  CLASS(timedisc_base), INTENT(INOUT) :: this
1833  CLASS(mesh_base), INTENT(IN) :: mesh
1834  CLASS(physics_base), INTENT(INOUT) :: physics
1835  CLASS(fluxes_base), INTENT(INOUT) :: fluxes
1836  CLASS(sources_base), POINTER :: sources
1837  REAL, OPTIONAL, INTENT(IN) :: dir_omega_(3)
1838  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VDIM), &
1839  OPTIONAL, INTENT(IN) :: accel_
1840  REAL, OPTIONAL, INTENT(IN) :: centrot(3)
1841  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VDIM) &
1842  :: velocity
1843  !------------------------------------------------------------------------!
1844  CLASS(marray_base), ALLOCATABLE &
1845  :: accel,dist_axis_projected,ephi_projected,eomega,tmpvec, &
1846  posvec,tmp
1847  REAL :: omega2,dir_omega(3)
1848  INTEGER :: k,l,err
1849  !------------------------------------------------------------------------!
1850  ALLOCATE(accel,dist_axis_projected,ephi_projected,posvec,eomega,tmpvec,tmp, &
1851  stat=err)
1852  IF (err.NE.0) CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
1853  "Unable to allocate memory.")
1854  ! initialize marrays
1855  tmp = marray_base()
1856  tmpvec = marray_base(3)
1857  eomega = marray_base(3)
1858  posvec = marray_base(3)
1859  dist_axis_projected = marray_base(physics%VDIM)
1860  ephi_projected = marray_base(physics%VDIM)
1861  accel = marray_base(physics%VDIM)
1862 
1863  ! do some sanity checks on the input data
1864  IF (PRESENT(dir_omega_)) THEN
1865  dir_omega(:) = dir_omega_(:)
1866  omega2 = sum(dir_omega(:)*dir_omega(:))
1867  ! omega must not be the zero vector
1868  IF (.NOT. (omega2.GT.0.0)) &
1869  CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
1870  "omega must not be the zero vector")
1871  ! normalize dir_omega
1872  dir_omega(:) = dir_omega(:) / sqrt(omega2)
1873  ELSE
1874  ! default is rotation in a plane perpendicular to ê_z
1875  dir_omega(1:3) = (/ 0.0, 0.0, 1.0 /)
1876  END IF
1877 
1878  ! compute (local) curvilinear vector components of the unit
1879  ! vectors pointing in the direction of the rotational axis
1880  tmpvec%data2d(:,1) = dir_omega(1)
1881  tmpvec%data2d(:,2) = dir_omega(2)
1882  tmpvec%data2d(:,3) = dir_omega(3)
1883  CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,tmpvec%data4d,eomega%data4d)
1884 
1885  ! 1. Get the curvilinear components of all vectors pointing
1886  ! from the center of rotation into each cell using cartesian
1887  ! coordinates.
1888  ! The cartesian coordinates are the components of the position
1889  ! vector with respect to the cartesian standard orthonormal
1890  ! basis {ê_x, ê_y, ê_z}.
1891  IF (PRESENT(centrot)) THEN
1892  IF ((mesh%rotsym.GT.0.0).AND.any(centrot(:).NE.0.0)) &
1893  CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
1894  "if rotational symmetry is enabled center of rotation must be the origin")
1895  ! Translate the position vector to the center of rotation.
1896  DO k=1,3
1897  tmpvec%data4d(:,:,:,k) = mesh%bccart(:,:,:,k) - centrot(k)
1898  END DO
1899  ! compute curvilinear components of translated position vectors
1900  CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,tmpvec%data4d,posvec%data4d)
1901  ELSE
1902  ! default center of rotation is the origin of the mesh
1903  ! -> use position vectors with respect to the origin
1904  posvec%data4d = mesh%posvec%bcenter
1905  END IF
1906 
1907  ! 2. Compute vectors pointing from axis of rotation into each cell
1908  ! and the azimuthal unit vector ephi
1909  tmp%data1d(:) = sum(posvec%data2d(:,:)*eomega%data2d(:,:),dim=2)
1910  posvec%data2d(:,1) = posvec%data2d(:,1) - tmp%data1d(:)*eomega%data2d(:,1)
1911  posvec%data2d(:,2) = posvec%data2d(:,2) - tmp%data1d(:)*eomega%data2d(:,2)
1912  posvec%data2d(:,3) = posvec%data2d(:,3) - tmp%data1d(:)*eomega%data2d(:,3)
1913 
1914  ! distance to axis
1915  tmp%data1d(:) = sqrt(sum(posvec%data2d(:,:)*posvec%data2d(:,:),dim=2))
1916  ! use tmpvec for temporary storage of ephi
1917  CALL cross_product(eomega%data2d(:,1),eomega%data2d(:,2),eomega%data2d(:,3), &
1918  posvec%data2d(:,1)/tmp%data1d(:),posvec%data2d(:,2)/tmp%data1d(:), &
1919  posvec%data2d(:,3)/tmp%data1d(:),tmpvec%data2d(:,1),tmpvec%data2d(:,2),tmpvec%data2d(:,3))
1920 
1921  ! project it onto the mesh
1922  SELECT TYPE(phys => physics)
1923  CLASS IS(physics_eulerisotherm)
1924  l = 1
1925  DO k=1,3
1926  ! check if vector component is used
1927  IF (btest(mesh%VECTOR_COMPONENTS,k-1)) THEN
1928  dist_axis_projected%data2d(:,l) = posvec%data2d(:,k)
1929  ephi_projected%data2d(:,l) = tmpvec%data2d(:,k)
1930  l = l + 1
1931  END IF
1932  END DO
1933  CLASS DEFAULT
1934  CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
1935  "physics not supported")
1936  END SELECT
1937 
1938  ! 3. determine the acceleration
1939  IF(PRESENT(accel_)) THEN
1940  accel%data4d(:,:,:,:) = accel_(:,:,:,:)
1941  ELSE
1942  ! Works only for centrot = (0,0,0)
1943  IF(PRESENT(centrot)) THEN
1944  IF(any(centrot(:).NE.0.0)) &
1945  CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
1946  "You are not allowed to define centrot without accel.")
1947  END IF
1948  ! This may not work for physics with unusual conservative variables.
1949  ! We assume: conservative momentum = density * velocity
1950  ! sanity check if someone derives new physics from the euler modules
1951  ! which does not have this property
1952  SELECT TYPE(phys => physics)
1953  TYPE IS(physics_eulerisotherm)
1954  ! supported -> do nothing
1955  TYPE IS(physics_euler)
1956  ! supported -> do nothing
1957  CLASS DEFAULT
1958  CALL this%Error("timedisc_base::GetCentrifugalVelocity", &
1959  "physics not supported")
1960  END SELECT
1961  SELECT TYPE(phys => physics)
1962  CLASS IS(physics_eulerisotherm)
1963  ! evaluate the right hand side for given (preliminary) initial data
1964  CALL this%ComputeRHS(mesh,phys,sources,fluxes,this%time,0.0, &
1965  this%pvar,this%cvar,this%checkdatabm,this%rhs)
1966  SELECT TYPE(pvar => this%pvar)
1967  CLASS IS(statevector_eulerisotherm)
1968  SELECT TYPE(rhs => this%rhs)
1969  CLASS IS(statevector_eulerisotherm)
1970  DO k=1,phys%VDIM
1971  accel%data2d(:,k) = -rhs%momentum%data2d(:,k) / pvar%density%data1d(:)
1972  END DO
1973  END SELECT
1974  END SELECT
1975  END SELECT
1976  END IF
1977 
1978  ! 4. compute absolute value of azimuthal velocity
1979  tmp%data1d(:) = sqrt(max(0.0,-sum(accel%data2d(:,:)*dist_axis_projected%data2d(:,:),dim=2)))
1980 
1981  ! 5. Compute components of the azimuthal velocity vector
1982  DO k=1,physics%VDIM
1983  velocity(:,:,:,k) = tmp%data3d(:,:,:) * ephi_projected%data4d(:,:,:,k)
1984  END DO
1985 
1986  CALL accel%Destroy()
1987  CALL dist_axis_projected%Destroy()
1988  CALL ephi_projected%Destroy()
1989  CALL posvec%Destroy()
1990  CALL eomega%Destroy()
1991  CALL tmpvec%Destroy()
1992  CALL tmp%Destroy()
1993  DEALLOCATE(accel,dist_axis_projected,ephi_projected,posvec,eomega,tmpvec,tmp)
1994 
1995  CONTAINS
1996 
1997  ELEMENTAL SUBROUTINE cross_product(ax,ay,az,bx,by,bz,aXbx,aXby,aXbz)
1998  IMPLICIT NONE
1999  !------------------------------------------------------------------------!
2000  REAL, INTENT(IN) :: ax,ay,az,bx,by,bz
2001  REAL, INTENT(OUT) :: axbx,axby,axbz
2002  !------------------------------------------------------------------------!
2003  axbx = ay*bz - az*by
2004  axby = az*bx - ax*bz
2005  axbz = ax*by - ay*bx
2006  END SUBROUTINE cross_product
2007 
2008  END FUNCTION getcentrifugalvelocity
2009 
2010  SUBROUTINE finalize_base(this)
2011  IMPLICIT NONE
2012  !------------------------------------------------------------------------!
2013  CLASS(timedisc_base) :: this
2014  !------------------------------------------------------------------------!
2015  IF (.NOT.this%Initialized()) &
2016  CALL this%Error("CloseTimedisc","not initialized")
2017  ! call boundary destructor
2018  CALL this%Boundary%Finalize()
2019 
2020  CALL this%pvar%Destroy()
2021  CALL this%ptmp%Destroy()
2022  CALL this%cvar%Destroy()
2023  CALL this%ctmp%Destroy()
2024  CALL this%cold%Destroy()
2025  CALL this%geo_src%Destroy()
2026  CALL this%src%Destroy()
2027  CALL this%rhs%Destroy()
2028 
2029  DEALLOCATE( &
2030  this%pvar,this%cvar,this%ptmp,this%ctmp, &
2031  this%geo_src,this%src,this%rhs, &
2032  this%xfluxdydz,this%yfluxdzdx,this%zfluxdxdy,this%amax,this%tol_abs,&
2033  this%dtmean,this%dtstddev,this%time)
2034 
2035  IF (ASSOCIATED(this%cerr)) THEN
2036  CALL this%cerr%Destroy()
2037  DEALLOCATE(this%cerr)
2038  END IF
2039  IF (ASSOCIATED(this%cerr_max)) THEN
2040  CALL this%cerr_max%Destroy()
2041  DEALLOCATE(this%cerr_max)
2042  END IF
2043  IF (ASSOCIATED(this%w)) DEALLOCATE(this%w)
2044  IF (ASSOCIATED(this%delxy))DEALLOCATE(this%delxy)
2045  IF (ASSOCIATED(this%shift))DEALLOCATE(this%shift)
2046 #ifdef PARALLEL
2047  IF(ASSOCIATED(this%buf)) DEALLOCATE(this%buf)
2048 #endif
2049  IF(ASSOCIATED(this%bflux)) DEALLOCATE(this%bflux)
2050  IF(ASSOCIATED(this%solution)) THEN
2051  CALL this%solution%Destroy()
2052  DEALLOCATE(this%solution)
2053  END IF
2054  END SUBROUTINE finalize_base
2055 
2056 
2057 END MODULE timedisc_base_mod
timedisc_base_mod::check_pmin
integer, parameter, public check_pmin
Definition: timedisc_base.f90:192
mesh_base
logging_base_mod::finalize
subroutine finalize(this)
Destructor of common class.
Definition: logging_base.f90:155
timedisc_base_mod::fargo_method
character(len=40), dimension(3), parameter fargo_method
Definition: timedisc_base.f90:197
timedisc_base_mod::check_all
integer, parameter, public check_all
Definition: timedisc_base.f90:190
sources_base_mod
generic source terms module providing functionaly common to all source terms
Definition: sources_base.f90:42
timedisc_base_mod::dtcause_cfl
integer, parameter, public dtcause_cfl
Definition: timedisc_base.f90:184
sources_base_mod::sources_base
Definition: sources_base.f90:57
logging_base_mod::default_mpi_real
integer, save default_mpi_real
default real type for MPI
Definition: logging_base.f90:80
timedisc_base_mod::finalize_base
subroutine finalize_base(this)
Definition: timedisc_base.f90:2011
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
boundary_generic_mod
Generic boundary module.
Definition: boundary_generic.f90:44
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
physics_eulerisotherm_mod::physics_eulerisotherm
Definition: physics_eulerisotherm.f90:57
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
logging_base_mod::logging_base
common data structure
Definition: logging_base.f90:56
timedisc_base_mod::computerhs
subroutine computerhs(this, Mesh, Physics, Sources, Fluxes, time, dt, pvar, cvar, checkdatabm, rhs)
compute the RHS of the spatially discretized PDE
Definition: timedisc_base.f90:993
timedisc_base_mod::timedisc_base
Definition: timedisc_base.f90:79
geometry_cartesian_mod::geometry_cartesian
Definition: geometry_cartesian.f90:41
geometry_logcylindrical_mod
defines properties of a 3D logcylindrical mesh
Definition: geometry_logcylindrical.f90:33
timedisc_base_mod::rejectsolution
subroutine rejectsolution(this, Mesh, Physics, Sources, Fluxes, time, dt)
Definition: timedisc_base.f90:906
timedisc_base_mod::getcentrifugalvelocity
real function, dimension(mesh%igmin:mesh%igmax, mesh%jgmin:mesh%jgmax, mesh%kgmin:mesh%kgmax, physics%vdim) getcentrifugalvelocity(this, Mesh, Physics, Fluxes, Sources, dir_omega_, accel_, centrot)
Compute velocity leading to a centrifugal acceleration with respect to some given axis of rotation wh...
Definition: timedisc_base.f90:1828
timedisc_base_mod::fargoadvectiony
subroutine fargoadvectiony(this, Fluxes, Mesh, Physics, Sources)
Calculates the linear transport step in Fargo Scheme along y-axis .
Definition: timedisc_base.f90:1612
geometry_cylindrical_mod
defines properties of a 3D cylindrical mesh
Definition: geometry_cylindrical.f90:36
timedisc_base_mod::fargoadvectionx
subroutine fargoadvectionx(this, Fluxes, Mesh, Physics, Sources)
Calculates the linear transport step in Fargo Scheme along x-axis .
Definition: timedisc_base.f90:1449
timedisc_base_mod::dtcause_fileio
integer, parameter, public dtcause_fileio
smallest ts due to fileio
Definition: timedisc_base.f90:187
timedisc_base_mod::integrationstep
subroutine integrationstep(this, Mesh, Physics, Sources, Fluxes, iter, config, IO)
Definition: timedisc_base.f90:630
timedisc_base_mod::setoutput
subroutine setoutput(this, Mesh, Physics, config, IO)
Definition: timedisc_base.f90:515
timedisc_base_mod::ssprk
integer, parameter, public ssprk
Definition: timedisc_base.f90:182
timedisc_base_mod::getcfl
pure real function getcfl(this)
Definition: timedisc_base.f90:1421
sources_gravity_mod
generic gravity terms module providing functionaly common to all gravity terms
Definition: sources_gravity.f90:41
timedisc_base_mod::calcbackgroundvelocity
subroutine calcbackgroundvelocity(this, Mesh, Physics, pvar, cvar, w)
Calculates new background velocity for fargo advection.
Definition: timedisc_base.f90:1761
timedisc_base_mod::computeerror
real function computeerror(this, Mesh, Physics, cvar_high, cvar_low)
Definition: timedisc_base.f90:937
timedisc_base_mod::check_nothing
integer, parameter, public check_nothing
Definition: timedisc_base.f90:189
reconstruction_base_mod
base module for reconstruction process
Definition: reconstruction_base.f90:35
timedisc_base_mod::dtcause_smallerr
integer, parameter, public dtcause_smallerr
Definition: timedisc_base.f90:186
physics_generic_mod
constructor for physics class
Definition: physics_generic.f90:35
timedisc_base_mod::check_csound
integer, parameter, public check_csound
Definition: timedisc_base.f90:191
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:64
timedisc_base_mod::check_tmin
integer, parameter, public check_tmin
Definition: timedisc_base.f90:195
physics_eulerisotherm_mod
physics module for 1D,2D and 3D isothermal Euler equations
Definition: physics_eulerisotherm.f90:44
timedisc_base_mod::acceptsolution
subroutine acceptsolution(this, Mesh, Physics, Sources, Fluxes, time, dt, iter)
Definition: timedisc_base.f90:858
physics_euler_mod::statevector_euler
Definition: physics_euler.f90:103
timedisc_base_mod::dtcause_erradj
integer, parameter, public dtcause_erradj
Definition: timedisc_base.f90:185
physics_base_mod::physics_base
named integer constants for flavour of state vectors
Definition: physics_base.f90:69
geometry_cartesian_mod
defines properties of a 3D cartesian mesh
Definition: geometry_cartesian.f90:36
timedisc_base_mod::checkdata
subroutine checkdata(this, Mesh, Physics, Fluxes, pvar, cvar, checkdatabm)
compute the RHS of the spatially discretized PDE
Definition: timedisc_base.f90:1334
timedisc_base_mod::modified_euler
integer, parameter, public modified_euler
Definition: timedisc_base.f90:178
cross_product
elemental subroutine cross_product(ax, ay, az, bx, by, bz, aXbx, aXby, aXbz)
Definition: timedisc_base.f90:1998
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
timedisc_base_mod::check_invalid
integer, parameter, public check_invalid
Definition: timedisc_base.f90:194
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
timedisc_base_mod::getorder
pure integer function getorder(this)
Definition: timedisc_base.f90:1412
physics_euler_mod
physics module for 1D,2D and 3D non-isothermal Euler equations
Definition: physics_euler.f90:41
geometry_logcylindrical_mod::geometry_logcylindrical
Definition: geometry_logcylindrical.f90:39
fluxes_base_mod::fluxes_base
Definition: fluxes_base.f90:61
timedisc_base_mod::adjusttimestep
real function adjusttimestep(this, maxerr, dtold)
adjust the time step
Definition: timedisc_base.f90:706
timedisc_base_mod::rk_fehlberg
integer, parameter, public rk_fehlberg
Definition: timedisc_base.f90:179
timedisc_base_mod::calctimestep
real function calctimestep(this, Mesh, Physics, Sources, Fluxes, time, dtcause)
Determines the CFL time step and time step limits due to source terms.
Definition: timedisc_base.f90:765
timedisc_base_mod::dormand_prince
integer, parameter, public dormand_prince
Definition: timedisc_base.f90:181
physics_euler_mod::physics_euler
Definition: physics_euler.f90:54
fluxes_base_mod
base module for numerical flux functions
Definition: fluxes_base.f90:39
timedisc_base_mod::cash_karp
integer, parameter, public cash_karp
Definition: timedisc_base.f90:180
physics_eulerisotherm_mod::statevector_eulerisotherm
Definition: physics_eulerisotherm.f90:121
timedisc_base_mod::inittimedisc
subroutine inittimedisc(this, Mesh, Physics, config, IO, ttype, tname)
Definition: timedisc_base.f90:216
geometry_cylindrical_mod::geometry_cylindrical
Definition: geometry_cylindrical.f90:41
sources_gravity_mod::sources_gravity
Definition: sources_gravity.f90:58
timedisc_base_mod
Definition: timedisc_base.f90:53
timedisc_base_mod::check_rhomin
integer, parameter, public check_rhomin
Definition: timedisc_base.f90:193
boundary_generic_mod::boundary_generic
Definition: boundary_generic.f90:66
Generated by   doxygen 1.8.14