88 :: pvar => null(),cvar => null(), &
89 ptmp => null(),ctmp => null(), &
91 src => null(),geo_src => null(), &
93 cerr => null(),cerr_max => null(), &
100 INTEGER :: dtcause,dtmincause
102 REAL,
POINTER :: dtmean, dtstddev
110 LOGICAL :: always_update_bccsound
117 INTEGER :: checkdatabm
118 REAL :: pmin,rhomin,tmin
123 REAL,
DIMENSION(:),
POINTER :: tol_abs
127 REAL,
DIMENSION(:,:,:,:),
POINTER :: phi,oldphi_s,&
129 REAL,
DIMENSION(:),
POINTER :: gamma
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()
161 procedure(solveode),
DEFERRED :: solveode
167 SUBROUTINE solveode(this,Mesh,Physics,Sources,Fluxes,time,dt,err)
173 CLASS(
sources_list),
ALLOCATABLE,
INTENT(INOUT) :: Sources
175 REAL,
INTENT(IN) :: time
176 REAL,
INTENT(INOUT) :: dt, err
207 "dynamic background velocity ", &
208 "user supplied fixed background velocity ", &
236 TYPE(
dict_typ),
POINTER :: config,IO
237 INTEGER,
INTENT(IN) :: ttype
238 CHARACTER(LEN=32),
INTENT(IN) :: tname
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
245 CALL this%InitLogging(ttype,tname)
247 IF (.NOT.physics%Initialized().OR..NOT.mesh%Initialized()) &
248 CALL this%Error(
"InitTimedisc",
"physics and/or mesh module uninitialized")
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, &
261 CALL this%Error(
"InitTimedisc",
"Unable to allocate memory.")
265 CALL physics%new_statevector(this%pvar,
primitive)
266 CALL physics%new_statevector(this%ptmp,
primitive)
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.
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
301 this%dtmincause = -99
305 CALL getattr(config,
"cfl", this%cfl, 0.4)
308 CALL getattr(config,
"dtlimit", this%dtlimit, epsilon(this%dtlimit)*this%stoptime)
311 CALL getattr(config,
"dtmax", this%dtmax, 5.0)
314 CALL getattr(config,
"maxiter", this%maxiter, 0)
317 CALL getattr(config,
"tol_rel", this%tol_rel, 0.01)
320 this%tol_abs(:) = 0.001
321 CALL getattr(config,
"tol_abs", this%tol_abs, this%tol_abs)
324 CALL getattr(config,
"rhstype", this%rhstype, 0)
330 CALL getattr(config,
"beta", this%beta, 0.08)
341 CALL getattr(config,
"always_update_bccsound", d, 1)
343 this%always_update_bccsound = .false.
345 this%always_update_bccsound = .true.
350 CALL getattr(config,
"checkdata", this%checkdatabm,
check_all)
353 CALL getattr(config,
"pmin", this%pmin, tiny(this%pmin))
356 CALL getattr(config,
"tmin", this%tmin, tiny(this%tmin))
359 CALL getattr(config,
"rhomin", this%rhomin, tiny(this%rhomin))
361 CALL getattr(config,
"output/"//trim(physics%pvarname(1)), d, 1)
363 CALL this%SetOutput(mesh,physics,config,io)
367 SELECT TYPE(phys => physics)
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!")
377 SELECT CASE(mesh%fargo%GetDirection())
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), &
384 this%buf(physics%VNUM+physics%PNUM,1:mesh%MININUM), &
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), &
393 this%buf(physics%VNUM+physics%PNUM,1:mesh%MINJNUM), &
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), &
402 this%buf(physics%VNUM+physics%PNUM,1:mesh%MINJNUM), &
409 IF (err.EQ.0.AND.mesh%fargo%GetType().EQ.1)
ALLOCATE(this%wfactor,mold=this%w,stat=err)
411 CALL this%Error(
"timedisc_base::InitTimedisc",
"Unable to allocate memory for fargo advection.")
414 SELECT CASE(mesh%fargo%GetType())
417 CALL this%Warning(
"timedisc_base::InitTimedisc",
"FARGO transport is experimental, use with care!")
421 this%shift(:,:) = 0.0
423 this%dq%data1d(:) = 0.0
425 this%dql%data1d(:) = 0.0
427 this%flux%data1d(:) = 0.0
428 IF (
ASSOCIATED(this%wfactor))
THEN
434 SELECT CASE(mesh%fargo%GetDirection())
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)
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)
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)
462 IF(mesh%shear_dir.EQ.2)
THEN
463 this%w(:,:) = -mesh%Q*mesh%omega*mesh%bcenter(:,mesh%JMIN,:,1)
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)
468 this%shift(:,:) = 0.0
469 shear_direction =
"south<->north"
472 this%dq%data1d(:) = 0.0
474 this%dql%data1d(:) = 0.0
476 this%flux%data1d(:) = 0.0
478 CALL this%Error(
"timedisc_base::InitTimedisc",
"unknown fargo advection scheme")
481 CALL this%Error(
"timedisc_base::InitTimedisc",
"selected physics doesn't support fargo transport")
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))
497 IF (this%tol_rel.LT.1.0.AND.this%order.GT.1)
THEN
500 this%cerr%data1d(:) = 0.
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)//
" %")
506 WRITE (info_str,
'(A)')
"disabled"
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)
515 IF (.NOT.
ASSOCIATED(this%w))
THEN
516 SELECT CASE(mesh%Geometry%GetAzimuthIndex())
518 ALLOCATE(this%w(mesh%IGMIN:mesh%IGMAX,mesh%KGMIN:mesh%KGMAX),stat = err)
520 ALLOCATE(this%w(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX),stat = err)
523 CALL this%Error(
"InitTimedisc",
"Unable to allocate memory for special RHS time stepping.")
527 CALL this%Info(
" special rhs: enabled")
529 CALL this%Error(
"InitTimedisc",
"unknown rhstype")
540 TYPE(
dict_typ),
POINTER :: config,IO
542 INTEGER :: valwrite,i
543 CHARACTER(LEN=128) :: key,pvar_key
544 LOGICAL :: writeSolution
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.
553 this%write_error = .false.
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)
562 NULLIFY(this%solution)
563 writesolution = .false.
566 CALL getattr(config,
"output/time", valwrite, 1)
568 CALL setattr(io,
"time", ref(this%time))
572 key = trim(physics%pvarname(i))
575 CALL getattr(config,
"output/" // trim(key), valwrite, 1)
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))
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))
588 key = trim(physics%cvarname(i))
589 IF(key.EQ.pvar_key) key = trim(key) //
"_cvar"
591 CALL getattr(config,
"output/" // trim(key), valwrite, 0)
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))
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))
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))
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))
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))
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))
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)
654 CLASS(
sources_list),
ALLOCATABLE,
INTENT(INOUT) :: Sources
657 TYPE(
dict_typ),
POINTER :: config,IO
659 REAL :: err,dtold,dt,time
661 INTENT(INOUT) :: iter
664 IF (mesh%fargo%GetType().EQ.1) &
665 CALL this%CalcBackgroundVelocity(mesh,physics,this%pvar,this%cvar)
667 SELECT CASE(mesh%fargo%GetDirection())
669 CALL physics%SubtractBackgroundVelocityX(mesh,this%w,this%pvar,this%cvar)
671 CALL physics%SubtractBackgroundVelocityY(mesh,this%w,this%pvar,this%cvar)
673 CALL physics%SubtractBackgroundVelocityZ(mesh,this%w,this%pvar,this%cvar)
682 this%dtmincause = this%dtcause
686 timestep:
DO WHILE (.NOT.this%break.AND.time+dt.LE.this%time+this%dt)
689 CALL this%SolveODE(mesh,physics,sources,fluxes,time,dt,err)
695 CALL this%Warning(
"timedisc_base::IntegrationStep",
"max error returned by SolveODE is NaN")
697 ELSE IF (err.GT.huge(err))
THEN
699 CALL this%Warning(
"timedisc_base::IntegrationStep",
"max error returned by SolveODE is Inf")
701 ELSE IF (err.GE.1.0)
THEN
703 CALL this%RejectSolution(mesh,physics,sources,fluxes,time,dt)
706 CALL this%AcceptSolution(mesh,physics,sources,fluxes,time,dtold,iter)
710 IF (dt.LT.this%dtlimit) this%break = .true.
714 IF (.NOT.this%break)
THEN
716 this%dt = time - this%time
722 IF (mesh%fargo%GetType().GT.0)
THEN
730 CALL this%boundary%CenterBoundary(mesh,physics,this%time,this%pvar,this%cvar)
731 SELECT CASE(mesh%fargo%GetDirection())
733 CALL this%FargoAdvectionX(fluxes,mesh,physics,sources)
735 CALL this%FargoAdvectionY(fluxes,mesh,physics,sources)
737 CALL this%FargoAdvectionZ(fluxes,mesh,physics,sources)
740 CALL physics%Convert2Primitive(this%cvar,this%pvar)
742 CALL this%ComputeRHS(mesh,physics,sources,fluxes,this%time,0.,this%pvar,this%cvar,&
743 this%checkdatabm,this%rhs)
745 CALL this%ComputeRHS(mesh,physics,sources,fluxes,this%time,dtold,this%pvar,this%cvar,&
746 this%checkdatabm,this%rhs)
748 this%cold = this%cvar
763 REAL :: maxerr,dtold,dtnew
767 INTENT(IN) :: maxerr,dtold
768 INTENT(INOUT) ::
this
771 IF (
this%tol_rel.GE.1.0)
THEN
780 IF(maxerr.GT.0.)
THEN
790 IF (
this%maxerrold.GT.0.)
THEN
791 dtnew = dttmp * exp(-
this%beta*log(
this%maxerrold/maxerr))
797 IF (maxerr.LT.1.0)
THEN
799 dtnew = min(dtnew,4.*dtold)
806 IF (dtnew.GT.dtold) dtnew = dttmp
808 dtnew = max(dtnew,0.25*dtold)
812 this%maxerrold = maxerr
825 CLASS(
sources_list),
ALLOCATABLE,
INTENT(INOUT) :: sources
827 REAL,
INTENT(IN) :: time
828 INTEGER,
INTENT(INOUT) :: dtcause
831 REAL :: dt_cfl, dt_src
832 REAL :: invdt_x,invdt_y,invdt_z
836 IF (.NOT.
this%always_update_bccsound)
THEN
837 SELECT TYPE(phys => physics)
839 SELECT TYPE(pvar =>
this%pvar)
841 CALL phys%UpdateSoundSpeed(pvar)
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)
853 SELECT TYPE(pvar =>
this%pvar)
855 CALL phys%CalculateWaveSpeeds(mesh,pvar,fluxes%minwav,fluxes%maxwav)
860 IF ((mesh%JNUM.EQ.1).AND.(mesh%KNUM.EQ.1))
THEN
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
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
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
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
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
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(:))
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)
896 dt_cfl =
this%cfl / invdt
901 IF (
ALLOCATED(sources))
THEN
904 CALL sources%CalcTimestep(mesh,physics,fluxes,
this%pvar,
this%cvar,time,dt_src,dtcause)
905 dt = min(dt_cfl,dt_src)
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
923 REAL,
DIMENSION(:),
POINTER :: bflux,bfold
925 INTENT(INOUT) :: 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
937 IF (mesh%INUM.GT.1)
THEN
941 bflux(1:
SIZE(fluxes%bxflux)) => fluxes%bxflux
942 bfold(1:
SIZE(fluxes%bxfold)) => fluxes%bxfold
946 IF (mesh%JNUM.GT.1)
THEN
947 bflux(1:
SIZE(fluxes%byflux)) => fluxes%byflux
948 bfold(1:
SIZE(fluxes%byfold)) => fluxes%byfold
952 IF (mesh%KNUM.GT.1)
THEN
953 bflux(1:
SIZE(fluxes%bzflux)) => fluxes%bzflux
954 bfold(1:
SIZE(fluxes%bzfold)) => fluxes%bzfold
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
976 this%cvar = this%cold
978 CALL this%ComputeRHS(mesh,physics,sources,fluxes,time,dt,this%pvar, &
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)
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)
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)
1002 this%n_adj = this%n_adj + 1
1005 IF (dt.LT.this%dtmin)
THEN
1007 this%dtmincause = this%dtcause
1016 CLASS(mesh_base),
INTENT(IN) :: mesh
1017 CLASS(physics_base),
INTENT(IN) :: physics
1018 CLASS(marray_compound),
INTENT(INOUT):: cvar_high,cvar_low
1025 REAL :: rel_err(physics%vnum)
1029 DO concurrent(i=1:
SIZE(this%cerr%data2d,dim=1))
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))
1035 rel_err(l) = maxval(this%cerr%data2d(:,l),mask=mesh%without_ghost_zones%mask1d(:))
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))
1045 maxerr = maxval(rel_err(:))
1049 CALL mpi_allreduce(mpi_in_place,maxerr,1,default_mpi_real,mpi_max,&
1050 mesh%comm_cart,ierror)
1073 SUBROUTINE computerhs(this,Mesh,Physics,Sources,Fluxes,time,dt,pvar,cvar,checkdatabm,rhs)
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
1090 LOGICAL :: have_potential
1093 REAL,
DIMENSION(:,:,:,:),
POINTER :: pot
1094 CLASS(sources_base),
POINTER :: sp => null()
1095 CLASS(sources_gravity),
POINTER :: grav => null()
1098 SELECT CASE(mesh%fargo%GetType())
1102 SELECT CASE(mesh%fargo%GetDirection())
1104 CALL physics%AddBackgroundVelocityX(mesh,this%w,pvar,cvar)
1106 CALL physics%AddBackgroundVelocityY(mesh,this%w,pvar,cvar)
1108 CALL physics%AddBackgroundVelocityZ(mesh,this%w,pvar,cvar)
1115 SELECT CASE(mesh%fargo%GetDirection())
1117 CALL physics%SubtractBackgroundVelocityX(mesh,this%w,pvar,cvar)
1119 CALL physics%SubtractBackgroundVelocityY(mesh,this%w,pvar,cvar)
1121 CALL physics%SubtractBackgroundVelocityZ(mesh,this%w,pvar,cvar)
1134 CALL this%boundary%CenterBoundary(mesh,physics,t,pvar,cvar)
1137 this%tmin.GT.1.e-10)
THEN
1138 SELECT TYPE(p => pvar)
1140 SELECT TYPE(c => cvar)
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)
1151 physics%PRESSURE.GT.0)
THEN
1154 SELECT TYPE(p => pvar)
1156 SELECT TYPE(c => cvar)
1159 p%pressure%data1d(:) &
1160 = max(p%pressure%data1d(:),this%pmin)
1161 CALL physics%Convert2Conservative(p,c)
1167 IF (this%always_update_bccsound)
THEN
1168 SELECT TYPE(phys => physics)
1170 SELECT TYPE(pvar => this%pvar)
1172 CALL phys%UpdateSoundSpeed(pvar)
1179 CALL this%CheckData(mesh,physics,fluxes,pvar,cvar,checkdatabm)
1182 CALL physics%GeometricalSources(mesh,pvar,cvar,this%geo_src)
1185 IF (
ALLOCATED(sources)) &
1186 CALL sources%ExternalSources(mesh,physics,fluxes,sources, &
1187 t,dt,pvar,cvar,this%src)
1192 SELECT CASE(mesh%fargo%GetType())
1197 IF (mesh%geometry%GetType().EQ.cartesian) this%geo_src%data1d(:) = 0.0
1198 SELECT CASE(mesh%fargo%GetDirection())
1200 CALL physics%AddFargoSourcesX(mesh,this%w,pvar,cvar,this%geo_src)
1201 CALL physics%SubtractBackgroundVelocityX(mesh,this%w,pvar,cvar)
1203 CALL physics%AddFargoSourcesY(mesh,this%w,pvar,cvar,this%geo_src)
1204 CALL physics%SubtractBackgroundVelocityY(mesh,this%w,pvar,cvar)
1206 CALL physics%AddFargoSourcesZ(mesh,this%w,pvar,cvar,this%geo_src)
1207 CALL physics%SubtractBackgroundVelocityZ(mesh,this%w,pvar,cvar)
1218 CALL fluxes%CalculateFluxes(mesh,physics,pvar,cvar,this%xfluxdydz, &
1219 this%yfluxdzdx,this%zfluxdxdy)
1224 rhsbdflux:
DO l=1,physics%VNUM+physics%PNUM
1226 DO k=mesh%KMIN,mesh%KMAX
1228 DO j=mesh%JMIN,mesh%JMAX
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)
1235 DO k=mesh%KMIN,mesh%KMAX
1237 DO i=mesh%IMIN,mesh%IMAX
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)
1244 DO j=mesh%JMIN,mesh%JMAX
1246 DO i=mesh%IMIN,mesh%IMAX
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)
1254 SELECT CASE(this%rhstype)
1257 DO l=1,physics%VNUM+physics%PNUM
1258 DO k=mesh%KMIN,mesh%KMAX
1259 DO j=mesh%JMIN,mesh%JMAX
1261 DO i=mesh%IMIN,mesh%IMAX
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)
1277 IF (
ALLOCATED(sources)) &
1278 sp => sources%GetSourcesPointer(gravity)
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.
1291 SELECT CASE(mesh%geometry%GetAzimuthIndex())
1293 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX
1294 DO j=mesh%JMIN-mesh%JP1,mesh%JMAX
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)
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)
1313 wp = this%w(i,k) + mesh%hy%bcenter(i,j,k)*mesh%OMEGA
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)
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)
1330 wp = this%w(i,k) + mesh%hy%bcenter(i,j,k)*mesh%OMEGA
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)
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)
1352 SELECT TYPE (phys => physics)
1354 DO k=mesh%KMIN,mesh%KMAX
1355 DO j=mesh%JMIN,mesh%JMAX
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
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))
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))
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) )
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) )
1396 this%geo_src%data4d(i,j,k,physics%ENERGY) = 0.
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) )
1418 CALL this%Error(
"timedisc_base::ComputeRHS",
"physics not supported with rhstype=1")
1423 DO l=1,physics%VNUM+physics%PNUM
1425 DO k=mesh%KMIN,mesh%KMAX
1426 DO j=mesh%JMIN,mesh%JMAX
1428 DO i=mesh%IMIN,mesh%IMAX
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)
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))
1446 DO k=mesh%KMIN,mesh%KMAX
1447 DO j=mesh%JMIN,mesh%JMAX
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)
1467 DO k=mesh%KMIN-mesh%KP1,mesh%KMAX
1468 DO j=mesh%JMIN-mesh%JP1,mesh%JMAX
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)
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)
1487 wp = this%w(i,j) + mesh%hz%bcenter(i,j,k)*mesh%OMEGA
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)
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)
1504 wp = this%w(i,j) + mesh%hz%bcenter(i,j,k)*mesh%OMEGA
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)
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)
1526 SELECT TYPE (phys => physics)
1528 DO k=mesh%KMIN,mesh%KMAX
1529 DO j=mesh%JMIN,mesh%JMAX
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
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))
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) )
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) )
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) )
1568 this%geo_src%data4d(i,j,k,physics%ENERGY) = 0.
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) )
1592 CALL this%Error(
"timedisc_base::ComputeRHS",
"physics not supported with rhstype=1")
1597 DO l=1,physics%VNUM+physics%PNUM
1599 DO k=mesh%KMIN,mesh%KMAX
1600 DO j=mesh%JMIN,mesh%JMAX
1602 DO i=mesh%IMIN,mesh%IMAX
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))
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))
1620 DO k=mesh%KMIN,mesh%KMAX
1621 DO j=mesh%JMIN,mesh%JMAX
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)
1641 CALL this%Error(
"timedisc_base::ComputeRHS",
"rhstype 1 is not supported for this geometry")
1647 DO l=1,physics%VNUM+physics%PNUM
1649 DO k=mesh%KMIN,mesh%KMAX
1650 DO j=mesh%JMIN,mesh%JMAX
1652 DO i=mesh%IMIN,mesh%IMAX
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)
1660 CALL this%Error(
"timedisc_base::ComputeRHS",
"only rhstype 0 or 1 supported")
1666 SUBROUTINE checkdata(this,Mesh,Physics,Fluxes,pvar,cvar,checkdatabm)
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
1685 SELECT TYPE(phys => physics)
1687 val = minval(phys%bccsound%data1d(:))
1690 CALL this%Warning(
"CheckData",
"Illegal speed of sound value less than 0.")
1694 CALL this%Warning(
"CheckData",
"check speed of sound selected, but bccsound not defined")
1699 SELECT TYPE(p => pvar)
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
1705 CALL this%Warning(
"CheckData",
"Pressure below allowed pmin value.")
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
1719 CALL this%Warning(
"CheckData",
"Density below allowed rhomin value.")
1726 IF(any((cvar%data1d.NE.cvar%data1d).OR.(pvar%data1d.NE.pvar%data1d)))
THEN
1728 CALL this%Warning(
"CheckData",
"Found NaN in pvar or cvar.")
1731 IF(any((cvar%data1d.GT.huge(cvar%data1d)).OR.pvar%data1d.GT.huge(pvar%data1d)))
THEN
1733 CALL this%Warning(
"CheckData",
"Found Infinity in pvar or cvar.")
1739 CALL mpi_allreduce(mpi_in_place, this%break, 1, mpi_logical, mpi_lor, &
1740 mesh%comm_cart, err)
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
1792 CHARACTER(LEN=80) :: str
1793 INTEGER :: status(MPI_STATUS_SIZE)
1795 REAL :: mpi_buf(2*Mesh%GNUM)
1796 REAL,
DIMENSION(Mesh%JMIN:Mesh%JMAX) :: buf
1801 this%delxy(:,:) = this%w(:,:) * this%dt / mesh%dlx%data3d(mesh%IMIN,:,:)
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)
1814 DO k = mesh%KGMIN,mesh%KGMAX
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))
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)
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)
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,:,:))
1832 this%dql%data3d(i,:,:) = 0.
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(:,:))
1840 this%flux%data3d(i,:,:) = this%cvar%data4d(i+1,:,:,l) - .5*this%dql%data3d(i+1,:,:)*(1. + this%delxy(:,:))
1842 this%cvar%data4d(i,:,:,l) = this%cvar%data4d(i,:,:,l) - this%delxy(:,:)*(this%flux%data3d(i,:,:) - this%flux%data3d(i-1,:,:))
1903 DO l=1,physics%VNUM+physics%PNUM
1905 DO k=mesh%KGMIN,mesh%KGMAX
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)
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
1946 CHARACTER(LEN=80) :: str
1947 INTEGER :: status(MPI_STATUS_SIZE)
1949 REAL :: mpi_buf(2*Mesh%GNUM)
1950 REAL,
DIMENSION(Mesh%JMIN:Mesh%JMAX) :: buf
1955 this%delxy(:,:) = this%w(:,:) * this%dt / mesh%dly%data3d(:,mesh%JMIN,:)
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)
1968 DO k = mesh%KGMIN,mesh%KGMAX
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))
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)
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)
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,:))
1986 this%dql%data3d(:,j,:) = 0.
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(:,:))
1994 this%flux%data3d(:,j,:) = this%cvar%data4d(:,j+1,:,l) - .5*this%dql%data3d(:,j+1,:)*(1. + this%delxy(:,:))
1996 this%cvar%data4d(:,j,:,l) = this%cvar%data4d(:,j,:,l) - this%delxy(:,:)*(this%flux%data3d(:,j,:) - this%flux%data3d(:,j-1,:))
2058 DO l=1,physics%VNUM+physics%PNUM
2060 DO k=mesh%KGMIN,mesh%KGMAX
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)
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
2100 CHARACTER(LEN=80) :: str
2101 INTEGER :: status(MPI_STATUS_SIZE)
2103 REAL :: mpi_buf(2*Mesh%GNUM)
2104 REAL,
DIMENSION(Mesh%JMIN:Mesh%JMAX) :: buf
2109 this%delxy(:,:) = this%w(:,:) * this%dt / mesh%dlz%data3d(:,:,mesh%KMIN)
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)
2122 DO j = mesh%JGMIN,mesh%JGMAX
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))
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)
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)
2138 WHERE (this%dq%data3d(:,:,k-1)*this%dq%data3d(:,:,k).LE.0.)
2140 this%dql%data3d(:,:,k) = 0.
2143 this%dql%data3d(:,:,k) = sign(min(abs(this%dq%data3d(:,:,k-1)),abs(this%dq%data3d(:,:,k))),this%dq%data3d(:,:,k))
2153 DO k=mesh%KMIN-1,mesh%KMAX
2157 WHERE(this%delxy(:,:).GT.0.)
2158 this%flux%data3d(:,:,k) = this%cvar%data4d(:,:,k,l) + .5 * this%dql%data3d(:,:,k) * (1. - this%delxy(:,:))
2160 this%flux%data3d(:,:,k) = this%cvar%data4d(:,:,k+1,l) - .5*this%dql%data3d(:,:,k+1)*(1. + this%delxy(:,:))
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))
2226 DO l=1,physics%VNUM+physics%PNUM
2228 DO j=mesh%JGMIN,mesh%JGMAX
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)
2247 CLASS(mesh_base),
INTENT(IN) :: Mesh
2248 CLASS(physics_base),
INTENT(INOUT) :: Physics
2249 CLASS(marray_compound),
INTENT(INOUT) :: pvar,cvar
2252 INTEGER :: i,j,k,v_idx
2259 SELECT TYPE(p => pvar)
2260 CLASS IS(statevector_eulerisotherm)
2261 v_idx = mesh%fargo%GetDirection()
2262 SELECT CASE(mesh%fargo%GetDirection())
2265 CALL physics%AddBackgroundVelocityX(mesh,this%w,pvar,cvar)
2266 DO k=mesh%KGMIN,mesh%KGMAX
2267 DO j=mesh%JGMIN,mesh%JGMAX
2269 wi = sum(p%velocity%data4d(mesh%IMIN:mesh%IMAX,j,k,v_idx)*this%wfactor(j,k))
2272 IF(mesh%dims(1).GT.1)
THEN
2273 CALL mpi_allreduce(mpi_in_place, wi, 1, default_mpi_real, mpi_sum, &
2279 this%w(j,k) = wi / mesh%INUM
2284 CALL physics%AddBackgroundVelocityY(mesh,this%w,pvar,cvar)
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
2290 wi = sum(p%velocity%data4d(i,mesh%JMIN:mesh%JMAX,k,v_idx)*this%wfactor(i,k))
2293 IF(mesh%dims(2).GT.1)
THEN
2294 CALL mpi_allreduce(mpi_in_place, wi, 1, default_mpi_real, mpi_sum, &
2298 this%w(i,k) = wi / mesh%JNUM
2303 CALL physics%AddBackgroundVelocityZ(mesh,this%w,pvar,cvar)
2305 v_idx = physics%VDIM
2306 DO j=mesh%JGMIN,mesh%JGMAX
2307 DO i=mesh%IGMIN,mesh%IGMAX
2309 wi = sum(p%velocity%data4d(i,j,mesh%KMIN:mesh%KMAX,v_idx)*this%wfactor(i,j))
2312 IF(mesh%dims(3).GT.1)
THEN
2313 CALL mpi_allreduce(mpi_in_place, wi, 1, default_mpi_real, mpi_sum, &
2319 this%w(i,j) = wi / mesh%KNUM
2324 CALL this%Error(
"timedisc_base::CalcBackgroundVelocity",
"physics currently not supported with fargo transport")
2356 dir_omega_,accel_,centrot)
RESULT(velocity)
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) &
2373 CLASS(marray_base),
ALLOCATABLE &
2374 :: accel,dist_axis_projected,ephi_projected,eomega,tmpvec, &
2376 REAL :: omega2,dir_omega(3)
2379 ALLOCATE(accel,dist_axis_projected,ephi_projected,posvec,eomega,tmpvec,tmp, &
2381 IF (err.NE.0)
CALL this%Error(
"timedisc_base::GetCentrifugalVelocity", &
2382 "Unable to allocate memory.")
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)
2393 IF (
PRESENT(dir_omega_))
THEN
2394 dir_omega(:) = dir_omega_(:)
2395 omega2 = sum(dir_omega(:)*dir_omega(:))
2397 IF (.NOT. (omega2.GT.0.0)) &
2398 CALL this%Error(
"timedisc_base::GetCentrifugalVelocity", &
2399 "omega must not be the zero vector")
2401 dir_omega(:) = dir_omega(:) / sqrt(omega2)
2404 dir_omega(1:3) = (/ 0.0, 0.0, 1.0 /)
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)
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")
2426 tmpvec%data4d(:,:,:,k) = mesh%bccart(:,:,:,k) - centrot(k)
2429 CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,tmpvec%data4d,posvec%data4d)
2433 posvec%data4d = mesh%posvec%bcenter
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)
2444 tmp%data1d(:) = sqrt(sum(posvec%data2d(:,:)*posvec%data2d(:,:),dim=2))
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))
2451 SELECT TYPE(phys => physics)
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)
2463 CALL this%Error(
"timedisc_base::GetCentrifugalVelocity", &
2464 "physics not supported")
2468 IF(
PRESENT(accel_))
THEN
2469 accel%data4d(:,:,:,:) = accel_(:,:,:,:)
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.")
2481 SELECT TYPE(phys => physics)
2487 CALL this%Error(
"timedisc_base::GetCentrifugalVelocity", &
2488 "physics not supported")
2490 SELECT TYPE(phys => physics)
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)
2500 accel%data2d(:,k) = -rhs%momentum%data2d(:,k) / pvar%density%data1d(:)
2508 tmp%data1d(:) = sqrt(max(0.0,-sum(accel%data2d(:,:)*dist_axis_projected%data2d(:,:),dim=2)))
2512 velocity(:,:,:,k) = tmp%data3d(:,:,:) * ephi_projected%data4d(:,:,:,k)
2515 DEALLOCATE(accel,dist_axis_projected,ephi_projected,posvec,eomega,tmpvec,tmp)
2522 REAL,
INTENT(IN) :: ax,ay,az,bx,by,bz
2523 REAL,
INTENT(OUT) :: axbx,axby,axbz
2525 axbx = ay*bz - az*by
2526 axby = az*bx - ax*bz
2527 axbz = ax*by - ay*bx
2537 IF (.NOT.this%Initialized()) &
2538 CALL this%Error(
"CloseTimedisc",
"not initialized")
2540 CALL this%Boundary%Finalize()
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)
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)
2555 IF(
ASSOCIATED(this%buf))
DEALLOCATE(this%buf)
2557 IF(
ASSOCIATED(this%bflux))
DEALLOCATE(this%bflux)
2558 IF(
ASSOCIATED(this%solution))
DEALLOCATE(this%solution)
elemental real function pc(tau, rho_s0, cs_inf, gamma)
Dictionary for generic data types.
type(logging_base), save this
base module for numerical flux functions
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
subroutine finalize(this)
Destructor of logging_base class.
base class for mesh arrays
derived class for compound of mesh arrays
integer, parameter east
named constant for eastern boundary
integer, parameter bottom
named constant for bottom boundary
integer, parameter south
named constant for southern boundary
integer, parameter top
named constant for top boundary
integer, parameter north
named constant for northern boundary
integer, parameter west
named constant for western boundary
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
container class to manage the list of source terms
elemental subroutine cross_product(ax, ay, az, bx, by, bz, aXbx, aXby, aXbz)