module mod_dt_srhd use mod_variables use mod_physics_vars implicit none private public :: setdt contains ! instantiate here for inlining in kernels pure subroutine to_primitive(u) use mod_con2prim real(dp), intent(inout) :: u(nw_phys) real(dp) :: rho,rhoh,pth,E real(dp) :: ssqr !! begin: call srhd_get_auxiliary(u) ssqr=(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2) if(srhd_eos)then call con2prim_eos(u(lfac_),u(xi_),u(iw_rho),ssqr,u(iw_e)) else call con2prim(u(lfac_),u(xi_),u(iw_rho),ssqr,u(iw_e)) endif !! end: call srhd_get_auxiliary(u) rho=u(iw_rho)/u(lfac_) rhoh=u(xi_)/u(lfac_)**2 !! begin: call srhd_get_pressure_eos(rho,rhoh,pth,E) if(srhd_eos)then E = (rhoh+dsqrt(rhoh**2+(srhd_gamma**2-1.0_dp)*rho**2)) & /(srhd_gamma+1.0_dp) pth = 0.5d0*gamma_1* (E-rho*(rho/E)) else pth =(rhoh-rho)/gamma_to_gamma_1 endif !! end: call srhd_get_pressure_eos(rho,rhoh,pth,E) u(iw_rho)=rho u(iw_mom(1))=u(lfac_)*u(iw_mom(1))/u(xi_) u(iw_mom(2))=u(lfac_)*u(iw_mom(2))/u(xi_) u(iw_mom(3))=u(lfac_)*u(iw_mom(3))/u(xi_) u(iw_e)=pth end subroutine to_primitive !> Returns maximum local signal speed from primitive state u in direction flux_dim; !> used in LLF/TVDLF flux estimation. pure real(dp) function get_cmax(u, x, flux_dim) result(wC) real(dp), intent(in) :: u(nw_phys) real(dp), intent(in) :: x(1:ndim) integer, intent(in) :: flux_dim real(dp) :: tmp1,tmp2,v2,csound2,vidim,cmax,cmin real(dp) :: rho,rhoh,pth,E ! input u is in primitive form ! auxiliaries are filled here rho=u(iw_rho) rhoh=u(xi_)/u(lfac_)**2 pth=u(iw_e) !!tmp1=u(iw_rho) !!tmp2=u(xi_)/u(lfac_)**2.0d0 !!v2=u(iw_e) !! begin: call srhd_get_csound2_prim_eos(tmp1,tmp2,v2,csound2) if(srhd_eos) then E = (rhoh+dsqrt(rhoh**2+(srhd_gamma**2-1.0_dp)& *rho**2))/(srhd_gamma+1.0_dp) csound2=(pth*((srhd_gamma+1.0_dp)& +gamma_1*(rho/E)**2))& /(2.0d0*rhoh) else csound2=srhd_gamma*pth/rhoh end if !! end call srhd_get_csound2_prim_eos(tmp1,tmp2,v2,csound2) v2 = 1.0d0-1.0d0/u(lfac_)**2 vidim = u(iw_mom(flux_dim))/u(lfac_) tmp2 = vidim**2 tmp1 = 1.0d0-v2*csound2-tmp2*(1.0d0-csound2) ! the radicand is positive analytically, as in estimate_speeds_minmax; clip ! it so an inadmissible state (c_s^2 > 1 from a collapsed rhoh) yields a ! degenerate speed rather than a NaN. get_cmax feeds setdt, where a single ! NaN survives the global MPI_ALLREDUCE and takes dt down for every rank. tmp2 = dsqrt(max(csound2*(1.0_dp-v2)*tmp1, 0.0_dp)) tmp1 = vidim*(1.0_dp-csound2) cmax = (tmp1+tmp2)/(1.0_dp-v2*csound2) cmin = (tmp1-tmp2)/(1.0_dp-v2*csound2) ! Limit by speed of light cmin = max(cmin, - 1.0d0) cmin = min(cmin, 1.0d0) cmax = max(cmax, - 1.0d0) cmax = min(cmax, 1.0d0) ! now take extremal value only for dt limit wC = max(dabs(cmax),dabs(cmin)) end function get_cmax subroutine phys_get_dt(w, x, dx, dtnew) real(dp), intent(in) :: w(nw_phys), x(1:ndim), dx(1:ndim) real(dp), intent(out) :: dtnew dtnew = huge(1.0d0) ! TODO: user-set time step limit, also in other physics modules!!! end subroutine phys_get_dt !>setdt - set dt for all levels between levmin and levmax. !> dtpar>0 --> use fixed dtpar for all level !> dtpar<=0 --> determine CFL limited timestep subroutine setdt() use mod_global_parameters integer :: iigrid, igrid, idims, ix1,ix2,ix3, ifile double precision :: dtmin_mype, factor, dx1,dx2,dx3, dtmax double precision :: dxinv(1:ndim), cmaxtot, cmax, u(1:nw_phys), cmax_mype double precision :: xloc(1:ndim), qdtnew if (dtpar<=zero) then dtmin_mype=bigdouble do iigrid=1,igridstail_active; igrid=igrids_active(iigrid) dx1=rnode(rpdx1_,igrid);dx2=rnode(rpdx2_,igrid) dx3=rnode(rpdx3_,igrid); dxinv(1)=one/dx1;dxinv(2)=one/dx2;dxinv(3)=one/dx3; do ix3=ixMlo3,ixMhi3 do ix2=ixMlo2,ixMhi2 do ix1=ixMlo1,ixMhi1 u = bg(1)%w(ix1, ix2, ix3, 1:nw_phys, igrid) call to_primitive(u) xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, igrid) cmaxtot = 0.0d0 do idims = 1, ndim cmax = get_cmax(u, xloc, idims) cmaxtot = cmaxtot + cmax * dxinv(idims) end do dtmin_mype = min( dtmin_mype, courantpar / cmaxtot ) u = bg(1)%w(ix1,ix2,ix3,1:nw_phys,igrid) xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, igrid) call phys_get_dt(u, xloc, one/dxinv(1:ndim), qdtnew) dtmin_mype = min( dtmin_mype, qdtnew ) end do end do end do end do else dtmin_mype=dtpar end if if (need_global_cmax) then cmax_mype=-bigdouble do iigrid=1,igridstail_active; igrid=igrids_active(iigrid) do ix3=ixMlo3,ixMhi3 do ix2=ixMlo2,ixMhi2 do ix1=ixMlo1,ixMhi1 u(1:nw_phys) = bg(1)%w(ix1, ix2, ix3, 1:nw_phys, igrid) call to_primitive(u) xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, igrid) do idims = 1, ndim cmax = get_cmax(u, xloc, idims) cmax_mype = max( cmax_mype, cmax ) end do end do end do end do end do end if if (dtmin_mype<dtmin) then write(unitterm,*)"Error: Time step too small!", dtmin_mype write(unitterm,*)"on processor:", mype, "at time:", global_time,& " step:", it write(unitterm,*)"Lower limit of time step:", dtmin crash=.true. end if if (slowsteps>it-it_init+1) then factor=one-(one-dble(it-it_init+1)/dble(slowsteps))**2 dtmin_mype=dtmin_mype*factor end if if(final_dt_reduction)then if(time_max-global_time<=dtmin) then final_dt_exit=.true. endif dtmin_mype=min(dtmin_mype,time_max-global_time) end if if (dtpar<=zero) then call MPI_ALLREDUCE(dtmin_mype,dt,1,MPI_DOUBLE_PRECISION,MPI_MIN, icomm,& ierrmpi) else dt=dtmin_mype end if if (need_global_cmax) then call MPI_ALLREDUCE(cmax_mype, cmax_global, 1, MPI_DOUBLE_PRECISION, MPI_MAX, icomm, & ierrmpi) end if if(any(dtsave(1:nfile)<bigdouble).or.any(tsave(isavet(1:nfile),& 1:nfile)<bigdouble))then dtmax = minval(ceiling(global_time/dtsave(1:nfile))*dtsave(1:nfile))-& global_time do ifile=1,nfile dtmax = min(tsave(isavet(ifile),ifile)-global_time,dtmax) end do if(dtmax > smalldouble)then dt=min(dt,dtmax) else ! dtmax=0 means dtsave is divisible by global_time dt=min(dt,minval(dtsave(1:nfile))) end if end if if(mype==0) then if(any(dtsave(1:nfile)<dt)) then write(unitterm,*) 'Warning: timesteps: ',dt,& ' exceeding output intervals ', dtsave(1:nfile) endif endif end subroutine setdt end module mod_dt_srhd