[Commits] [svn:einsteintoolkit] GRHydro/branches/hot_and_MHD_temp_dev/src/ (Rev. 168)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Mon Nov 1 00:49:09 CDT 2010
User: cott
Date: 2010/11/01 12:49 AM
Modified:
/branches/hot_and_MHD_temp_dev/src/
GRHydro_Eigenproblem_Marquina.F90, GRHydro_Marquina.F90
Log:
* add finite-T to Marquina
File Changes:
Directory: /branches/hot_and_MHD_temp_dev/src/
==============================================
File [modified]: GRHydro_Eigenproblem_Marquina.F90
Delta lines: +740 -0
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_Eigenproblem_Marquina.F90 2010-11-01 05:33:26 UTC (rev 167)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_Eigenproblem_Marquina.F90 2010-11-01 05:49:09 UTC (rev 168)
@@ -785,6 +785,746 @@
end subroutine eigenproblem_marquina
+subroutine eigenproblem_marquina_hot(handle,rhor,velxr,velyr,&
+ velzr,epsr,rhol,velxl,velyl,velzl,epsl,temp,ye_l,ye_r,gxx,gxy,gxz,&
+ gyy,gyz,gzz,u,det,alp,beta,densl,sxl,syl,szl,taul,&
+ densr,sxr,syr,szr,taur,flux1,flux2,flux3,flux4,flux5)
+
+ USE GRHydro_Scalars
+
+ implicit none
+
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_REAL rhor,velxr,velyr,velzr,epsr,w_lorentzr
+ CCTK_REAL rhol,velxl,velyl,velzl,epsl,w_lorentzl
+ CCTK_REAL densl,sxl,syl,szl,taul
+ CCTK_REAL densr,sxr,syr,szr,taur
+ CCTK_REAL temp,ye_l,ye_r
+
+ CCTK_REAL lam(5),p(5,5),q(5,5),dw(5)
+ CCTK_REAL rfluxr(5),rfluxl(5)
+ CCTK_REAL gxx,gxy,gxz,gyy,gyz,gzz,u,det
+ CCTK_REAL alp,beta,flux1r,flux2r,flux3r,flux4r,flux5r
+ CCTK_REAL flux1l,flux2l,flux3l,flux4l,flux5l
+ CCTK_REAL flux1,flux2,flux3,flux4,flux5
+
+!!$ LOCAL VARS
+
+ CCTK_REAL du(5),aa(5,5),qdiff1,qdiff2,qdiff3,qdiff4,qdiff5
+ integer i,j,k,l
+ CCTK_REAL paug(5,6),tmp1,tmp2,sump,summ,f_du(5)
+ CCTK_REAL leivec1l(5),leivec2l(5),leivec3l(5),leivecpl(5),leivecml(5)
+ CCTK_REAL leivec1r(5),leivec2r(5),leivec3r(5),leivecpr(5),leivecmr(5)
+ CCTK_REAL reivec1l(5),reivec2l(5),reivec3l(5),reivecpl(5),reivecml(5)
+ CCTK_REAL reivec1r(5),reivec2r(5),reivec3r(5),reivecpr(5),reivecmr(5)
+ CCTK_REAL lam1l,lam2l,lam3l,lamml,lampl,lamm_nobetal,lamp_nobetal
+ CCTK_REAL lam1r,lam2r,lam3r,lammr,lampr,lamm_nobetar,lamp_nobetar
+ CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamp_tmp
+
+ CCTK_REAL cs2l,cs2r,one,two
+ CCTK_REAL vlowxr,vlowyr,vlowzr,v2r,wr
+ CCTK_REAL vlowxl,vlowyl,vlowzl,v2l,wl
+
+ CCTK_REAL lamp_nobeta,lamm_nobeta
+
+ CCTK_REAL pressl,dpdrhol,dpdepsl,enthalpyl,kappal
+ CCTK_REAL pressr,dpdrhor,dpdepsr,enthalpyr,kappar
+ CCTK_REAL axpl,axml,vxpl,vxml,xsil,dltl
+ CCTK_REAL axpr,axmr,vxpr,vxmr,xsir,dltr
+ CCTK_REAL cxx,cxy,cxz,gam,vxa,vxb
+ CCTK_INT handle
+ CCTK_REAL invrhol, betainvalp, invfacl
+ CCTK_REAL invrhor, invfacr
+ CCTK_REAL invulpfac, invulmfac, invurpfac, invurmfac
+
+!!$ Warning, warning. Nasty hack follows
+
+#if !USE_EOS_OMNI
+#ifdef _EOS_BASE_INC_
+#undef _EOS_BASE_INC_
+#endif
+#include "EOS_Base.inc"
+#endif
+
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+ real*8 :: xtemp = 0.0d0
+ real*8 :: xye = 0.0d0
+! end EOS Omni vars
+#endif
+
+ one = 1.0d0
+ two = 2.0d0
+
+!!$ LEFT
+
+!!$ Set required fluid quantities
+ invrhol = one / rhol
+#if USE_EOS_OMNI
+ call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rhol,epsl,temp,ye_l,pressl,keyerr,anyerr)
+ call EOS_Omni_cs2(handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rhol,epsl,temp,ye_l,cs2l,keyerr,anyerr)
+ call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rhol,epsl,temp,ye_l,dpdepsl,keyerr,anyerr)
+#else
+ pressl = EOS_Pressure(handle,rhol,epsl)
+ dpdrhol = EOS_DPressByDRho(handle,rhol,epsl)
+ dpdepsl = EOS_DPressByDEps(handle,rhol,epsl)
+ cs2l = (dpdrhol + pressl * dpdepsl * invrhol**2)/ &
+ (1.0d0 + epsl + pressl*invrhol)
+#endif
+
+ enthalpyl = one + epsl + pressl * invrhol
+
+ vlowxl = gxx*velxl + gxy*velyl + gxz*velzl
+ vlowyl = gxy*velxl + gyy*velyl + gyz*velzl
+ vlowzl = gxz*velxl + gyz*velyl + gzz*velzl
+ v2l = vlowxl*velxl + vlowyl*velyl + vlowzl*velzl
+
+ if (v2l .ge. 1.d0-MARQUINA_SMALL) then
+ v2l = 1.d0-MARQUINA_SMALL
+ endif
+!!$ Assume consistent primitive data
+
+ wl = one / sqrt(one - v2l)
+
+!!$ EIGENVALUES
+
+ betainvalp = beta / alp
+ lam1l = velxl - betainvalp
+ lam2l = lam1l
+ lam3l = lam1l
+ lamp_tmp = cs2l*(one-v2l)*(u*(one-v2l*cs2l) - velxl**2*(one-cs2l))
+ if (lamp_tmp .le. MARQUINA_SMALL) then
+ lamp_tmp = MARQUINA_SMALL;
+ endif
+ invfacl = one/(one-v2l*cs2l)
+ lamp_nobetal = (velxl*(one-cs2l) + sqrt(lamp_tmp))*invfacl
+ lamm_nobetal = (velxl*(one-cs2l) - sqrt(lamp_tmp))*invfacl
+
+ lampl = lamp_nobetal - betainvalp
+ lamml = lamm_nobetal - betainvalp
+
+!!$ RIGHT
+
+!!$ Set required fluid quantities
+ invrhor = one / rhor
+#if USE_EOS_OMNI
+ call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rhor,epsr,temp,ye_r,pressr,keyerr,anyerr)
+ call EOS_Omni_cs2(handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rhor,epsr,temp,ye_r,cs2r,keyerr,anyerr)
+ call EOS_Omni_DPressByDEps(handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rhor,epsr,temp,ye_r,dpdepsr,keyerr,anyerr)
+#else
+ pressr = EOS_Pressure(handle,rhor,epsr)
+ dpdrhor = EOS_DPressByDRho(handle,rhor,epsr)
+ dpdepsr = EOS_DPressByDEps(handle,rhor,epsr)
+ cs2r = (dpdrhor + pressr * dpdepsr * invrhor**2)/ &
+ (1.0d0 + epsr + pressr*invrhor)
+#endif
+ enthalpyr = one + epsr + pressr * invrhor
+
+ vlowxr = gxx*velxr + gxy*velyr + gxz*velzr
+ vlowyr = gxy*velxr + gyy*velyr + gyz*velzr
+ vlowzr = gxz*velxr + gyz*velyr + gzz*velzr
+ v2r = vlowxr*velxr + vlowyr*velyr + vlowzr*velzr
+
+ if (v2r .ge. 1.d0-MARQUINA_SMALL) then
+ v2r = 1.d0-MARQUINA_SMALL
+ endif
+!!$ Assume consistent primitive data
+ wr = one / sqrt(one - v2r)
+
+!!$ EIGENVALUES
+
+ lam1r = velxr - betainvalp
+ lam2r = lam1r
+ lam3r = lam1r
+ lamp_tmp = cs2r*(one-v2r)*(u*(one-v2r*cs2r) - velxr**2*(one-cs2r))
+ if (lamp_tmp .le. MARQUINA_SMALL) then
+ lamp_tmp = MARQUINA_SMALL;
+ endif
+ invfacr = one/(one-v2r*cs2r)
+ lamp_nobetar = (velxr*(one-cs2r) + sqrt(lamp_tmp))*invfacr
+ lamm_nobetar = (velxr*(one-cs2r) - sqrt(lamp_tmp))*invfacr
+
+ lampr = lamp_nobetar - betainvalp
+ lammr = lamm_nobetar - betainvalp
+
+!!$ FINAL
+
+ lam1 = dmax1(dabs(lam1l),dabs(lam1r))
+ lam2 = lam1
+ lam3 = lam1
+ lamp = dmax1(dabs(lampl),dabs(lampr))
+ lamm = dmax1(dabs(lamml),dabs(lammr))
+
+!!$ lam(1) = lamm
+!!$ lam(2) = lam1
+!!$ lam(3) = lam2
+!!$ lam(4) = lam3
+!!$ lam(5) = lamp
+
+ lam(1) = lamm
+ lam(5) = lam1
+ lam(3) = lam2
+ lam(4) = lam3
+ lam(2) = lamp
+
+!!$ LEFT
+
+!!$ Compute some auxiliary quantities
+
+ invulpfac = one / (u - velxl*lamp_nobetal)
+ invulmfac = one / (u - velxl*lamm_nobetal)
+ axpl = (u - velxl*velxl)*invulpfac
+ axml = (u - velxl*velxl)*invulmfac
+ vxpl = (velxl - lamp_nobetal)*invulpfac
+ vxml = (velxl - lamm_nobetal)*invulmfac
+
+!!$ Calculate associated right eigenvectors
+
+ kappal = dpdepsl / (dpdepsl - rhol * cs2l)
+
+
+!!$ Right eigenvector # 1
+
+ reivec1l(1) = kappal / (enthalpyl * wl)
+ reivec1l(2) = vlowxl
+ reivec1l(3) = vlowyl
+ reivec1l(4) = vlowzl
+ reivec1l(5) = one - reivec1l(1)
+
+!!$ Right eigenvector # 2
+
+ reivec2l(1) = wl * vlowyl
+ reivec2l(2) = enthalpyl * (gxy + two * wl * wl * vlowxl * vlowyl)
+ reivec2l(3) = enthalpyl * (gyy + two * wl * wl * vlowyl * vlowyl)
+ reivec2l(4) = enthalpyl * (gyz + two * wl * wl * vlowyl * vlowzl)
+ reivec2l(5) = vlowyl * wl * (two * wl * enthalpyl - one)
+
+!!$ Right eigenvector # 3
+
+ reivec3l(1) = wl * vlowzl
+ reivec3l(2) = enthalpyl * (gxz + two * wl * wl * vlowxl * vlowzl)
+ reivec3l(3) = enthalpyl * (gyz + two * wl * wl * vlowyl * vlowzl)
+ reivec3l(4) = enthalpyl * (gzz + two * wl * wl * vlowzl * vlowzl)
+ reivec3l(5) = vlowzl * wl * (two * wl * enthalpyl - one)
+
+!!$ Right + eigenvector
+
+ reivecpl(1) = one
+ reivecpl(2) = enthalpyl * wl * (vlowxl - vxpl)
+ reivecpl(3) = enthalpyl * wl * vlowyl
+ reivecpl(4) = enthalpyl * wl * vlowzl
+ reivecpl(5) = enthalpyl * wl * axpl - one
+
+!!$ Right - eigenvector
+
+ reivecml(1) = one
+ reivecml(2) = enthalpyl * wl * (vlowxl - vxml)
+ reivecml(3) = enthalpyl * wl * vlowyl
+ reivecml(4) = enthalpyl * wl * vlowzl
+ reivecml(5) = enthalpyl * wl * axml - one
+
+!!$ RIGHT
+
+!!$ Compute some auxiliary quantities
+
+ invurpfac = one / (u - velxr*lamp_nobetar)
+ invurmfac = one / (u - velxr*lamm_nobetar)
+ axpr = (u - velxr*velxr)*invurpfac
+ axmr = (u - velxr*velxr)*invurmfac
+ vxpr = (velxr - lamp_nobetar)*invurpfac
+ vxmr = (velxr - lamm_nobetar)*invurmfac
+
+!!$ Calculate associated right eigenvectors
+
+ kappar = dpdepsr / (dpdepsr - rhor * cs2r)
+
+!!$ Right eigenvector # 1
+
+ reivec1r(1) = kappar / (enthalpyr * wr)
+ reivec1r(2) = vlowxr
+ reivec1r(3) = vlowyr
+ reivec1r(4) = vlowzr
+ reivec1r(5) = one - reivec1r(1)
+
+!!$ Right eigenvector # 2
+
+ reivec2r(1) = wr * vlowyr
+ reivec2r(2) = enthalpyr * (gxy + two * wr * wr * vlowxr * vlowyr)
+ reivec2r(3) = enthalpyr * (gyy + two * wr * wr * vlowyr * vlowyr)
+ reivec2r(4) = enthalpyr * (gyz + two * wr * wr * vlowyr * vlowzr)
+ reivec2r(5) = vlowyr * wr * (two * wr * enthalpyr - one)
+
+!!$ Right eigenvector # 3
+
+ reivec3r(1) = wr * vlowzr
+ reivec3r(2) = enthalpyr * (gxz + two * wr * wr * vlowxr * vlowzr)
+ reivec3r(3) = enthalpyr * (gyz + two * wr * wr * vlowyr * vlowzr)
+ reivec3r(4) = enthalpyr * (gzz + two * wr * wr * vlowzr * vlowzr)
+ reivec3r(5) = vlowzr * wr * (two * wr * enthalpyr - one)
+
+!!$ Right + eigenvector
+
+ reivecpr(1) = one
+ reivecpr(2) = enthalpyr * wr * (vlowxr - vxpr)
+ reivecpr(3) = enthalpyr * wr * vlowyr
+ reivecpr(4) = enthalpyr * wr * vlowzr
+ reivecpr(5) = enthalpyr * wr * axpr - one
+
+!!$ Right - eigenvector
+
+ reivecmr(1) = one
+ reivecmr(2) = enthalpyr * wr * (vlowxr - vxmr)
+ reivecmr(3) = enthalpyr * wr * vlowyr
+ reivecmr(4) = enthalpyr * wr * vlowzr
+ reivecmr(5) = enthalpyr * wr * axmr - one
+
+!!$ Calculate associated left eigenvectors if requested
+
+ if (ANALYTICAL) then
+
+ cxx = gyy * gzz - gyz * gyz
+ cxy = gxz * gyz - gxy * gzz
+ cxz = gxy * gyz - gxz * gyy
+ gam = gxx * cxx + gxy * cxy + gxz * cxz
+
+!!$ LEFT
+
+ xsil = cxx - gam * velxl * velxl
+ dltl = enthalpyl**3 * wl * (kappal - one) * (vxml - vxpl) * xsil
+
+!!$ Left eigenvector # 1
+
+ tmp1 = wl / (kappal - one)
+
+ leivec1l(1) = tmp1 * (enthalpyl - wl)
+ leivec1l(2) = tmp1 * wl * velxl
+ leivec1l(3) = tmp1 * wl * velyl
+ leivec1l(4) = tmp1 * wl * velzl
+ leivec1l(5) =-tmp1 * wl
+
+!!$ Left eigenvector # 2
+
+ tmp1 = one / (xsil * enthalpyl)
+
+ leivec2l(1) = (gyz * vlowzl - gzz * vlowyl) * tmp1
+ leivec2l(2) = (gzz * vlowyl - gyz * vlowzl) * tmp1 * velxl
+ leivec2l(3) = (gzz * (one - velxl * vlowxl) + gxz * vlowzl * velxl ) * tmp1
+ leivec2l(4) = (gyz * (velxl * vlowxl - one) - gxz * velxl * vlowyl) * tmp1
+ leivec2l(5) = (gyz * vlowzl - gzz * vlowyl) * tmp1
+
+!!$ Left eigenvector # 3
+
+ leivec3l(1) = (gyz * vlowyl - gyy * vlowzl) * tmp1
+ leivec3l(2) = (gyy * vlowzl - gyz * vlowyl) * tmp1 * velxl
+ leivec3l(3) = (gyz * (velxl * vlowxl - one) - gxy * velxl * vlowzl) * tmp1
+ leivec3l(4) = (gyy * (one - velxl * vlowxl) + gxy * velxl * vlowyl) * tmp1
+ leivec3l(5) = (gyz * vlowyl - gyy * vlowzl) * tmp1
+
+!!$ Left + eigenvector
+
+ tmp1 = enthalpyl * enthalpyl / dltl
+ tmp2 = wl * wl * xsil
+
+ leivecpl(1) = - (enthalpyl * wl * vxml * xsil + (one - kappal) * &
+ (vxml * (tmp2 - cxx) - gam * velxl) - kappal * tmp2 * vxml) * tmp1
+ leivecpl(2) = - (cxx * (one - kappal * axml) + (two * kappal - one) * &
+ vxml * (tmp2 * velxl - cxx * velxl)) * tmp1
+ leivecpl(3) = - (cxy * (one - kappal * axml) + (two * kappal - one) * &
+ vxml * (tmp2 * velyl - cxy * velxl)) * tmp1
+ leivecpl(4) = - (cxz * (one - kappal * axml) + (two * kappal - one) * &
+ vxml * (tmp2 * velzl - cxz * velxl)) * tmp1
+ leivecpl(5) = - ((one - kappal) * (vxml * (tmp2 - cxx) - gam * velxl) - &
+ kappal * tmp2 * vxml) * tmp1
+
+!!$ Left - eigenvector
+
+ leivecml(1) = (enthalpyl * wl * vxpl * xsil + (one - kappal) * &
+ (vxpl * (tmp2 - cxx) - gam * velxl) - kappal * tmp2 * vxpl) * tmp1
+ leivecml(2) = (cxx * (one - kappal * axpl) + (two * kappal - one) * &
+ vxpl * (tmp2 * velxl - cxx * velxl)) * tmp1
+ leivecml(3) = (cxy * (one - kappal * axpl) + (two * kappal - one) * &
+ vxpl * (tmp2 * velyl - cxy * velxl)) * tmp1
+ leivecml(4) = (cxz * (one - kappal * axpl) + (two * kappal - one) * &
+ vxpl * (tmp2 * velzl - cxz * velxl)) * tmp1
+ leivecml(5) = ((one - kappal) * (vxpl * (tmp2 - cxx) - gam * &
+ velxl) - kappal * tmp2 * vxpl) * tmp1
+
+!!$ RIGHT
+
+ xsir = cxx - gam * velxr * velxr
+ dltr = enthalpyr**3 * wr * (kappar - one) * (vxmr - vxpr) * xsir
+
+!!$ Left eigenvector # 1
+
+ tmp1 = wr / (kappar - one)
+
+ leivec1r(1) = tmp1 * (enthalpyr - wr)
+ leivec1r(2) = tmp1 * wr * velxr
+ leivec1r(3) = tmp1 * wr * velyr
+ leivec1r(4) = tmp1 * wr * velzr
+ leivec1r(5) =-tmp1 * wr
+
+!!$ Left eigenvector # 2
+
+ tmp1 = one / (xsir * enthalpyr)
+
+ leivec2r(1) = (gyz * vlowzr - gzz * vlowyr) * tmp1
+ leivec2r(2) = (gzz * vlowyr - gyz * vlowzr) * tmp1 * velxr
+ leivec2r(3) = (gzz * (one - velxr * vlowxr) + gxz * vlowzr * velxr) * tmp1
+ leivec2r(4) = (gyz * (velxr * vlowxr - one) - gxz * velxr * vlowyr) * tmp1
+ leivec2r(5) = (gyz * vlowzr - gzz * vlowyr) * tmp1
+
+!!$ Left eigenvector # 3
+
+ leivec3r(1) = (gyz * vlowyr - gyy * vlowzr) * tmp1
+ leivec3r(2) = (gyy * vlowzr - gyz * vlowyr) * tmp1 * velxr
+ leivec3r(3) = (gyz * (velxr * vlowxr - one) - gxy * velxr * vlowzr) * tmp1
+ leivec3r(4) = (gyy * (one - velxr * vlowxr) + gxy * velxr * vlowyr) * tmp1
+ leivec3r(5) = (gyz * vlowyr - gyy * vlowzr) * tmp1
+
+!!$ Left + eigenvector
+
+ tmp1 = enthalpyr * enthalpyr / dltr
+ tmp2 = wr * wr * xsir
+
+ leivecpr(1) = - (enthalpyr * wr * vxmr * xsir + (one - kappar) * &
+ (vxmr * (tmp2 - cxx) - gam * velxr) - kappar * tmp2 * vxmr) * tmp1
+ leivecpr(2) = - (cxx * (one - kappar * axmr) + (two * kappar - one) * &
+ vxmr * (tmp2 * velxr - cxx * velxr)) * tmp1
+ leivecpr(3) = - (cxy * (one - kappar * axmr) + (two * kappar - one) * &
+ vxmr * (tmp2 * velyr - cxy * velxr)) * tmp1
+ leivecpr(4) = - (cxz * (one - kappar * axmr) + (two * kappar - one) * &
+ vxmr * (tmp2 * velzr - cxz * velxr)) * tmp1
+ leivecpr(5) = - ((one - kappar) * (vxmr * (tmp2 - cxx) - gam * velxr) - &
+ kappar * tmp2 * vxmr) * tmp1
+
+!!$ Left - eigenvector
+
+ leivecmr(1) = (enthalpyr * wr * vxpr * xsir + (one - kappar) * &
+ (vxpr * (tmp2 - cxx) - gam * velxr) - kappar * tmp2 * vxpr) * tmp1
+ leivecmr(2) = (cxx * (one - kappar * axpr) + (two * kappar - one) * &
+ vxpr * (tmp2 * velxr - cxx * velxr)) * tmp1
+ leivecmr(3) = (cxy * (one - kappar * axpr) + (two * kappar - one) * &
+ vxpr * (tmp2 * velyr - cxy * velxr)) * tmp1
+ leivecmr(4) = (cxz * (one - kappar * axpr) + (two * kappar - one) * &
+ vxpr * (tmp2 * velzr - cxz * velxr)) * tmp1
+ leivecmr(5) = ((one - kappar) * (vxpr * (tmp2 - cxx) - gam * &
+ velxr) - kappar * tmp2 * vxpr) * tmp1
+
+ endif
+
+!!$ LEFT
+!!$ PUT RIGHT EIGENVECTORS IN THE P MATRIX
+
+!!$ p(:,1) = reivecml(:)
+!!$ p(:,2) = reivec1l(:)
+!!$ p(:,3) = reivec2l(:)
+!!$ p(:,4) = reivec3l(:)
+!!$ p(:,5) = reivecpl(:)
+
+ p(:,1) = reivecml(:)
+ p(:,2) = reivecpl(:)
+ p(:,3) = reivec2l(:)
+ p(:,4) = reivec3l(:)
+ p(:,5) = reivec1l(:)
+
+!!$ Calculate change in u:
+
+ du(1) = densl
+ du(2) = sxl
+ du(3) = syl
+ du(4) = szl
+ du(5) = taul
+
+ if (ANALYTICAL) then
+
+ if (FAST) then
+
+ sump = 0.0d0
+ summ = 0.0d0
+
+ do i=1,5
+ sump = sump + (lamp - lam1) * leivecpl(i) * du(i)
+ summ = summ + (lamm - lam1) * leivecml(i) * du(i)
+ enddo
+
+ vxa = sump + summ
+ vxb =-(sump * vxpl + summ * vxml)
+
+ rfluxl(1) = lam1 * du(1) + vxa
+ rfluxl(2) = lam1 * du(2) + enthalpyl * wl * (vlowxl * vxa + vxb)
+ rfluxl(3) = lam1 * du(3) + enthalpyl * wl * (vlowyl * vxa)
+ rfluxl(4) = lam1 * du(4) + enthalpyl * wl * (vlowzl * vxa)
+ rfluxl(5) = lam1 * du(5) + enthalpyl * wl * (velxl * vxb + vxa) - vxa
+
+ else
+
+!!$ PUT LEFT EIGENVECTORS IN THE Q MATRIX
+
+ q(1,:) = leivecml(:)
+ q(2,:) = leivecpl(:)
+ q(3,:) = leivec2l(:)
+ q(4,:) = leivec3l(:)
+ q(5,:) = leivec1l(:)
+
+ do i=1,5
+ dw(i) = 0.0d0
+ do j=1,5
+ dw(i) = dw(i) + q(i,j) * du(j)
+ enddo
+ enddo
+
+ do i = 1, 5
+ rfluxl(i) = 0.d0
+ do j = 1, 5
+ rfluxl(i) = rfluxl(i) + p(i,j) * lam(j) * dw(j)
+ end do
+ end do
+
+ endif
+
+ else
+
+!!$ Solve for characteristic variable change, dw
+
+ dw=du
+ aa = p
+
+ do i=1,5
+ paug(:,i) = p(:,i)
+ enddo
+
+!same, but in old F77 style!
+!!$ do i=1,5
+!!$ dw(i) = du(i)
+!!$ do j=1,5
+!!$ aa(i,j) = p(i,j)
+!!$ end do
+!!$ enddo
+!!$
+!!$ do i=1,5
+!!$ paug(i,1) = p(i,1)
+!!$ paug(i,2) = p(i,2)
+!!$ paug(i,3) = p(i,3)
+!!$ paug(i,4) = p(i,4)
+!!$ paug(i,5) = p(i,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)
+
+!!$ Get lower left triangle to be all zeros
+ do i=1,5
+!!$ First, make diagonal element 1
+ tmp1 = one / paug(i,i)
+ do j=i,6
+ paug(i,j) = paug(i,j)*tmp1
+ enddo
+!!$ Now, get rid of everything below that diagonal
+ do j=i+1,5
+ tmp1 = - (paug(j,i))
+ do k=i,6
+ paug(j,k) = paug(j,k) + tmp1*paug(i,k)
+ enddo
+ enddo
+ enddo
+!!$ Back substitute
+ f_du(5) = paug(5,6)
+ do i=4,1,-1
+ f_du(i) = paug(i,6)
+ do j=i+1,5
+ f_du(i) = f_du(i) - paug(i,j)*f_du(j)
+ enddo
+ enddo
+
+ dw(1) = f_du(1)
+ dw(2) = f_du(2)
+ dw(3) = f_du(3)
+ dw(4) = f_du(4)
+ dw(5) = f_du(5)
+
+ do i = 1, 5
+ rfluxl(i) = 0.d0
+ do j = 1, 5
+ rfluxl(i) = rfluxl(i) + p(i,j) * lam(j) * dw(j)
+ end do
+ end do
+
+ endif
+
+ flux1l = rfluxl(1)
+ flux2l = rfluxl(2)
+ flux3l = rfluxl(3)
+ flux4l = rfluxl(4)
+ flux5l = rfluxl(5)
+
+!!$ RIGHT
+
+!!$ PUT RIGHT EIGENVECTORS IN THE P MATRIX
+
+!!$ p(:,1) = reivecmr(:)
+!!$ p(:,2) = reivec1r(:)
+!!$ p(:,3) = reivec2r(:)
+!!$ p(:,4) = reivec3r(:)
+!!$ p(:,5) = reivecpr(:)
+
+ p(:,1) = reivecmr(:)
+ p(:,2) = reivecpr(:)
+ p(:,3) = reivec2r(:)
+ p(:,4) = reivec3r(:)
+ p(:,5) = reivec1r(:)
+
+ du(1) = densr
+ du(2) = sxr
+ du(3) = syr
+ du(4) = szr
+ du(5) = taur
+
+ if (ANALYTICAL) then
+
+ if (FAST) then
+
+ sump = 0.0d0
+ summ = 0.0d0
+
+ do i=1,5
+ sump = sump + (lamp - lam1) * leivecpr(i) * du(i)
+ summ = summ + (lamm - lam1) * leivecmr(i) * du(i)
+ enddo
+
+ vxa = sump + summ
+ vxb =-(sump * vxpr + summ * vxmr)
+
+ rfluxr(1) = lam1 * du(1) + vxa
+ rfluxr(2) = lam1 * du(2) + enthalpyr * wr * (vlowxr * vxa + vxb)
+ rfluxr(3) = lam1 * du(3) + enthalpyr * wr * (vlowyr * vxa)
+ rfluxr(4) = lam1 * du(4) + enthalpyr * wr * (vlowzr * vxa)
+ rfluxr(5) = lam1 * du(5) + enthalpyr * wr * (velxr * vxb + vxa) - vxa
+
+ else
+
+!!$ PUT LEFT EIGENVECTORS IN THE Q MATRIX
+
+ q(1,:) = leivecmr(:)
+ q(2,:) = leivecpr(:)
+ q(3,:) = leivec2r(:)
+ q(4,:) = leivec3r(:)
+ q(5,:) = leivec1r(:)
+
+ do i=1,5
+ dw(i) = 0.0d0
+ do j=1,5
+ dw(i) = dw(i) + q(i,j) * du(j)
+ enddo
+ enddo
+
+ do i = 1, 5
+ rfluxr(i) = 0.d0
+ do j = 1, 5
+ rfluxr(i) = rfluxr(i) + p(i,j) * lam(j) * dw(j)
+ end do
+ end do
+
+ endif
+
+ else
+
+!!$ Solve for characteristic variable change, dw
+
+ do i=1,5
+ dw(i) = du(i)
+ do j=1,5
+ aa(i,j) = p(i,j)
+ end do
+ enddo
+
+ do i=1,5
+ paug(i,1) = p(i,1)
+ paug(i,2) = p(i,2)
+ paug(i,3) = p(i,3)
+ paug(i,4) = p(i,4)
+ paug(i,5) = p(i,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)
+
+!!$ Get lower left triangle to be all zeros
+
+ do i=1,5
+!!$ First, make diagonal element 1
+ tmp1 = one / paug(i,i)
+ do j=i,6
+ paug(i,j) = paug(i,j)*tmp1
+ enddo
+!!$ Now, get rid of everything below that diagonal
+ do j=i+1,5
+ tmp1 = - (paug(j,i))
+ do k=i,6
+ paug(j,k) = paug(j,k) + tmp1*paug(i,k)
+ enddo
+ enddo
+ enddo
+!!$ Back substitute
+
+ f_du(5) = paug(5,6)
+ do i=4,1,-1
+ f_du(i) = paug(i,6)
+ do j=i+1,5
+ f_du(i) = f_du(i) - paug(i,j)*f_du(j)
+ enddo
+ enddo
+
+ dw(1) = f_du(1)
+ dw(2) = f_du(2)
+ dw(3) = f_du(3)
+ dw(4) = f_du(4)
+ dw(5) = f_du(5)
+
+ do i = 1, 5
+ rfluxr(i) = 0.d0
+ do j = 1, 5
+ rfluxr(i) = rfluxr(i) + p(i,j) * lam(j) * dw(j)
+ end do
+ end do
+
+ endif
+
+ flux1r = rfluxr(1)
+ flux2r = rfluxr(2)
+ flux3r = rfluxr(3)
+ flux4r = rfluxr(4)
+ flux5r = rfluxr(5)
+
+ flux1 = flux1r - flux1l
+ flux2 = flux2r - flux2l
+ flux3 = flux3r - flux3l
+ flux4 = flux4r - flux4l
+ flux5 = flux5r - flux5l
+
+ return
+
+end subroutine eigenproblem_marquina_hot
+
/*@@
@routine eigenproblem_marquina
@date Wed Feb 13 12:27:59 2002
File [modified]: GRHydro_Marquina.F90
Delta lines: +80 -33
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_Marquina.F90 2010-11-01 05:33:26 UTC (rev 167)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_Marquina.F90 2010-11-01 05:49:09 UTC (rev 168)
@@ -47,6 +47,7 @@
avg_det,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,&
pressp,pressm_i, &
tmp_w_lorentzp, tmp_w_lorentzm_i, w_lorentzp,w_lorentzm_i, usendh, psi4h
+ CCTK_REAL :: xtemp
integer :: m
integer :: i,j,k
@@ -212,44 +213,90 @@
!!$eigenvalues and right eigenvectors
if (flux_direction == 1) then
+
+ if(evolve_temp.eq.0) then
+ call eigenproblem_marquina(GRHydro_eos_handle,&
+ primm_i(1),primm_i(2), &
+ primm_i(3),primm_i(4),primm_i(5),primp(1), &
+ primp(2),primp(3),primp(4),primp(5), &
+ gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
+ usendh,avg_det,avg_alp,avg_beta,consp(1),consp(2),&
+ consp(3), consp(4), consp(5),consm_i(1),consm_i(2), &
+ consm_i(3),consm_i(4),consm_i(5),marquinaflux(1), &
+ marquinaflux(2),marquinaflux(3),marquinaflux(4), &
+ marquinaflux(5))
+ else
+ xtemp = temperature(i,j,k)
+ call eigenproblem_marquina_hot(GRHydro_eos_handle,&
+ primm_i(1),primm_i(2), &
+ primm_i(3),primm_i(4),primm_i(5),primp(1), &
+ primp(2),primp(3),primp(4),primp(5), &
+ xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),y_e_plus(i,j,k),&
+ gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
+ usendh,avg_det,avg_alp,avg_beta,consp(1),consp(2),&
+ consp(3), consp(4), consp(5),consm_i(1),consm_i(2), &
+ consm_i(3),consm_i(4),consm_i(5),marquinaflux(1), &
+ marquinaflux(2),marquinaflux(3),marquinaflux(4), &
+ marquinaflux(5))
+
+ endif
- call eigenproblem_marquina(GRHydro_eos_handle,&
- primm_i(1),primm_i(2), &
- primm_i(3),primm_i(4),primm_i(5),primp(1), &
- primp(2),primp(3),primp(4),primp(5), &
- gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
- usendh,avg_det,avg_alp,avg_beta,consp(1),consp(2),&
- consp(3), consp(4), consp(5),consm_i(1),consm_i(2), &
- consm_i(3),consm_i(4),consm_i(5),marquinaflux(1), &
- marquinaflux(2),marquinaflux(3),marquinaflux(4), &
- marquinaflux(5))
-
else if (flux_direction == 2) then
+
+ if(evolve_temp.eq.0) then
+ call eigenproblem_marquina(GRHydro_eos_handle,&
+ primm_i(1),primm_i(3), &
+ primm_i(4),primm_i(2),primm_i(5),primp(1), &
+ primp(3),primp(4),primp(2),primp(5), &
+ gyyh,gyzh,gxyh,gzzh,gxzh,gxxh, &
+ usendh,avg_det,avg_alp,avg_beta,consp(1),consp(3),&
+ consp(4), consp(2), consp(5),consm_i(1),consm_i(3), &
+ consm_i(4),consm_i(2),consm_i(5),marquinaflux(1), &
+ marquinaflux(3),marquinaflux(4),marquinaflux(2), &
+ marquinaflux(5))
+ else
+ xtemp = temperature(i,j,k)
+ call eigenproblem_marquina(GRHydro_eos_handle,&
+ primm_i(1),primm_i(3), &
+ primm_i(4),primm_i(2),primm_i(5),primp(1), &
+ primp(3),primp(4),primp(2),primp(5), &
+ xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),y_e_plus(i,j,k),&
+ gyyh,gyzh,gxyh,gzzh,gxzh,gxxh, &
+ usendh,avg_det,avg_alp,avg_beta,consp(1),consp(3),&
+ consp(4), consp(2), consp(5),consm_i(1),consm_i(3), &
+ consm_i(4),consm_i(2),consm_i(5),marquinaflux(1), &
+ marquinaflux(3),marquinaflux(4),marquinaflux(2), &
+ marquinaflux(5))
+
+ endif
- call eigenproblem_marquina(GRHydro_eos_handle,&
- primm_i(1),primm_i(3), &
- primm_i(4),primm_i(2),primm_i(5),primp(1), &
- primp(3),primp(4),primp(2),primp(5), &
- gyyh,gyzh,gxyh,gzzh,gxzh,gxxh, &
- usendh,avg_det,avg_alp,avg_beta,consp(1),consp(3),&
- consp(4), consp(2), consp(5),consm_i(1),consm_i(3), &
- consm_i(4),consm_i(2),consm_i(5),marquinaflux(1), &
- marquinaflux(3),marquinaflux(4),marquinaflux(2), &
- marquinaflux(5))
-
else
-
- call eigenproblem_marquina(GRHydro_eos_handle,&
- primm_i(1),primm_i(4), &
- primm_i(2),primm_i(3),primm_i(5),primp(1), &
- primp(4),primp(2),primp(3),primp(5), &
- gzzh,gxzh,gyzh,gxxh,gxyh,gyyh, &
- usendh,avg_det,avg_alp,avg_beta,consp(1),consp(4),&
- consp(2), consp(3), consp(5),consm_i(1),consm_i(4), &
- consm_i(2),consm_i(3),consm_i(5),marquinaflux(1), &
- marquinaflux(4),marquinaflux(2),marquinaflux(3), &
- marquinaflux(5))
+ if(evolve_temp.eq.0) then
+ call eigenproblem_marquina(GRHydro_eos_handle,&
+ primm_i(1),primm_i(4), &
+ primm_i(2),primm_i(3),primm_i(5),primp(1), &
+ primp(4),primp(2),primp(3),primp(5), &
+ gzzh,gxzh,gyzh,gxxh,gxyh,gyyh, &
+ usendh,avg_det,avg_alp,avg_beta,consp(1),consp(4),&
+ consp(2), consp(3), consp(5),consm_i(1),consm_i(4), &
+ consm_i(2),consm_i(3),consm_i(5),marquinaflux(1), &
+ marquinaflux(4),marquinaflux(2),marquinaflux(3), &
+ marquinaflux(5))
+ else
+ xtemp = temperature(i,j,k)
+ call eigenproblem_marquina(GRHydro_eos_handle,&
+ primm_i(1),primm_i(4), &
+ primm_i(2),primm_i(3),primm_i(5),primp(1), &
+ primp(4),primp(2),primp(3),primp(5), &
+ xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),y_e_plus(i,j,k),&
+ gzzh,gxzh,gyzh,gxxh,gxyh,gyyh, &
+ usendh,avg_det,avg_alp,avg_beta,consp(1),consp(4),&
+ consp(2), consp(3), consp(5),consm_i(1),consm_i(4), &
+ consm_i(2),consm_i(3),consm_i(5),marquinaflux(1), &
+ marquinaflux(4),marquinaflux(2),marquinaflux(3), &
+ marquinaflux(5))
+ endif
end if
fplus = 0.d0
More information about the Commits
mailing list