!> Module for handling split source terms (split from the fluxes) module mod_source_mhd 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: rc_fl, radiative_cooling_add_source use mod_global_parameters, only:cmax_global 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 if (.not. qsourcesplit) then !--------------------------------- ! unsplit sources !--------------------------------- do idim = 1, ndim field = gravity_field(wCT, x, idim) wnew(iw_mom(idim)) = wnew(iw_mom(idim)) + qdt * field * wCT(iw_rho) wnew(iw_e) = wnew(iw_e) + qdt * field * wCT(iw_mom(idim)) 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 !--------------------------------- wnew(psi_)=wnew(psi_)*dexp(-qdt*cmax_global*mhd_glm_alpha/minval(dr)) end if end subroutine addsource_local subroutine addsource_nonlocal(qdt, dtfactor, qtC, wCTprim, qt, wnew, x, dx, idir, & qsourcesplit) 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 end subroutine addsource_nonlocal subroutine addsource_compact(qdt, dtfactor, qtC, wCTprim1, wCTprim2, wCTprim3, qt, wnew, x, dx, & qsourcesplit) use mod_global_parameters, only: dt, ndir, smalldouble, courantpar 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 ! .. local .. real(dp) :: laplb_cd2 real(dp) :: Jdir1,Jdir2,Jdir3 integer :: idir real(dp) :: Te_c, rho_c, pth_c real(dp) :: gradT(3), Bmag2, Bmag, Bmag2_safe, bgradT real(dp) :: bhat(3), gradT_perp(3), gradTperp_mag real(dp) :: chi, sig_par, sig_perp real(dp) :: q_sat, f_sat_par, f_sat_perp, tau_par, tau_perp real(dp) :: cf2_tc, cmax2_tc real(dp) :: Qpar_proj, Qperp_proj_k integer :: k_tc if (.not. qsourcesplit) then !--------------------------------- ! unsplit sources !--------------------------------- ! using the compact stencil formulation and adopting constant eta ! > eta*(laplacian of B)_idir added to B_idir ! > B_idir*[eta*(laplacian of B)_idir] added to e_ do idir = 1,ndim laplb_cd2 = & ( & wCTprim1(iw_mag(1)-1+idir,3) & - 2*wCTprim1(iw_mag(1)-1+idir,2) & + wCTprim1(iw_mag(1)-1+idir,1) & ) & / dx(1)**2 + & ( & wCTprim2(iw_mag(1)-1+idir,3) & - 2*wCTprim2(iw_mag(1)-1+idir,2) & + wCTprim2(iw_mag(1)-1+idir,1) & ) & / dx(2)**2 + & ( & wCTprim3(iw_mag(1)-1+idir,3) & - 2*wCTprim3(iw_mag(1)-1+idir,2) & + wCTprim3(iw_mag(1)-1+idir,1) & ) & / dx(3)**2 wnew(iw_mag(1)-1+idir) = wnew(iw_mag(1)-1+idir) + qdt*mhd_eta*laplb_cd2 wnew(iw_e) = wnew(iw_e) + qdt * mhd_eta * laplb_cd2 * wCTprim1(iw_mag(1)-1+idir,2) enddo ! > eta* J**2 added to e Jdir1 = (wCTprim2(iw_mag(3),3) - wCTprim2(iw_mag(3),1)) & / 2.0_dp/dx(2) & - (wCTprim3(iw_mag(2),3) - wCTprim3(iw_mag(2),1)) & / 2.0_dp/dx(3) Jdir2 = (wCTprim3(iw_mag(1),3) - wCTprim3(iw_mag(1),1)) & /2.0_dp/dx(3) & - (wCTprim1(iw_mag(3),3) - wCTprim1(iw_mag(3),1)) & / 2.0_dp/dx(1) Jdir3 = (wCTprim1(iw_mag(2),3) - wCTprim1(iw_mag(2),1)) & /2.0_dp/dx(1) & - (wCTprim2(iw_mag(1),3) - wCTprim2(iw_mag(1),1)) & / 2.0_dp/dx(2) wnew(iw_e) = wnew(iw_e) + qdt*mhd_eta*(Jdir1**2+Jdir2**2+Jdir3**2) ! ---- 1. Local thermodynamic state ---- Te_c = wCTprim1(iw_e,2) / wCTprim1(iw_rho,2) rho_c = wCTprim1(iw_rho,2) pth_c = wCTprim1(iw_e,2) ! ---- 2. Temperature gradient ---- gradT(1) = (wCTprim1(iw_e,3)/wCTprim1(iw_rho,3) - wCTprim1(iw_e,1)/wCTprim1(iw_rho,1)) / (2.d0*dx(1)) gradT(2) = (wCTprim2(iw_e,3)/wCTprim2(iw_rho,3) - wCTprim2(iw_e,1)/wCTprim2(iw_rho,1)) / (2.d0*dx(2)) gradT(3) = (wCTprim3(iw_e,3)/wCTprim3(iw_rho,3) - wCTprim3(iw_e,1)/wCTprim3(iw_rho,1)) / (2.d0*dx(3)) !--------------------------------- ! Anisotropic TC code here !--------------------------------- ! ---- 3. Field direction ---- Bmag2 = 0.0d0 do k_tc = 1, ndir Bmag2 = Bmag2 + wCTprim1(iw_mag(k_tc),2)**2 end do Bmag = sqrt(Bmag2) Bmag2_safe = max(Bmag2, smalldouble**2) do k_tc = 1, ndir bhat(k_tc) = wCTprim1(iw_mag(k_tc),2) * Bmag / Bmag2_safe end do ! ---- 4. Decompose gradT into parallel scalar + perpendicular vector ---- bgradT = bhat(1)*gradT(1) + bhat(2)*gradT(2) + bhat(3)*gradT(3) do k_tc = 1, ndir gradT_perp(k_tc) = gradT(k_tc) - bgradT * bhat(k_tc) end do gradTperp_mag = sqrt(gradT_perp(1)**2 + gradT_perp(2)**2 + gradT_perp(3)**2) ! ---- 5. Conductivities ---- if (tc_kappa_par > 0.0d0) then sig_par = tc_kappa_par else sig_par = tc_kappa0_par * Te_c**2.5d0 end if if (tc_kappa_perp > 0.0d0) then sig_perp = tc_kappa_perp else chi = htc_Cchi * Bmag * Te_c**1.5d0 / rho_c sig_perp = sig_par / (1.0d0 + chi**2) end if ! ---- 6. Saturation and relaxation times (per channel) ---- cf2_tc = (Bmag2 + mhd_gamma*pth_c) / rho_c cmax2_tc = 0.0d0 do k_tc = 1, ndim cmax2_tc = max(cmax2_tc, 0.5d0*(cf2_tc + sqrt(abs( & cf2_tc**2 - 4.0d0*mhd_gamma*pth_c*wCTprim1(iw_mag(k_tc),2)**2/rho_c**2)))) end do q_sat = 1.5d0 * rho_c * (pth_c / rho_c)**1.5d0 f_sat_par = 1.0d0 / (1.0d0 + abs(sig_par * bgradT ) / q_sat) f_sat_perp = 1.0d0 / (1.0d0 + abs(sig_perp * gradTperp_mag) / q_sat) tau_par = max(4.d0*dt, f_sat_par *sig_par *Te_c*courantpar**2*(mhd_gamma-1.0d0) / (pth_c*cmax2_tc)) tau_perp = max(4.d0*dt, f_sat_perp*sig_perp*Te_c*courantpar**2*(mhd_gamma-1.0d0) / (pth_c*cmax2_tc)) ! ---- 7. Decompose current q onto current field direction ---- Qpar_proj = wCTprim1(q_(1),2)*bhat(1) + wCTprim1(q_(2),2)*bhat(2) + wCTprim1(q_(3),2)*bhat(3) ! ---- 8. Relax parallel and perpendicular parts, each at its own rate ---- do k_tc = 1, ndir Qperp_proj_k = wCTprim1(q_(k_tc),2) - Qpar_proj * bhat(k_tc) wnew(q_(k_tc)) = wnew(q_(k_tc)) & - qdt*(f_sat_par *sig_par *bgradT + Qpar_proj)*bhat(k_tc)/tau_par & - qdt*(f_sat_perp*sig_perp*gradT_perp(k_tc) + Qperp_proj_k)/tau_perp end do !--------------------------------- ! Field-aligned TC code here !--------------------------------- else !--------------------------------- ! split sources !--------------------------------- ! Not yet implemented end if 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_mom(2))=u(iw_mom(2))/u(iw_rho) u(iw_mom(3))=u(iw_mom(3))/u(iw_rho) u(iw_e)=gamma_1*(u(iw_e)-0.5_dp*& (u(iw_rho)*(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2)+& u(iw_mag(1))**2+u(iw_mag(2))**2+u(iw_mag(3))**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_mhd