!!****************************************************************************** !! !! This file is part of the AMUN source code, a program to perform !! Newtonian or relativistic magnetohydrodynamical simulations on uniform or !! adaptive mesh. !! !! Copyright (C) 2008-2024 Grzegorz Kowal !! !! This program is free software: you can redistribute it and/or modify !! it under the terms of the GNU General Public License as published by !! the Free Software Foundation, either version 3 of the License, or !! (at your option) any later version. !! !! This program is distributed in the hope that it will be useful, !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! GNU General Public License for more details. !! !! You should have received a copy of the GNU General Public License !! along with this program. If not, see . !! !!****************************************************************************** !! !! module: SCHEMES !! !! The module provides and interface to numerical schemes, subroutines to !! calculate variable increment and one dimensional Riemann solvers. !! !! If you implement a new set of equations, you have to add at least one !! corresponding update_flux subroutine, and one Riemann solver. !! !!****************************************************************************** ! module schemes implicit none ! Rieman solver and state vectors names ! character(len=255) , save :: name_sol = "" character(len=255) , save :: name_sts = "primitive" ! KEPES solver ! logical , save :: kepes = .false. ! 4-vector reconstruction flag ! logical , save :: states_4vec = .false. ! high order fluxes ! logical , save :: high_order_flux = .false. ! interfaces for procedure pointers ! abstract interface subroutine update_flux_iface(dx, q, f) real(kind=8), dimension(:) , intent(in) :: dx real(kind=8), dimension(:,:,:,:) , intent(in) :: q real(kind=8), dimension(:,:,:,:,:), intent(out) :: f end subroutine subroutine riemann_solver_iface(ql, qr, fi) real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi end subroutine end interface ! pointer to the Riemann solver ! procedure(riemann_solver_iface), pointer, save :: riemann_solver => null() ! by default everything is private ! private ! declare public subroutines ! public :: initialize_schemes, finalize_schemes, print_schemes public :: update_flux !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! ! subroutine INITIALIZE_SCHEMES: ! ----------------------------- ! ! Subroutine initiate the module by setting module parameters and subroutine ! pointers. ! ! Arguments: ! ! verbose - a logical flag turning the information printing; ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine initialize_schemes(verbose, status) use equations , only : eqsys, eos use helpers , only : print_message, uppercase use parameters, only : get_parameter implicit none logical, intent(in) :: verbose integer, intent(out) :: status character(len=64) :: solver = "HLL" character(len=64) :: statev = "primitive" character(len=64) :: flux = "off" character(len=256) :: msg character(len=*), parameter :: fmt = "('The selected Riemann solver " // & "is not implemented for the ',a,': ',a,'.'" // & ",a,4x,'Available Riemann solvers are: ', a, '.')" character(len=*), parameter :: loc = 'SCHEMES::initialize_schemes()' !------------------------------------------------------------------------------- ! status = 0 call get_parameter("riemann_solver" , solver) call get_parameter("state_variables", statev) ! depending on the system of equations initialize the module variables ! select case(trim(eqsys)) !--- HYDRODYNAMICS --- ! case("hd", "HD", "hydro", "HYDRO", "hydrodynamic", "HYDRODYNAMIC") select case(trim(eos)) case("iso", "ISO", "isothermal", "ISOTHERMAL") select case(trim(solver)) case("hll", "HLL") name_sol = "HLL" riemann_solver => riemann_hll case("roe", "ROE") name_sol = "ROE" riemann_solver => riemann_hd_iso_roe case("kepes", "KEPES") name_sol = "KEPES" kepes = .true. riemann_solver => riemann_hd_iso_kepes case default if (verbose) then write(msg,fmt) 'isothermal hydrodynamics', & "'" // uppercase(trim(solver)) // "'", & new_line('a'), & "'HLL', 'ROE', 'KEPES'" call print_message(loc, msg) end if status = 1 end select case("adi", "ADI", "adiabatic", "ADIABATIC") select case(trim(solver)) case("hll", "HLL") name_sol = "HLL" riemann_solver => riemann_hll case("hllc", "HLLC") name_sol = "HLLC" riemann_solver => riemann_hd_hllc case("roe", "ROE") name_sol = "ROE" riemann_solver => riemann_hd_adi_roe case("kepes", "KEPES") name_sol = "KEPES" kepes = .true. riemann_solver => riemann_hd_adi_kepes case default if (verbose) then write(msg,fmt) 'adiabatic hydrodynamics', & "'" // uppercase(trim(solver)) // "'", & new_line('a'), & "'HLL', 'HLLC', 'ROE', 'KEPES'" call print_message(loc, msg) end if status = 1 end select end select !--- MAGNETOHYDRODYNAMICS --- ! case("mhd", "MHD", "magnetohydrodynamic", "MAGNETOHYDRODYNAMIC") select case(trim(eos)) case("iso", "ISO", "isothermal", "ISOTHERMAL") select case(trim(solver)) case("hll", "HLL") name_sol = "HLL" riemann_solver => riemann_hll case ("hlld", "HLLD") name_sol = "HLLD" riemann_solver => riemann_mhd_iso_hlld case("roe", "ROE") name_sol = "ROE" riemann_solver => riemann_mhd_iso_roe case("kepes", "KEPES") name_sol = "KEPES" kepes = .true. riemann_solver => riemann_mhd_iso_kepes case default if (verbose) then write(msg,fmt) 'isothermal magnetohydrodynamics', & "'" // uppercase(trim(solver)) // "'", & new_line('a'), & "'HLL', 'HLLD', 'ROE', 'KEPES'" call print_message(loc, msg) end if status = 1 end select case("adi", "ADI", "adiabatic", "ADIABATIC") select case(trim(solver)) case ("hll", "HLL") name_sol = "HLL" riemann_solver => riemann_hll case ("hllc", "HLLC") name_sol = "HLLC" riemann_solver => riemann_mhd_hllc case ("hlld", "HLLD") name_sol = "HLLD" riemann_solver => riemann_mhd_adi_hlld case("roe", "ROE") name_sol = "ROE" riemann_solver => riemann_mhd_adi_roe case("kepes", "KEPES") name_sol = "KEPES" kepes = .true. riemann_solver => riemann_mhd_adi_kepes case default if (verbose) then write(msg,fmt) 'adiabatic magnetohydrodynamics', & "'" // uppercase(trim(solver)) // "'", & new_line('a'), & "'HLL', 'HLLC', 'HLLD', 'ROE', 'KEPES'" call print_message(loc, msg) end if status = 1 end select end select !--- SPECIAL RELATIVITY HYDRODYNAMICS --- ! case("srhd", "SRHD") select case(trim(eos)) case("adi", "ADI", "adiabatic", "ADIABATIC") select case(trim(statev)) case("4vec", "4-vector", "4VEC", "4-VECTOR") name_sts = "4-vector" states_4vec = .true. case default name_sts = "primitive" end select select case(trim(solver)) case ("hll", "HLL") name_sol = "HLL" riemann_solver => riemann_hll case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M") name_sol = "HLLC (Mignone & Bodo 2005)" riemann_solver => riemann_srhd_hllc case default if (verbose) then write(msg,fmt) 'adiabatic special relativity hydrodynamics', & "'" // uppercase(trim(solver)) // "'", & new_line('a'), & "'HLL', 'HLLC'" call print_message(loc, msg) end if status = 1 end select end select !--- SPECIAL RELATIVITY MAGNETOHYDRODYNAMICS --- ! case("srmhd", "SRMHD") select case(trim(eos)) case("adi", "ADI", "adiabatic", "ADIABATIC") select case(trim(statev)) case("4vec", "4-vector", "4VEC", "4-VECTOR") name_sts = "4-vector" states_4vec = .true. case default name_sts = "primitive" end select select case(trim(solver)) case ("hll", "HLL") name_sol = "HLL" riemann_solver => riemann_hll case("hllc", "HLLC", "hllcm", "HLLCM", "hllc-m", "HLLC-M") name_sol = "HLLC (Mignone & Bodo 2006)" riemann_solver => riemann_srmhd_hllc case default if (verbose) then write(msg,fmt) 'adiabatic special relativity magnetohydrodynamics',& "'" // uppercase(trim(solver)) // "'", & new_line('a'), & "'HLL', 'HLLC'" call print_message(loc, msg) end if status = 1 end select end select end select ! flag for higher order flux correction ! call get_parameter("high_order_flux", flux) select case(trim(flux)) case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") high_order_flux = .true. case default high_order_flux = .false. end select !------------------------------------------------------------------------------- ! end subroutine initialize_schemes ! !=============================================================================== ! ! subroutine FINALIZE_SCHEMES: ! --------------------------- ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine finalize_schemes(status) implicit none integer, intent(out) :: status !------------------------------------------------------------------------------- ! status = 0 nullify(riemann_solver) !------------------------------------------------------------------------------- ! end subroutine finalize_schemes ! !=============================================================================== ! ! subroutine PRINT_SCHEMES: ! ------------------------ ! ! Subroutine prints module parameters and settings. ! ! Arguments: ! ! verbose - a logical flag turning the information printing; ! !=============================================================================== ! subroutine print_schemes(verbose) use helpers, only : print_section, print_parameter implicit none logical, intent(in) :: verbose !------------------------------------------------------------------------------- ! if (verbose) then call print_section(verbose, "Schemes") call print_parameter(verbose, "Riemann solver" , name_sol) call print_parameter(verbose, "state variables", name_sts) if (high_order_flux) then call print_parameter(verbose, "high order flux correction", "on" ) else call print_parameter(verbose, "high order flux correction", "off") end if end if !------------------------------------------------------------------------------- ! end subroutine print_schemes ! !=============================================================================== ! ! subroutine UPDATE_FLUX: ! ---------------------- ! ! Subroutine solves the Riemann problem along each direction and calculates ! the numerical fluxes, which are used later to calculate the conserved ! variable increments. ! ! Arguments: ! ! dx - the spatial step; ! q - the array of primitive variables; ! s - the state arrays for primitive variables; ! f - the array of numerical fluxes; ! !=============================================================================== ! subroutine update_flux(dx, q, s, f) use coordinates, only : nn => bcells, nbl, neu use equations , only : relativistic use equations , only : nf, ivars, ivx, ivz implicit none real(kind=8), dimension(:) , intent(in) :: dx real(kind=8), dimension(:,:,:,:) , intent(in) :: q real(kind=8), dimension(:,:,:,:,:,:), intent(inout) :: s real(kind=8), dimension(:,:,:,:,:) , intent(out) :: f integer :: n, i, j, k, l, p real(kind=8) :: vm real(kind=8), dimension(nf,nn,2) :: qi real(kind=8), dimension(nf,nn) :: fi !------------------------------------------------------------------------------- ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ ! in the relativistic case, apply the reconstruction on variables ! using the four-velocities if requested ! if (relativistic .and. states_4vec) then f(:,:,:,:,1) = q(:,:,:,:) #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn do i = 1, nn vm = sqrt(1.0d+00 - sum(q(ivx:ivz,i,j,k)**2)) f(ivx:ivz,i,j,k,1) = f(ivx:ivz,i,j,k,1) / vm end do ! i = 1, nn end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ call reconstruct_interfaces(dx(:), f(:,:,:,:,1), s(:,:,:,:,:,:)) ! convert the states' four-velocities back to velocities ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn do i = 1, nn do l = 1, 2 do p = 1, NDIMS vm = sqrt(1.0d+00 + sum(s(ivx:ivz,i,j,k,l,p)**2)) s(ivx:ivz,i,j,k,l,p) = s(ivx:ivz,i,j,k,l,p) / vm end do ! p = 1, ndims end do ! l = 1, 2 end do ! i = 1, nn end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ else call reconstruct_interfaces(dx(:), q(:,:,:,:), s(:,:,:,:,:,:)) end if f(:,:,:,:,:) = 0.0d+00 ! calculate the flux along the X-direction ! #if NDIMS == 3 do k = nbl, neu #endif /* NDIMS == 3 */ do j = nbl, neu ! copy states to directional lines with proper vector component ordering ! do n = 1, nf qi(n,:,1:2) = s(ivars(1,n),:,j,k,1:2,1) end do ! call one dimensional Riemann solver in order to obtain numerical fluxes ! call numerical_flux(qi(:,:,:), fi(:,:)) ! update the array of fluxes ! do n = 1, nf f(ivars(1,n),:,j,k,1) = fi(n,:) end do end do ! j = nbl, neu ! calculate the flux along the Y direction ! do i = nbl, neu ! copy directional variable vectors to pass to the one dimensional solver ! do n = 1, nf qi(n,:,1:2) = s(ivars(2,n),i,:,k,1:2,2) end do ! call one dimensional Riemann solver in order to obtain numerical fluxes ! call numerical_flux(qi(:,:,:), fi(:,:)) ! update the array of fluxes ! do n = 1, nf f(ivars(2,n),i,:,k,2) = fi(n,:) end do end do ! i = nbl, neu #if NDIMS == 3 end do ! k = nbl, neu #endif /* NDIMS == 3 */ #if NDIMS == 3 ! calculate the flux along the Z direction ! do j = nbl, neu do i = nbl, neu ! copy directional variable vectors to pass to the one dimensional solver ! do n = 1, nf qi(n,:,1:2) = s(ivars(3,n),i,j,:,1:2,3) end do ! call one dimensional Riemann solver in order to obtain numerical fluxes ! call numerical_flux(qi(:,:,:), fi(:,:)) ! update the array of fluxes ! do n = 1, nf f(ivars(3,n),i,j,:,3) = fi(n,:) end do end do ! i = nbl, neu end do ! j = nbl, neu #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine update_flux ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine RECONSTRUCT_INTERFACES: ! --------------------------------- ! ! Subroutine reconstructs the Riemann states along all directions. ! ! Arguments: ! ! dx - the spatial step; ! q - the array of primitive variables; ! qi - the array of reconstructed states (2 in each direction); ! !=============================================================================== ! subroutine reconstruct_interfaces(dx, q, qi) ! include external procedures ! use equations , only : nf use equations , only : positive use interpolations, only : interfaces ! local variables are not implicit by default ! implicit none ! subroutine arguments ! real(kind=8), dimension(:) , intent(in) :: dx real(kind=8), dimension(:,:,:,:) , intent(in) :: q real(kind=8), dimension(:,:,:,:,:,:), intent(out) :: qi ! local variables ! integer :: p ! !------------------------------------------------------------------------------- ! ! interpolate interfaces for all flux variables ! do p = 1, nf call interfaces(positive(p), dx(1:NDIMS), q (p,:,:,:), & qi(p,:,:,:,1:2,1:NDIMS)) end do !------------------------------------------------------------------------------- ! end subroutine reconstruct_interfaces ! !=============================================================================== ! ! subroutine NUMERICAL_FLUX: ! ------------------------- ! ! Subroutine prepares Riemann states and calls the selected Riemann solver ! in order to get the numerical flux. If requested, the resulting flux ! is corrected by higher order correction terms. ! ! Arguments: ! ! q - the array of primitive variables at the Riemann states; ! f - the output array of fluxes; ! !=============================================================================== ! subroutine numerical_flux(q, f) use coordinates, only : nn => bcells use equations , only : nf, ibx, ibp use equations , only : magnetized, cmax implicit none real(kind=8), dimension(:,:,:), intent(inout) :: q real(kind=8), dimension(:,:) , intent(out) :: f integer :: i real(kind=8) :: bx, bp real(kind=8), dimension(nf,nn) :: ql, qr !------------------------------------------------------------------------------- ! if (kepes) then call riemann_solver(q(:,:,1), q(:,:,2), f(:,:)) else ! copy the state vectors ! ql(:,:) = q(:,:,1) qr(:,:) = q(:,:,2) ! obtain the state values for Bx and Psi for the GLM-MHD equations ! if (magnetized) then do i = 1, size(q, 2) bx = 5.0d-01 * ((qr(ibx,i) + ql(ibx,i)) & - (qr(ibp,i) - ql(ibp,i)) / cmax) bp = 5.0d-01 * ((qr(ibp,i) + ql(ibp,i)) & - (qr(ibx,i) - ql(ibx,i)) * cmax) ql(ibx,i) = bx qr(ibx,i) = bx ql(ibp,i) = bp qr(ibp,i) = bp end do end if ! get the numerical fluxes ! call riemann_solver(ql(:,:), qr(:,:), f(:,:)) ! higher order flux corrections ! call higher_order_flux_correction(f) end if !------------------------------------------------------------------------------- ! end subroutine numerical_flux ! !******************************************************************************* ! ! GENERIC RIEMANN SOLVERS ! !******************************************************************************* ! !=============================================================================== ! ! subroutine RIEMANN_HLL: ! ---------------------- ! ! Subroutine solves one dimensional general Riemann problem using ! the Harten-Lax-van Leer (HLL) method. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Harten, A., Lax, P. D., Van Leer, B. ! "On Upstream Differencing and Godunov-Type Schemes for Hyperbolic ! Conservation Laws", ! SIAM Review, 1983, Volume 25, Number 1, pp. 35-61, ! https://doi.org/10.1137/1025002 ! !=============================================================================== ! subroutine riemann_hll(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: sl, sr real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension( 2,nn) :: cl, cr real(kind=8), dimension(nf ) :: wl, wr !------------------------------------------------------------------------------- ! call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl, cl) call fluxspeed(qr, ur, fr, cr) do i = 1, nn sl = min(cl(1,i), cr(1,i)) sr = max(cl(2,i), cr(2,i)) if (sl >= 0.0d+00) then fi(:,i) = fl(:,i) else if (sr <= 0.0d+00) then fi(:,i) = fr(:,i) else ! sl < 0 < sr wl = sl * ul(:,i) - fl(:,i) wr = sr * ur(:,i) - fr(:,i) fi(:,i) = (sl * wr - sr * wl) / (sr - sl) end if ! sl < 0 < sr end do !------------------------------------------------------------------------------- ! end subroutine riemann_hll ! !******************************************************************************* ! ! ISOTHERMAL HYDRODYNAMIC RIEMANN SOLVERS ! !******************************************************************************* ! !=============================================================================== ! ! subroutine RIEMANN_HD_ISO_ROE: ! ----------------------------- ! ! Subroutine solves one dimensional Riemann problem using ! the Roe's method. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Roe, P. L. ! "Approximate Riemann Solvers, Parameter Vectors, and Difference ! Schemes", ! Journal of Computational Physics, 1981, 43, pp. 357-372 ! [2] Toro, E. F., ! "Riemann Solvers and Numerical Methods for Fluid dynamics", ! Springer-Verlag Berlin Heidelberg, 2009 ! !=============================================================================== ! subroutine riemann_hd_iso_roe(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, csnd use equations , only : idn, ivx, ivy, ivz use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: sdl, sdr, sds, cs, ch real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension(nf ) :: qs, lm, al, df logical , save :: first = .true. real(kind=8) , save :: chi real(kind=8), dimension(4,4), save :: rvec, lvec !$omp threadprivate(first, chi, rvec, lvec) !------------------------------------------------------------------------------- ! if (first) then chi = 5.0d-01 / csnd rvec(:,1) = [ 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ] rvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] rvec(:,3) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ] rvec(:,4) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ] lvec(:,1) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] lvec(:,2) = [ chi , 0.0d+00, 0.0d+00, - chi ] lvec(:,3) = [ 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ] lvec(:,4) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ] first = .false. end if call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl) call fluxspeed(qr, ur, fr) do i = 1, nn sdl = sqrt(ql(idn,i)) sdr = sqrt(qr(idn,i)) sds = sdl + sdr qs(idn) = sdl * sdr qs(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds cs = sign(csnd, qs(ivx)) ch = sign(chi , qs(ivx)) lm(1) = qs(ivx) + cs lm(2) = qs(ivx) lm(3) = qs(ivx) lm(4) = qs(ivx) - cs rvec(1,2) = lm(1) rvec(4,2) = lm(4) rvec(1,3) = qs(ivy) rvec(4,3) = qs(ivy) rvec(1,4) = qs(ivz) rvec(4,4) = qs(ivz) lvec(1,1) = - ch * lm(4) lvec(2,1) = - qs(ivy) lvec(3,1) = - qs(ivz) lvec(4,1) = ch * lm(1) lvec(1,2) = ch lvec(4,2) = - ch al = abs(lm) * matmul(lvec, ur(:,i) - ul(:,i)) df = matmul(al, rvec) fi(:,i) = 5.0d-01 * ((fl(:,i) + fr(:,i)) - df) end do !------------------------------------------------------------------------------- ! end subroutine riemann_hd_iso_roe ! !=============================================================================== ! ! subroutine RIEMANN_HD_ISO_KEPES: ! ------------------------------- ! ! Subroutine solves one dimensional isothermal hydrodynamic Riemann problem ! using the entropy stable KEPES method. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Winters, A. R., Czernik, C., Schily, M. B., Gassner, G. J., ! "Entropy stable numerical approximations for the isothermal and ! polytropic Euler equations", ! BIT Numerical Mathematics (2020) 60:791–824, ! https://doi.org/10.1007/s10543-019-00789-w ! !=============================================================================== ! subroutine riemann_hd_iso_kepes(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, csnd, csnd2 use equations , only : idn, ivx, ivy, ivz, imx, imy, imz implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: v2l, v2r, cs real(kind=8) :: dnl, dna, vxa, vya, vza real(kind=8), dimension(nf) :: v, lm, tm logical , save :: first = .true. real(kind=8), dimension(4,4), save :: rm !$omp threadprivate(first, rm) !------------------------------------------------------------------------------- ! if (first) then rm(:,1) = [ 1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] rm(:,2) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ] rm(:,3) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ] rm(:,4) = [ 1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] first = .false. end if do i = 1, nn v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i)) v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i)) dnl = lmean(ql(idn,i), qr(idn,i)) dna = amean(ql(idn,i), qr(idn,i)) vxa = amean(ql(ivx,i), qr(ivx,i)) vya = amean(ql(ivy,i), qr(ivy,i)) vza = amean(ql(ivz,i), qr(ivz,i)) v(idn) = 5.0d-01 * (v2l - v2r) if (qr(idn,i) > ql(idn,i)) then v(idn) = v(idn) + csnd2 * log(qr(idn,i) / ql(idn,i)) else if (ql(idn,i) > qr(idn,i)) then v(idn) = v(idn) - csnd2 * log(ql(idn,i) / qr(idn,i)) end if v(ivx) = qr(ivx,i) - ql(ivx,i) v(ivy) = qr(ivy,i) - ql(ivy,i) v(ivz) = qr(ivz,i) - ql(ivz,i) cs = sign(csnd, vxa) lm(:) = [ vxa + cs, vxa, vxa, vxa - cs ] tm(1) = 5.0d-01 * dnl / csnd2 tm(2) = dnl tm(3) = dnl tm(4) = tm(1) rm(2:4,1) = [ lm(1), vya, vza ] rm(2:4,4) = [ lm(4), vya, vza ] fi(idn,i) = dnl * vxa fi(imx,i) = fi(idn,i) * vxa + csnd2 * dna fi(imy,i) = fi(idn,i) * vya fi(imz,i) = fi(idn,i) * vza fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm)) end do !------------------------------------------------------------------------------- ! end subroutine riemann_hd_iso_kepes ! !******************************************************************************* ! ! ADIABATIC HYDRODYNAMIC RIEMANN SOLVERS ! !******************************************************************************* ! !=============================================================================== ! ! subroutine RIEMANN_HD_HLLC: ! -------------------------- ! ! Subroutine solves one dimensional Riemann problem using the HLLC ! method by Gurski or Li. The HLLC and HLLC-C differ by definitions of ! the tangential components of the velocity and magnetic field. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Toro, E. F., ! "The HLLC Riemann solver", ! Shock Waves, 2019, 29, 8, 1065-1082, ! https://doi.org/10.1007/s00193-019-00912-4 ! !=============================================================================== ! subroutine riemann_hd_hllc(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf use equations , only : idn, ivy, ivz, imx, imy, imz, ien use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: dn, pr real(kind=8) :: sl, sr, sm, slmm, srmm real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension( 2,nn) :: cl, cr real(kind=8), dimension(nf ) :: wl, wr, ui !------------------------------------------------------------------------------- ! call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl, cl) call fluxspeed(qr, ur, fr, cr) do i = 1, nn sl = min(cl(1,i), cr(1,i)) sr = max(cl(2,i), cr(2,i)) if (sl >= 0.0d+00) then fi(:,i) = fl(:,i) else if (sr <= 0.0d+00) then fi(:,i) = fr(:,i) else ! sl < 0 < sr wl = sl * ul(:,i) - fl(:,i) wr = sr * ur(:,i) - fr(:,i) dn = wr(idn) - wl(idn) sm = (wr(imx) - wl(imx)) / dn pr = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn if (sm > 0.0d+00) then slmm = sl - sm ui(idn) = wl(idn) / slmm ui(imx) = ui(idn) * sm ui(imy) = ui(idn) * ql(ivy,i) ui(imz) = ui(idn) * ql(ivz,i) ui(ien) = (wl(ien) + sm * pr) / slmm fi(:,i) = sl * ui - wl else if (sm < 0.0d+00) then srmm = sr - sm ui(idn) = wr(idn) / srmm ui(imx) = ui(idn) * sm ui(imy) = ui(idn) * qr(ivy,i) ui(imz) = ui(idn) * qr(ivz,i) ui(ien) = (wr(ien) + sm * pr) / srmm fi(:,i) = sr * ui - wr else fi(:,i) = (sl * wr - sr * wl) / (sr - sl) end if end if end do !------------------------------------------------------------------------------- ! end subroutine riemann_hd_hllc ! !=============================================================================== ! ! subroutine RIEMANN_HD_ADI_ROE: ! ----------------------------- ! ! Subroutine solves one dimensional Riemann problem using ! the Roe's method. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Roe, P. L. ! "Approximate Riemann Solvers, Parameter Vectors, and Difference ! Schemes", ! Journal of Computational Physics, 1981, 43, pp. 357-372 ! https://doi.org/10.1016/0021-9991(81)90128-5 ! [2] Toro, E. F., ! "Riemann Solvers and Numerical Methods for Fluid dynamics", ! Springer-Verlag Berlin Heidelberg, 2009 ! !=============================================================================== ! subroutine riemann_hd_adi_roe(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, adiabatic_index use equations , only : idn, ivx, ivy, ivz, ipr, ien use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: sdl, sdr, sds, vv, vh, c2, cs, vc, fa, fb, fc, f1, f2, f3 real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension(nf ) :: qs, lm, al, df logical , save :: first = .true. real(kind=8) , save :: adi_m1, adi_m1x real(kind=8), dimension(5,5), save :: rvec, lvec !$omp threadprivate(first, adi_m1, adi_m1x, rvec, lvec) !------------------------------------------------------------------------------- ! if (first) then adi_m1 = adiabatic_index - 1.0d+00 adi_m1x = adiabatic_index / adi_m1 rvec(:,1) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ] rvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] rvec(:,3) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ] rvec(:,4) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ] rvec(:,5) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] lvec(:,1) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] lvec(:,2) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] lvec(:,3) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ] lvec(:,4) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ] lvec(:,5) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] first = .false. end if call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl) call fluxspeed(qr, ur, fr) do i = 1, nn sdl = sqrt(ql(idn,i)) sdr = sqrt(qr(idn,i)) sds = sdl + sdr qs(idn) = sdl * sdr qs(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds qs(ien) = ((ul(ien,i) + ql(ipr,i)) / sdl & + (ur(ien,i) + qr(ipr,i)) / sdr) / sds ! enthalpy vv = sum(qs(ivx:ivz) * qs(ivx:ivz)) vh = 5.0d-01 * vv c2 = adi_m1 * (qs(ien) - vh) cs = sign(sqrt(c2), qs(ivx)) fa = adi_m1 / c2 fb = 5.0d-01 * fa fc = 5.0d-01 / cs vc = cs * qs(ivx) f1 = fb * vh f2 = fc * qs(ivx) f3 = - fb * qs(ivx) lm(1) = qs(ivx) + cs lm(2) = qs(ivx) lm(3) = qs(ivx) lm(4) = qs(ivx) lm(5) = qs(ivx) - cs rvec(1,2) = lm(1) rvec(1,3) = qs(ivy) rvec(1,4) = qs(ivz) rvec(1,5) = qs(ien) + vc rvec(2,2) = qs(ivx) rvec(2,3) = qs(ivy) rvec(2,4) = qs(ivz) rvec(2,5) = vh rvec(3,5) = qs(ivy) rvec(4,5) = qs(ivz) rvec(5,2) = lm(5) rvec(5,3) = qs(ivy) rvec(5,4) = qs(ivz) rvec(5,5) = qs(ien) - vc lvec(1,1) = f1 - f2 lvec(1,2) = f3 + fc lvec(1,3) = - fb * qs(ivy) lvec(1,4) = - fb * qs(ivz) lvec(1,5) = fb lvec(2,1) = - fa * vh + 1.0d+00 lvec(2,2) = fa * qs(ivx) lvec(2,3) = fa * qs(ivy) lvec(2,4) = fa * qs(ivz) lvec(2,5) = - fa lvec(3,1) = - qs(ivy) lvec(4,1) = - qs(ivz) lvec(5,1) = f1 + f2 lvec(5,2) = f3 - fc lvec(5,3) = lvec(1,3) lvec(5,4) = lvec(1,4) lvec(5,5) = fb al = abs(lm) * matmul(lvec, ur(:,i) - ul(:,i)) df = matmul(al, rvec) fi(:,i) = 5.0d-01 * ((fl(:,i) + fr(:,i)) - df) end do !------------------------------------------------------------------------------- ! end subroutine riemann_hd_adi_roe ! !=============================================================================== ! ! subroutine RIEMANN_HD_ADI_KEPES: ! ------------------------------- ! ! Subroutine solves one dimensional adiabatic hydrodynamic Riemann problem ! using the entropy stable KEPES method. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Winters, A. R., Derigs, D., Gassner, G. J., Walch, S. ! "A uniquely defined entropy stable matrix dissipation operator ! for high Mach number ideal MHD and compressible Euler ! simulations", ! Journal of Computational Physics, 2017, 332, pp. 274-289 ! https://doi.org/10.1016/j.jcp.2016.12.006 ! !=============================================================================== ! subroutine riemann_hd_adi_kepes(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, adiabatic_index use equations , only : idn, ivx, ivy, ivz, ipr, imx, imy, imz, ien implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: bl, br, v2l, v2r, cs, vc real(kind=8) :: btl, dnl, prl, bta, dna, vxa, vya, vza, vva, pra, ent real(kind=8), dimension(nf) :: v, lm, tm logical , save :: first = .true. real(kind=8) , save :: adi_m1, adi_m1x real(kind=8), dimension(5,5), save :: rm !$omp threadprivate(first, adi_m1, adi_m1x, rm) !------------------------------------------------------------------------------- ! if (first) then adi_m1 = adiabatic_index - 1.0d+00 adi_m1x = adi_m1 / adiabatic_index rm(1,:) = [ 1.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00 ] rm(2,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] rm(3,:) = [ 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00 ] rm(4,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00 ] rm(5,:) = [ 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 0.0d+00 ] first = .false. end if do i = 1, nn bl = ql(idn,i) / ql(ipr,i) br = qr(idn,i) / qr(ipr,i) v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i)) v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i)) btl = lmean(bl, br) bta = amean(bl, br) dnl = lmean(ql(idn,i), qr(idn,i)) dna = amean(ql(idn,i), qr(idn,i)) vxa = amean(ql(ivx,i), qr(ivx,i)) vya = amean(ql(ivy,i), qr(ivy,i)) vza = amean(ql(ivz,i), qr(ivz,i)) pra = dna / bta prl = dnl / btl vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i)) ent = 1.0d+00 / adi_m1x / btl + vva v(idn) = 5.0d-01 * (bl * v2l - br * v2r) if (qr(idn,i) > ql(idn,i)) then v(idn) = v(idn) + log(qr(idn,i) / ql(idn,i)) else if (ql(idn,i) > qr(idn,i)) then v(idn) = v(idn) - log(ql(idn,i) / qr(idn,i)) end if if (br > bl) then v(idn) = v(idn) + log(br / bl) / adi_m1 else if (bl > br) then v(idn) = v(idn) - log(bl / br) / adi_m1 end if v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i) v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i) v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i) v(ipr) = bl - br cs = sign(sqrt(adiabatic_index * pra / dnl), vxa) vc = vxa * cs lm(1) = vxa + cs lm(2) = vxa lm(3) = vxa lm(4) = vxa lm(5) = vxa - cs tm(1) = 5.0d-01 * dnl / adiabatic_index tm(2) = dnl * adi_m1x tm(3) = pra tm(4) = pra tm(5) = tm(1) rm(2,1) = lm(1) rm(3,1) = vya rm(4,1) = vza rm(2,2) = vxa rm(3,2) = vya rm(4,2) = vza rm(2,5) = lm(5) rm(3,5) = vya rm(4,5) = vza rm(5,:) = [ ent + vc, vva, vya, vza, ent - vc ] fi(idn,i) = dnl * vxa fi(imx,i) = fi(idn,i) * vxa + pra fi(imy,i) = fi(idn,i) * vya fi(imz,i) = fi(idn,i) * vza fi(ien,i) = (prl / adi_m1 + pra + dnl * vva) * vxa fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm)) end do !------------------------------------------------------------------------------- ! end subroutine riemann_hd_adi_kepes ! !******************************************************************************* ! ! RELATIVISTIC ADIABATIC HYDRODYNAMIC RIEMANN SOLVERS ! !******************************************************************************* ! !=============================================================================== ! ! subroutine RIEMANN_SRHD_HLLC: ! ---------------------------- ! ! Subroutine solves one dimensional Riemann problem using ! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC) ! by Mignone & Bodo. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Mignone, A. & Bodo, G. ! "An HLLC Riemann solver for relativistic flows - I. Hydrodynamics", ! Monthly Notices of the Royal Astronomical Society, ! 2005, Volume 364, Pages 126-136 ! !=============================================================================== ! subroutine riemann_srhd_hllc(ql, qr, fi) use algebra , only : quadratic use coordinates, only : nn => bcells use equations , only : nf, idn, imx, imy, imz, ien use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i, nr real(kind=8) :: sl, sr, srml, sm real(kind=8) :: pr, dv real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension( 2,nn) :: cl, cr real(kind=8), dimension(nf ) :: wl, wr, uh, us, fh real(kind=8), dimension(3) :: a real(kind=8), dimension(2) :: x !------------------------------------------------------------------------------- ! call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl, cl) call fluxspeed(qr, ur, fr, cr) do i = 1, nn sl = min(cl(1,i), cr(1,i)) sr = max(cl(2,i), cr(2,i)) if (sl >= 0.0d+00) then fi(:,i) = fl(:,i) else if (sr <= 0.0d+00) then fi(:,i) = fr(:,i) else ! sl < 0 < sr srml = sr - sl wl(:) = sl * ul(:,i) - fl(:,i) wr(:) = sr * ur(:,i) - fr(:,i) uh(:) = ( wr(:) - wl(:)) / srml fh(:) = (sl * wr(:) - sr * wl(:)) / srml wl(ien) = wl(ien) + wl(idn) wr(ien) = wr(ien) + wr(idn) ! prepare the quadratic coefficients (eq. 18 in [1]) ! a(1) = uh(imx) a(2) = - (fh(imx) + uh(ien) + uh(idn)) a(3) = fh(ien) + fh(idn) ! solve the quadratic equation ! nr = quadratic(a(1:3), x(1:2)) ! if Δ < 0, just use the HLL flux ! if (nr < 1) then fi(:,i) = fh(:) else ! get the contact dicontinuity speed ! sm = x(1) ! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux ! if ((sm <= sl) .or. (sm >= sr)) then fi(:,i) = fh(:) else ! calculate total pressure (eq. 17 in [1]) ! pr = fh(imx) - (fh(ien) + fh(idn)) * sm ! if the pressure is negative, use the HLL flux ! if (pr <= 0.0d+00) then fi(:,i) = fh(:) else ! depending in the sign of the contact dicontinuity speed, calculate the proper ! state and corresponding flux ! if (sm > 0.0d+00) then ! calculate the conserved variable vector (eqs. 16 in [1]) ! dv = sl - sm us(idn) = wl(idn) / dv us(imy) = wl(imy) / dv us(imz) = wl(imz) / dv us(ien) = (wl(ien) + pr * sm) / dv us(imx) = (us(ien) + pr) * sm us(ien) = us(ien) - us(idn) ! calculate the flux (eq. 14 in [1]) ! fi(:,i) = fl(:,i) + sl * (us(:) - ul(:,i)) else if (sm < 0.0d+00) then ! calculate the conserved variable vector (eqs. 16 in [1]) ! dv = sr - sm us(idn) = wr(idn) / dv us(imy) = wr(imy) / dv us(imz) = wr(imz) / dv us(ien) = (wr(ien) + pr * sm) / dv us(imx) = (us(ien) + pr) * sm us(ien) = us(ien) - us(idn) ! calculate the flux (eq. 14 in [1]) ! fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i)) else ! intermediate flux is constant across the contact discontinuity and all fluxes ! except the parallel momentum one are zero ! fi(idn,i) = 0.0d+00 fi(imx,i) = pr fi(imy,i) = 0.0d+00 fi(imz,i) = 0.0d+00 fi(ien,i) = 0.0d+00 end if ! sm == 0 end if ! p* < 0 end if ! sl < sm < sr end if ! nr < 1 end if ! sl < 0 < sr end do !------------------------------------------------------------------------------- ! end subroutine riemann_srhd_hllc ! !******************************************************************************* ! ! ISOTHERMAL MAGNETOHYDRODYNAMIC RIEMANN SOLVERS ! !******************************************************************************* ! !=============================================================================== ! ! subroutine RIEMANN_ISO_MHD_HLLD: ! ------------------------------- ! ! Subroutine solves one dimensional Riemann problem using the modified ! isothermal HLLD method by Mignone. The difference is in allowing a density ! jumps across the Alfvén waves, due to the inner intermediate states being ! an average of the states between two Alfvén waves, which also include ! two slow magnetosonic waves. Moreover, due to the different left and right ! intermediate state densities the left and right Alfvén speeds might be ! different too. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Mignone, A., ! "A simple and accurate Riemann solver for isothermal MHD", ! Journal of Computational Physics, 2007, 225, pp. 1427-1441, ! https://doi.org/10.1016/j.jcp.2007.01.033 ! !=============================================================================== ! subroutine riemann_mhd_iso_hlld(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, idn, imx, imy, imz, ibx, iby, ibz, ibp use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm real(kind=8) :: bx, b2, dvl, dvr, dnl, dnr real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension( 2,nn) :: cl, cr real(kind=8), dimension(nf ) :: wl, wr, wcl, wcr, uil, uir !------------------------------------------------------------------------------- ! call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl, cl) call fluxspeed(qr, ur, fr, cr) do i = 1, nn sl = min(cl(1,i), cr(1,i)) sr = max(cl(2,i), cr(2,i)) if (sl >= 0.0d+00) then fi(:,i) = fl(:,i) else if (sr <= 0.0d+00) then fi(:,i) = fr(:,i) else bx = ql(ibx,i) b2 = ql(ibx,i) * qr(ibx,i) wl(:) = sl * ul(:,i) - fl(:,i) wr(:) = sr * ur(:,i) - fr(:,i) sm = (wr(imx) - wl(imx)) / (wr(idn) - wl(idn)) srml = sr - sl slmm = sl - sm srmm = sr - sm uil(idn) = wl(idn) / slmm uir(idn) = wr(idn) / srmm uil(imx) = uil(idn) * sm uir(imx) = uir(idn) * sm uil(ibx) = bx uir(ibx) = bx uil(ibp) = ql(ibp,i) uir(ibp) = qr(ibp,i) sml = sm - abs(bx) / sqrt(uil(idn)) smr = sm + abs(bx) / sqrt(uir(idn)) dvl = slmm * wl(idn) - b2 dvr = srmm * wr(idn) - b2 dnr = uir(idn) / dvr dnl = uil(idn) / dvl ! check degeneracy Sl* -> Sl or Sr* -> Sr ! if (sml > sl .and. smr < sr .and. min(dvl, dvr) > 0.0d+00) then if (sml > 0.0d+00) then ! sl* > 0 uil(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl fi(:,i) = sl * uil(:) - wl(:) else if (smr < 0.0d+00) then ! sr* < 0 uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr fi(:,i) = sr * uir(:) - wr(:) else ! sl* <= 0 <= sr* if (smr > sml) then uil(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr wcl(:) = (sml - sl) * uil(:) + wl(:) wcr(:) = (smr - sr) * uir(:) + wr(:) dvl = smr - sml fi(:,i) = (sml * wcr(:) - smr * wcl(:)) / dvl else ! Bx = 0 -> Sₘ = 0 fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if end if ! sl* <= 0 <= sr* else if (sml > sl .and. dvl > 0.0d+00) then ! sr* > sr uil(imy) = dnl * (slmm * wl(imy) - bx * wl(iby)) uil(imz) = dnl * (slmm * wl(imz) - bx * wl(ibz)) uil(iby) = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl uil(ibz) = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl if (sml > 0.0d+00) then ! sl* > 0 fi(:,i) = sl * uil(:) - wl(:) else if (sr > sml) then wcl(:) = (sml - sl) * uil(:) + wl(:) dvl = sr - sml fi(:,i) = (sml * wr(:) - sr * wcl(:)) / dvl else ! Bx = 0 -> Sₘ = 0 fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if ! sl* < 0 else if (smr < sr .and. dvr > 0.0d+00) then ! sl* < sl uir(imy) = dnr * (srmm * wr(imy) - bx * wr(iby)) uir(imz) = dnr * (srmm * wr(imz) - bx * wr(ibz)) uir(iby) = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr uir(ibz) = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr if (smr < 0.0d+00) then ! sr* < 0 fi(:,i) = sr * uir(:) - wr(:) else if (smr > sl) then wcr(:) = (smr - sr) * uir(:) + wr(:) dvr = smr - sl fi(:,i) = (sl * wcr(:) - smr * wl(:)) / dvr else ! Bx = 0 -> Sₘ = 0 fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if ! sr* > 0 else ! sl* < sl & sr* > sr fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if end if end do !------------------------------------------------------------------------------- ! end subroutine riemann_mhd_iso_hlld ! !=============================================================================== ! ! subroutine RIEMANN_MHD_ISO_ROE: ! ------------------------------ ! ! Subroutine solves one dimensional Riemann problem using ! the Roe's method. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Stone, J. M. & Gardiner, T. A., ! "ATHENA: A New Code for Astrophysical MHD", ! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177 ! [2] Toro, E. F., ! "Riemann Solvers and Numerical Methods for Fluid dynamics", ! Springer-Verlag Berlin Heidelberg, 2009 ! !=============================================================================== ! subroutine riemann_mhd_iso_roe(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, csnd2 use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, & ibx, iby, ibz, ibp use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: sdl, sdr, sds, xx, yy, bty, btz, br, br2 real(kind=8) :: cc, ca, cs, cf, cc2, ca2, cs2, cf2, alf, als real(kind=8) :: vqstr, sqrtd, qs, qf, norm, css, cff real(kind=8) :: as_prime, af_prime, as, af, aspb, afpb real(kind=8) :: f1, f2, f3, f4 real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension(nf ) :: qi real(kind=8), dimension(6 ) :: lm, al, df logical , save :: first = .true. integer , dimension(6) , save :: ivs real(kind=8), dimension(6,6), save :: rvec, lvec !$omp threadprivate(first, ivs, rvec, lvec) !------------------------------------------------------------------------------- ! if (first) then rvec(:,:) = 0.0d+00 lvec(:,:) = 0.0d+00 ivs = [ idn, imx, imy, imz, iby, ibz ] first = .false. end if call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl) call fluxspeed(qr, ur, fr) do i = 1, nn sdl = sqrt(ql(idn,i)) sdr = sqrt(qr(idn,i)) sds = sdl + sdr qi(idn) = sdl * sdr qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds qi(ibx) = ql(ibx,i) qi(ibp) = ql(ibp,i) f1 = ql(iby,i) - qr(iby,i) f2 = ql(ibz,i) - qr(ibz,i) xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds) yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn) br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz) br = sqrt(br2) if (br2 > 0.0d+00) then bty = qi(iby) / br btz = qi(ibz) / br else bty = 1.0d+00 btz = 0.0d+00 end if cc2 = csnd2 + xx ca2 = qi(ibx) * qi(ibx) / qi(idn) f1 = br2 / qi(idn) f2 = ca2 + f1 f3 = f2 - cc2 f4 = sqrt(f3 * f3 + 4.0d+00 * cc2 * f1) cf2 = 5.0d-01 * (f2 + cc2 + f4) cs2 = cc2 * ca2 / cf2 if (cs2 < cc2 .and. cc2 < cf2) then f1 = (cc2 - cs2) / (cf2 - cs2) alf = sqrt(f1) als = sqrt(1.0d+00 - f1) else if (cc2 >= cf2) then alf = 1.0d+00 als = 0.0d+00 else alf = 0.0d+00 als = 1.0d+00 end if cc = sqrt(cc2) ca = sqrt(ca2) cf = sign(sqrt(cf2), qi(ivx)) cs = sign(sqrt(cs2), qi(ivx)) lm(1) = qi(ivx) + cf lm(2) = qi(ivx) + ca lm(3) = qi(ivx) + cs lm(4) = qi(ivx) - cs lm(5) = qi(ivx) - ca lm(6) = qi(ivx) - cf sqrtd = sqrt(qi(idn)) f1 = sqrt(yy) qf = cf * sign(alf, qi(ibx)) / f1 qs = cs * sign(als, qi(ibx)) / f1 f1 = cc / (f1 * sqrtd) af_prime = alf * f1 as_prime = als * f1 rvec(1,1) = alf rvec(1,2) = alf * lm(1) rvec(1,3) = alf * qi(ivy) - qs * bty rvec(1,4) = alf * qi(ivz) - qs * btz rvec(1,5) = as_prime * bty rvec(1,6) = as_prime * btz rvec(6,1) = alf rvec(6,2) = alf * lm(6) rvec(6,3) = alf * qi(ivy) + qs * bty rvec(6,4) = alf * qi(ivz) + qs * btz rvec(6,5) = rvec(1,5) rvec(6,6) = rvec(1,6) rvec(3,1) = als rvec(3,2) = als * lm(3) rvec(3,3) = als * qi(ivy) + qf * bty rvec(3,4) = als * qi(ivz) + qf * btz rvec(3,5) = - af_prime * bty rvec(3,6) = - af_prime * btz rvec(4,1) = als rvec(4,2) = als * lm(4) rvec(4,3) = als * qi(ivy) - qf * bty rvec(4,4) = als * qi(ivz) - qf * btz rvec(4,5) = rvec(3,5) rvec(4,6) = rvec(3,6) f1 = sign(1.0d+00 / sqrtd, qi(ibx)) rvec(2,3) = btz rvec(2,4) = - bty rvec(2,5) = - btz * f1 rvec(2,6) = bty * f1 rvec(5,3) = - rvec(2,3) rvec(5,4) = - rvec(2,4) rvec(5,5) = rvec(2,5) rvec(5,6) = rvec(2,6) norm = 2.0d+00 * cc2 cff = alf * cf / norm css = als * cs / norm norm = norm / yy qf = qf / norm qs = qs / norm f1 = qi(idn) / norm af = af_prime * f1 as = as_prime * f1 f1 = br / norm afpb = af_prime * f1 aspb = as_prime * f1 vqstr = (qi(ivy) * bty + qi(ivz) * btz) f1 = qs * vqstr lvec(1,1) = - cff * lm(6) + f1 - aspb lvec(1,2) = cff lvec(1,3) = - qs * bty lvec(1,4) = - qs * btz lvec(1,5) = as * bty lvec(1,6) = as * btz lvec(6,1) = cff * lm(1) - f1 - aspb lvec(6,2) = - lvec(1,2) lvec(6,3) = - lvec(1,3) lvec(6,4) = - lvec(1,4) lvec(6,5) = lvec(1,5) lvec(6,6) = lvec(1,6) f1 = qf * vqstr lvec(3,1) = - css * lm(4) - f1 + afpb lvec(3,2) = css lvec(3,3) = qf * bty lvec(3,4) = qf * btz lvec(3,5) = - af * bty lvec(3,6) = - af * btz lvec(4,1) = css * lm(3) + f1 + afpb lvec(4,2) = - lvec(3,2) lvec(4,3) = - lvec(3,3) lvec(4,4) = - lvec(3,4) lvec(4,5) = lvec(3,5) lvec(4,6) = lvec(3,6) f1 = sign(5.0d-01 * sqrtd, qi(ibx)) lvec(2,1) = - 5.0d-01 * (qi(ivy) * btz - qi(ivz) * bty) lvec(2,3) = 5.0d-01 * btz lvec(2,4) = - 5.0d-01 * bty lvec(2,5) = - f1 * btz lvec(2,6) = f1 * bty lvec(5,1) = - lvec(2,1) lvec(5,3) = - lvec(2,3) lvec(5,4) = - lvec(2,4) lvec(5,5) = lvec(2,5) lvec(5,6) = lvec(2,6) al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i)) df = matmul(al, rvec) fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:)) fi(ibx,i) = fl(ibx,i) fi(ibp,i) = fl(ibp,i) end do !------------------------------------------------------------------------------- ! end subroutine riemann_mhd_iso_roe ! !=============================================================================== ! ! subroutine RIEMANN_MHD_ISO_KEPES: ! -------------------------------- ! ! Subroutine solves one dimensional isothermal magnetohydrodynamic ! Riemann problem using the entropy stable KEPES method. The method is ! a modification of the method described in [1] for the isothermal equation ! of state. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M., ! "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic ! field divergence diminishing ideal magnetohydrodynamics equations", ! Journal of Computational Physics, 2018, 364, pp. 420-467, ! https://doi.org/10.1016/j.jcp.2018.03.002 ! [2] Winters, A. R., Czernik, C., Schily, M. B., Gassner, G. J., ! "Entropy stable numerical approximations for the isothermal and ! polytropic Euler equations", ! BIT Numerical Mathematics (2020) 60:791–824, ! https://doi.org/10.1007/s10543-019-00789-w ! !=============================================================================== ! subroutine riemann_mhd_iso_kepes(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, csnd, csnd2, ch => cglm use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, & ibx, iby, ibz, ibp implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: dna, pta, vxa, vya, vza, bxa, bya, bza, bpa real(kind=8) :: dnl, v2l, v2r, b2l, b2r, bp2, das, daf, csq real(kind=8) :: ca, cs, cf, ca2, cs2, cf2, x2, x3 real(kind=8) :: alf, als, f1, f2, f3, f4 real(kind=8), dimension(nf) :: v, lm, tm logical , save :: first = .true. real(kind=8), dimension(8,8), save :: rm !$omp threadprivate(first, rm) !------------------------------------------------------------------------------- ! if (first) then rm(:,:) = 0.0d+00 rm(5,4) = 1.0d+00 rm(8,4) = 1.0d+00 rm(5,5) = 1.0d+00 rm(8,5) = - 1.0d+00 first = .false. end if do i = 1, nn v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i)) v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i)) b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i)) b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i)) dnl = lmean(ql(idn,i), qr(idn,i)) dna = amean(ql(idn,i), qr(idn,i)) vxa = amean(ql(ivx,i), qr(ivx,i)) vya = amean(ql(ivy,i), qr(ivy,i)) vza = amean(ql(ivz,i), qr(ivz,i)) bxa = amean(ql(ibx,i), qr(ibx,i)) bya = amean(ql(iby,i), qr(iby,i)) bza = amean(ql(ibz,i), qr(ibz,i)) bpa = amean(ql(ibp,i), qr(ibp,i)) pta = 5.0d-01 * amean(b2l, b2r) + csnd2 * dna v(idn) = 5.0d-01 * (v2l - v2r) if (qr(idn,i) > ql(idn,i)) then v(idn) = v(idn) + csnd2 * log(qr(idn,i) / ql(idn,i)) else if (ql(idn,i) > qr(idn,i)) then v(idn) = v(idn) - csnd2 * log(ql(idn,i) / qr(idn,i)) end if v(ivx) = qr(ivx,i) - ql(ivx,i) v(ivy) = qr(ivy,i) - ql(ivy,i) v(ivz) = qr(ivz,i) - ql(ivz,i) v(ibx) = qr(ibx,i) - ql(ibx,i) v(iby) = qr(iby,i) - ql(iby,i) v(ibz) = qr(ibz,i) - ql(ibz,i) v(ibp) = qr(ibp,i) - ql(ibp,i) bp2 = bya * bya + bza * bza if (bp2 > 0.0d+00) then f1 = sqrt(bp2) x2 = bya / f1 x3 = bza / f1 else x2 = 0.0d+00 x3 = 0.0d+00 end if ca2 = bxa * bxa / dnl f1 = bp2 / dnl f2 = ca2 + f1 f3 = f2 - csnd2 f4 = sqrt(f3 * f3 + 4.0d+00 * csnd2 * f1) cf2 = 5.0d-01 * (f2 + csnd2 + f4) cs2 = csnd2 * ca2 / cf2 if (cs2 < csnd2 .and. csnd2 < cf2) then f1 = (csnd2 - cs2) / (cf2 - cs2) alf = sqrt(f1) als = sqrt(1.0d+00 - f1) else if (csnd2 >= cf2) then alf = 1.0d+00 als = 0.0d+00 else alf = 0.0d+00 als = 1.0d+00 end if cf = sign(sqrt(cf2), vxa) ca = sqrt(ca2) cs = sign(sqrt(cs2), vxa) lm(1) = vxa + cf lm(2) = vxa + ca lm(3) = vxa + cs lm(4) = vxa + ch lm(5) = vxa - ch lm(6) = vxa - cs lm(7) = vxa - ca lm(8) = vxa - cf tm(1) = 5.0d-01 / (dnl * csnd2) tm(2) = 5.0d-01 / dnl**2 tm(3) = tm(1) tm(4) = 5.0d-01 tm(5) = tm(4) tm(6) = tm(1) tm(7) = tm(2) tm(8) = tm(1) das = dnl * als daf = dnl * alf csq = csnd * sqrt(dnl) f1 = daf f2 = sign(das, bxa) * cs f3 = als * csq rm(1,1) = f1 rm(2,1) = f1 * (vxa + cf) rm(3,1) = f1 * vya - f2 * x2 rm(4,1) = f1 * vza - f2 * x3 rm(6,1) = f3 * x2 rm(7,1) = f3 * x3 rm(1,8) = f1 rm(2,8) = f1 * (vxa - cf) rm(3,8) = f1 * vya + f2 * x2 rm(4,8) = f1 * vza + f2 * x3 rm(6,8) = rm(6,1) rm(7,8) = rm(7,1) f1 = dnl * sqrt(dna) rm(3,2) = f1 * x3 rm(4,2) = - f1 * x2 rm(6,2) = - dnl * x3 rm(7,2) = dnl * x2 rm(3,7) = - rm(3,2) rm(4,7) = - rm(4,2) rm(6,7) = rm(6,2) rm(7,7) = rm(7,2) f1 = das f2 = sign(daf, bxa) * cf f3 = - alf * csq rm(1,3) = f1 rm(2,3) = f1 * (vxa + cs) rm(3,3) = f1 * vya + f2 * x2 rm(4,3) = f1 * vza + f2 * x3 rm(6,3) = f3 * x2 rm(7,3) = f3 * x3 rm(1,6) = f1 rm(2,6) = f1 * (vxa - cs) rm(3,6) = f1 * vya - f2 * x2 rm(4,6) = f1 * vza - f2 * x3 rm(6,6) = rm(6,3) rm(7,6) = rm(7,3) fi(idn,i) = dnl * vxa fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta fi(imy,i) = fi(idn,i) * vya - bxa * bya fi(imz,i) = fi(idn,i) * vza - bxa * bza fi(ibx,i) = ch * bpa fi(iby,i) = vxa * bya - bxa * vya fi(ibz,i) = vxa * bza - bxa * vza fi(ibp,i) = ch * bxa fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm)) end do !------------------------------------------------------------------------------- ! end subroutine riemann_mhd_iso_kepes ! !******************************************************************************* ! ! ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS ! !******************************************************************************* ! !=============================================================================== ! ! subroutine RIEMANN_MHD_HLLC: ! --------------------------- ! ! Subroutine solves one dimensional Riemann problem using the HLLC method, ! by Toro. In the HLLC method the tangential components of the velocity are ! discontinuous actoss the contact dicontinuity. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Toro, E. F., Spruce, M., & Speares, W. ! "Restoration of the contact surface in the HLL-Riemann solver", ! Shock Waves, 1994, Volume 4, Issue 1, pp. 25-34 ! !=============================================================================== ! subroutine riemann_mhd_hllc(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, idn, ivy, ivz, imx, imy, imz, ien, & ibx, iby, ibz, ibp use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: dn, pt, vy, vz, bx, by, bz, vb, b2 real(kind=8) :: sl, sr, sm, srml, slmm, srmm real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension( 2,nn) :: cl, cr real(kind=8), dimension(nf ) :: wl, wr, ui !------------------------------------------------------------------------------- ! call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl, cl) call fluxspeed(qr, ur, fr, cr) do i = 1, nn sl = min(cl(1,i), cr(1,i)) sr = max(cl(2,i), cr(2,i)) if (sl >= 0.0d+00) then fi(:,i) = fl(:,i) else if (sr <= 0.0d+00) then fi(:,i) = fr(:,i) else ! sl < 0 < sr bx = ql(ibx,i) b2 = ql(ibx,i) * qr(ibx,i) srml = sr - sl wl(:) = sl * ul(:,i) - fl(:,i) wr(:) = sr * ur(:,i) - fr(:,i) dn = wr(idn) - wl(idn) sm = (wr(imx) - wl(imx)) / dn pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2 if (b2 > 0.0d+00) then vy = (wr(imy) - wl(imy)) / dn vz = (wr(imz) - wl(imz)) / dn by = (wr(iby) - wl(iby)) / srml bz = (wr(ibz) - wl(ibz)) / srml vb = sm * bx + vy * by + vz * bz if (sm > 0.0d+00) then ! sm > 0 slmm = sl - sm ui(idn) = wl(idn) / slmm ui(imx) = ui(idn) * sm ui(imy) = ui(idn) * vy ui(imz) = ui(idn) * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ql(ibp,i) ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm fi(:,i) = sl * ui(:) - wl(:) else if (sm < 0.0d+00) then ! sm < 0 srmm = sr - sm ui(idn) = wr(idn) / srmm ui(imx) = ui(idn) * sm ui(imy) = ui(idn) * vy ui(imz) = ui(idn) * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = qr(ibp,i) ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm fi(:,i) = sr * ui(:) - wr(:) else ! sm = 0 ! when Sₘ = 0 all variables are continuous, therefore the flux reduces ! to the HLL one ! fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if ! sm = 0 else ! Bₓ = 0 if (sm > 0.0d+00) then ! sm > 0 slmm = sl - sm ui(idn) = wl(idn) / slmm ui(imx) = ui(idn) * sm ui(imy) = ui(idn) * ql(ivy,i) ui(imz) = ui(idn) * ql(ivz,i) ui(ibx) = 0.0d+00 ui(iby) = wl(iby) / slmm ui(ibz) = wl(ibz) / slmm ui(ibp) = ql(ibp,i) ui(ien) = (wl(ien) + sm * pt) / slmm fi(:,i) = sl * ui(:) - wl(:) else if (sm < 0.0d+00) then ! sm < 0 srmm = sr - sm ui(idn) = wr(idn) / srmm ui(imx) = ui(idn) * sm ui(imy) = ui(idn) * qr(ivy,i) ui(imz) = ui(idn) * qr(ivz,i) ui(ibx) = 0.0d+00 ui(iby) = wr(iby) / srmm ui(ibz) = wr(ibz) / srmm ui(ibp) = qr(ibp,i) ui(ien) = (wr(ien) + sm * pt) / srmm fi(:,i) = sr * ui(:) - wr(:) else ! sm = 0 ! when Sₘ = 0 all variables are continuous, therefore the flux reduces ! to the HLL one ! fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if ! sm = 0 end if ! Bₓ = 0 end if ! sl < 0 < sr end do !------------------------------------------------------------------------------- ! end subroutine riemann_mhd_hllc ! !=============================================================================== ! ! subroutine RIEMANN_ADI_MHD_HLLD: ! ------------------------------- ! ! Subroutine solves one dimensional Riemann problem using the adiabatic HLLD ! method described by Miyoshi & Kusano. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Miyoshi, T. & Kusano, K., ! "A multi-state HLL approximate Riemann solver for ideal ! magnetohydrodynamics", ! Journal of Computational Physics, 2005, 208, pp. 315-344 ! !=============================================================================== ! subroutine riemann_mhd_adi_hlld(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, idn, imx, imy, imz, ien, ibx, iby, ibz, ibp use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: sl, sr, sm, srml, slmm, srmm real(kind=8) :: dn, bx, b2, pt, vy, vz, by, bz, vb real(kind=8) :: dnl, dnr, cal, car, sml, smr real(kind=8) :: dv, dvl, dvr real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension( 2,nn) :: cl, cr real(kind=8), dimension(nf ) :: wl, wr, wcl, wcr, ui !------------------------------------------------------------------------------- ! call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl, cl) call fluxspeed(qr, ur, fr, cr) do i = 1, nn sl = min(cl(1,i), cr(1,i)) sr = max(cl(2,i), cr(2,i)) if (sl >= 0.0d+00) then ! sl ≥ 0 fi(:,i) = fl(:,i) else if (sr <= 0.0d+00) then ! sr ≤ 0 fi(:,i) = fr(:,i) else ! sl < 0 < sr bx = ql(ibx,i) b2 = ql(ibx,i) * qr(ibx,i) srml = sr - sl wl(:) = sl * ul(:,i) - fl(:,i) wr(:) = sr * ur(:,i) - fr(:,i) dn = wr(idn) - wl(idn) sm = (wr(imx) - wl(imx)) / dn slmm = sl - sm srmm = sr - sm dnl = wl(idn) / slmm dnr = wr(idn) / srmm pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2 cal = abs(bx) / sqrt(dnl) car = abs(bx) / sqrt(dnr) sml = sm - cal smr = sm + car ! calculate division factors (also used to determine degeneracies) ! dvl = slmm * wl(idn) - b2 dvr = srmm * wr(idn) - b2 ! check degeneracy Sl* -> Sl or Sr* -> Sr ! if (sml > sl .and. smr < sr .and. min(dvl, dvr) > 0.0d+00) then ! no degeneracy ! choose the correct state depending on the speed signs ! if (sml > 0.0d+00) then ! sl* > 0 ! primitive variables for the outer left intermediate state ! vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl vb = sm * bx + vy * by + vz * bz ! conservative variables for the outer left intermediate state ! ui(idn) = dnl ui(imx) = dnl * sm ui(imy) = dnl * vy ui(imz) = dnl * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ul(ibp,i) ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm ! the outer left intermediate flux ! fi(:,i) = sl * ui(:) - wl(:) else if (smr < 0.0d+00) then ! sr* < 0 ! primitive variables for the outer right intermediate state ! vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr vb = sm * bx + vy * by + vz * bz ! conservative variables for the outer right intermediate state ! ui(idn) = dnr ui(imx) = dnr * sm ui(imy) = dnr * vy ui(imz) = dnr * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ur(ibp,i) ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm ! the outer right intermediate flux ! fi(:,i) = sr * ui(:) - wr(:) else ! sl* < 0 < sr* ! separate cases with non-zero and zero Bx ! if (b2 > 0.0d+00) then ! primitive variables for the outer left intermediate state ! vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl vb = sm * bx + vy * by + vz * bz ! conservative variables for the outer left intermediate state ! ui(idn) = dnl ui(imx) = dnl * sm ui(imy) = dnl * vy ui(imz) = dnl * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ul(ibp,i) ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm ! vector of the left-going Alfvén wave ! wcl(:) = (sml - sl) * ui(:) + wl(:) ! primitive variables for the outer right intermediate state ! vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr vb = sm * bx + vy * by + vz * bz ! conservative variables for the outer right intermediate state ! ui(idn) = dnr ui(imx) = dnr * sm ui(imy) = dnr * vy ui(imz) = dnr * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ur(ibp,i) ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm ! vector of the right-going Alfvén wave ! wcr(:) = (smr - sr) * ui(:) + wr(:) ! prepare constant primitive variables of the intermediate states ! dv = abs(bx) * (sqrt(dnr) + sqrt(dnl)) vy = (wcr(imy) - wcl(imy)) / dv vz = (wcr(imz) - wcl(imz)) / dv dv = car + cal by = (wcr(iby) - wcl(iby)) / dv bz = (wcr(ibz) - wcl(ibz)) / dv vb = sm * bx + vy * by + vz * bz ! choose the correct state depending on the sign of contact discontinuity ! advection speed ! if (sm > 0.0d+00) then ! sm > 0 ! conservative variables for the inmost left intermediate state ! ui(idn) = dnl ui(imx) = dnl * sm ui(imy) = dnl * vy ui(imz) = dnl * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ul(ibp,i) ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal ! the inmost left intermediate flux ! fi(:,i) = sml * ui(:) - wcl(:) else if (sm < 0.0d+00) then ! sm < 0 ! conservative variables for the inmost right intermediate state ! ui(idn) = dnr ui(imx) = dnr * sm ui(imy) = dnr * vy ui(imz) = dnr * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ur(ibp,i) ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car ! the inmost right intermediate flux ! fi(:,i) = smr * ui(:) - wcr(:) else ! sm = 0 ! in the case when Sₘ = 0 and Bₓ² > 0, all variables are continuous, therefore ! the flux can be averaged from the Alfvén waves using a simple HLL formula; ! fi(:,i) = (sml * wcr(:) - smr * wcl(:)) / dv end if ! sm = 0 else ! Bx = 0 if (sm > 0.0d+00) then ! primitive variables for the outer left intermediate state ! vy = slmm * wl(imy) / dvl vz = slmm * wl(imz) / dvl by = wl(idn) * wl(iby) / dvl bz = wl(idn) * wl(ibz) / dvl vb = vy * by + vz * bz ! conservative variables for the outer left intermediate state ! ui(idn) = dnl ui(imx) = dnl * sm ui(imy) = dnl * vy ui(imz) = dnl * vz ui(ibx) = 0.0d+00 ui(iby) = by ui(ibz) = bz ui(ibp) = ul(ibp,i) ui(ien) = (wl(ien) + sm * pt) / slmm ! the outer left intermediate flux ! fi(:,i) = sl * ui(:) - wl(:) else if (sm < 0.0d+00) then ! primitive variables for the outer right intermediate state ! vy = ( srmm * wr(imy)) / dvr vz = ( srmm * wr(imz)) / dvr by = (wr(idn) * wr(iby)) / dvr bz = (wr(idn) * wr(ibz)) / dvr vb = vy * by + vz * bz ! conservative variables for the outer right intermediate state ! ui(idn) = dnr ui(imx) = dnr * sm ui(imy) = dnr * vy ui(imz) = dnr * vz ui(ibx) = 0.0d+00 ui(iby) = by ui(ibz) = bz ui(ibp) = ur(ibp,i) ui(ien) = (wr(ien) + sm * pt) / srmm ! the outer right intermediate flux ! fi(:,i) = sr * ui(:) - wr(:) else ! Sm = 0 ! since Bx = 0 and Sm = 0, then revert to the HLL flux ! fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if end if ! Bx = 0 end if ! sl* < 0 < sr* else ! some degeneracies ! separate degeneracies ! if (sml > sl .and. dvl > 0.0d+00) then ! sr* > sr ! primitive variables for the outer left intermediate state ! vy = ( slmm * wl(imy) - bx * wl(iby)) / dvl vz = ( slmm * wl(imz) - bx * wl(ibz)) / dvl by = (wl(idn) * wl(iby) - bx * wl(imy)) / dvl bz = (wl(idn) * wl(ibz) - bx * wl(imz)) / dvl vb = sm * bx + vy * by + vz * bz ! conservative variables for the outer left intermediate state ! ui(idn) = dnl ui(imx) = dnl * sm ui(imy) = dnl * vy ui(imz) = dnl * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ul(ibp,i) ui(ien) = (wl(ien) + sm * pt - bx * vb) / slmm ! choose the correct state depending on the speed signs ! if (sml > 0.0d+00) then ! sl* > 0 ! the outer left intermediate flux ! fi(:,i) = sl * ui(:) - wl(:) else ! sl* <= 0 ! separate cases with non-zero and zero Bx ! if (b2 > 0.0d+00) then ! vector of the left-going Alfvén wave ! wcl(:) = (sml - sl) * ui(:) + wl(:) ! primitive variables for the inner left intermediate state ! dv = srmm * dnr + cal * dnl vy = (wr(imy) - wcl(imy)) / dv vz = (wr(imz) - wcl(imz)) / dv dv = sr - sml by = (wr(iby) - wcl(iby)) / dv bz = (wr(ibz) - wcl(ibz)) / dv vb = sm * bx + vy * by + vz * bz ! conservative variables for the inner left intermediate state ! ui(idn) = dnl ui(imx) = dnl * sm ui(imy) = dnl * vy ui(imz) = dnl * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ul(ibp,i) ui(ien) = - (wcl(ien) + sm * pt - bx * vb) / cal ! choose the correct state depending on the sign of contact discontinuity ! advection speed ! if (sm >= 0.0d+00) then ! sm >= 0 ! the inner left intermediate flux ! fi(:,i) = sml * ui(:) - wcl(:) else ! sm < 0 ! vector of the left-going Alfvén wave ! wcr(:) = (sm - sml) * ui(:) + wcl(:) ! calculate the average flux over the right inner intermediate state ! fi(:,i) = (sm * wr(:) - sr * wcr(:)) / srmm end if ! sm < 0 else ! bx = 0 ! no Alfvén wave, so revert to the HLL flux ! fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if end if ! sl* < 0 else if (smr < sr .and. dvr > 0.0d+00) then ! sl* < sl ! primitive variables for the outer right intermediate state ! vy = ( srmm * wr(imy) - bx * wr(iby)) / dvr vz = ( srmm * wr(imz) - bx * wr(ibz)) / dvr by = (wr(idn) * wr(iby) - bx * wr(imy)) / dvr bz = (wr(idn) * wr(ibz) - bx * wr(imz)) / dvr vb = sm * bx + vy * by + vz * bz ! conservative variables for the outer right intermediate state ! ui(idn) = dnr ui(imx) = dnr * sm ui(imy) = dnr * vy ui(imz) = dnr * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ur(ibp,i) ui(ien) = (wr(ien) + sm * pt - bx * vb) / srmm ! choose the correct state depending on the speed signs ! if (smr < 0.0d+00) then ! sr* < 0 ! the outer right intermediate flux ! fi(:,i) = sr * ui(:) - wr(:) else ! sr* > 0 ! separate cases with non-zero and zero Bx ! if (b2 > 0.0d+00) then ! vector of the right-going Alfvén wave ! wcr(:) = (smr - sr) * ui(:) + wr(:) ! primitive variables for the inner right intermediate state ! dv = slmm * dnl - car * dnr vy = (wl(imy) - wcr(imy)) / dv vz = (wl(imz) - wcr(imz)) / dv dv = sl - smr by = (wl(iby) - wcr(iby)) / dv bz = (wl(ibz) - wcr(ibz)) / dv vb = sm * bx + vy * by + vz * bz ! conservative variables for the inner left intermediate state ! ui(idn) = dnr ui(imx) = dnr * sm ui(imy) = dnr * vy ui(imz) = dnr * vz ui(ibx) = bx ui(iby) = by ui(ibz) = bz ui(ibp) = ur(ibp,i) ui(ien) = (wcr(ien) + sm * pt - bx * vb) / car ! choose the correct state depending on the sign of contact discontinuity ! advection speed ! if (sm <= 0.0d+00) then ! sm <= 0 ! the inner right intermediate flux ! fi(:,i) = smr * ui(:) - wcr(:) else ! sm > 0 ! vector of the right-going Alfvén wave ! wcl(:) = (sm - smr) * ui(:) + wcr(:) ! calculate the average flux over the left inner intermediate state ! fi(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm end if ! sm > 0 else ! bx = 0 ! no Alfvén wave, so revert to the HLL flux ! fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if end if ! sr* > 0 else ! sl* < sl & sr* > sr ! so far we revert to HLL flux in the case of degeneracies ! fi(:,i) = (sl * wr(:) - sr * wl(:)) / srml end if ! sl* < sl & sr* > sr end if ! deneneracies end if ! sl < 0 < sr end do !------------------------------------------------------------------------------- ! end subroutine riemann_mhd_adi_hlld ! !=============================================================================== ! ! subroutine RIEMANN_MHD_ADI_ROE: ! ------------------------------ ! ! Subroutine solves one dimensional Riemann problem using ! the Roe's method. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Stone, J. M. & Gardiner, T. A., ! "ATHENA: A New Code for Astrophysical MHD", ! The Astrophysical Journal Suplement Series, 2008, 178, pp. 137-177 ! [2] Toro, E. F., ! "Riemann Solvers and Numerical Methods for Fluid dynamics", ! Springer-Verlag Berlin Heidelberg, 2009 ! !=============================================================================== ! subroutine riemann_mhd_adi_roe(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, adiabatic_index use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien, & ibx, iby, ibz, ibp use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: sdl, sdr, sds real(kind=8) :: xx, yy, pml, pmr real(kind=8) :: cc2, ca2, cf2, cs2, cc, ca, cf, cs real(kind=8) :: v2, v2h, br2, br, hp, ayy, sqty real(kind=8) :: bty, btz, qf, qs, sqrtd, vqstr, vbet, norm real(kind=8) :: alf, als, af_prime, as_prime, afpbb, aspbb, afpb, aspb real(kind=8) :: f1, f2, f3, f4, f5 real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension(nf ) :: qi real(kind=8), dimension(7 ) :: lm, al, df logical , save :: first = .true. integer , dimension(7) , save :: ivs real(kind=8) , save :: adi_m1, adi_m2, adi_m2d1 real(kind=8), dimension(7,7), save :: rvec, lvec !$omp threadprivate(first, ivs, adi_m1, adi_m2, adi_m2d1, rvec, lvec) !------------------------------------------------------------------------------- ! if (first) then adi_m1 = adiabatic_index - 1.0d+00 adi_m2 = adiabatic_index - 2.0d+00 adi_m2d1 = adi_m2 / adi_m1 rvec(:,:) = 0.0d+00 lvec(:,:) = 0.0d+00 rvec(4,1) = 1.0d+00 ivs = [ idn, imx, imy, imz, ipr, iby, ibz ] first = .false. end if call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl) call fluxspeed(qr, ur, fr) do i = 1, nn pml = 5.0d-01 * sum(ql(ibx:ibz,i) * ql(ibx:ibz,i)) pmr = 5.0d-01 * sum(qr(ibx:ibz,i) * qr(ibx:ibz,i)) sdl = sqrt(ql(idn,i)) sdr = sqrt(qr(idn,i)) sds = sdl + sdr qi(idn) = sdl * sdr qi(ivx:ivz) = (sdl * ql(ivx:ivz,i) + sdr * qr(ivx:ivz,i)) / sds qi(ipr) = ((ul(ien,i) + ql(ipr,i) + pml) / sdl & + (ur(ien,i) + qr(ipr,i) + pmr) / sdr) / sds qi(iby:ibz) = (sdr * ql(iby:ibz,i) + sdl * qr(iby:ibz,i)) / sds qi(ibx) = ql(ibx,i) qi(ibp) = ql(ibp,i) f1 = qr(iby,i) - ql(iby,i) f2 = qr(ibz,i) - ql(ibz,i) xx = 5.0d-01 * (f1 * f1 + f2 * f2) / (sds * sds) yy = 5.0d-01 * (ql(idn,i) + qr(idn,i)) / qi(idn) br2 = qi(iby) * qi(iby) + qi(ibz) * qi(ibz) if (br2 > 0.0d+00) then br = sqrt(br2) bty = qi(iby) / br btz = qi(ibz) / br else br = 0.0d+00 bty = 1.0d+00 btz = 0.0d+00 end if v2 = sum(qi(ivx:ivz) * qi(ivx:ivz)) v2h = 5.0d-01 * v2 ayy = adi_m1 - adi_m2 * yy sqty = sqrt(ayy) ca2 = qi(ibx) * qi(ibx) hp = qi(ien) - (ca2 + br2) / qi(idn) ca2 = ca2 / qi(idn) f1 = adi_m1 * (hp - v2h) f2 = adi_m2 * xx if (f1 > f2) then cc2 = f1 - f2 cc = sqrt(cc2) f1 = ayy * br2 / qi(idn) f2 = ca2 + f1 f3 = f2 + cc2 f4 = f2 - cc2 f5 = sqrt(f4 * f4 + 4.0d+00 * cc2 * f1) cf2 = 5.0d-01 * (f3 + f5) cs2 = cc2 * ca2 / cf2 else cf2 = ca2 + ayy * br2 / qi(idn) cc2 = 0.0d+00 cc = 0.0d+00 cs2 = 0.0d+00 end if if (cs2 < cc2 .and. cc2 < cf2) then f1 = (cc2 - cs2) / (cf2 - cs2) alf = sqrt(f1) als = sqrt(1.0d+00 - f1) else if (cc2 >= cf2) then alf = 1.0d+00 als = 0.0d+00 else alf = 0.0d+00 als = 1.0d+00 end if cf = sign(sqrt(cf2), qi(ivx)) ca = sign(sqrt(ca2), qi(ivx)) cs = sign(sqrt(cs2), qi(ivx)) lm(1) = qi(ivx) + cf lm(2) = qi(ivx) + ca lm(3) = qi(ivx) + cs lm(4) = qi(ivx) lm(5) = qi(ivx) - cs lm(6) = qi(ivx) - ca lm(7) = qi(ivx) - cf vbet = (qi(ivy) * bty + qi(ivz) * btz) / sqty sqrtd = sqrt(qi(idn)) qf = cf * sign(alf, qi(ibx)) qs = cs * sign(als, qi(ibx)) f1 = cc / sqrtd af_prime = f1 * alf as_prime = f1 * als f1 = br / sqty afpbb = af_prime * f1 aspbb = as_prime * f1 f1 = qs / sqty f2 = qs * vbet f3 = as_prime / sqty rvec(1,1) = alf rvec(1,2) = alf * lm(1) rvec(1,3) = alf * qi(ivy) - f1 * bty rvec(1,4) = alf * qi(ivz) - f1 * btz rvec(1,5) = alf * (hp + qi(ivx) * cf) - f2 + aspbb rvec(1,6) = f3 * bty rvec(1,7) = f3 * btz rvec(7,1) = alf rvec(7,2) = alf * lm(7) rvec(7,3) = alf * qi(ivy) + f1 * bty rvec(7,4) = alf * qi(ivz) + f1 * btz rvec(7,5) = alf * (hp - qi(ivx) * cf) + f2 + aspbb rvec(7,6) = rvec(1,6) rvec(7,7) = rvec(1,7) f1 = qf / sqty f2 = qf * vbet f3 = - af_prime / sqty rvec(3,1) = als rvec(3,2) = als * lm(3) rvec(3,3) = als * qi(ivy) + f1 * bty rvec(3,4) = als * qi(ivz) + f1 * btz rvec(3,5) = als * (hp + qi(ivx) * cs) + f2 - afpbb rvec(3,6) = f3 * bty rvec(3,7) = f3 * btz rvec(5,1) = als rvec(5,2) = als * lm(5) rvec(5,3) = als * qi(ivy) - f1 * bty rvec(5,4) = als * qi(ivz) - f1 * btz rvec(5,5) = als * (hp - qi(ivx) * cs) - f2 - afpbb rvec(5,6) = rvec(3,6) rvec(5,7) = rvec(3,7) f1 = sign(1.0d+00, qi(ivx)) f2 = sign(1.0d+00 / sqrtd, qi(ibx)) rvec(2,3) = f1 * btz rvec(2,4) = - f1 * bty rvec(2,5) = f1 * (qi(ivy) * btz - qi(ivz) * bty) rvec(2,6) = - f2 * btz rvec(2,7) = f2 * bty rvec(6,3) = - rvec(2,3) rvec(6,4) = - rvec(2,4) rvec(6,5) = - rvec(2,5) rvec(6,6) = rvec(2,6) rvec(6,7) = rvec(2,7) rvec(4,2) = qi(ivx) rvec(4,3) = qi(ivy) rvec(4,4) = qi(ivz) rvec(4,5) = v2h + adi_m2d1 * xx norm = 2.0d+00 * cc2 f1 = sqty * br / norm afpb = af_prime * f1 aspb = as_prime * f1 sqty = sqty / norm vqstr = (qi(ivy) * bty + qi(ivz) * btz) * sqty f1 = adi_m1 * alf / norm f2 = qs * sqty f3 = as_prime * qi(idn) * sqty f4 = alf * cf / norm f5 = qs * vqstr lvec(1,1) = f1 * (v2 - hp) - f4 * lm(7) + f5 - aspb lvec(1,2) = - f1 * qi(ivx) + f4 lvec(1,3) = - f1 * qi(ivy) - f2 * bty lvec(1,4) = - f1 * qi(ivz) - f2 * btz lvec(1,5) = f1 lvec(1,6) = f3 * bty - f1 * qi(iby) lvec(1,7) = f3 * btz - f1 * qi(ibz) lvec(7,1) = f1 * (v2 - hp) + f4 * lm(1) - f5 - aspb lvec(7,2) = - f1 * qi(ivx) - f4 lvec(7,3) = - f1 * qi(ivy) + f2 * bty lvec(7,4) = - f1 * qi(ivz) + f2 * btz lvec(7,5) = f1 lvec(7,6) = lvec(1,6) lvec(7,7) = lvec(1,7) f1 = adi_m1 * als / norm f2 = qf * sqty f3 = af_prime * qi(idn) * sqty f4 = als * cs / norm f5 = qf * vqstr lvec(3,1) = f1 * (v2 - hp) - f4 * lm(5) - f5 + afpb lvec(3,2) = - f1 * qi(ivx) + f4 lvec(3,3) = - f1 * qi(ivy) + f2 * bty lvec(3,4) = - f1 * qi(ivz) + f2 * btz lvec(3,5) = f1 lvec(3,6) = - f3 * bty - f1 * qi(iby) lvec(3,7) = - f3 * btz - f1 * qi(ibz) lvec(5,1) = f1 * (v2 - hp) + f4 * lm(3) + f5 + afpb lvec(5,2) = - f1 * qi(ivx) - f4 lvec(5,3) = - f1 * qi(ivy) - f2 * bty lvec(5,4) = - f1 * qi(ivz) - f2 * btz lvec(5,5) = f1 lvec(5,6) = lvec(3,6) lvec(5,7) = lvec(3,7) f1 = sign(5.0d-01, qi(ivx)) * bty f2 = sign(5.0d-01, qi(ivx)) * btz f3 = sign(1.0d+00, qi(ivx)) * sign(sqrtd, qi(ibx)) lvec(2,1) = qi(ivz) * f1 - qi(ivy) * f2 lvec(2,3) = f2 lvec(2,4) = - f1 lvec(2,6) = - f2 * f3 lvec(2,7) = f1 * f3 lvec(6,1) = - lvec(2,1) lvec(6,3) = - lvec(2,3) lvec(6,4) = - lvec(2,4) lvec(6,6) = lvec(2,6) lvec(6,7) = lvec(2,7) f1 = 1.0d+00 / cc2 lvec(4,1) = 1.0d+00 - f1 * (v2h - adi_m2d1 * xx) lvec(4,2) = f1 * qi(ivx) lvec(4,3) = f1 * qi(ivy) lvec(4,4) = f1 * qi(ivz) lvec(4,5) = - f1 lvec(4,6) = f1 * qi(iby) lvec(4,7) = f1 * qi(ibz) al = abs(lm) * matmul(lvec, ur(ivs,i) - ul(ivs,i)) df = matmul(al, rvec) fi(ivs,i) = 5.0d-01 * ((fl(ivs,i) + fr(ivs,i)) - df(:)) fi(ibx,i) = fl(ibx,i) fi(ibp,i) = fl(ibp,i) end do !------------------------------------------------------------------------------- ! end subroutine riemann_mhd_adi_roe ! !=============================================================================== ! ! subroutine RIEMANN_MHD_ADI_KEPES: ! -------------------------------- ! ! Subroutine solves one dimensional adiabatic magnetohydrodynamic ! Riemann problem using the entropy stable KEPES method. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., Bohm, M., ! "Ideal GLM-MHD: About the entropy consistent nine-wave magnetic ! field divergence diminishing ideal magnetohydrodynamics equations", ! Journal of Computational Physics, 2018, 364, pp. 420-467, ! https://doi.org/10.1016/j.jcp.2018.03.002 ! !=============================================================================== ! subroutine riemann_mhd_adi_kepes(ql, qr, fi) use coordinates, only : nn => bcells use equations , only : nf, adiabatic_index, cglm use equations , only : idn, imx, imy, imz, ivx, ivy, ivz, & ibx, iby, ibz, ibp, ipr, ien implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i real(kind=8) :: dna, pra, vxa, vya, vza, bxa, bya, bza, bpa, dnl, prl, pta real(kind=8) :: v2l, v2r, b2l, b2r, bl, br, bta, btl, bp2, eka, ema, vva real(kind=8) :: ch, ca, cs, cf, ca2, cs2, cf2, aa2, x2, x3, ub2, uba real(kind=8) :: alf, als, das, daf, csq, wps, wms, wpf, wmf real(kind=8) :: f1, f2, f3, f4, f5 real(kind=8), dimension(nf) :: v, lm, tm logical , save :: first = .true. real(kind=8) , save :: adi_m1, adi_m1x, adi_m1xi real(kind=8), dimension(9,9), save :: rm !$omp threadprivate(first, adi_m1, adi_m1x, adi_m1xi, rm) !------------------------------------------------------------------------------- ! if (first) then adi_m1 = adiabatic_index - 1.0d+00 adi_m1x = adi_m1 / adiabatic_index adi_m1xi = adiabatic_index / adi_m1 rm(:,:) = 0.0d+00 rm(1,5) = 1.0d+00 rm(6,4) = 1.0d+00 rm(6,6) = 1.0d+00 first = .false. end if do i = 1, nn v2l = sum(ql(ivx:ivz,i) * ql(ivx:ivz,i)) v2r = sum(qr(ivx:ivz,i) * qr(ivx:ivz,i)) b2l = sum(ql(ibx:ibz,i) * ql(ibx:ibz,i)) b2r = sum(qr(ibx:ibz,i) * qr(ibx:ibz,i)) bl = ql(idn,i) / ql(ipr,i) br = qr(idn,i) / qr(ipr,i) btl = lmean(bl, br) bta = amean(bl, br) dnl = lmean(ql(idn,i), qr(idn,i)) dna = amean(ql(idn,i), qr(idn,i)) vxa = amean(ql(ivx,i), qr(ivx,i)) vya = amean(ql(ivy,i), qr(ivy,i)) vza = amean(ql(ivz,i), qr(ivz,i)) bxa = amean(ql(ibx,i), qr(ibx,i)) bya = amean(ql(iby,i), qr(iby,i)) bza = amean(ql(ibz,i), qr(ibz,i)) bpa = amean(ql(ibp,i), qr(ibp,i)) prl = lmean(ql(ipr,i), qr(ipr,i)) pra = dna / bta eka = 5.0d-01 * amean(v2l, v2r) ema = 5.0d-01 * amean(b2l, b2r) pta = pra + ema vva = 5.0d-01 * sum(ql(ivx:ivz,i) * qr(ivx:ivz,i)) ub2 = 5.0d-01 * amean(ql(ivx,i) * b2l, qr(ivx,i) * b2r) uba = amean(sum(ql(ivx:ivz,i) * ql(ibx:ibz,i)), & sum(qr(ivx:ivz,i) * qr(ibx:ibz,i))) v(idn) = 5.0d-01 * (bl * v2l - br * v2r) if (qr(idn,i) > ql(idn,i)) then v(idn) = v(idn) + log(qr(idn,i) / ql(idn,i)) else if (ql(idn,i) > qr(idn,i)) then v(idn) = v(idn) - log(ql(idn,i) / qr(idn,i)) end if if (br > bl) then v(idn) = v(idn) + log(br / bl) / adi_m1 else if (bl > br) then v(idn) = v(idn) - log(bl / br) / adi_m1 end if v(ivx) = br * qr(ivx,i) - bl * ql(ivx,i) v(ivy) = br * qr(ivy,i) - bl * ql(ivy,i) v(ivz) = br * qr(ivz,i) - bl * ql(ivz,i) v(ipr) = bl - br v(ibx) = br * qr(ibx,i) - bl * ql(ibx,i) v(iby) = br * qr(iby,i) - bl * ql(iby,i) v(ibz) = br * qr(ibz,i) - bl * ql(ibz,i) v(ibp) = br * qr(ibp,i) - bl * ql(ibp,i) bp2 = bya * bya + bza * bza if (bp2 > 0.0d+00) then f1 = sqrt(bp2) x2 = bya / f1 x3 = bza / f1 else x2 = 0.0d+00 x3 = 0.0d+00 end if aa2 = adiabatic_index * pra / dnl ca2 = bxa * bxa / dnl f1 = bp2 / dnl f2 = ca2 + f1 f3 = f2 - aa2 f4 = sqrt(f3 * f3 + 4.0d+00 * aa2 * f1) cf2 = 5.0d-01 * (f2 + aa2 + f4) cs2 = aa2 * ca2 / cf2 if (cs2 < aa2 .and. aa2 < cf2) then f1 = (aa2 - cs2) / (cf2 - cs2) alf = sqrt(f1) als = sqrt(1.0d+00 - f1) else if (aa2 >= cf2) then alf = 1.0d+00 als = 0.0d+00 else alf = 0.0d+00 als = 1.0d+00 end if cf = sign(sqrt(cf2), vxa) cs = sign(sqrt(cs2), vxa) ca = sign(sqrt(ca2), vxa) ch = sign( cglm, vxa) lm(1) = vxa + cf lm(2) = vxa + ca lm(3) = vxa + cs lm(4) = vxa + ch lm(5) = vxa lm(6) = vxa - ch lm(7) = vxa - cs lm(8) = vxa - ca lm(9) = vxa - cf f1 = 5.0d-01 / bta f2 = adiabatic_index * dnl tm(1) = 5.0d-01 / f2 tm(2) = f1 / dnl**2 tm(3) = tm(1) tm(4) = f1 tm(5) = dnl * adi_m1x tm(6) = tm(4) tm(7) = tm(1) tm(8) = tm(2) tm(9) = tm(1) das = dnl * als daf = dnl * alf csq = sqrt(f2 / bta) f3 = vva + adi_m1xi * prl / dnl f4 = vya * x2 + vza * x3 f5 = sqrt(adiabatic_index * bp2 / (dnl * bta)) f1 = dnl * (als * f3 - alf * f5) f2 = dnl * (als * cs * vxa + sign(alf, bxa) * cf * f4) wps = f1 + f2 wms = f1 - f2 f1 = dnl * (alf * f3 + als * f5) f2 = dnl * (alf * cf * vxa - sign(als, bxa) * cs * f4) wpf = f1 + f2 wmf = f1 - f2 f1 = daf f2 = sign(das, bxa) * cs f3 = als * csq rm(1,1) = f1 rm(2,1) = f1 * lm(1) rm(3,1) = f1 * vya - f2 * x2 rm(4,1) = f1 * vza - f2 * x3 rm(5,1) = wpf rm(7,1) = f3 * x2 rm(8,1) = f3 * x3 rm(1,9) = f1 rm(2,9) = f1 * lm(9) rm(3,9) = f1 * vya + f2 * x2 rm(4,9) = f1 * vza + f2 * x3 rm(5,9) = wmf rm(7,9) = rm(7,1) rm(8,9) = rm(8,1) f1 = sign(dnl * sqrt(dna), vxa) rm(3,2) = f1 * x3 rm(4,2) = - f1 * x2 rm(5,2) = - f1 * (x2 * vza - x3 * vya) rm(7,2) = - dnl * x3 rm(8,2) = dnl * x2 rm(3,8) = - rm(3,2) rm(4,8) = - rm(4,2) rm(5,8) = - rm(5,2) rm(7,8) = rm(7,2) rm(8,8) = rm(8,2) f1 = das f2 = sign(daf, bxa) * cf f3 = - alf * csq rm(1,3) = f1 rm(2,3) = f1 * lm(3) rm(3,3) = f1 * vya + f2 * x2 rm(4,3) = f1 * vza + f2 * x3 rm(5,3) = wps rm(7,3) = f3 * x2 rm(8,3) = f3 * x3 rm(1,7) = f1 rm(2,7) = f1 * lm(7) rm(3,7) = f1 * vya - f2 * x2 rm(4,7) = f1 * vza - f2 * x3 rm(5,7) = wms rm(7,7) = rm(7,3) rm(8,7) = rm(8,3) f1 = sign(1.0d+00, vxa) rm(5,4) = bxa + bpa * f1 rm(5,6) = bxa - bpa * f1 rm(9,4) = f1 rm(9,6) = - f1 rm(2,5) = vxa rm(3,5) = vya rm(4,5) = vza rm(5,5) = vva fi(idn,i) = dnl * vxa fi(imx,i) = fi(idn,i) * vxa - bxa * bxa + pta fi(imy,i) = fi(idn,i) * vya - bxa * bya fi(imz,i) = fi(idn,i) * vza - bxa * bza fi(ibx,i) = cglm * bpa fi(iby,i) = vxa * bya - bxa * vya fi(ibz,i) = vxa * bza - bxa * vza fi(ibp,i) = cglm * bxa fi(ien,i) = fi(idn,i) * (1.0d+00 / (adi_m1 * btl) - eka) & + (fi(imx,i) * vxa + fi(imy,i) * vya + fi(imz,i) * vza) & + (fi(ibx,i) * bxa + fi(iby,i) * bya + fi(ibz,i) * bza) & + fi(ibp,i) * bpa - ub2 + bxa * uba & - cglm * amean(ql(ibx,i) * ql(ibp,i), qr(ibx,i) * qr(ibp,i)) fi(:,i) = fi(:,i) - 5.0d-01 * matmul(rm, (abs(lm) * tm) * matmul(v, rm)) end do !------------------------------------------------------------------------------- ! end subroutine riemann_mhd_adi_kepes ! !******************************************************************************* ! ! RELATIVISTIC ADIABATIC MAGNETOHYDRODYNAMIC RIEMANN SOLVERS ! !******************************************************************************* ! !=============================================================================== ! ! subroutine RIEMANN_SRMHD_HLLC: ! ----------------------------- ! ! Subroutine solves one dimensional Riemann problem using ! the Harten-Lax-van Leer method with contact discontinuity resolution (HLLC) ! by Mignone & Bodo. ! ! Arguments: ! ! ql, qr - the primitive variables of the left and right Riemann states; ! fi - the numerical flux at the cell interface; ! ! References: ! ! [1] Mignone, A. & Bodo, G. ! "An HLLC Riemann solver for relativistic flows - II. ! Magnetohydrodynamics", ! Monthly Notices of the Royal Astronomical Society, ! 2006, Volume 368, Pages 1040-1054 ! !=============================================================================== ! subroutine riemann_srmhd_hllc(ql, qr, fi) use algebra , only : quadratic use coordinates, only : nn => bcells use equations , only : nf, idn, ivx, imx, imy, imz, ien, ibx, iby, ibz, ibp use equations , only : prim2cons, fluxspeed implicit none real(kind=8), dimension(:,:), intent(in) :: ql, qr real(kind=8), dimension(:,:), intent(out) :: fi integer :: i, nr real(kind=8) :: sl, sr, srml, sm real(kind=8) :: bx2 real(kind=8) :: pt, dv, fc, uu, ff, uf real(kind=8) :: vv, vb, gi real(kind=8), dimension(nf,nn) :: ul, ur, fl, fr real(kind=8), dimension( 2,nn) :: cl, cr real(kind=8), dimension(nf ) :: wl, wr, uh, us, fh real(kind=8), dimension(3) :: a, vs real(kind=8), dimension(2) :: x !------------------------------------------------------------------------------- ! call prim2cons(ql, ul) call prim2cons(qr, ur) call fluxspeed(ql, ul, fl, cl) call fluxspeed(qr, ur, fr, cr) do i = 1, nn sl = min(cl(1,i), cr(1,i)) sr = max(cl(2,i), cr(2,i)) if (sl >= 0.0d+00) then fi(:,i) = fl(:,i) else if (sr <= 0.0d+00) then fi(:,i) = fr(:,i) else ! sl < 0 < sr srml = sr - sl wl(:) = sl * ul(:,i) - fl(:,i) wr(:) = sr * ur(:,i) - fr(:,i) uh(:) = ( wr(:) - wl(:)) / srml fh(:) = (sl * wr(:) - sr * wl(:)) / srml wl(ien) = wl(ien) + wl(idn) wr(ien) = wr(ien) + wr(idn) ! calculate Bₓ² ! bx2 = ql(ibx,i) * qr(ibx,i) ! calculate the contact dicontinuity speed solving eq. (41) ! if (bx2 >= 1.0d-16) then ! prepare the quadratic coefficients, (eq. 42 in [1]) ! uu = sum(uh(iby:ibz) * uh(iby:ibz)) ff = sum(fh(iby:ibz) * fh(iby:ibz)) uf = sum(uh(iby:ibz) * fh(iby:ibz)) a(1) = uh(imx) - uf a(2) = uu + ff - (fh(imx) + uh(ien) + uh(idn)) a(3) = fh(ien) + fh(idn) - uf ! solve the quadratic equation ! nr = quadratic(a(1:3), x(1:2)) ! if Δ < 0, just use the HLL flux ! if (nr < 1) then fi(:,i) = fh(:) else ! get the contact dicontinuity speed ! sm = x(1) ! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux ! if ((sm <= sl) .or. (sm >= sr)) then fi(:,i) = fh(:) else ! substitute magnetic field components from the HLL state (eq. 37 in [1]) ! us(ibx) = ql(ibx,i) us(iby) = uh(iby) us(ibz) = uh(ibz) us(ibp) = ql(ibp,i) ! calculate velocity components without Bₓ normalization (eq. 38 in [1]) ! vs(1) = sm vs(2) = (us(iby) * sm - fh(iby)) / us(ibx) vs(3) = (us(ibz) * sm - fh(ibz)) / us(ibx) ! calculate v² and v.B ! vv = sum(vs(1:3) * vs( 1: 3)) vb = sum(vs(1:3) * us(ibx:ibz)) ! calculate inverse of Lorentz factor ! gi = 1.0d+00 - vv ! calculate total pressure (eq. 40 in [1]) ! pt = fh(imx) + gi * bx2 - (fh(ien) + fh(idn) - us(ibx) * vb) * sm ! if the pressure is negative, use the HLL flux ! if (pt <= 0.0d+00) then fi(:,i) = fh(:) else ! depending in the sign of the contact dicontinuity speed, calculate the proper ! state and corresponding flux ! if (sm > 0.0d+00) then ! calculate the conserved variable vector (eqs. 43-46 in [1]) ! dv = sl - sm fc = (sl - ql(ivx,i)) / dv us(idn) = fc * ul(idn,i) us(imy) = (sl * ul(imy,i) - fl(imy,i) & - us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv us(imz) = (sl * ul(imz,i) - fl(imz,i) & - us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv us(ien) = (sl * (ul(ien,i) + ul(idn,i)) & - (fl(ien,i) + fl(idn,i)) & + pt * sm - us(ibx) * vb) / dv - us(idn) ! calculate normal component of momentum ! us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb ! calculate the flux (eq. 14 in [1]) ! fi(:,i) = fl(:,i) + sl * (us(:) - ul(:,i)) else if (sm < 0.0d+00) then ! calculate the conserved variable vector (eqs. 43-46 in [1]) ! dv = sr - sm fc = (sr - qr(ivx,i)) / dv us(idn) = fc * ur(idn,i) us(imy) = (sr * ur(imy,i) - fr(imy,i) & - us(ibx) * (gi * us(iby) + vs(2) * vb)) / dv us(imz) = (sr * ur(imz,i) - fr(imz,i) & - us(ibx) * (gi * us(ibz) + vs(3) * vb)) / dv us(ien) = (sr * (ur(ien,i) + ur(idn,i)) & - (fr(ien,i) + fr(idn,i)) & + pt * sm - us(ibx) * vb) / dv - us(idn) ! calculate normal component of momentum ! us(imx) = (us(ien) + us(idn) + pt) * sm - us(ibx) * vb ! calculate the flux (eq. 14 in [1]) ! fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i)) else ! intermediate flux is constant across the contact discontinuity so use HLL flux ! fi(:,i) = fh(:) end if ! sm == 0 end if ! p* < 0 end if ! sl < sm < sr end if ! nr < 1 else ! Bx² > 0 ! prepare the quadratic coefficients (eq. 47 in [1]) ! a(1) = uh(imx) a(2) = - (fh(imx) + uh(ien) + uh(idn)) a(3) = fh(ien) + fh(idn) ! solve the quadratic equation ! nr = quadratic(a(1:3), x(1:2)) ! if Δ < 0, just use the HLL flux ! if (nr < 1) then fi(:,i) = fh(:) else ! get the contact dicontinuity speed ! sm = x(1) ! if the contact discontinuity speed exceeds the sonic speeds, use the HLL flux ! if ((sm <= sl) .or. (sm >= sr)) then fi(:,i) = fh(:) else ! calculate total pressure (eq. 48 in [1]) ! pt = fh(imx) - (fh(ien) + fh(idn)) * sm ! if the pressure is negative, use the HLL flux ! if (pt <= 0.0d+00) then fi(:,i) = fh(:) else ! depending in the sign of the contact dicontinuity speed, calculate the proper ! state and corresponding flux ! if (sm > 0.0d+00) then ! calculate the conserved variable vector (eqs. 43-46 in [1]) ! dv = sl - sm us(idn) = wl(idn) / dv us(imy) = wl(imy) / dv us(imz) = wl(imz) / dv us(ien) = (wl(ien) + pt * sm) / dv us(imx) = (us(ien) + pt) * sm us(ien) = us(ien) - us(idn) us(ibx) = 0.0d+00 us(iby) = wl(iby) / dv us(ibz) = wl(ibz) / dv us(ibp) = ql(ibp,i) ! calculate the flux (eq. 27 in [1]) ! fi(:,i) = fl(:,i) + sl * (us(:) - ul(:,i)) else if (sm < 0.0d+00) then ! calculate the conserved variable vector (eqs. 43-46 in [1]) ! dv = sr - sm us(idn) = wr(idn) / dv us(imy) = wr(imy) / dv us(imz) = wr(imz) / dv us(ien) = (wr(ien) + pt * sm) / dv us(imx) = (us(ien) + pt) * sm us(ien) = us(ien) - us(idn) us(ibx) = 0.0d+00 us(iby) = wr(iby) / dv us(ibz) = wr(ibz) / dv us(ibp) = qr(ibp,i) ! calculate the flux (eq. 27 in [1]) ! fi(:,i) = fr(:,i) + sr * (us(:) - ur(:,i)) else ! intermediate flux is constant across the contact discontinuity so use HLL flux ! fi(:,i) = fh(:) end if ! sm == 0 end if ! p* < 0 end if ! sl < sm < sr end if ! nr < 1 end if ! Bx² == 0 or > 0 end if ! sl < 0 < sr end do !------------------------------------------------------------------------------- ! end subroutine riemann_srmhd_hllc ! !******************************************************************************* ! ! SUPPORTING SUBROUTINES/FUNCTIONS ! !******************************************************************************* ! !=============================================================================== ! ! subroutine HIGHER_ORDER_FLUX_CORRECTION: ! --------------------------------------- ! ! Subroutine adds higher order corrections to the numerical fluxes. ! ! Arguments: ! ! f - the vector of numerical fluxes; ! !=============================================================================== ! subroutine higher_order_flux_correction(f) ! include external procedures ! use interpolations, only : order ! local variables are not implicit by default ! implicit none ! subroutine arguments ! real(kind=8), dimension(:,:), intent(inout) :: f ! local variables ! integer :: n ! local arrays to store the states ! real(kind=8), dimension(size(f,1),size(f,2)) :: f2, f4 ! !------------------------------------------------------------------------------- ! ! depending on the scheme order calculate and add higher order corrections, ! if required ! if (.not. high_order_flux .or. order <= 2) return ! get the length of vector ! n = size(f, 2) ! 2nd order correction ! f2(:,2:n-1) = (2.0d+00 * f(:,2:n-1) - (f(:,3:n) + f(:,1:n-2))) / 2.4d+01 f2(:,1 ) = 0.0d+00 f2(:, n ) = 0.0d+00 ! 4th order correction ! if (order > 4) then f4(:,3:n-2) = 3.0d+00 * (f(:,1:n-4) + f(:,5:n ) & - 4.0d+00 * (f(:,2:n-3) + f(:,4:n-1)) & + 6.0d+00 * f(:,3:n-2)) / 6.4d+02 f4(:,1 ) = 0.0d+00 f4(:,2 ) = 0.0d+00 f4(:, n-1) = 0.0d+00 f4(:, n ) = 0.0d+00 ! correct the flux for 4th order ! f(:,:) = f(:,:) + f4(:,:) end if ! correct the flux for 2nd order ! f(:,:) = f(:,:) + f2(:,:) !------------------------------------------------------------------------------- ! end subroutine higher_order_flux_correction ! !=============================================================================== ! ! function LMEAN: ! -------------- ! ! Function calculates the logarithmic mean using the optimized algorithm ! by Ranocha et al. ! ! Arguments: ! ! l, r - the values to calculate the mean of; ! ! References: ! ! [1] Ranocha et al., ! "Efficient implementation of modern entropy stable and kinetic energy ! preserving discontinuous Galerkin methods for conservation laws", ! http://arxiv.org/abs/2112.10517v1 ! !=============================================================================== ! real(kind=8) function lmean(l, r) implicit none real(kind=8), intent(in) :: l, r real(kind=8) :: u, d, s real(kind=8), parameter :: c1 = 2.0d+00, c2 = 2.0d+00 / 3.0d+00, & c3 = 4.0d-01, c4 = 2.0d+00 / 7.0d+00 !------------------------------------------------------------------------------- ! d = r - l s = r + l u = (d / s)**2 if (u < 1.0d-04) then lmean = s / (c1 + u * (c2 + u * (c3 + u * c4))) else if (d >= 0.0d+00) then s = log(r / l) else s = log(l / r) end if lmean = abs(d) / s end if !------------------------------------------------------------------------------- ! end function lmean ! !=============================================================================== ! ! function AMEAN: ! -------------- ! ! Function calculate the arithmetic mean. ! ! Arguments: ! ! l, r - the values to calculate the mean of; ! !=============================================================================== ! real(kind=8) function amean(l, r) implicit none real(kind=8), intent(in) :: l, r !------------------------------------------------------------------------------- ! amean = 5.0d-01 * (l + r) !------------------------------------------------------------------------------- ! end function amean !=============================================================================== ! end module schemes