!!****************************************************************************** !! !! module: scheme - handling the actual solver of the set of equations !! !! Copyright (C) 2008-2010 Grzegorz Kowal !! !!****************************************************************************** !! !! This file is part of Godunov-AMR. !! !! Godunov-AMR 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. !! !! Godunov-AMR 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 scheme implicit none ! the maximal speed in the system ! real, save :: cmax contains ! !=============================================================================== ! ! update: subroutine sweeps over all directions and integrates the directional ! derivatives of the flux in order to get the increment of solution ! !=============================================================================== ! subroutine update(u, du, dxi, dyi, dzi) use config , only : im, jm, km use variables , only : nvr, nqt, nfl use variables , only : idn, imx, imy, imz #ifdef ADI use variables , only : ien #endif /* ADI */ #ifdef MHD use variables , only : ibx, iby, ibz #ifdef GLM use variables , only : iph #endif /* GLM */ #endif /* MHD */ implicit none ! input arguments ! real, dimension(nqt,im,jm,km) , intent(in) :: u real, dimension(nqt,im,jm,km) , intent(out) :: du real , intent(in) :: dxi, dyi, dzi ! local variables ! integer :: i, j, k, im1, jm1, km1, ip1, jp1, kp1 real :: dx, dy, dz ! local temporary arrays ! real, dimension(nvr,im) :: ux real, dimension(nqt,im) :: fx real, dimension(nvr,jm) :: uy real, dimension(nqt,jm) :: fy #if NDIMS == 3 real, dimension(nvr,km) :: uz real, dimension(nqt,km) :: fz #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! ! reset the increment array ! du(:,:,:,:) = 0.0 ! prepare the spacial increment ! dx = 1.0d0 / dxi dy = 1.0d0 / dyi #if NDIMS == 3 dz = 1.0d0 / dzi #endif /* NDIMS == 3 */ ! update along X-direction ! do k = 1, km #if NDIMS == 3 #ifdef MHD km1 = max(k - 1, 1) kp1 = min(k + 1,km) #endif /* MHD */ #endif /* NDIMS == 3 */ do j = 1, jm #ifdef MHD jm1 = max(j - 1, 1) jp1 = min(j + 1,jm) #endif /* MHD */ ! copy directional vectors of variables for the one dimensional solver ! do i = 1, im ux(idn,i) = u(idn,i,j,k) ux(imx,i) = u(imx,i,j,k) ux(imy,i) = u(imy,i,j,k) ux(imz,i) = u(imz,i,j,k) #ifdef ADI ux(ien,i) = u(ien,i,j,k) #endif /* ADI */ #ifdef MHD ux(ibx,i) = u(ibx,i,j,k) ux(iby,i) = u(iby,i,j,k) ux(ibz,i) = u(ibz,i,j,k) #ifdef GLM ux(iph,i) = u(iph,i,j,k) #endif /* GLM */ #endif /* MHD */ end do ! execute solver (returns fluxes for the update) ! #ifdef HLL call hll (im, dx, ux(:,:), fx(:,:)) #endif /* HLL */ #ifdef HLLC call hllc(im, dx, ux(:,:), fx(:,:)) #endif /* HLLC */ #ifdef HLLD call hlld(im, dx, ux(:,:), fx(:,:)) #endif /* HLLD */ ! update the arrays of increments ! do i = 1, im du(idn,i,j,k) = du(idn,i,j,k) + dxi * fx(idn,i) du(imx,i,j,k) = du(imx,i,j,k) + dxi * fx(imx,i) du(imy,i,j,k) = du(imy,i,j,k) + dxi * fx(imy,i) du(imz,i,j,k) = du(imz,i,j,k) + dxi * fx(imz,i) #ifdef ADI du(ien,i,j,k) = du(ien,i,j,k) + dxi * fx(ien,i) #endif /* ADI */ #ifdef MHD ! update magnetic variables ! #ifdef GLM du(ibx,i,j,k) = du(ibx,i,j,k) + dxi * fx(ibx,i) du(iby,i,j,k) = du(iby,i,j,k) + dxi * fx(iby,i) du(ibz,i,j,k) = du(ibz,i,j,k) + dxi * fx(ibz,i) du(iph,i,j,k) = du(iph,i,j,k) + dxi * fx(iph,i) #endif /* GLM */ #endif /* MHD */ end do end do end do ! update along Y-direction ! do k = 1, km #if NDIMS == 3 #ifdef MHD km1 = max(k - 1, 1) kp1 = min(k + 1,km) #endif /* MHD */ #endif /* NDIMS == 3 */ do i = 1, im #ifdef MHD im1 = max(i - 1, 1) ip1 = min(i + 1,im) #endif /* MHD */ ! copy directional vectors of variables for the one dimensional solver ! do j = 1, jm uy(idn,j) = u(idn,i,j,k) uy(imx,j) = u(imy,i,j,k) uy(imy,j) = u(imz,i,j,k) uy(imz,j) = u(imx,i,j,k) #ifdef ADI uy(ien,j) = u(ien,i,j,k) #endif /* ADI */ #ifdef MHD uy(ibx,j) = u(iby,i,j,k) uy(iby,j) = u(ibz,i,j,k) uy(ibz,j) = u(ibx,i,j,k) #ifdef GLM uy(iph,j) = u(iph,i,j,k) #endif /* GLM */ #endif /* MHD */ end do ! execute solver (returns fluxes for the update) ! #ifdef HLL call hll (jm, dy, uy(:,:), fy(:,:)) #endif /* HLL */ #ifdef HLLC call hllc(jm, dy, uy(:,:), fy(:,:)) #endif /* HLLC */ #ifdef HLLD call hlld(jm, dy, uy(:,:), fy(:,:)) #endif /* HLLD */ ! update the arrays of increments ! do j = 1, jm du(idn,i,j,k) = du(idn,i,j,k) + dyi * fy(idn,j) du(imx,i,j,k) = du(imx,i,j,k) + dyi * fy(imz,j) du(imy,i,j,k) = du(imy,i,j,k) + dyi * fy(imx,j) du(imz,i,j,k) = du(imz,i,j,k) + dyi * fy(imy,j) #ifdef ADI du(ien,i,j,k) = du(ien,i,j,k) + dyi * fy(ien,j) #endif /* ADI */ #ifdef MHD ! update magnetic variables ! #ifdef GLM du(ibx,i,j,k) = du(ibx,i,j,k) + dyi * fy(ibz,j) du(iby,i,j,k) = du(iby,i,j,k) + dyi * fy(ibx,j) du(ibz,i,j,k) = du(ibz,i,j,k) + dyi * fy(iby,j) du(iph,i,j,k) = du(iph,i,j,k) + dyi * fy(iph,j) #endif /* GLM */ #endif /* MHD */ end do end do end do #if NDIMS == 3 ! update along Z-direction ! do j = 1, jm #ifdef MHD jm1 = max(j - 1, 1) jp1 = min(j + 1,jm) #endif /* MHD */ do i = 1, im #ifdef MHD im1 = max(i - 1, 1) ip1 = min(i + 1,im) #endif /* MHD */ ! copy directional vectors of variables for the one dimensional solver ! do k = 1, km uz(idn,k) = u(idn,i,j,k) uz(imx,k) = u(imz,i,j,k) uz(imy,k) = u(imx,i,j,k) uz(imz,k) = u(imy,i,j,k) #ifdef ADI uz(ien,k) = u(ien,i,j,k) #endif /* ADI */ #ifdef MHD uz(ibx,k) = u(ibz,i,j,k) uz(iby,k) = u(ibx,i,j,k) uz(ibz,k) = u(iby,i,j,k) #ifdef GLM uz(iph,k) = u(iph,i,j,k) #endif /* GLM */ #endif /* MHD */ end do ! execute solver (returns fluxes for the update) ! #ifdef HLL call hll (km, dz, uz(:,:), fz(:,:)) #endif /* HLL */ #ifdef HLLC call hllc(km, dz, uz(:,:), fz(:,:)) #endif /* HLLC */ #ifdef HLLD call hlld(km, dz, uz(:,:), fz(:,:)) #endif /* HLLD */ ! update the arrays of increments ! do k = 1, km du(idn,i,j,k) = du(idn,i,j,k) + dzi * fz(idn,k) du(imx,i,j,k) = du(imx,i,j,k) + dzi * fz(imy,k) du(imy,i,j,k) = du(imy,i,j,k) + dzi * fz(imz,k) du(imz,i,j,k) = du(imz,i,j,k) + dzi * fz(imx,k) #ifdef ADI du(ien,i,j,k) = du(ien,i,j,k) + dzi * fz(ien,k) #endif /* ADI */ #ifdef MHD ! update magnetic variables ! #ifdef GLM du(ibx,i,j,k) = du(ibx,i,j,k) + dzi * fz(iby,k) du(iby,i,j,k) = du(iby,i,j,k) + dzi * fz(ibz,k) du(ibz,i,j,k) = du(ibz,i,j,k) + dzi * fz(ibx,k) du(iph,i,j,k) = du(iph,i,j,k) + dzi * fz(iph,k) #endif /* GLM */ #endif /* MHD */ end do end do end do #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! end subroutine update #ifdef HLL ! !=============================================================================== ! ! hll: subroutine computes the approximated flux using the HLL method ! !=============================================================================== ! subroutine hll(n, h, u, f) use interpolation, only : reconstruct use variables , only : nvr, nfl, nqt use variables , only : ivx, ivz #ifdef MHD use variables , only : ibx, iby, ibz #ifdef GLM use variables , only : iph #endif /* GLM */ #endif /* MHD */ implicit none ! input/output arguments ! integer , intent(in) :: n real , intent(in) :: h real, dimension(nvr,n), intent(in) :: u real, dimension(nqt,n), intent(out) :: f ! local variables ! integer :: p, i real, dimension(nvr,n) :: q, ql, qr, ul, ur real, dimension(nqt,n) :: fl, fr, fn real, dimension(n) :: cl, cr real :: al, ar, ap, div ! !------------------------------------------------------------------------------- ! ! calculate the primitive variables ! call cons2prim(n, u, q) ! reconstruct the left and right states of the primitive variables ! do p = 1, nfl call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) end do #ifdef MHD ! reconstruct the left and right states of the magnetic field components ! do p = ibx, ibz call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) end do #ifdef GLM ! reconstruct the left and right states of the scalar potential ! call reconstruct(n, h, q(iph,:), ql(iph,:), qr(iph,:)) ! obtain the state values for Bx and Psi for the GLM-MHD equations ! cl(:) = 0.5d0 * ((qr(ibx,:) + ql(ibx,:)) - (qr(iph,:) - ql(iph,:)) / cmax) cr(:) = 0.5d0 * ((qr(iph,:) + ql(iph,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) ql(ibx,:) = cl(:) qr(ibx,:) = cl(:) ql(iph,:) = cr(:) qr(iph,:) = cr(:) #endif /* GLM */ #endif /* MHD */ ! calculate conservative variables at states ! call prim2cons(n, ql, ul) call prim2cons(n, qr, ur) ! calculate fluxes and speeds ! call fluxspeed(n, ql, ul, fl, cl) call fluxspeed(n, qr, ur, fr, cr) ! iterate over all points ! do i = 1, n ! calculate min and max and intermediate speeds: eq. (67) ! al = min(ql(ivx,i) - cl(i), qr(ivx,i) - cr(i)) ar = max(ql(ivx,i) + cl(i), qr(ivx,i) + cr(i)) ! calculate HLL flux ! if (al .ge. 0.0) then fn(:,i) = fl(:,i) else if (ar .le. 0.0) then fn(:,i) = fr(:,i) else ap = ar * al div = 1.0 / (ar - al) fn(:,i) = div * (ar * fl(:,i) - al * fr(:,i) + ap * (ur(:,i) - ul(:,i))) end if end do ! calculate numerical flux ! f( 1:nfl,2:n) = - fn( 1:nfl,2:n) + fn( 1:nfl,1:n-1) #ifdef MHD #ifdef GLM f(ibx:ibz,2:n) = - fn(ibx:ibz,2:n) + fn(ibx:ibz,1:n-1) f(iph ,2:n) = - fn(iph ,2:n) + fn(iph ,1:n-1) #endif /* GLM */ #endif /* MHD */ !------------------------------------------------------------------------------- ! end subroutine hll #endif /* HLL */ #ifdef HLLC !=============================================================================== ! ! hllc: subroutine to compute flux approximated by HLLC method (HYDRO only) ! ([1] Batten et al., 1997, JSC, 18, 6, 1553) ! ([2] Miyoshi & Kusano, 2005, JCP, 208, 315) ! !=============================================================================== ! subroutine hllc(n, h, u, f) use interpolation, only : reconstruct use variables , only : nvr, nfl, nqt use variables , only : idn, imx, imy, imz, ien, ivx, ivy, ivz, ipr implicit none ! input/output arguments ! integer , intent(in) :: n real , intent(in) :: h real, dimension(nvr,n), intent(in) :: u real, dimension(nqt,n), intent(out) :: f ! local variables ! integer :: p, i real, dimension(nvr,n) :: ql, qr, q, ul, ur real, dimension(nfl,n) :: fl, fr, fn real, dimension(n) :: cl, cr, cm real :: sl, sr, sm, sml, smr, srmv, slmv, srmm, slmm & , smvl, smvr, div, pt real, dimension(nvr) :: q1l, q1r, u1l, u1r ! !------------------------------------------------------------------------------- ! ! obtain the primitive variables ! call cons2prim(n, u(:,:), q(:,:)) ! reconstruct the left and right states of the primitive variables ! do p = 1, nfl call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) end do ! obtain the conservative variables at both states ! call prim2cons(n, ql(:,:), ul(:,:)) call prim2cons(n, qr(:,:), ur(:,:)) ! calculate the physical fluxes and speeds ! call fluxspeed(n, ql(:,:), ul(:,:), fl(:,:), cl(:)) call fluxspeed(n, qr(:,:), ur(:,:), fr(:,:), cr(:)) ! iterate along the direction ! do i = 1, n ! calculate the minimum and maxximum speeds ! sl = min(ql(ivx,i) - cl(i), qr(ivx,i) - cr(i)) sr = max(ql(ivx,i) + cl(i), qr(ivx,i) + cr(i)) ! all speeds >= 0, left side flux ! if (sl .ge. 0.0) then fn(:,i) = fl(:,i) ! all speeds <= 0, right side flux ! else if (sr .le. 0.0) then fn(:,i) = fr(:,i) ! intermediate states ! else ! sl < 0 & sr > 0 ! useful differences ! slmv = sl - ql(ivx,i) srmv = sr - qr(ivx,i) ! speed of the contact discontinuity (eq. 34 [1], 14 [2]) ! div = srmv * qr(idn,i) - slmv * ql(idn,i) sml = (srmv * ur(imx,i) - slmv * ul(imx,i) & - qr(ipr,i) + ql(ipr,i)) / div div = slmv * ql(idn,i) - srmv * qr(idn,i) smr = (slmv * ul(imx,i) - srmv * ur(imx,i) & - ql(ipr,i) + qr(ipr,i)) / div sm = 0.5d0 * (sml + smr) if (sm .eq. 0.0d0) then ! calculate the left intermediate state ! pt = ql(ipr,i) - ul(imx,i) * slmv u1l(idn) = ql(idn,i) * slmv / sl u1l(imx) = 0.0d0 u1l(imy) = u1l(idn) * ql(ivy,i) u1l(imz) = u1l(idn) * ql(ivz,i) if (sl .eq. 0.0d0) then u1l(ien) = ul(ien,i) else u1l(ien) = (slmv * ul(ien,i) - ql(ipr,i) * ql(ivx,i)) / sl end if ! calculate right intermediate state ! pt = qr(ipr,i) - ur(imx,i) * srmv u1r(idn) = qr(idn,i) * srmv / sr u1r(imx) = 0.0d0 u1r(imy) = u1r(idn) * qr(ivy,i) u1r(imz) = u1r(idn) * qr(ivz,i) if (sr .eq. 0.0d0) then u1r(ien) = ur(ien,i) else u1r(ien) = (srmv * ur(ien,i) - qr(ipr,i) * qr(ivx,i)) / sr endif ! calculate intermediate flux ! fn(:,i) = 0.5d0 * (fl(:,i) + sl * (u1l(:) - ul(:,i)) & + fr(:,i) + sr * (u1r(:) - ur(:,i))) else ! useful differences ! slmm = sl - sm srmm = sr - sm smvl = sm - ql(ivx,i) smvr = sm - qr(ivx,i) ! intermediate discontinuities ! if (sm .gt. 0.0d0) then ! pressure of intermediate states (eq. 36 [1], 16 [2]) ! pt = ql(ipr,i) + ql(idn,i) * slmv * smvl ! calculate the left intermediate state ! u1l(idn) = ql(idn,i) * slmv / slmm u1l(imx) = u1l(idn) * sm u1l(imy) = u1l(idn) * ql(ivy,i) u1l(imz) = u1l(idn) * ql(ivz,i) if (slmm .eq. 0.0d0) then u1l(ien) = ul(ien,i) else u1l(ien) = (slmv * ul(ien,i) - ql(ipr,i) * ql(ivx,i) & + pt * sm) / slmm end if ! calculate the left intermediate flux ! fn(:,i) = fl(:,i) + sl * (u1l(:) - ul(:,i)) else if (sm .lt. 0.0) then ! pressure of intermediate states (eq. 36 [1], 16 [2]) ! pt = qr(ipr,i) + qr(idn,i) * srmv * smvr ! calculate the right intermediate state ! u1r(idn) = qr(idn,i) * srmv / srmm u1r(imx) = u1r(idn) * sm u1r(imy) = u1r(idn) * qr(ivy,i) u1r(imz) = u1r(idn) * qr(ivz,i) if (srmm .eq. 0.0d0) then u1r(ien) = ur(ien,i) else u1r(ien) = (srmv * ur(ien,i) - qr(ipr,i) * qr(ivx,i) & + pt * sm) / srmm end if ! calculate the right intermediate flux ! fn(:,i) = fr(:,i) + sr * (u1r(:) - ur(:,i)) end if end if end if end do ! calculate numerical flux ! f(:,2:n) = - fn(:,2:n) + fn(:,1:n-1) !------------------------------------------------------------------------------- ! end subroutine hllc #endif /* HLLC */ #ifdef MHD #ifdef HLLD #ifdef ISO ! !=============================================================================== ! ! hlld: subroutine computes the approximated flux using the HLLD method ! for the isothermal equation of state ! !=============================================================================== ! subroutine hlld(n, h, u, f) use interpolation, only : reconstruct use variables , only : nvr, nfl, nqt use variables , only : idn, imx, imy, imz, ivx, ivy, ivz use variables , only : ibx, iby, ibz #ifdef GLM use variables , only : iph #endif /* GLM */ implicit none ! input/output arguments ! integer , intent(in) :: n real , intent(in) :: h real, dimension(nvr,n), intent(in) :: u real, dimension(nqt,n), intent(out) :: f ! local variables ! integer :: p, i real, dimension(nvr,n) :: q, ql, qr, ul, ur real, dimension(nqt,n) :: fl, fr, fn real, dimension(n) :: cl, cr real, dimension(nvr) :: u1l, u1r, u2 real :: sl, sr, srl, srml, sm, sml, smr real :: dnm, mxm, sqd, div, fac, bxs ! !------------------------------------------------------------------------------- ! ! calculate the primitive variables ! call cons2prim(n, u, q) ! reconstruct the left and right states of the primitive variables ! do p = 1, nfl call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) end do ! reconstruct the left and right states of the magnetic field components ! do p = ibx, ibz call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) end do #ifdef GLM ! reconstruct the left and right states of the scalar potential ! call reconstruct(n, h, q(iph,:), ql(iph,:), qr(iph,:)) ! obtain the state values for Bx and Psi for the GLM-MHD equations ! cl(:) = 0.5d0 * ((qr(ibx,:) + ql(ibx,:)) - (qr(iph,:) - ql(iph,:)) / cmax) cr(:) = 0.5d0 * ((qr(iph,:) + ql(iph,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) ql(ibx,:) = cl(:) qr(ibx,:) = cl(:) ql(iph,:) = cr(:) qr(iph,:) = cr(:) #endif /* GLM */ ! calculate conservative variables at states ! call prim2cons(n, ql, ul) call prim2cons(n, qr, ur) ! calculate fluxes and speeds ! call fluxspeed(n, ql, ul, fl, cl) call fluxspeed(n, qr, ur, fr, cr) ! iterate over all points and calculate the HLLD flux ! do i = 1, n ! calculate min and max and intermediate speeds: eq. (67) ! sl = min(ql(ivx,i) - cl(i), qr(ivx,i) - cr(i)) sr = max(ql(ivx,i) + cl(i), qr(ivx,i) + cr(i)) ! all speeds > 0, left side flux ! if (sl .ge. 0.0) then fn(:,i) = fl(:,i) ! all speeds < 0, right side flux ! else if (sr .le. 0.0) then fn(:,i) = fr(:,i) ! intermediate states ! else ! sl < 0 & sr > 0 ! product and difference of speeds ! srl = sr * sl srml = sr - sl ! density of the intermediate state (eq. 20 and 21) ! dnm = (sr * ur(idn,i) - sl * ul(idn,i) - fr(idn,i) + fl(idn,i)) / srml mxm = (sr * ur(imx,i) - sl * ul(imx,i) - fr(imx,i) + fl(imx,i)) / srml sqd = sqrt(dnm) ! fluxes for density and x-momentum are the same for all intermediate states (eq. 22 and 23) ! fn(idn,i) = (sr * fl(idn,i) - sl * fr(idn,i) & + srl * (ur(idn,i) - ul(idn,i))) / srml fn(imx,i) = (sr * fl(imx,i) - sl * fr(imx,i) & + srl * (ur(imx,i) - ul(imx,i))) / srml #ifdef GLM ! fluxes for parallel magnetic component and the scalar potential is the same ! as well fn(ibx,i) = (sr * fl(ibx,i) - sl * fr(ibx,i) & + srl * (ur(ibx,i) - ul(ibx,i))) / srml fn(iph,i) = (sr * fl(iph,i) - sl * fr(iph,i) & + srl * (ur(iph,i) - ul(iph,i))) / srml #endif /* GLM */ ! the speed of contact discontinuity (from eq. 15 and eq. 17) ! sm = fn(idn,i) / dnm ! Alfven speeds (eq. 29) ! sml = sm - abs(ql(ibx,i)) / sqd smr = sm + abs(qr(ibx,i)) / sqd ! calculate the left intermediate state ! u1l(idn) = dnm u1l(imx) = mxm div = (sl - sml) * (sl - smr) if (sm .eq. ql(ivx,i) .or. div .eq. 0.0 .or. ql(ibx,i) .eq. 0.0) then u1l(imy) = dnm * ql(ivy,i) u1l(imz) = dnm * ql(ivz,i) u1l(iby) = ql(iby,i) u1l(ibz) = ql(ibz,i) else fac = ql(ibx,i) * (sm - ql(ivx,i)) / div u1l(imy) = dnm * ql(ivy,i) - ql(iby,i) * fac u1l(imz) = dnm * ql(ivz,i) - ql(ibz,i) * fac fac = (ql(idn,i) * (sl - ql(ivx,i))**2 - ql(ibx,i)**2) & / (dnm * div) u1l(iby) = ql(iby,i) * fac u1l(ibz) = ql(ibz,i) * fac end if ! calculate the right intermediate state ! u1r(idn) = dnm u1r(imx) = mxm div = (sr - sml) * (sr - smr) if (sm .eq. qr(ivx,i) .or. div .eq. 0.0 .or. qr(ibx,i) .eq. 0.0) then u1r(imy) = dnm * qr(ivy,i) u1r(imz) = dnm * qr(ivz,i) u1r(iby) = qr(iby,i) u1r(ibz) = qr(ibz,i) else fac = qr(ibx,i) * (sm - qr(ivx,i)) / div u1r(imy) = dnm * qr(ivy,i) - qr(iby,i) * fac u1r(imz) = dnm * qr(ivz,i) - qr(ibz,i) * fac fac = (qr(idn,i) * (sr - qr(ivx,i))**2 - qr(ibx,i)**2) & / (dnm * div) u1r(iby) = qr(iby,i) * fac u1r(ibz) = qr(ibz,i) * fac end if ! intermediate discontinuities ! if (sml .ge. 0.0) then ! calculate the left intermediate flux ! fn(imy,i) = fl(imy,i) + sl * (u1l(imy) - ul(imy,i)) ! eq. (38) fn(imz,i) = fl(imz,i) + sl * (u1l(imz) - ul(imz,i)) fn(iby,i) = fl(iby,i) + sl * (u1l(iby) - ul(iby,i)) fn(ibz,i) = fl(ibz,i) + sl * (u1l(ibz) - ul(ibz,i)) else if (smr .le. 0.0) then ! calculate right intermediate flux ! fn(imy,i) = fr(imy,i) + sr * (u1r(imy) - ur(imy,i)) ! eq. (38) fn(imz,i) = fr(imz,i) + sr * (u1r(imz) - ur(imz,i)) fn(iby,i) = fr(iby,i) + sr * (u1r(iby) - ur(iby,i)) fn(ibz,i) = fr(ibz,i) + sr * (u1r(ibz) - ur(ibz,i)) else ! sml < 0 & smr > 0 ! normal component of magnetic field multiplied by sqrt(dnm) ! if (ql(ibx,i) .ge. 0.0) then bxs = sqd else bxs = - sqd end if ! calculate the intermediate state (eq. 34-37) ! u2(imy) = 0.5d0 * (u1r(imy) + u1l(imy) + bxs * (u1r(iby) - u1l(iby))) u2(imz) = 0.5d0 * (u1r(imz) + u1l(imz) + bxs * (u1r(ibz) - u1l(ibz))) u2(iby) = 0.5d0 * (u1r(iby) + u1l(iby) + (u1r(imy) - u1l(imy)) / bxs) u2(ibz) = 0.5d0 * (u1r(ibz) + u1l(ibz) + (u1r(imz) - u1l(imz)) / bxs) ! calculate the intermediate flux (eq. 24) ! fn(imy,i) = sm * u2(imy) - ql(ibx,i) * u2(iby) fn(imz,i) = sm * u2(imz) - ql(ibx,i) * u2(ibz) fn(iby,i) = sm * u2(iby) - ql(ibx,i) * u2(imy) / dnm fn(ibz,i) = sm * u2(ibz) - ql(ibx,i) * u2(imz) / dnm end if end if end do ! calculate numerical flux ! f( 1:nfl,2:n) = - fn( 1:nfl,2:n) + fn( 1:nfl,1:n-1) f(ibx:ibz,2:n) = - fn(ibx:ibz,2:n) + fn(ibx:ibz,1:n-1) #ifdef GLM f(iph ,2:n) = - fn(iph ,2:n) + fn(iph ,1:n-1) #endif /* GLM */ !------------------------------------------------------------------------------- ! end subroutine hlld #endif /* ISO */ #ifdef ADI ! !=============================================================================== ! ! hlld: subroutine computes the approximated flux using the HLLD method ! for the adiabatic equation of state ! !=============================================================================== ! subroutine hlld(n, h, u, f) use config , only : gamma use interpolation, only : reconstruct use variables , only : nvr, nfl, nqt use variables , only : idn, imx, imy, imz, ien, ivx, ivy, ivz, ipr use variables , only : ibx, iby, ibz #ifdef GLM use variables , only : iph #endif /* GLM */ implicit none ! input/output arguments ! integer , intent(in) :: n real , intent(in) :: h real, dimension(nvr,n), intent(in) :: u real, dimension(nqt,n), intent(out) :: f ! local variables ! integer :: p, i real, dimension(nvr,n) :: q, ql, qr, ul, ur real, dimension(nqt,n) :: fl, fr, fn real, dimension(n) :: cl, cr real, dimension(nvr) :: u1l, u1r, u2, q1l, q1r, q2 real :: sl, sr, slmv, srmv, slmm, srmm, sm, smvl, smvr & , sml, smr real :: ptl, ptr, pt, bx2, div, fac, bxs, dlsq, drsq ! !------------------------------------------------------------------------------- ! ! calculate the primitive variables ! call cons2prim(n, u, q) ! reconstruct the left and right states of the primitive variables ! do p = 1, nfl call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) end do ! reconstruct the left and right states of the magnetic field components ! do p = ibx, ibz call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:)) end do #ifdef GLM ! reconstruct the left and right states of the scalar potential ! call reconstruct(n, h, q(iph,:), ql(iph,:), qr(iph,:)) ! obtain the state values for Bx and Psi for the GLM-MHD equations ! cl(:) = 0.5d0 * ((qr(ibx,:) + ql(ibx,:)) - (qr(iph,:) - ql(iph,:)) / cmax) cr(:) = 0.5d0 * ((qr(iph,:) + ql(iph,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax) ql(ibx,:) = cl(:) qr(ibx,:) = cl(:) ql(iph,:) = cr(:) qr(iph,:) = cr(:) #endif /* GLM */ ! calculate conservative variables at states ! call prim2cons(n, ql, ul) call prim2cons(n, qr, ur) ! calculate fluxes and speeds ! call fluxspeed(n, ql, ul, fl, cl) call fluxspeed(n, qr, ur, fr, cr) ! iterate over all points and calculate the HLLD flux ! do i = 1, n ! calculate min and max and intermediate speeds: eq. (67) ! sl = min(ql(ivx,i) - cl(i), qr(ivx,i) - cr(i)) sr = max(ql(ivx,i) + cl(i), qr(ivx,i) + cr(i)) ! all speeds > 0, left side flux ! if (sl .ge. 0.0) then fn(:,i) = fl(:,i) ! all speeds < 0, right side flux ! else if (sr .le. 0.0) then fn(:,i) = fr(:,i) ! intermediate states ! else ! sl < 0 & sr > 0 ! calculate the total left and right pressures ! ptl = ql(ipr,i) + 0.5d0 * sum(ql(ibx:ibz,i)**2) ptr = qr(ipr,i) + 0.5d0 * sum(qr(ibx:ibz,i)**2) ! useful speed differences ! slmv = sl - ql(ivx,i) srmv = sr - qr(ivx,i) ! the speed of contact discontinuity (eq. 38, average from the both states) ! div = slmv * ql(idn,i) - srmv * qr(idn,i) slmm = (slmv * ul(imx,i) - srmv * ur(imx,i) - ptl + ptr) / div div = srmv * qr(idn,i) - slmv * ql(idn,i) srmm = (srmv * ur(imx,i) - slmv * ul(imx,i) - ptr + ptl) / div sm = 0.5d0 * (slmm + srmm) ! more useful speed differences ! slmm = sl - sm srmm = sr - sm smvl = sm - ql(ivx,i) smvr = sm - qr(ivx,i) bx2 = ql(ibx,i) * qr(ibx,i) ! pressure of the intermediate states (eq. 41) ! pt = 0.5d0 * (ptl + ptr + ql(idn,i) * slmv * smvl & + qr(idn,i) * srmv * smvr) ! calculate the left intermediate state variables ! q1l(idn) = ql(idn,i) * slmv / slmm q1l(ivx) = sm q1l(ibx) = ql(ibx,i) div = ql(idn,i) * slmv * slmm - bx2 if ((sm .eq. ql(ivx,i)) .or. (div .eq. 0.0) & .or. (bx2 .ge. gamma * ql(ipr,i)) & .or. (sl .eq. (ql(ivx,i) + cl(i))) & .or. (sl .eq. (ql(ivx,i) - cl(i)))) then q1l(ivy) = ql(ivy,i) q1l(ivz) = ql(ivz,i) q1l(iby) = ql(iby,i) q1l(ibz) = ql(ibz,i) else fac = ql(ibx,i) * smvl / div q1l(ivy) = ql(ivy,i) - ql(iby,i) * fac q1l(ivz) = ql(ivz,i) - ql(ibz,i) * fac fac = (ql(idn,i) * slmv**2 - bx2) / div q1l(iby) = ql(iby,i) * fac q1l(ibz) = ql(ibz,i) * fac end if ! convert the left intermediate state to the conservative form ! u1l(idn) = q1l(idn) u1l(imx) = q1l(idn) * q1l(ivx) u1l(imy) = q1l(idn) * q1l(ivy) u1l(imz) = q1l(idn) * q1l(ivz) if (slmm .ne. 0.0) then u1l(ien) = (slmv * ul(ien,i) - ptl * ql(ivx,i) + pt * sm & + ql(ibx,i) * (sum(ql(ivx:ivz,i) * ql(ibx:ibz,i)) & - sum(q1l(ivx:ivz) * q1l(ibx:ibz)))) / slmm else u1l(ien) = ul(ien,i) end if u1l(ibx) = q1l(ibx) u1l(iby) = q1l(iby) u1l(ibz) = q1l(ibz) #ifdef GLM u1l(iph) = ul(iph,i) #endif /* GLM */ ! calculate the right intermediate state variables ! q1r(idn) = qr(idn,i) * srmv / srmm q1r(ivx) = sm q1r(ibx) = qr(ibx,i) div = qr(idn,i) * srmv * srmm - bx2 if ((sm .eq. qr(ivx,i)) .or. (div .eq. 0.0) & .or. (bx2 .ge. gamma * qr(ipr,i)) & .or. (sr .eq. (qr(ivx,i) + cr(i))) & .or. (sr .eq. (qr(ivx,i) - cr(i)))) then q1r(ivy) = qr(ivy,i) q1r(ivz) = qr(ivz,i) q1r(iby) = qr(iby,i) q1r(ibz) = qr(ibz,i) else fac = qr(ibx,i) * smvr / div q1r(ivy) = qr(ivy,i) - qr(iby,i) * fac q1r(ivz) = qr(ivz,i) - qr(ibz,i) * fac fac = (qr(idn,i) * srmv**2 - bx2) / div q1r(iby) = qr(iby,i) * fac q1r(ibz) = qr(ibz,i) * fac end if ! convert the right intermediate state to the conservative form ! u1r(idn) = q1r(idn) u1r(imx) = q1r(idn) * q1r(ivx) u1r(imy) = q1r(idn) * q1r(ivy) u1r(imz) = q1r(idn) * q1r(ivz) if (srmm .ne. 0.0) then u1r(ien) = (srmv * ur(ien,i) - ptr * qr(ivx,i) + pt * sm & + qr(ibx,i) * (sum(qr(ivx:ivz,i) * qr(ibx:ibz,i)) & - sum(q1r(ivx:ivz) * q1r(ibx:ibz)))) / srmm else u1r(ien) = ur(ien,i) end if u1r(ibx) = q1r(ibx) u1r(iby) = q1r(iby) u1r(ibz) = q1r(ibz) #ifdef GLM u1r(iph) = ur(iph,i) #endif /* GLM */ ! Alfven speeds (eq. 51) ! sml = sm - abs(ql(ibx,i)) / sqrt(q1l(idn)) smr = sm + abs(qr(ibx,i)) / sqrt(q1r(idn)) ! intermediate discontinuities ! if (sml .ge. 0.0d0) then ! calculate the left intermediate flux ! fn(:,i) = fl(:,i) + sl * (u1l(:) - ul(:,i)) else if (smr .le. 0.0d0) then ! calculate the right intermediate flux ! fn(:,i) = fr(:,i) + sr * (u1r(:) - ur(:,i)) else ! sml < 0 & smr > 0 ! obtain the normal component of magnetic field ! if (ql(ibx,i) .gt. 0.0d0) then bxs = 1.0d0 else if (ql(ibx,i) .lt. 0.0d0) then bxs = -1.0d0 else bxs = 0.0d0 end if ! compute the density root squares ! dlsq = sqrt(q1l(idn)) drsq = sqrt(q1r(idn)) div = dlsq + drsq ! calculate the velocity components ! q2(ivx) = sm q2(ivy) = (dlsq * q1l(ivy) + drsq * q1r(ivy) & + (q1r(iby) - q1l(iby)) * bxs) / div q2(ivz) = (dlsq * q1l(ivz) + drsq * q1r(ivz) & + (q1r(ibz) - q1l(ibz)) * bxs) / div ! calculate the magnetic field components ! q2(ibx) = ql(ibx,i) q2(iby) = (dlsq * q1r(iby) + drsq * q1l(iby) & + dlsq * drsq * (q1r(ivy) - q1l(ivy)) * bxs) / div q2(ibz) = (dlsq * q1r(ibz) + drsq * q1l(ibz) & + dlsq * drsq * (q1r(ivz) - q1l(ivz)) * bxs) / div if (sm .ge. 0.0) then ! convert the left Alfven intermediate state to the conservative form ! u2(idn) = u1l(idn) u2(imx) = u1l(idn) * q2(ivx) u2(imy) = u1l(idn) * q2(ivy) u2(imz) = u1l(idn) * q2(ivz) u2(ien) = u1l(ien) - dlsq * (sum(q1l(ivx:ivz) * q1l(ibx:ibz)) & - sum(q2 (ivx:ivz) * q2 (ibx:ibz))) * bxs u2(ibx) = u1l(ibx) u2(iby) = q2(iby) u2(ibz) = q2(ibz) #ifdef GLM u2(iph) = u1l(iph) #endif /* GLM */ ! calculate the numerical flux ! fn(:,i) = fl(:,i) + sml * u2(:) - (sml - sl) * u1l(:) - sl * ul(:,i) else ! sm < 0 ! convert the right Alfven intermediate state to the conservative form ! u2(idn) = u1r(idn) u2(imx) = u1r(idn) * q2(ivx) u2(imy) = u1r(idn) * q2(ivy) u2(imz) = u1r(idn) * q2(ivz) u2(ien) = u1r(ien) + drsq * (sum(q1r(ivx:ivz) * q1r(ibx:ibz)) & - sum(q2 (ivx:ivz) * q2 (ibx:ibz))) * bxs u2(ibx) = u1r(ibx) u2(iby) = q2(iby) u2(ibz) = q2(ibz) #ifdef GLM u2(iph) = u1r(iph) #endif /* GLM */ ! calculate the numerical flux ! fn(:,i) = fr(:,i) + smr * u2(:) - (smr - sr) * u1r(:) - sr * ur(:,i) end if end if end if end do ! calculate the numerical flux derivative ! f( 1:nfl,2:n) = - fn( 1:nfl,2:n) + fn( 1:nfl,1:n-1) f(ibx:ibz,2:n) = - fn(ibx:ibz,2:n) + fn(ibx:ibz,1:n-1) #ifdef GLM f(iph ,2:n) = - fn(iph ,2:n) + fn(iph ,1:n-1) #endif /* GLM */ !------------------------------------------------------------------------------- ! end subroutine hlld #endif /* ADI */ #endif /* HLLD */ #endif /* MHD */ ! !=============================================================================== ! ! fluxspeed: subroutine computes fluxes and speeds for a given set of equations ! !=============================================================================== ! subroutine fluxspeed(n, q, u, f, c) use config , only : gamma, csnd, csnd2 use variables, only : nvr, nqt use variables, only : idn, imx, imy, imz, ivx, ivy, ivz #ifdef ADI use variables, only : ipr, ien #endif /* ADI */ #ifdef MHD use variables, only : ibx, iby, ibz #ifdef GLM use variables, only : iph #endif /* GLM */ #endif /* MHD */ implicit none ! input/output arguments ! integer , intent(in) :: n real, dimension(nvr,n), intent(in) :: q, u real, dimension(nqt,n), intent(out) :: f real, dimension(n) , intent(out) :: c ! local variables ! integer :: i real :: bb, pm, vb, cs, cb, ca ! !------------------------------------------------------------------------------- ! ! sweep over all points ! do i = 1, n ! compute fluxes ! f(idn,i) = u(imx,i) #ifdef ADI f(imx,i) = q(ivx,i) * u(imx,i) + q(ipr,i) #endif /* ADI */ #ifdef ISO f(imx,i) = q(ivx,i) * u(imx,i) + q(idn,i) * csnd2 #endif /* ISO */ f(imy,i) = q(ivx,i) * u(imy,i) f(imz,i) = q(ivx,i) * u(imz,i) #ifdef ADI f(ien,i) = q(ivx,i) * (u(ien,i) + q(ipr,i)) #endif /* ADI */ #ifdef MHD bb = sum(q(ibx:ibz,i) * q(ibx:ibz,i)) pm = 0.5 * bb vb = sum(q(ivx:ivz,i) * q(ibx:ibz,i)) f(imx,i) = f(imx,i) - q(ibx,i) * q(ibx,i) + pm f(imy,i) = f(imy,i) - q(ibx,i) * q(iby,i) f(imz,i) = f(imz,i) - q(ibx,i) * q(ibz,i) #ifdef ADI f(ien,i) = f(ien,i) + q(ivx,i) * pm - q(ibx,i) * vb #endif /* ADI */ #ifdef GLM f(ibx,i) = q(iph,i) f(iby,i) = q(ivx,i) * q(iby,i) - q(ibx,i) * q(ivy,i) f(ibz,i) = q(ivx,i) * q(ibz,i) - q(ibx,i) * q(ivz,i) f(iph,i) = q(ibx,i) #endif /* GLM */ #endif /* MHD */ ! compute speeds ! #ifdef MHD #ifdef ADI cs = gamma * q(ipr,i) #endif /* ADI */ #ifdef ISO cs = csnd2 * q(idn,i) #endif /* ISO */ cb = cs + bb ca = q(ibx,i) * q(ibx,i) c(i) = sqrt(0.5 * (cb + sqrt(max(0.0, cb * cb - 4.0 * cs * ca))) / q(idn,i)) #else /* MHD */ #ifdef ADI c(i) = sqrt(gamma * q(ipr,i) / q(idn,i)) #endif /* ADI */ #ifdef ISO c(i) = csnd #endif /* ISO */ #endif /* MHD */ end do !------------------------------------------------------------------------------- ! end subroutine fluxspeed #ifdef MHD ! !=============================================================================== ! ! emf: subroutine computes magnetic fluxes (electromotive force) ! !=============================================================================== ! subroutine emf(n, v, b, f) implicit none ! input/output arguments ! integer , intent(in) :: n real, dimension(3,n), intent(in) :: v, b real, dimension(3,n), intent(out) :: f ! local variables ! integer :: i ! !------------------------------------------------------------------------------- ! ! sweep over all points ! do i = 1, n f(1,i) = 0.0 f(2,i) = v(1,i) * b(2,i) - b(1,i) * v(2,i) f(3,i) = v(1,i) * b(3,i) - b(1,i) * v(3,i) end do ! !------------------------------------------------------------------------------- ! end subroutine emf #endif /* MHD */ ! !=============================================================================== ! ! cons2prim: subroutine converts primitive variables to conservative ! !=============================================================================== ! subroutine cons2prim(n, u, q) use config , only : gammam1 use variables, only : nvr use variables, only : idn, imx, imy, imz, ivx, ivy, ivz #ifdef ADI use variables, only : ipr, ien #endif /* ADI */ #ifdef MHD use variables, only : ibx, iby, ibz #ifdef GLM use variables, only : iph #endif /* GLM */ #endif /* MHD */ implicit none ! input/output arguments ! integer , intent(in) :: n real, dimension(nvr,n), intent(in) :: u real, dimension(nvr,n), intent(out) :: q ! local variables ! integer :: i real :: dni, ei, ek, em ! !------------------------------------------------------------------------------- ! do i = 1, n dni = 1.0 / u(idn,i) q(idn,i) = u(idn,i) q(ivx,i) = dni * u(imx,i) q(ivy,i) = dni * u(imy,i) q(ivz,i) = dni * u(imz,i) #ifdef ADI ek = 0.5 * sum(u(imx:imz,i) * q(ivx:ivz,i)) ei = u(ien,i) - ek #ifdef MHD em = 0.5 * sum(u(ibx:ibz,i) * u(ibx:ibz,i)) ei = ei - em #endif /* MHD */ q(ipr,i) = gammam1 * ei #endif /* ADI */ #ifdef MHD q(ibx,i) = u(ibx,i) q(iby,i) = u(iby,i) q(ibz,i) = u(ibz,i) #ifdef GLM q(iph,i) = u(iph,i) #endif /* GLM */ #endif /* MHD */ end do !------------------------------------------------------------------------------- ! end subroutine cons2prim ! !=============================================================================== ! ! prim2cons: subroutine converts primitive variables to conservative ! !=============================================================================== ! subroutine prim2cons(n, q, u) use config , only : gammam1i use variables, only : nvr use variables, only : idn, imx, imy, imz, ivx, ivy, ivz #ifdef ADI use variables, only : ipr, ien #endif /* ADI */ #ifdef MHD use variables, only : ibx, iby, ibz #ifdef GLM use variables, only : iph #endif /* GLM */ #endif /* MHD */ implicit none ! input/output arguments ! integer , intent(in) :: n real, dimension(nvr,n), intent(in) :: q real, dimension(nvr,n), intent(out) :: u ! local variables ! integer :: i real :: ei, ek, em ! !------------------------------------------------------------------------------- ! do i = 1, n u(idn,i) = q(idn,i) u(imx,i) = q(idn,i) * q(ivx,i) u(imy,i) = q(idn,i) * q(ivy,i) u(imz,i) = q(idn,i) * q(ivz,i) #ifdef ADI ei = gammam1i * q(ipr,i) ek = 0.5 * sum(u(imx:imz,i) * q(ivx:ivz,i)) u(ien,i) = ei + ek #endif /* ADI */ #ifdef MHD #ifdef ADI em = 0.5 * sum(q(ibx:ibz,i) * q(ibx:ibz,i)) u(ien,i) = u(ien,i) + em #endif /* ADI */ u(ibx,i) = q(ibx,i) u(iby,i) = q(iby,i) u(ibz,i) = q(ibz,i) #ifdef GLM u(iph,i) = q(iph,i) #endif /* GLM */ #endif /* MHD */ end do !------------------------------------------------------------------------------- ! end subroutine prim2cons ! !=============================================================================== ! ! maxspeed: function to calculate maximum speed in the system ! !=============================================================================== ! function maxspeed(u) use config , only : im, jm, km, ib, ie, jb, je, kb, ke #ifdef ADI use config , only : gamma #endif /* ADI */ #ifdef ISO use config , only : csnd, csnd2 #endif /* ISO */ use variables , only : nvr, nqt use variables , only : idn, ivx, ivz #ifdef ADI use variables , only : ipr #endif /* ADI */ #ifdef MHD use variables , only : ibx, iby, ibz #endif /* MHD */ implicit none ! input arguments ! real, dimension(nqt,im,jm,km), intent(in) :: u ! local variables ! integer :: i, j, k real :: vv, v, c #ifdef MHD real :: bb #endif /* MHD */ real :: maxspeed ! local arrays ! real, dimension(nvr,im) :: q ! !------------------------------------------------------------------------------- ! maxspeed = 0.0 ! iterate over all points and calculate maximum speed ! do k = kb, ke do j = jb, je call cons2prim(im, u(1:nqt,1:im,j,k), q(1:nqt,1:im)) do i = ib, ie ! calculate the velocity ! vv = sum(q(ivx:ivz,i)**2) v = sqrt(vv) #ifdef MHD bb = sum(q(ibx:ibz,i)**2) #endif /* MHD */ ! calculate the maximum characteristic speed ! #ifdef MHD #ifdef ADI c = sqrt((gamma * q(ipr,i) + bb) / q(idn,i)) #endif /* ADI */ #ifdef ISO c = sqrt(csnd2 + bb / q(idn,i)) #endif /* ISO */ #else /* MHD */ #ifdef ADI c = sqrt(gamma * q(ipr,i) / q(idn,i)) #endif /* ADI */ #ifdef ISO c = csnd #endif /* ISO */ #endif /* MHD */ ! calculate maximum of the speed ! maxspeed = max(maxspeed, v + c) end do end do end do ! !------------------------------------------------------------------------------- ! end function maxspeed !=============================================================================== ! end module