!> @file euler_riemann_2d.f90 !> @brief 4-component 2D Euler face fluxes (Rusanov, HLLC) with face-normal rotation. !! !! Q = (rho, rho*u, rho*v, E). For a face with normal in direction dir (1=x, 2=y), !! normal-momentum index mn and tangential mt: !! dir=1 -> mn=2 (rho*u), mt=3 (rho*v); dir=2 -> mn=3 (rho*v), mt=2 (rho*u). !! Pressure uses BOTH velocity components (full KE) -- this is why the 3-component !! 1D solvers cannot be reused directly. module euler_riemann_2d use precision, only: wp use solver_state_2d, only: neq2d use option_registry, only: flux_rusanov implicit none private public :: face_flux_2d, face_flux_2d_normal, wall_flux_2d_normal contains !> Inviscid solid-wall (no-penetration) face flux through unit normal nrm. !! At a solid wall v.n = 0, so the Euler flux collapses to the pressure term: !! F = [ 0, p*nx, p*ny, 0 ], !! with p the interior-side pressure and (nx,ny) the boundary-face normal. !! Mass and energy flux are EXACTLY zero by construction (a "weak"/pressure !! wall), so no spurious mass crosses the wall and no spurious entropy is !! generated by an acoustic over-pressure -- unlike a ghost-reflection wall, !! whose Riemann star pressure carries an O(rho*c*|v_n|) error wherever the !! reconstructed interior normal velocity is non-zero (always, on a curved !! wall). `q` is the interior-side (reconstructed) conserved state; `gam` the !! ratio of specific heats. nrm need not be normalised in sign: the residual's !! +xi/+eta divergence orientation makes the resulting wall pressure force act !! into the domain at both min and max edges. pure subroutine wall_flux_2d_normal(q, nrm, gam, flux) real(wp), intent(in) :: q(neq2d), nrm(2), gam real(wp), intent(out) :: flux(neq2d) real(wp) :: p p = (gam - 1.0_wp) * (q(4) - 0.5_wp * (q(2)**2 + q(3)**2) / q(1)) p = max(p, 0.0_wp) ! defense-in-depth: a non-physical face must not pull inward flux(1) = 0.0_wp flux(2) = p * nrm(1) flux(3) = p * nrm(2) flux(4) = 0.0_wp end subroutine wall_flux_2d_normal pure subroutine prim_split(q, dir, gam, rho, un, ut, p, c, mn, mt) real(wp), intent(in) :: q(neq2d), gam integer, intent(in) :: dir real(wp), intent(out) :: rho, un, ut, p, c integer, intent(out) :: mn, mt real(wp) :: u, v rho = q(1) u = q(2) / rho v = q(3) / rho p = (gam - 1.0_wp) * (q(4) - 0.5_wp * rho * (u * u + v * v)) ! Clamp the sound-speed radicand: a non-physical (negative-pressure) ! reconstructed face would otherwise give c = NaN that poisons the whole ! HLLC flux (sL/sR/star region). Mirrors the c_roe clamp in hllc_flux_2d; ! production paths already abort upstream on non-positive p (FVS/FDS guards), ! so this is defense-in-depth that keeps the flux finite for any input. c = sqrt(max(gam * p / rho, 0.0_wp)) if (dir == 1) then un = u; ut = v; mn = 2; mt = 3 else un = v; ut = u; mn = 3; mt = 2 end if end subroutine prim_split pure subroutine normal_phys_flux(q, dir, gam, f) real(wp), intent(in) :: q(neq2d), gam integer, intent(in) :: dir real(wp), intent(out) :: f(neq2d) real(wp) :: rho, un, ut, p, c integer :: mn, mt call prim_split(q, dir, gam, rho, un, ut, p, c, mn, mt) f(1) = rho * un f(mn) = rho * un * un + p f(mt) = rho * un * ut f(4) = un * (q(4) + p) end subroutine normal_phys_flux subroutine face_flux_2d(qL, qR, dir, scheme, gam, flux) real(wp), intent(in) :: qL(neq2d), qR(neq2d), gam integer, intent(in) :: dir character(len=*), intent(in) :: scheme real(wp), intent(out) :: flux(neq2d) if (dir /= 1 .and. dir /= 2) error stop 'euler_riemann_2d: dir must be 1 or 2' select case (trim(scheme)) case (flux_rusanov) call rusanov_flux_2d(qL, qR, dir, gam, flux) case default call hllc_flux_2d(qL, qR, dir, gam, flux) end select end subroutine face_flux_2d !> Face flux through an arbitrary UNIT normal nrm=(nx,ny). Rotates the !! conserved momentum into the face frame, reuses the axis solver (dir=1), !! then rotates the momentum flux back to (x,y). nrm must be a unit vector. subroutine face_flux_2d_normal(qL, qR, nrm, scheme, gam, flux) real(wp), intent(in) :: qL(neq2d), qR(neq2d), nrm(2), gam character(len=*), intent(in) :: scheme real(wp), intent(out) :: flux(neq2d) real(wp) :: qLr(neq2d), qRr(neq2d), fr(neq2d), nx, ny nx = nrm(1); ny = nrm(2) ! Rotate conserved momentum into (normal, tangential). qLr(1) = qL(1) qLr(2) = qL(2) * nx + qL(3) * ny qLr(3) = -qL(2) * ny + qL(3) * nx qLr(4) = qL(4) qRr(1) = qR(1) qRr(2) = qR(2) * nx + qR(3) * ny qRr(3) = -qR(2) * ny + qR(3) * nx qRr(4) = qR(4) ! Proven axis solver treats index 2 as normal-momentum, 3 as tangential. call face_flux_2d(qLr, qRr, 1, scheme, gam, fr) ! Rotate the momentum flux back to (x, y). flux(1) = fr(1) flux(2) = fr(2) * nx - fr(3) * ny flux(3) = fr(2) * ny + fr(3) * nx flux(4) = fr(4) end subroutine face_flux_2d_normal pure subroutine rusanov_flux_2d(qL, qR, dir, gam, flux) real(wp), intent(in) :: qL(neq2d), qR(neq2d), gam integer, intent(in) :: dir real(wp), intent(out) :: flux(neq2d) real(wp) :: fL(neq2d), fR(neq2d) real(wp) :: rhoL, unL, utL, pL, cL, rhoR, unR, utR, pR, cR, smax integer :: mn, mt call prim_split(qL, dir, gam, rhoL, unL, utL, pL, cL, mn, mt) call prim_split(qR, dir, gam, rhoR, unR, utR, pR, cR, mn, mt) call normal_phys_flux(qL, dir, gam, fL) call normal_phys_flux(qR, dir, gam, fR) smax = max(abs(unL) + cL, abs(unR) + cR) flux = 0.5_wp * (fL + fR) - 0.5_wp * smax * (qR - qL) end subroutine rusanov_flux_2d pure subroutine hllc_flux_2d(qL, qR, dir, gam, flux) real(wp), intent(in) :: qL(neq2d), qR(neq2d), gam integer, intent(in) :: dir real(wp), intent(out) :: flux(neq2d) real(wp) :: rhoL, unL, utL, pL, cL, rhoR, unR, utR, pR, cR real(wp) :: fL(neq2d), fR(neq2d) real(wp) :: sL, sR, sStar, factor, uStar(neq2d) real(wp) :: rt, un_roe, ut_roe, c_roe, hL, hR, h_roe integer :: mn, mt call prim_split(qL, dir, gam, rhoL, unL, utL, pL, cL, mn, mt) call prim_split(qR, dir, gam, rhoR, unR, utR, pR, cR, mn, mt) call normal_phys_flux(qL, dir, gam, fL) call normal_phys_flux(qR, dir, gam, fR) hL = (qL(4) + pL) / rhoL hR = (qR(4) + pR) / rhoR rt = sqrt(rhoR / rhoL) un_roe = (unL + rt * unR) / (1.0_wp + rt) ut_roe = (utL + rt * utR) / (1.0_wp + rt) h_roe = (hL + rt * hR) / (1.0_wp + rt) ! Clamp the Roe radicand at 0: for two physical (p>0) states Roe's identity ! keeps it >= 0, but a non-physical reconstructed face (negative pressure) ! can drive it negative, giving c_roe = NaN that poisons sL/sR and the whole ! flux. max(...,0) is a no-op for valid states and lets a bad face yield a ! finite flux (the upstream FDS positivity guard then aborts cleanly). c_roe = sqrt(max((gam - 1.0_wp) * (h_roe - 0.5_wp * (un_roe * un_roe + ut_roe * ut_roe)), 0.0_wp)) sL = min(unL - cL, un_roe - c_roe) sR = max(unR + cR, un_roe + c_roe) if (sL >= 0.0_wp) then flux = fL; return end if if (sR <= 0.0_wp) then flux = fR; return end if sStar = (pR - pL + rhoL * unL * (sL - unL) - rhoR * unR * (sR - unR)) / & (rhoL * (sL - unL) - rhoR * (sR - unR)) if (sStar >= 0.0_wp) then factor = rhoL * (sL - unL) / (sL - sStar) uStar(1) = factor uStar(mn) = factor * sStar uStar(mt) = factor * utL uStar(4) = factor * (qL(4) / rhoL + (sStar - unL) * (sStar + pL / (rhoL * (sL - unL)))) flux = fL + sL * (uStar - qL) else factor = rhoR * (sR - unR) / (sR - sStar) uStar(1) = factor uStar(mn) = factor * sStar uStar(mt) = factor * utR uStar(4) = factor * (qR(4) / rhoR + (sStar - unR) * (sStar + pR / (rhoR * (sR - unR)))) flux = fR + sR * (uStar - qR) end if end subroutine hllc_flux_2d end module euler_riemann_2d