!> Module for handling split source terms (split from the fluxes) module mod_source_ffhd use mod_variables use mod_physics_vars use mod_global_parameters, only: ndim implicit none public !> How to apply dimensional splitting to the source terms, see !> @ref discretization.md !> defaulting to sfs integer :: sourcesplit = 0 integer, parameter :: sourcesplit_sfs = 0 integer, parameter :: sourcesplit_sf = 1 integer, parameter :: sourcesplit_ssf = 2 integer, parameter :: sourcesplit_ssfss = 3 public :: add_split_source contains ! instantiate the templated functions here for inlining: subroutine addsource_local(qdt, dtfactor, qtC, wCT, wCTprim, qt, wnew, x, dr, & qsourcesplit) use mod_usr, only: addsource_usr use mod_usr, only: gravity_field use mod_radiative_cooling, only: radiative_cooling_add_source real(dp), intent(in) :: qdt, dtfactor, qtC, qt real(dp), intent(in) :: wCT(nw_phys), wCTprim(nw_phys) real(dp), intent(in) :: x(1:ndim), dr(ndim) real(dp), intent(inout) :: wnew(nw_phys) logical, intent(in) :: qsourcesplit ! .. local .. integer :: idim real(dp) :: field, mag, divb if (.not. qsourcesplit) then !--------------------------------- ! unsplit sources !--------------------------------- do idim = 1, ndim field = gravity_field(wCT, x, idim) mag = wCTprim(iw_b1-1+idim) wnew(iw_mom(1)) = wnew(iw_mom(1)) + qdt * field * wCT(iw_rho) * mag wnew(iw_e) = wnew(iw_e) + qdt * field * wCT(iw_mom(1)) * mag end do call radiative_cooling_add_source(qdt,wCT,wCTprim,wnew,x) call addsource_usr(qdt, qt, wCT, wCTprim, wnew, x, .false.) else !--------------------------------- ! split sources !--------------------------------- ! Not yet implemented end if end subroutine addsource_local subroutine addsource_nonlocal(qdt, dtfactor, qtC, wCTprim, qt, wnew, x, dx, idir, & qsourcesplit) use mod_global_parameters, only: dt, cmax_global, courantpar, third real(dp), intent(in) :: qdt, dtfactor, qtC, qt real(dp), intent(in) :: wCTprim(nw_phys,5) real(dp), intent(in) :: x(1:ndim), dx(1:ndim) real(dp), intent(inout) :: wnew(nw_phys) integer, intent(in) :: idir logical, intent(in) :: qsourcesplit ! .. local .. real(dp) :: tau, htc_qrsc, sigT real(dp) :: Te(1:5), gradT real(dp) :: mag5(1:5), divb, mag if (.not. qsourcesplit) then !--------------------------------- ! unsplit sources !--------------------------------- ! > p*divb mag5(1:5) = wCTprim(iw_b1-1+idir,1:5) divb = (8*mag5(4)-8*mag5(2)-mag5(5)+mag5(1))/12.0_dp/dx(idir) wnew(iw_mom(1)) = wnew(iw_mom(1)) + qdt*wCTprim(iw_e,3)*divb !> gradient of temperature: Te(1:5) = wCTprim(iw_e,1:5) / wCTprim(iw_rho,1:5) mag = wCTprim(iw_b1-1+idir,3) gradT = (8.d0*(Te(4)-Te(2))-Te(5)+Te(1))/(12.d0*dx(idir)) sigT = hypertc_kappa * sqrt(Te(3)**5) tau = max(4.d0*dt, sigT*Te(3)*courantpar**2*gamma_1/& (wCTprim(iw_e,3)*cmax_global**2)) htc_qrsc = sigT * mag * gradT wnew(iw_q) = wnew(iw_q) - qdt * (htc_qrsc + wCTprim(iw_q,3)*third) / tau else !--------------------------------- ! split sources !--------------------------------- ! Not yet implemented end if end subroutine addsource_nonlocal subroutine addsource_compact(qdt, dtfactor, qtC, wCTprim1, wCTprim2, wCTprim3, qt, wnew, x, dx, & qsourcesplit) real(dp), intent(in) :: qdt, dtfactor, qtC, qt real(dp), intent(in) :: wCTprim1(nw_phys,3),wCTprim2(nw_phys,3),wCTprim3(nw_phys,3) real(dp), intent(in) :: x(1:ndim), dx(1:ndim) real(dp), intent(inout) :: wnew(nw_phys) logical, intent(in) :: qsourcesplit end subroutine addsource_compact pure subroutine to_primitive(u) real(dp), intent(inout) :: u(nw_phys) u(iw_mom(1)) = u(iw_mom(1))/u(iw_rho) u(iw_e) = gamma_1 * (u(iw_e) - 0.5_dp * u(iw_rho) * & u(iw_mom(1))**2 ) end subroutine to_primitive subroutine add_split_source(prior) use mod_global_parameters use mod_ghostcells_update use mod_physics, only: phys_req_diagonal, phys_global_source_after use mod_comm_lib, only: mpistop logical, intent(in) :: prior real(dp) :: qt real(dp) :: dr(ndim), xloc(ndim) ! Physical cell size handed to the source term. dr is the *logical* ! spacing, which is not a length in a curvilinear system and, under ! LOG_RADIUS, is not even proportional to one. real(dp) :: dloc(ndim) integer :: iigrid, n, ix1, ix2, ix3 integer :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3 real(dp) :: wprim(nw_phys), wnew(nw_phys), wCT(nw_phys) logical :: src_active src_active = .false. if ((.not.prior).and.(sourcesplit==sourcesplit_sf .or. & sourcesplit==sourcesplit_ssf)) return if (prior) then qt=global_time else qt=global_time+dt end if if(any_source_split) then ! add split source terms ixOmin1=ixGlo1+nghostcells;ixOmin2=ixGlo2+nghostcells ixOmin3=ixGlo3+nghostcells;ixOmax1=ixGhi1-nghostcells ixOmax2=ixGhi2-nghostcells;ixOmax3=ixGhi3-nghostcells; select case (sourcesplit) case (sourcesplit_sfs) do iigrid=1,igridstail_active n = igrids_active(iigrid) dr = rnode(rpdx1_:rnodehi, n) do ix3=ixOmin3,ixOmax3 do ix2=ixOmin2,ixOmax2 do ix1=ixOmin1,ixOmax1 xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, n) dloc = dr wCT = bg(1)%w(ix1,ix2,ix3, 1:nw_phys, n) wnew = wCT wprim = wCT call to_primitive( wprim ) call addsource_local(0.5d0*dt,& 0.5d0, qt, wCT,& wprim, qt, wnew, xloc, dloc, .true. ) bg(1)%w(ix1,ix2,ix3, 1:nw_phys, n) = wnew(1:nw_phys) ! TBD end do end do end do end do src_active = .true. case default write(unitterm,*)'No such type/not yet implemented sourcesplit=',sourcesplit call mpistop("Error: Unknown type of sourcesplit!") end select if (.not. prior .and. associated(phys_global_source_after)) then call phys_global_source_after(dt, qt, src_active) end if if (src_active) then call getbc(qt,0.d0,ps,iwstart,nwgc,phys_req_diagonal) end if end if end subroutine add_split_source end module mod_source_ffhd