Commits November 2010

commits@lists.einsteintoolkit.org

9 participants
258 discussions

[svn:einsteintoolkit] www/release-info/ (Rev. 359)
by knarf＠cct.lsu.edu 01 Nov '10

01 Nov '10

User: knarf Date: 2010/11/01 08:41 AM Modified: /release-info/ parse_testsuite_results.cgi Log: more color File Changes: Directory: /release-info/ ========================= File [modified]: parse_testsuite_results.cgi Delta lines: +4 -4 =================================================================== --- release-info/parse_testsuite_results.cgi 2010-11-01 13:38:02 UTC (rev 358) +++ release-info/parse_testsuite_results.cgi 2010-11-01 13:41:22 UTC (rev 359) @@ -135,16 +135,16 @@ $allhosts{$host}{$nprocs}{"ran"}++; if ($results{$thorn}{$test}{$host}{$nprocs}{"result"} =~ /Success/) { $allhosts{$host}{$nprocs}{"success"}++; - $table_entries .= "<td$color align='center'><font color='00FF00'>â</font></td>"; + $table_entries .= "<td$color align='center'><font color='#00FF00'>â</font></td>"; } else { $allhosts{$host}{$nprocs}{"fail"}++; my $details = $results{$thorn}{$test}{$host}{$nprocs}{"details"}; my $entry = "<b>$details</b>"; if ($details =~ /(\d+) files compared, (\d+) differ, (\d+) differ significantly/) { - $entry = "<font color='FF0000'>$3</font>/$1"; + $entry = "<font color='#FF0000'>$3</font>/$1"; } if ($details =~ /(\d+) files missing, (\d+) files compared, (\d+) differ/) { - $entry = "<font color='FF0000'>$1</font> missing, <font color='FF0000'>$3</font>/$2"; + $entry = "<font color='#FF0000'>$1</font> missing, <font color='#FF0000'>$3</font>/$2"; } $table_entries .= "<td$color align='center'>$entry</td>"; } @@ -221,6 +221,6 @@ $table_entries = $table_entries_h . $table_entries; -print '<style type="text/css"> th, td {padding: 0em;}</style>\n'; +print '<style type="text/css"> th, td {padding: 0em;}</style>'."\n"; print "<table border='1'>$table_entries</table>\n";

1 0

[svn:einsteintoolkit] www/release-info/ (Rev. 358)
by knarf＠cct.lsu.edu 01 Nov '10

01 Nov '10

User: knarf Date: 2010/11/01 08:38 AM Modified: /release-info/ parse_testsuite_results.cgi Log: more color File Changes: Directory: /release-info/ ========================= File [modified]: parse_testsuite_results.cgi Delta lines: +3 -3 =================================================================== --- release-info/parse_testsuite_results.cgi 2010-11-01 13:33:44 UTC (rev 357) +++ release-info/parse_testsuite_results.cgi 2010-11-01 13:38:02 UTC (rev 358) @@ -135,16 +135,16 @@ $allhosts{$host}{$nprocs}{"ran"}++; if ($results{$thorn}{$test}{$host}{$nprocs}{"result"} =~ /Success/) { $allhosts{$host}{$nprocs}{"success"}++; - $table_entries .= "<td$color align='center'>â</td>"; + $table_entries .= "<td$color align='center'><font color='00FF00'>â</font></td>"; } else { $allhosts{$host}{$nprocs}{"fail"}++; my $details = $results{$thorn}{$test}{$host}{$nprocs}{"details"}; my $entry = "<b>$details</b>"; if ($details =~ /(\d+) files compared, (\d+) differ, (\d+) differ significantly/) { - $entry = "<b>$3</b>/$1"; + $entry = "<font color='FF0000'>$3</font>/$1"; } if ($details =~ /(\d+) files missing, (\d+) files compared, (\d+) differ/) { - $entry = "<b>$1</b> missing, <b>$3</b>/$2"; + $entry = "<font color='FF0000'>$1</font> missing, <font color='FF0000'>$3</font>/$2"; } $table_entries .= "<td$color align='center'>$entry</td>"; }

1 0

[svn:einsteintoolkit] www/release-info/ (Rev. 357)
by knarf＠cct.lsu.edu 01 Nov '10

01 Nov '10

User: knarf Date: 2010/11/01 08:33 AM Modified: /release-info/ parse_testsuite_results.cgi Log: no padding File Changes: Directory: /release-info/ ========================= File [modified]: parse_testsuite_results.cgi Delta lines: +2 -0 =================================================================== --- release-info/parse_testsuite_results.cgi 2010-11-01 13:29:35 UTC (rev 356) +++ release-info/parse_testsuite_results.cgi 2010-11-01 13:33:44 UTC (rev 357) @@ -221,4 +221,6 @@ $table_entries = $table_entries_h . $table_entries; +print '<style type="text/css"> th, td {padding: 0em;}</style>\n'; + print "<table border='1'>$table_entries</table>\n";

1 0

[svn:einsteintoolkit] www/release-info/ (Rev. 356)
by knarf＠cct.lsu.edu 01 Nov '10

01 Nov '10

User: knarf Date: 2010/11/01 08:29 AM Modified: /release-info/ numrel__1.log, numrel__2.log Log: update results File Changes: Directory: /release-info/ ========================= File [modified]: numrel__1.log Delta lines: +7 -9 =================================================================== --- release-info/numrel__1.log 2010-10-31 17:58:33 UTC (rev 355) +++ release-info/numrel__1.log 2010-11-01 13:29:35 UTC (rev 356) @@ -868,9 +868,9 @@ Hydro_InitExcision: x_pugh_ppm Success: 72 files compared, 34 differ in the last digits Hydro_InitExcision: x_pugh_tvd - Failure: 72 files missing, 0 files compared, 0 differ + Success: 72 files compared, 18 differ in the last digits IDAnalyticBH: kerr - Failure: 45 files missing, 0 files compared, 0 differ + Success: 45 files compared, 1 differ in the last digits IDAnalyticBH: test_bl Success: 8 files identical IDAnalyticBH: test_misner @@ -880,7 +880,7 @@ IDAxiBrillBH: test_axibrill_nostagger Success: 96 files compared, 69 differ in the last digits IDAxiOddBrillBH: E2xeon_test_axioddbh - Success: 88 files compared, 85 differ in the last digits + Success: 88 files compared, 87 differ in the last digits IDAxiOddBrillBH: test_axioddbh Success: 120 files compared, 72 differ in the last digits IDConstraintViolate: ConstraintViolate @@ -1025,7 +1025,7 @@ Summary for configuration sim - Time -> Sun Oct 31 11:27:04 CDT 2010 + Time -> Sun Oct 31 17:22:11 CDT 2010 Host -> numrel07.cct.lsu.edu User -> knarf @@ -1034,15 +1034,13 @@ Runnable tests -> 132 Total number of thorns -> 157 Number of tested thorns -> 39 - Number of tests passed -> 128 + Number of tests passed -> 130 Number passed only to - set tolerance -> 76 - Number failed -> 4 + set tolerance -> 78 + Number failed -> 2 Tests failed: - x_pugh_tvd (from Hydro_InitExcision) - kerr (from IDAnalyticBH) test_waveell (from IDScalarWaveElliptic) test_binary_1 (from WaveBinarySource) File [modified]: numrel__2.log Delta lines: +3 -3 =================================================================== --- release-info/numrel__2.log 2010-10-31 17:58:33 UTC (rev 355) +++ release-info/numrel__2.log 2010-11-01 13:29:35 UTC (rev 356) @@ -875,7 +875,7 @@ IDAxiBrillBH: test_axibrill_nostagger Success: 96 files compared, 69 differ in the last digits IDAxiOddBrillBH: E2xeon_test_axioddbh - Success: 88 files compared, 87 differ in the last digits + Success: 88 files compared, 85 differ in the last digits IDAxiOddBrillBH: test_axioddbh Success: 120 files compared, 74 differ in the last digits IDConstraintViolate: ConstraintViolate @@ -1016,8 +1016,8 @@ Summary for configuration sim - Time -> Sun Oct 31 11:47:32 CDT 2010 - Host -> numrel07.cct.lsu.edu + Time -> + Host -> User -> knarf Total available tests -> 174

1 0

[svn:einsteintoolkit] GRHydro/branches/hot_and_MHD_temp_dev/src/ (Rev. 169)
by cott＠tapir.caltech.edu 31 Oct '10

31 Oct '10

User: cott Date: 2010/11/01 01:43 AM Modified: /branches/hot_and_MHD_temp_dev/src/ GRHydro_Con2Prim.F90, GRHydro_Eigenproblem.F90, GRHydro_HLLEPoly.F90, GRHydro_ParamCheck.F90, GRHydro_Prim2Con.F90, GRHydro_RiemannSolve.F90, GRHydro_UpdateMask.F90 Log: * bunch of small bug fixes in the 'hot' part of the code File Changes: Directory: /branches/hot_and_MHD_temp_dev/src/ ============================================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +2 -3 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_Con2Prim.F90 2010-11-01 05:49:09 UTC (rev 168) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_Con2Prim.F90 2010-11-01 06:43:05 UTC (rev 169) @@ -580,9 +580,9 @@ #if USE_EOS_OMNI call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,temp,ye,xpress,keyerr,anyerr) - pold = max(1.d-10,xpress) + pold = max(1.d-15,xpress) #else - pold = max(1.d-10,EOS_Pressure(handle, rho, epsilon)) + pold = max(1.d-15,EOS_Pressure(handle, rho, epsilon)) #endif @@ -615,7 +615,6 @@ #if USE_EOS_OMNI call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,temp,ye,xpress,keyerr,anyerr) - f = pold - xpress #else f = pold - EOS_Pressure(handle, rho, epsilon) File [modified]: GRHydro_Eigenproblem.F90 Delta lines: +73 -0 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_Eigenproblem.F90 2010-11-01 05:49:09 UTC (rev 168) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_Eigenproblem.F90 2010-11-01 06:43:05 UTC (rev 169) @@ -136,6 +136,79 @@ end subroutine eigenvalues + +subroutine eigenvalues_hot(handle,rho,velx,vely,velz,eps, & + temp,ye,w_lorentz,lam,gxx,gxy,gxz,gyy,gyz,gzz,u,alp,beta) + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL rho,velx,vely,velz,eps,w_lorentz + CCTK_REAL lam(5) + CCTK_REAL gxx,gxy,gxz,gyy,gyz,gzz + CCTK_REAL alp,beta,u + CCTK_REAL temp,ye + + CCTK_REAL cs2,one,two + CCTK_REAL vlowx,vlowy,vlowz,v2,w + CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamm_nobeta,lamp_nobeta + CCTK_INT handle + CCTK_REAL dpdrho,dpdeps,press + +#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 +! end EOS Omni vars +#endif + + one = 1.0d0 + two = 2.0d0 + +!!$ Set required fluid quantities + call EOS_Omni_cs2(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,eps,temp,ye,cs2,keyerr,anyerr) + + 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 = w_lorentz + +!!$ Calculate eigenvalues + + lam1 = velx - beta/alp + lam2 = velx - beta/alp + lam3 = velx - beta/alp + lamp_nobeta = (velx*(one-cs2) + sqrt(cs2*(one-v2)*& + (u*(one-v2*cs2) - velx**2*(one-cs2))))/(one-v2*cs2) + lamm_nobeta = (velx*(one-cs2) - sqrt(cs2*(one-v2)*& + (u*(one-v2*cs2) - velx**2*(one-cs2))))/(one-v2*cs2) + + lamp = lamp_nobeta - beta/alp + lamm = lamm_nobeta - beta/alp + + lam(1) = lamm + lam(2) = lam1 + lam(3) = lam2 + lam(4) = lam3 + lam(5) = lamp + + +end subroutine eigenvalues_hot + /*@@ @routine eigenproblem @date Sat Jan 26 01:27:59 2002 File [modified]: GRHydro_HLLEPoly.F90 Delta lines: +123 -48 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_HLLEPoly.F90 2010-11-01 05:49:09 UTC (rev 168) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_HLLEPoly.F90 2010-11-01 06:43:05 UTC (rev 169) @@ -46,6 +46,7 @@ CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r + CCTK_REAL :: xtemp CCTK_INT :: type_bits, trivial, not_trivial @@ -179,22 +180,47 @@ !!$ Eigenvalues and fluxes either side of the cell interface if (flux_direction == 1) then - call eigenvalues(GRHydro_eos_handle,& - rhominus(i+xoffset,j+yoffset,k+zoffset),& - velxminus(i+xoffset,j+yoffset,k+zoffset),& - velyminus(i+xoffset,j+yoffset,k+zoffset),& - velzminus(i+xoffset,j+yoffset,k+zoffset),& - epsminus(i+xoffset,j+yoffset,k+zoffset),& - w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& - lamminus,gxxh,gxyh,gxzh,gyyh,& - gyzh,gzzh,& - usendh,avg_alp,avg_beta) - call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& - velxplus(i,j,k),velyplus(i,j,k),& - velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& - lamplus,gxxh,gxyh,gxzh,gyyh,& - gyzh,gzzh,& - usendh,avg_alp,avg_beta) + if(evolve_temp.ne.1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + lamplus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + xtemp,& + y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k), & + xtemp,y_e_plus(i,j,k),& + w_lorentzplus(i,j,k),& + lamplus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + endif call num_x_flux(consp(1),consp(2),consp(3),consp(4),consp(5),& fplus(1),fplus(2),fplus(3),fplus(4),& fplus(5),velxplus(i,j,k),pressplus(i,j,k),& @@ -206,22 +232,47 @@ pressminus(i+xoffset,j+yoffset,k+zoffset),& avg_det,avg_alp,avg_beta) else if (flux_direction == 2) then - call eigenvalues(GRHydro_eos_handle,& - rhominus(i+xoffset,j+yoffset,k+zoffset),& - velyminus(i+xoffset,j+yoffset,k+zoffset),& - velzminus(i+xoffset,j+yoffset,k+zoffset),& - velxminus(i+xoffset,j+yoffset,k+zoffset),& - epsminus(i+xoffset,j+yoffset,k+zoffset),& - w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& - lamminus,gyyh,gyzh,gxyh,gzzh,& - gxzh,gxxh,& - usendh,avg_alp,avg_beta) - call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& - velyplus(i,j,k),velzplus(i,j,k),& - velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& - lamplus,gyyh,gyzh,gxyh,gzzh,& - gxzh,gxxh,& - usendh,avg_alp,avg_beta) + if(evolve_temp.ne.1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velyplus(i,j,k),velzplus(i,j,k),& + velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + lamplus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + xtemp,& + y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,rhoplus(i,j,k),& + velyplus(i,j,k),velzplus(i,j,k),& + velxplus(i,j,k),epsplus(i,j,k),& + xtemp,y_e_plus(i,j,k),& + w_lorentzplus(i,j,k),& + lamplus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + endif call num_x_flux(consp(1),consp(3),consp(4),consp(2),consp(5),& fplus(1),fplus(3),fplus(4),fplus(2),& fplus(5),velyplus(i,j,k),pressplus(i,j,k),& @@ -233,22 +284,46 @@ pressminus(i+xoffset,j+yoffset,k+zoffset),& avg_det,avg_alp,avg_beta) else if (flux_direction == 3) then - call eigenvalues(GRHydro_eos_handle,& - rhominus(i+xoffset,j+yoffset,k+zoffset),& - velzminus(i+xoffset,j+yoffset,k+zoffset),& - velxminus(i+xoffset,j+yoffset,k+zoffset),& - velyminus(i+xoffset,j+yoffset,k+zoffset),& - epsminus(i+xoffset,j+yoffset,k+zoffset),& - w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& - lamminus,gzzh,gxzh,gyzh,& - gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) - call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& - velzplus(i,j,k),velxplus(i,j,k),& - velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& - lamplus,gzzh,gxzh,gyzh,& - gxxh,gxyh,gyyh,& - usendh,avg_alp,avg_beta) + if(evolve_temp.ne.1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velzplus(i,j,k),velxplus(i,j,k),& + velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + lamplus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvalues_hot(GRHydro_eos_handle,rhoplus(i,j,k),& + velzplus(i,j,k),velxplus(i,j,k),& + velyplus(i,j,k),epsplus(i,j,k),& + xtemp,y_e_plus(i,j,k),& + w_lorentzplus(i,j,k),& + lamplus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + endif call num_x_flux(consp(1),consp(4),consp(2),consp(3),consp(5),& fplus(1),fplus(4),fplus(2),fplus(3),& fplus(5),velzplus(i,j,k),pressplus(i,j,k),avg_det,& File [modified]: GRHydro_ParamCheck.F90 Delta lines: +1 -1 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_ParamCheck.F90 2010-11-01 05:49:09 UTC (rev 168) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_ParamCheck.F90 2010-11-01 06:43:05 UTC (rev 169) @@ -98,7 +98,7 @@ evolve_Y_e = 0 endif - if (CCTK_EQUALS(temperature_and_entropy_evolution_method,"GRHydro")) then + if (CCTK_EQUALS(temp_and_ent_evolution_method,"GRHydro")) then evolve_temp = 1 else evolve_temp = 0 File [modified]: GRHydro_Prim2Con.F90 Delta lines: +0 -3 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-01 05:49:09 UTC (rev 168) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-01 06:43:05 UTC (rev 169) @@ -141,9 +141,6 @@ endif - - - end subroutine primitive2conservative /*@@ File [modified]: GRHydro_RiemannSolve.F90 Delta lines: +1 -1 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_RiemannSolve.F90 2010-11-01 05:49:09 UTC (rev 168) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_RiemannSolve.F90 2010-11-01 06:43:05 UTC (rev 169) @@ -39,7 +39,7 @@ CCTK_INT :: i,j,k if (CCTK_EQUALS(riemann_solver,"HLLE")) then - + call GRHydro_HLLE(CCTK_PASS_FTOF) if (evolve_tracer .ne. 0) then File [modified]: GRHydro_UpdateMask.F90 Delta lines: +8 -8 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_UpdateMask.F90 2010-11-01 05:49:09 UTC (rev 168) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_UpdateMask.F90 2010-11-01 06:43:05 UTC (rev 169) @@ -297,10 +297,10 @@ vely(i,j,k) = 0.0d0 velz(i,j,k) = 0.0d0 #if USE_EOS_OMNI - call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) - call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& - rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) #else press(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps(i,j,k)) eps(i,j,k) = EOS_SpecificIntEnergy(eos_handle, rho_min, press(i,j,k)) @@ -321,9 +321,9 @@ vely_p(i,j,k) = 0.0d0 velz_p(i,j,k) = 0.0d0 #if USE_EOS_OMNI - call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr) - call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr) #else press_p(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps_p(i,j,k)) @@ -346,9 +346,9 @@ vely_p_p(i,j,k) = 0.0d0 velz_p_p(i,j,k) = 0.0d0 #if USE_EOS_OMNI - call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr) - call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr) #else press_p_p(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps_p_p(i,j,k))

1 0

[svn:einsteintoolkit] HydroBase/trunk/ (Rev. 32)
by cott＠tapir.caltech.edu 31 Oct '10

31 Oct '10

User: cott Date: 2010/11/01 12:58 AM Modified: /trunk/ param.ccl Log: * change parameter name from: temperature_and_entropy_evolution_method to temp_and_ent_evolution_method to make Intel compiler happy. File Changes: Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +1 -1 =================================================================== --- trunk/param.ccl 2010-10-31 16:14:06 UTC (rev 31) +++ trunk/param.ccl 2010-11-01 05:58:34 UTC (rev 32) @@ -41,7 +41,7 @@ "none" :: "Evolution for Y_e is disabled" } "none" -KEYWORD temperature_and_entropy_evolution_method "Evolution method for temperature and entropy" +KEYWORD temp_and_ent_evolution_method "Evolution method for temperature and entropy" { "none" :: "Evolution for temperature and entropy is disabled" } "none"

1 0

[svn:einsteintoolkit] GRHydro/branches/hot_and_MHD_temp_dev/src/ (Rev. 168)
by cott＠tapir.caltech.edu 31 Oct '10

31 Oct '10

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

1 0

[svn:einsteintoolkit] GRHydro/branches/hot_and_MHD_temp_dev/src/ (Rev. 167)
by cott＠tapir.caltech.edu 31 Oct '10

31 Oct '10

User: cott Date: 2010/11/01 12:33 AM Modified: /branches/hot_and_MHD_temp_dev/src/ GRHydro_Prim2Con.F90 Log: * add hot p2c routines File Changes: Directory: /branches/hot_and_MHD_temp_dev/src/ ============================================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +188 -53 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-10-31 16:19:38 UTC (rev 166) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-01 05:33:26 UTC (rev 167) @@ -46,44 +46,104 @@ integer :: i, j, k CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr - - do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 - do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 - do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + CCTK_REAL :: xtemp - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + if(evolve_temp.ne.1) then + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + + gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + + call prim2con(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,& + gyzl,gzzl, & + avg_detl,densminus(i,j,k),sxminus(i,j,k),& + syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& + rhominus(i,j,k), & + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& + epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k)) + + call prim2con(GRHydro_eos_handle, gxxr,gxyr,gxzr,& + gyyr,gyzr,gzzr, & + avg_detr, densplus(i,j,k),sxplus(i,j,k),& + syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& + rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + w_lorentzplus(i,j,k)) + + end do + end do + end do + else + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + + gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) - call prim2con(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& - syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), & - velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& - epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k)) - - call prim2con(GRHydro_eos_handle, gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& - syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& - rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& - velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& - w_lorentzplus(i,j,k)) + ! we do not have plus/minus vars for temperature since + ! eps is reconstructed. Hence, we do not update the + ! variable 'temperature' in prim2con at the interfaces + xtemp = temperature(i,j,k) + call prim2con_hot(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,& + gyzl,gzzl, & + avg_detl,densminus(i,j,k),sxminus(i,j,k),& + syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& + rhominus(i,j,k), & + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& + epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k),& + xtemp,y_e_minus(i,j,k)) - end do - end do - end do + ! we do not have plus/minus vars for temperature since + ! eps is reconstructed. Hence, we do not update the + ! variable 'temperature' in prim2con at the interfaces + xtemp = temperature(i,j,k) + call prim2con_hot(GRHydro_eos_handle, gxxr,gxyr,gxzr,& + gyyr,gyzr,gzzr, & + avg_detr, densplus(i,j,k),sxplus(i,j,k),& + syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& + rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + w_lorentzplus(i,j,k),xtemp, & + y_e_plus(i,j,k)) + + end do + end do + end do + endif + + + + end subroutine primitive2conservative /*@@ @@ -154,6 +214,59 @@ end subroutine prim2con +subroutine prim2con_hot(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & + dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w, & + temp,ye) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,& + deps, dpress, w, vlowx, vlowy, vlowz + CCTK_REAL :: temp,ye + CCTK_INT :: handle + +#if !USE_EOS_OMNI +#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 +! end EOS Omni vars +#endif + + w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz & + *dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz& + *dvely*dvelz)) + +#if USE_EOS_OMNI + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) +#else + dpress = EOS_Pressure(handle, drho, deps) +#endif + + vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz + vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz + vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz + + ddens = sqrt(det) * drho * w + dsx = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowx + dsy = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowy + dsz = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowz + dtau = sqrt(det) * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens + +end subroutine prim2con_hot + + /*@@ @routine Primitive2ConservativeCells @date Sun Mar 10 21:16:20 2002 @@ -180,26 +293,7 @@ CCTK_INT :: i, j, k CCTK_REAL :: det - - !$OMP PARALLEL DO PRIVATE(i, j) - do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 - do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 - do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) - call prim2con(GRHydro_eos_handle,gxx(i,j,k),& - gxy(i,j,k),gxz(i,j,k),& - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& - tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& - eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) - end do - end do - end do - !$OMP END PARALLEL DO - if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 @@ -212,7 +306,48 @@ !$OMP END PARALLEL DO endif + if(evolve_temp.ne.1) then + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + + det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + call prim2con(GRHydro_eos_handle,gxx(i,j,k),& + gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end do + end do + end do + !$OMP END PARALLEL DO + else + !$OMP PARALLEL DO PRIVATE(i, j) + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + + det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + call prim2con_hot(GRHydro_eos_handle,gxx(i,j,k),& + gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & + temperature(i,j,k),y_e(i,j,k)) + end do + end do + end do + !$OMP END PARALLEL DO + endif + + end subroutine Primitive2ConservativeCells

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits November 2010