!------------------------------------------------------- ! Compute the Eigenvalues for the M1 radiative ! transfer HLL Riemann solver. ! ! Original file sent by Matthias González, was ! used in HERACLES code ! (ui.adsabs.harvard.edu/abs/2007A%26A...464..429G) ! under the CeCILL licence, which is compatible with ! the GNU GPL licence. !------------------------------------------------------- program compute_RT_riemann_solver_eigenvalues implicit none integer, parameter :: dp = SELECTED_REAL_KIND(15) integer, parameter :: ndim = 3 real(dp), parameter :: pi = 4.D0*DATAN(1.d0) real(dp), parameter :: c = 1.0 ! How many points to go through? ! NOTE: we always also include zero, so the resulting ! arrays will have size n_points + 1 integer, parameter :: n_points = 99 ! Make this array allocatable because it might be too ! large for the stack. real(dp), dimension(:,:,:), allocatable :: eigenvals ! eigenvalues. allocate(eigenvals(0:n_points, 0:n_points, 1:4)) call compute_eigenvals() call write_eigenvals() contains ! Subroutine COMPUTE_eigenvals ! !> Computes the tabulated eigenvalues for the M1 radiative !! transfer solver. !< subroutine compute_eigenvals implicit none integer :: i,j real(dp) :: normef ! |F| real(dp) :: theta ! Arguments needed for external eigenvalue calculation routine real(dp), dimension(4) :: WR real(dp), dimension(4) :: WI integer,parameter :: LDVL=1 integer,parameter :: LDVR=1 real(dp), dimension(LDVL,4) :: VL real(dp), dimension(LDVR,4) :: VR integer :: info integer :: lda=4 integer,parameter :: lwork=20 real(dp),dimension(lwork) :: work real(dp) :: E ! Radiation energy real(dp),dimension(1:3 ) :: F ! Radiation flux real(dp),dimension(3,3 ) :: Dedd ! Eddington tensor real(dp),dimension(3,3 ) :: Dedd_dE ! dDedd/dE * E real(dp),dimension(3,3,1:3) :: Dedd_dF ! dDedd/dF * E real(dp), dimension(4,4) :: mat ! Jacobi matrix E=1. ! Photon Energy F=0. ! Photon Flux do i=0,n_points ! loop over ||f|| normef=float(i)/n_points do j=0,n_points ! loop over 0 < theta < pi: angle w.r.t. interface theta=j*pi/n_points F(1)=c*E*cos(theta)*normef F(2)=c*E*sin(theta)*normef ! Compute Eddington tensor and its derivatives ! (derivatives are multiplied by Ei) call cal_Dedd(E,F,Dedd,Dedd_dE,Dedd_dF) ! Get Jacobi Matrix (for (E, Fx) only, but using ! gradients in up to all 3 dimensions) mat=0. mat(1,1) = 0. mat(1,2) = 1. mat(2,1) = c**2*(Dedd(1,1)+Dedd_dE(1,1 )) mat(2,2) = c**2* Dedd_dF(1,1,1) if(ndim.gt.1) then mat(1,3) = 0. mat(2,3) = c**2* Dedd_dF(1,1,2) mat(3,1) = c**2*(Dedd(1,2)+Dedd_dE(1,2 )) mat(3,2) = c**2* Dedd_dF(1,2,1) mat(3,3) = c**2* Dedd_dF(1,2,2) endif if(ndim.gt.2) then mat(1,4) = 0. mat(2,4) = c**2* Dedd_dF(1,1,3) mat(3,4) = c**2* Dedd_dF(1,2,3) mat(4,1) = c**2*(Dedd(1,3)+Dedd_dE(1,3 )) mat(4,2) = c**2* Dedd_dF(1,3,1) mat(4,3) = c**2* Dedd_dF(1,3,2) mat(4,4) = c**2* Dedd_dF(1,3,3) endif ! normalized the matrix to get eigenvalues in units of c mat(2:4,1)=mat(2:4,1)/c**2 mat(2:4,2:4)=mat(2:4,2:4)/c ! DGEEV computes the eigenvalues and, optionally, the left ! and/or right eigenvectors for GE matrices. It's a function ! from LAPACK ! Warning : dsyev works only for symmetric matrix -> use dgeev instead call dgeev('N', & ! Don't compute left eigenvectors 'N', & ! Don't compute right eigenvectors 4, & ! Order of matrix mat, & ! Matrix for which to calculate the Eigenvalues lda, & ! leading dimension of array mat. (which index to start with) WR, & ! Real parts of eigenvalues WI, & ! Imaginary parts of eigenvalues VL, & ! array for left eigenvectors (unused) LDVL, & ! leading dimension for VL array VR, & ! right eigenvectors (unused) LDVR, & ! leading dimension for VR array WORK, & ! workspace for this routine to work with. Should be empty array LWORK,& ! Dimension of array WORK. LWORK >= max(1, 3*N) info & ! error code. Success if info = 0 ) ! sort the real part of the eigenvalues in ascending order ! ! warning : the eigenvalues may be not real but complex... ! in particular it happens when eps->0 and f->1 ! note for Mladen: I took out the eps part call bubble_sort(4,WR) eigenvals(i,j,1:4)=WR(1:4) enddo enddo return end subroutine compute_eigenvals !########################################################### !########################################################### !########################################################### !########################################################### subroutine bubble_sort(n,arr) implicit none integer :: n,i,j real(dp) :: a real(dp), dimension(n) :: arr do j = 1,n do i = 2,n if(arr(i) < arr(i-1))then a = arr(i ) arr(i ) = arr(i-1) arr(i-1) = a endif enddo enddo return end subroutine bubble_sort !########################################################### !########################################################### !########################################################### !########################################################### ! Subroutine CAL_DEDD ! !> Computes the Eddington tensor Dedd and its derivatives !! with respect to Er and Fr, Dedd_dE and Dedd_dF, !! respectively. !! We'll need them to compute the Jacobi-Matrix dF/dU, !! where !! U = (E_photon, F_photon); !! F = (F_photon, c^2 P) !! P = D E_photon. !! The derivatives of P will be done by part, e.g. !! dP/dE_photon = D + dD/dE * E !! which is why we directly multiply the derivatives by E. !< subroutine cal_Dedd(E,F,Dedd,Dedd_dE,Dedd_dF) implicit none real(dp) , intent(in ) :: E ! Radiation energy real(dp),dimension(1:ndim ), intent(in ) :: F ! Radiation flux real(dp),dimension(3,3 ), intent(out) :: Dedd ! Eddington tensor real(dp),dimension(3,3 ), intent(out) :: Dedd_dE ! dD / dE * Ei real(dp),dimension(3,3,1:ndim), intent(out) :: Dedd_dF ! dD / dF * Ei integer :: i real(dp):: normeF ! || F || real(dp):: chi ! chi parameter for Eddington tensor real(dp):: chiprime ! dchi/df real(dp):: g ! diagonal part of Edd. tensor real(dp):: h ! anisotropic part of Edd. tensor real(dp):: gprime ! dg / df real(dp):: hprime ! dh / df real(dp):: nx ! x component of unit vector real(dp):: ny ! y component of unit vector real(dp):: nz ! z component of unit vector real(dp):: fx ! reduced flux x component real(dp):: fy ! reduced flux y component real(dp):: fz ! reduced flux z component real(dp):: cE ! c * E Dedd = 0.d0 ; Dedd_dE = 0.d0 ; Dedd_dF = 0.d0 fx = 0.d0 ; fy = 0.d0 ; fz = 0.d0 cE = c * E fx = F(1)/cE if(ndim.gt.1) fy = F(2)/cE if(ndim.gt.2) fz = F(3)/cE normef = 0. do i=1,ndim normef = normef + (F(i)/cE)**2 enddo normef = sqrt(normef) if(normef.gt.1) normef = 1.d0 ! For chi, f <= 1 is required chi = (3.d0+4.d0*normef**2) / (5.d0+2.d0*sqrt(4.d0-3.d0*normef**2)) chiprime = 2.d0*normef / sqrt(4.d0 - 3.d0 * normef**2) ! dchi/df g = 0.5d0 * (1.d0 - chi) ! diagonal part of tensor h = 0.5d0 * (3.d0 * chi-1.d0) ! anisotropic part of tensor gprime = -chiprime * 0.5d0 ! dg/df hprime = 3.d0 * chiprime * 0.5d0 ! dh/df nx = 0.d0 ; ny = 0.d0 ; nz = 0.d0 ! unit vector elements ! n is a unit vector normef = 0.d0 do i=1,ndim normef = normef + (F(i)/cE)**2 enddo normef = sqrt(normef) if(normef.gt.1.e-8) then nx = fx/normef ny = fy/normef nz = fz/normef endif ! Eddington tensor Dedd(1,1) = g + h*nx*nx ; Dedd(1,2) = h*nx*ny ; Dedd(1,3) = h*nx*nz Dedd(2,1) = h*ny*nx ; Dedd(2,2) = g + h*ny*ny ; Dedd(2,3) = h*ny*nz Dedd(3,1) = h*nz*nx ; Dedd(3,2) = h*nz*ny ; Dedd(3,3) = g + h*nz*nz ! (Eddington tensor derivative with respect to Ei) * Ei ! remember dchi/dE = dchi/df df/dE = - dchi/df |F|/(cE^2) = -dchi/df * f / E Dedd_dE(1,1) = - normef * (gprime + hprime*nx*nx) Dedd_dE(1,2) = - normef * ( hprime*nx*ny) Dedd_dE(1,3) = - normef * ( hprime*nx*nz) Dedd_dE(2,1) = - normef * ( hprime*ny*nx) Dedd_dE(2,2) = - normef * (gprime + hprime*ny*ny) Dedd_dE(2,3) = - normef * ( hprime*ny*nz) Dedd_dE(3,1) = - normef * ( hprime*nz*nx) Dedd_dE(3,2) = - normef * ( hprime*nz*ny) Dedd_dE(3,3) = - normef * (gprime + hprime*nz*nz) if(normef.gt.1.e-8) then ! (Eddington tensor derivative with respect to Fx) * Ei ! remember ! dchi/dFx = dchi/df df/dFx ! df/dFx = nx / (c E) ! dfx/dFx = 1 / (c E) ! dg/dFx = dg/df df/dFx ! d(h nx^2)/dFx = dh/df df/dFx nx^2 + h d(nx^2) / dFx ! d(nx^2)/dFx = 2 nx dnx/dfx dfx/dFx ! dnx / dfx = 1/f (1 - nx^2) ! dny / dFx = - nx ny / (f c E) ! so that ! d/dFx (g + h nx^2) = 1/(c E) [ dg/df nx + dh/df nx^3 + 2 h nx (1 - nx^2) / f ] ! d/dFx (h nx ny) = 1/(c E) [ dh/df nx^2 ny + h ny (1 - 2 nx^2) / f ] ! d/dFx (g + h ny ny) = 1/(c E) [ dg/df nx + dh/df nx ny^2 - 2 h nx ny^2 / f ] Dedd_dF(1,1,1) = ( gprime*nx + hprime*nx*nx*nx + h*(-2.d0*nx*nx*nx /normef+2.d0*nx/normef) )/c Dedd_dF(1,2,1) = ( hprime*nx*nx*ny + h*(-2.d0*nx*nx*ny /normef+ ny/normef) )/c Dedd_dF(1,3,1) = ( hprime*nx*nx*nz + h*(-2.d0*nx*nx*nz /normef+ nz/normef) )/c Dedd_dF(2,1,1) = ( hprime*nx*ny*nx + h*(-2.d0*nx*ny*nx /normef+ ny/normef) )/c Dedd_dF(2,2,1) = ( gprime*nx + hprime*nx*ny*ny + h*(-2.d0*nx*ny*ny /normef ) )/c Dedd_dF(2,3,1) = ( hprime*nx*ny*nz + h*(-2.d0*nx*ny*nz /normef ) )/c Dedd_dF(3,1,1) = ( hprime*nx*nz*nx + h*(-2.d0*nx*nz*nx /normef+ nz/normef) )/c Dedd_dF(3,2,1) = ( hprime*nx*nz*ny + h*(-2.d0*nx*nz*ny /normef ) )/c Dedd_dF(3,3,1) = ( gprime*nx + hprime*nx*nz*nz + h*(-2.d0*nx*nz*nz /normef ) )/c if(ndim.gt.1) then ! (Eddington tensor derivative with respect to Fy) * Ei Dedd_dF(1,1,2) = ( gprime*ny + hprime*ny*nx*nx + h*(-2.d0*ny*nx*nx /normef ) )/c Dedd_dF(1,2,2) = ( hprime*ny*nx*ny + h*(-2.d0*ny*nx*ny /normef+ nx/normef) )/c Dedd_dF(1,3,2) = ( hprime*ny*nx*nz + h*(-2.d0*ny*nx*nz /normef ) )/c Dedd_dF(2,1,2) = ( hprime*ny*ny*nx + h*(-2.d0*ny*ny*nx /normef+ nx/normef) )/c Dedd_dF(2,2,2) = ( gprime*ny + hprime*ny*ny*ny + h*(-2.d0*ny*ny*ny /normef+2.d0*ny/normef) )/c Dedd_dF(2,3,2) = ( hprime*ny*ny*nz + h*(-2.d0*ny*ny*nz /normef+ nz/normef) )/c Dedd_dF(3,1,2) = ( hprime*ny*nz*nx + h*(-2.d0*ny*nz*nx /normef ) )/c Dedd_dF(3,2,2) = ( hprime*ny*nz*ny + h*(-2.d0*ny*nz*ny /normef+ nz/normef) )/c Dedd_dF(3,3,2) = ( gprime*ny + hprime*ny*nz*nz + h*(-2.d0*ny*nz*nz /normef ) )/c if(ndim.gt.2) then ! (Eddington tensor derivative with respect to Fz) * Ei Dedd_dF(1,1,3) = ( gprime*nz + hprime*nz*nx*nx + h*(-2.d0*nz*nx*nx /normef ) )/c Dedd_dF(1,2,3) = ( hprime*nz*nx*ny + h*(-2.d0*nz*nx*ny /normef ) )/c Dedd_dF(1,3,3) = ( hprime*nz*nx*nz + h*(-2.d0*nz*nx*nz /normef+ nx/normef) )/c Dedd_dF(2,1,3) = ( hprime*nz*ny*nx + h*(-2.d0*nz*ny*nx /normef ) )/c Dedd_dF(2,2,3) = ( gprime*nz + hprime*nz*ny*ny + h*(-2.d0*nz*ny*ny /normef ) )/c Dedd_dF(2,3,3) = ( hprime*nz*ny*nz + h*(-2.d0*nz*ny*nz /normef+ ny/normef) )/c Dedd_dF(3,1,3) = ( hprime*nz*nz*nx + h*(-2.d0*nz*nz*nx /normef+ nx/normef) )/c Dedd_dF(3,2,3) = ( hprime*nz*nz*ny + h*(-2.d0*nz*nz*ny /normef+ ny/normef) )/c Dedd_dF(3,3,3) = ( gprime*nz + hprime*nz*nz*nz + h*(-2.d0*nz*nz*nz /normef+2.d0*nz/normef) )/c endif endif endif return end subroutine cal_Dedd !########################################################### !########################################################### !########################################################### !########################################################### ! Write the results. subroutine write_eigenvals implicit none integer :: i, j ! Write for C input !-------------------- ! Note that the eigenvalues are sorted by increasing value, ! and we only need the min and the max. So we only write ! the indices 1 and 4. ! You can use this output to directly copy-paste eigenvalues ! into /src/rt/GEAR/rt_riemann_HLL_eigenvalues.h ! ! write(*, "(A2)") "{ " ! do i = 0, n_points ! write(*, "(A4)") " { " ! do j = 0, n_points ! write(*, "(A5, E12.6, A3, E12.6, A3)") " {", eigenvals(i,j,1), " , ", eigenvals(i,j,4), " }," ! enddo ! write(*, "(A4)") " }," ! enddo ! write(*, "(A2)") "};" ! Write for yourself !-------------------- ! do i = 0, n_points ! do j = 0, n_points ! write(*, "(I3,x,I3,x,4(E12.6,x))") i, j, eigenvals(i, j, :) ! enddo ! enddo ! Write for python3 plotting !-------------------------- open(unit=1, file="eigenvals.txt") write(1, "(A2,I5)") "# ", n_points + 1 do i = 0, n_points do j = 0, n_points write(1, "(4(E12.6,x))") eigenvals(i, j, :) enddo enddo end subroutine write_eigenvals end program compute_RT_riemann_solver_eigenvalues