c Mahc_Evolve: Numerical Evolution of the General Relativistic Hydro Equations c Copyright (C) 2001 Mark Miller /*@@ @file vis_contrib.F @date December 1999 @author Mark Miller @desc Viscocity part for HRSC hydro update. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" /*@@ @routine vis_contrib @date December 1999 @author Mark Miller @desc Viscocity part for HRSC hydro update. @enddesc @calls @calledby @history @endhistory @@*/ ccc#define DEBUG_EIGEN #define FAST_MATRIX_INVERT subroutine vis_contrib(dens_vis,sx_vis,sy_vis,sz_vis,tau_vis, . dens,sx,sy,sz,tau,rho,velx,vely,velz, . eps,press,gxx,gxy,gxz,gyy,gyz,gzz, . uxx,uxy,uxz,uyy,uyz,uzz,det,dpdrho,dpdeps,alp, . dens_r,sx_r,sy_r,sz_r,tau_r,dens_l,sx_l,sy_l,sz_l,tau_l, . betax,betay,betaz,fmi,some_hydro_dust,some_hydro_mahc) implicit none CCTK_REAL dens_vis,sx_vis,sy_vis,sz_vis,tau_vis, . dens,sx,sy,sz,tau,rho,velx,vely,velz, . eps,press,gxx,gxy,gxz,gyy,gyz,gzz, . uxx,uxy,uxz,uyy,uyz,uzz,det,alp,dpdrho,dpdeps, . dens_r,sx_r,sy_r,sz_r,tau_r,dens_l,sx_l,sy_l,sz_l,tau_l, . betax,betay,betaz,bxx,bxy,bxz,byx,byy,byz, . bzx,bzy,bzz integer fmi logical some_hydro_dust,some_hydro_mahc DECLARE_CCTK_PARAMETERS c LOCAL VARS CCTK_REAL lam1,lam2,lam3,lamp,lamm,lam(5) CCTK_REAL eivec1(5),eivec2(5),eivec3(5),eivecp(5),eivecm(5) CCTK_REAL cs2,w,one,v2,two,vlowx,vlowy,vlowz CCTK_REAL p(5,5),dw(5),du(5),vis_test(5),vis(5),aa(5,5),dd integer i,j,k,l,indx(5) CCTK_REAL lamp_nobeta,lamm_nobeta,enthalpy,maxtestval #ifdef DEBUG_EIGEN CCTK_REAL dfxdy(5,5),df0dy(5,5),dfxdf0(5,5),test(5,5) CCTK_REAL vv(5,5),ww(5),maxvis #endif CCTK_REAL paug(5,6),temp,f_du(5) integer ii,jj,kk one = 1.0d0 two = 2.0d0 cs2 = (dpdrho + press * dpdeps / (rho**2))/ . (1.0d0 + eps + press/rho) vlowx = gxx*velx + gxy*vely + gxz*velz vlowy = gxy*velx + gyy*vely + gyz*velz vlowz = gxz*velx + gyz*vely + gzz*velz v2 = vlowx*velx + vlowy*vely + vlowz*velz w = one/sqrt(one - v2) c PERFECT FLUID if(some_hydro_mahc) then c EIGENVALUES lam1 = velx - betax/alp lam2 = velx - betax/alp lam3 = velx - betax/alp lamp_nobeta = (velx*(one-cs2) + sqrt(cs2*(one-v2)* . (uxx*(one-v2*cs2) - velx*velx*(one-cs2))))/(one-v2*cs2) lamm_nobeta = (velx*(one-cs2) - sqrt(cs2*(one-v2)* . (uxx*(one-v2*cs2) - velx*velx*(one-cs2))))/(one-v2*cs2) lamp = lamp_nobeta - betax/alp lamm = lamm_nobeta - betax/alp lam(1) = lam1 lam(2) = lam2 lam(3) = lam3 lam(4) = lamp lam(5) = lamm enthalpy = one + eps + press/rho c EIGENVECTORS c eigenvector # 1 eivec1(1) = dpdeps / ( (one + eps)*w*dpdeps - rho*w*dpdrho) eivec1(2) = vlowx eivec1(3) = vlowy eivec1(4) = vlowz eivec1(5) = one - eivec1(1) c someones wrong eigenvector 2 c eivec2(1) = w * vely c eivec2(2) = two * enthalpy * w * w * vlowx * vely c eivec2(3) = enthalpy * (one + two * w * w * vlowy * vely) c eivec2(4) = two * enthalpy * w * w * vlowz * vely c eivec2(5) = w * vely * (two * enthalpy * w - one) c MARKS COOL EIGENVECTOR 2 eivec2(1) = w * vlowy eivec2(2) = (enthalpy) * (gxy + two * w * w * vlowx * vlowy) eivec2(3) = (enthalpy) * (gyy + two * w * w * vlowy * vlowy) eivec2(4) = (enthalpy) * (gyz + two * w * w * vlowy * vlowz) eivec2(5) = vlowy * w * (two * w * (enthalpy) - one) c someones wrong eigenvector 3 c eivec3(1) = w * velz c eivec3(2) = two * enthalpy * w * w * vlowx * velz c eivec3(3) = two * enthalpy * w * w * vlowy * velz c eivec3(4) = enthalpy * (one + two * w * w * vlowz * velz) c eivec3(5) = w * velz * (two * enthalpy * w - one) c MARKS COOL EIGENVECTOR 3 eivec3(1) = w * vlowz eivec3(2) = (enthalpy) * (gxz + two * w * w * vlowx * vlowz) eivec3(3) = (enthalpy) * (gyz + two * w * w * vlowy * vlowz) eivec3(4) = (enthalpy) * (gzz + two * w * w * vlowz * vlowz) eivec3(5) = vlowz * w * (two * w * (enthalpy) - one) c + eigenvector eivecp(1) = one eivecp(2) = (enthalpy) * w * ( vlowx - . (velx - lamp_nobeta)/(uxx - velx * lamp_nobeta)) eivecp(3) = (enthalpy) * w * vlowy eivecp(4) = (enthalpy) * w * vlowz eivecp(5) = (enthalpy) * w * (uxx - velx*velx)/ . (uxx - velx*lamp_nobeta) - one c - eigenvector eivecm(1) = one eivecm(2) = (enthalpy) * w * ( vlowx - . (velx - lamm_nobeta)/(uxx - velx * lamm_nobeta)) eivecm(3) = (enthalpy) * w * vlowy eivecm(4) = (enthalpy) * w * vlowz eivecm(5) = (enthalpy) * w * (uxx - velx*velx)/ . (uxx - velx*lamm_nobeta) - one endif c DUST if(some_hydro_dust) then c EIGENVALUES lam1 = velx - betax/alp lam2 = velx - betax/alp lam3 = velx - betax/alp lamp_nobeta = velx lamm_nobeta = velx lamp = lamp_nobeta - betax/alp lamm = lamm_nobeta - betax/alp lam(1) = lam1 lam(2) = lam2 lam(3) = lam3 lam(4) = lamp lam(5) = lamm c EIGENVECTORS c eigenvector # 1 eivec1(1) = 1.0d0 eivec1(2) = 0.0d0 eivec1(3) = 0.0d0 eivec1(4) = 0.0d0 eivec1(5) = 0.0d0 c eigenvector # 2 eivec2(1) = 0.0d0 eivec2(2) = 0.0d0 eivec2(3) = 1.0d0 eivec2(4) = 0.0d0 eivec2(5) = 0.0d0 c eigenvector # 3 eivec3(1) = 0.0d0 eivec3(2) = 0.0d0 eivec3(3) = 0.0d0 eivec3(4) = 1.0d0 eivec3(5) = 0.0d0 c + eigenvector eivecp(1) = 0.0d0 eivecp(2) = 1.0d0 eivecp(3) = 0.0d0 eivecp(4) = 0.0d0 eivecp(5) = 0.0d0 c - eigenvector eivecm(1) = 0.0d0 eivecm(2) = 0.0d0 eivecm(3) = 0.0d0 eivecm(4) = 0.0d0 eivecm(5) = 1.0d0 endif c PUT EIGENVECTORS IN THE P MATRIX do i=1,5 p(i,1) = eivec1(i) p(i,2) = eivec2(i) p(i,3) = eivec3(i) p(i,4) = eivecp(i) p(i,5) = eivecm(i) enddo c here, check that lam and p actually solve the eigensystem c defined by the matrix dF^x/dF^0. Numerically construct c this matrix, and check: #ifdef DEBUG_EIGEN dfxdy(1,1) = (-(betax/alp)+velx)*w dfxdy(1,2) = rho*(w+(-(betax/alp)+velx)* . (gxx*velx+gxy*vely+gxz*velz)*w**3) dfxdy(1,3) = rho*(-(betax/alp)+velx)* . (gxy*velx+gyy*vely+gyz*velz)*w**3 dfxdy(1,4) = rho*(-(betax/alp)+velx)* . (gxz*velx+gyz*vely+gzz*velz)*w**3 dfxdy(1,5) = 0.0d0 dfxdy(2,1) = dpdrho+((1.000d0+dpdrho+eps)* & (-betax+alp*velx)*(gxx*velx+gxy*v & ely+gxz*velz)*w**2)/alp dfxdy(2,2) = ((press+rho+eps*rho)*w**2* & (-(betax*gxx)+alp*(2.000d0*gxx*velx+ & gxy*vely+gxz*velz)-2.000d0*(betax-alp*velx)*(gxx*velx+gxy*vely+g & xz*velz)**2*w**2))/alp dfxdy(2,3) = ((press+rho+eps*rho)*(-betax+alp*velx)* & w*(gxy*w+2.000d0*(gxx*v & elx+gxy*vely+gxz*velz)*(gxy*velx+gyy*vely+gyz*velz)*w**3))/alp dfxdy(2,4) = ((press+rho+eps*rho)*(-betax+alp*velx)* & w*(gxz*w+2.000d0*(gxx*v & elx+gxy*vely+gxz*velz)*(gxz*velx+gyz*vely+gzz*velz)*w**3))/alp dfxdy(2,5) = dpdeps+((dpdeps+rho)* & (-betax+alp*velx)*(gxx*velx+gxy*vely+gxz* & velz)*w**2)/alp dfxdy(3,1) = ((1.000d0+dpdrho+eps)*(-betax+alp*velx)* & (gxy*velx+gyy*vely+gyz & *velz)*w**2)/alp dfxdy(3,2) = ((press+rho+eps*rho)*w**2* & (-(betax*gxy)+alp*(2.000d0*gxy*velx+ & gyy*vely+gyz*velz)-2.000d0*(betax-alp*velx)*(gxx*velx+gxy*vely+g & xz*velz)*(gxy*velx+gyy*vely+gyz*velz)*w**2))/alp dfxdy(3,3) = ((press+rho+eps*rho)*(-betax+alp*velx)*w* & (gyy*w+2.000d0*(gxy*v & elx+gyy*vely+gyz*velz)**2*w**3))/alp dfxdy(3,4) = ((press+rho+eps*rho)*(-betax+alp*velx)*w* & (gyz*w+2.000d0*(gxy*v & elx+gyy*vely+gyz*velz)*(gxz*velx+gyz*vely+gzz*velz)*w**3))/alp dfxdy(3,5) = ((dpdeps+rho)*(-betax+alp*velx)* & (gxy*velx+gyy*vely+gyz*velz)*w & **2)/alp dfxdy(4,1) = ((1.000d0+dpdrho+eps)*(-betax+alp*velx)* & (gxz*velx+gyz*vely+gzz & *velz)*w**2)/alp dfxdy(4,2) = ((press+rho+eps*rho)*w**2* & (-(betax*gxz)+alp*(2.000d0*gxz*velx+ & gyz*vely+gzz*velz)-2.000d0*(betax-alp*velx)*(gxx*velx+gxy*vely+g & xz*velz)*(gxz*velx+gyz*vely+gzz*velz)*w**2))/alp dfxdy(4,3) = ((press+rho+eps*rho)*(-betax+alp*velx)*w* & (gyz*w+2.000d0*(gxy*v & elx+gyy*vely+gyz*velz)*(gxz*velx+gyz*vely+gzz*velz)*w**3))/alp dfxdy(4,4) = ((press+rho+eps*rho)*(-betax+alp*velx)*w* & (gzz*w+2.000d0*(gxz*v & elx+gyz*vely+gzz*velz)**2*w**3))/alp dfxdy(4,5) = ((dpdeps+rho)*(-betax+alp*velx)* & (gxz*velx+gyz*vely+gzz*velz)*w & **2)/alp dfxdy(5,1) = (betax*dpdrho+(betax-alp*velx)*w- & (1.000d0+dpdrho+eps)*(betax-a & lp*velx)*w**2)/alp dfxdy(5,2) = (w*(-(alp*rho)+alp*(press+rho+eps*rho)*w+ & rho*(betax-alp*velx)* & (gxx*velx+gxy*vely+gxz*velz)*w**2-2.000d0*(press+rho+eps*rho)*(b & etax-alp*velx)*(gxx*velx+gxy*vely+gxz*velz)*w**3))/alp dfxdy(5,3) = ((-betax+alp*velx)* & (gxy*velx+gyy*vely+gyz*velz)*w**3*(-rho+2.0 & 00d0*(press+rho+eps*rho)*w))/alp dfxdy(5,4) = ((-betax+alp*velx)* & (gxz*velx+gyz*vely+gzz*velz)*w**3*(-rho+2.0 & 00d0*(press+rho+eps*rho)*w))/alp dfxdy(5,5) = (betax*dpdeps-(dpdeps+rho)*(betax-alp*velx)*w**2)/alp df0dy(1,1) = w df0dy(1,2) = rho*(gxx*velx+gxy*vely+gxz*velz)*w**3 df0dy(1,3) = rho*(gxy*velx+gyy*vely+gyz*velz)*w**3 df0dy(1,4) = rho*(gxz*velx+gyz*vely+gzz*velz)*w**3 df0dy(1,5) = 0.0d0 df0dy(2,1) = (1.000d0+dpdrho+eps)* & (gxx*velx+gxy*vely+gxz*velz)*w**2 df0dy(2,2) = (press+rho+eps*rho)*w* & (gxx*w+2.000d0*(gxx*velx+gxy*vely+gxz*ve & lz)**2*w**3) df0dy(2,3) = (press+rho+eps*rho)*w* & (gxy*w+2.000d0*(gxx*velx+gxy*vely+gxz*ve & lz)*(gxy*velx+gyy*vely+gyz*velz)*w**3) df0dy(2,4) = (press+rho+eps*rho)*w* & (gxz*w+2.000d0*(gxx*velx+gxy*vely+gxz*ve & lz)*(gxz*velx+gyz*vely+gzz*velz)*w**3) df0dy(2,5) = (dpdeps+rho)*(gxx*velx+gxy*vely+gxz*velz)*w**2 df0dy(3,1) = (1.000d0+dpdrho+eps)* & (gxy*velx+gyy*vely+gyz*velz)*w**2 df0dy(3,2) = (press+rho+eps*rho)*w* & (gxy*w+2.000d0*(gxx*velx+gxy*vely+gxz*ve & lz)*(gxy*velx+gyy*vely+gyz*velz)*w**3) df0dy(3,3) = (press+rho+eps*rho)*w* & (gyy*w+2.000d0*(gxy*velx+gyy*vely+gyz*ve & lz)**2*w**3) df0dy(3,4) = (press+rho+eps*rho)*w* & (gyz*w+2.000d0*(gxy*velx+gyy*vely+gyz*ve & lz)*(gxz*velx+gyz*vely+gzz*velz)*w**3) df0dy(3,5) = (dpdeps+rho)*(gxy*velx+gyy*vely+gyz*velz)*w**2 df0dy(4,1) = (1.000d0+dpdrho+eps)* & (gxz*velx+gyz*vely+gzz*velz)*w**2 df0dy(4,2) = (press+rho+eps*rho)*w* & (gxz*w+2.000d0*(gxx*velx+gxy*vely+gxz*ve & lz)*(gxz*velx+gyz*vely+gzz*velz)*w**3) df0dy(4,3) = (press+rho+eps*rho)*w* & (gyz*w+2.000d0*(gxy*velx+gyy*vely+gyz*ve & lz)*(gxz*velx+gyz*vely+gzz*velz)*w**3) df0dy(4,4) = (press+rho+eps*rho)*w* & (gzz*w+2.000d0*(gxz*velx+gyz*vely+gzz*ve & lz)**2*w**3) df0dy(4,5) = (dpdeps+rho)*(gxz*velx+gyz*vely+gzz*velz)*w**2 df0dy(5,1) = w*(-1.000d0+w+eps*w)+dpdrho*(-1.000d0+w**2) df0dy(5,2) = (gxx*velx+gxy*vely+gxz*velz)*w**3* & (-rho+2.000d0*(press+rho+eps & *rho)*w) df0dy(5,3) = (gxy*velx+gyy*vely+gyz*velz)*w**3* & (-rho+2.000d0*(press+rho+eps & *rho)*w) df0dy(5,4) = (gxz*velx+gyz*vely+gzz*velz)*w**3* & (-rho+2.000d0*(press+rho+eps & *rho)*w) df0dy(5,5) = rho*w**2+dpdeps*(-1.000d0+w**2) c invert df0dy do i=1,5 do j=1,5 aa(i,j) = df0dy(i,j) enddo enddo call mahc_svdcmp(aa,5,5,5,5,ww,vv) do i=1,5 do j=1,5 dfxdf0(i,j) = 0.0d0 do k = 1,5 do l = 1,5 dfxdf0(i,j) = dfxdf0(i,j) + . dfxdy(i,k) * vv(k,l) * aa(j,l) / ww(l) enddo enddo enddo enddo maxtestval = 0.0 do i=1,5 do j=1,5 test(i,j) = 0.0d0 do k=1,5 test(i,j) = test(i,j) + dfxdf0(i,k)*p(k,j) enddo if(abs(test(i,j)) .gt. maxtestval) then maxtestval = abs(test(i,j)) endif enddo enddo maxtestval = maxtestval + 1.0e-25 do j=1,5 write(*,*) 'eigenvector ',j do i=1,5 write(*,*) ' ',test(i,j),p(i,j)*lam(j) test(i,j) = test(i,j) - p(i,j) * lam(j) enddo write(*,*) ' ' enddo if(maxval(abs(test))/maxtestval .gt. 1.0d-9) then write(*,*) 'vis_contrib: eigenproblem problem',maxval(abs(test)) call CCTK_WARN(0,"vis_contrib") endif c now, compute inverse of P: do i=1,5 do j=1,5 aa(i,j) = p(i,j) enddo enddo call mahc_svdcmp(aa,5,5,5,5,ww,vv) c calculate the test viscosity components du(1) = dens_r - dens_l du(2) = sx_r - sx_l du(3) = sy_r - sy_l du(4) = sz_r - sz_l du(5) = tau_r - tau_l do i=1,5 vis_test(i) = 0.0d0 if(abs(ww(i)) .le. 1.0d-12) then write(*,*) 'vis_contrib: just thought ya might want to know: ' write(*,*) ' p matrix almost singular! ' call CCTK_WARN(0,"vis_contrib") endif do j=1,5 do k=1,5 do l=1,5 vis_test(i) = vis_test(i) - 0.5d0 * p(i,j) * . abs(lam(j)) * vv(j,l) * aa(k,l) / ww(l) * du(k) enddo enddo enddo write(*,*) 'debug art vis: ',vis_test(i) enddo #endif c calculate change in u: du(1) = dens_r - dens_l du(2) = sx_r - sx_l du(3) = sy_r - sy_l du(4) = sz_r - sz_l du(5) = tau_r - tau_l c solve for characteristic variable change, dw do i=1,5 dw(i) = du(i) do j=1,5 aa(i,j) = p(i,j) enddo enddo #ifdef FAST_MATRIX_INVERT c switch columns to be in the order (1 4 2 3 5) do ii=1,5 paug(ii,1) = p(ii,1) paug(ii,2) = p(ii,4) paug(ii,3) = p(ii,2) paug(ii,4) = p(ii,3) paug(ii,5) = p(ii,5) enddo paug(1,6) = du(1) paug(2,6) = du(2) paug(3,6) = du(3) paug(4,6) = du(4) paug(5,6) = du(5) c switch equation order if fmi is on if(fmi .eq. 1) then do ii=1,6 temp = paug(1,ii) paug(1,ii) = paug(5,ii) paug(5,ii) = temp enddo endif c get lower left triangle to be all zeros do ii=1,5 c first, make diagonal element 1 temp = paug(ii,ii) do jj=ii,6 paug(ii,jj) = paug(ii,jj)/temp enddo c now, get rid of everything below that diagonal do jj=ii+1,5 temp = - (paug(jj,ii)) do kk=ii,6 paug(jj,kk) = paug(jj,kk) + temp*paug(ii,kk) enddo enddo enddo c back substitute f_du(5) = paug(5,6) do ii=4,1,-1 f_du(ii) = paug(ii,6) do jj=ii+1,5 f_du(ii) = f_du(ii) - paug(ii,jj)*f_du(jj) enddo enddo dw(1) = f_du(1) dw(2) = f_du(3) dw(3) = f_du(4) dw(4) = f_du(2) dw(5) = f_du(5) #else call mahc_ludcmp(aa,5,5,indx,dd) call mahc_lubksb(aa,5,5,indx,dw) call mahc_mprove(p,aa,5,5,indx,du,dw) #endif do i=1,5 vis(i) = 0.0d0 do j=1,5 vis(i) = vis(i) - 0.5d0 * p(i,j) * abs(lam(j)) * dw(j) enddo #ifdef DEBUG_EIGEN write(*,*) 'coded art vis : ',vis(i) #endif enddo #ifdef DEBUG_EIGEN maxvis = 0.5d0 * (maxval(abs(vis)) + maxval(abs(vis_test))) + . 1.0d-29 if(maxval(abs(vis - vis_test))/maxvis .gt. 1.0d-8) then write(*,*) 'vis_contrib: precision problem? ', . maxval(abs(vis - vis_test))/maxvis call CCTK_WARN(0,"vis_contrib") endif #endif dens_vis = vis(1) sx_vis = vis(2) sy_vis = vis(3) sz_vis = vis(4) tau_vis = vis(5) return end