Commits November 2012

commits@lists.einsteintoolkit.org

9 participants
90 discussions

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 179)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_BondiM_new.F90, GRHydro_Bondi_new.F90 Log: GRHydro_InitData: port smoother velocity matching to mag. bondi routine. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM_new.F90 Delta lines: +41 -12 =================================================================== --- trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:06 UTC (rev 178) +++ trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:09 UTC (rev 179) @@ -64,7 +64,7 @@ CCTK_REAL :: r_bondi(N_points), logr_bondi(N_points), rho_bondi(N_points), u_bondi(N_points), v_bondi(N_points) CCTK_REAL :: drhodr, det, rhocheck, rhocheck2, riso, rnew, rsch, ucheck CCTK_REAL :: uiso, uisocheck, vcheck, ucheck2, vcheck2, xhat,yhat, zhat, xp, yp, zp - CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm + CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm,dudr,uisocheck2,auiso,buiso character(400) :: debug_message @@ -118,7 +118,8 @@ endif enddo - rhotmp = -1. !!$ start with guess +!!$ rhotmp = -1. !!$ start with guess + rhotmp=rhos !!$ start with value at sonic point! do i=imin,N_points rspher = r_bondi(i) @@ -132,7 +133,8 @@ v_bondi(i) = vtmp end do - rhotmp = -1. +!!$ rhotmp = -1. + rhotmp=rhos !!$ start with value at sonic point! do i=imin-1,1,-1 rspher = r_bondi(i) @@ -165,27 +167,42 @@ !!$ write(debug_message,'(a,4f22.14)') "i=1500:",r_bondi(1500),rho_bondi(1500),& !!$ u_bondi(1500),v_bondi(1500) !!$ call CCTK_INFO(debug_message) - -!!$ // find the derivative near r=M +!!$ // find the derivative near r=M in isotropic coords = r=9/4M in schwarzschild; rnew = 2.25 * M - j = floor ( 0.5 + (log10(rnew) - logrmin) / dlogr ) + j = floor ((log10(rnew) - logrmin) / dlogr + 1.0) +!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck = rho_bondi(j) call find_bondi_solution(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot ) uisocheck = 4.0*vcheck/3.0 +!!$ the previous point was r=M in isotropic coords = r=9/4M in schwarzschild; this one is r=1.01M in isotropic rnew = 0.25 * 3.02**2 * M/1.01 - j = floor( 0.5 + (log10(rnew) - logrmin) / dlogr ) + j = floor((log10(rnew) - logrmin) / dlogr + 1.0) +!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck2 = rho_bondi(j) call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot ) + uisocheck2 = vcheck2 / (1.0 - 1.0/2.02) / (1.0+ 1.0/2.02) + drhodr = 100.0*(rhocheck2-rhocheck)/M +!!$ Don't divide by M here, to simplify the math + dudr = 100.0*(uisocheck2-uisocheck) + write(debug_message,'(a,3f22.14)') 'Rhocheck:',rhocheck,rhocheck2,drhodr call CCTK_INFO(debug_message) + write(debug_message,'(a,3f22.14)') 'Ucheck:',uisocheck,uisocheck2,dudr + call CCTK_INFO(debug_message) nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) + + + write(debug_message,'(a,3f22.14)') 'Lower bound coordinates:',x(1,1,1),y(1,1,1),z(1,1,1) + call CCTK_INFO(debug_message) + write(debug_message,'(a,3f22.14)') 'Upper bound coordinates:',x(nx,ny,nz),y(nx,ny,nz),z(nx,ny,nz) + call CCTK_INFO(debug_message) do i=1,nx do j=1,ny @@ -203,12 +220,17 @@ if(roverm > ONEmTINY) then rsch = 0.25 * ( 2.*riso + M)**2 / riso - jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) +!!$ jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) + jb = floor( (log10(rsch) - logrmin) / dlogr + 1.0) +!!$ jb = NINT((log10(rsch) - logrmin) / dlogr + 1.0) - if(jb > N_points)jb = N_points + if(jb >= N_points) jb = N_points-1 - rhotmp = rho_bondi(jb) - call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) + rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*& + (rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb)) +!!$ call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, & +!!$ rho_bondi(jb),rho_bondi(jb+1)) + call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) rho(i,j,k) = rhotmp uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso) else @@ -218,7 +240,14 @@ rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0 endif utmp = Kval * rho(i,j,k)**( gam ) / (gam - 1.) - uiso = uisocheck * riso / M + +!!$ match to uiso and dudr at roverm=1 +!!$ a R + b R^3 ---> a+b = uisocheck; a+3b = dudr +!!$ b = (dudr-uisocheck)/2; a=3*uisocheck-dudr)/2 + + auiso = 1.5*uisocheck-0.5*dudr + buiso = 0.5*dudr-0.5*uisocheck + uiso = roverm*(auiso+buiso*roverm**2) endif eps(i,j,k) = utmp/rhotmp File [modified]: GRHydro_Bondi_new.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:06 UTC (rev 178) +++ trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:09 UTC (rev 179) @@ -148,7 +148,7 @@ u_bondi(i),v_bondi(i) end do close(47) - + !!$ write(debug_message,'(a,4f22.14)') "i=1:",r_bondi(1),rho_bondi(1),& !!$ u_bondi(1),v_bondi(1) !!$ call CCTK_INFO(debug_message)

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 178)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_Bondi_new.F90 Log: GRHydro_InitData: exchange find_bondi_solution_bracket for find_bondi_solution and improve smoothing of v field. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Bondi_new.F90 Delta lines: +10 -7 =================================================================== --- trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:02 UTC (rev 177) +++ trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:06 UTC (rev 178) @@ -162,13 +162,14 @@ !!$ u_bondi(1500),v_bondi(1500) !!$ call CCTK_INFO(debug_message) -!!$ // find the derivative near r=M +!!$ // find the derivative near r=M in isotropic coords = r=9/4M in schwarzschild; rnew = 2.25 * M j = floor ((log10(rnew) - logrmin) / dlogr + 1.0) !!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck = rho_bondi(j) - call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, & - rho_bondi(j),rho_bondi(j+1)) +!!$ call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, & +!!$ rho_bondi(j),rho_bondi(j+1)) + call find_bondi_solution(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot) uisocheck = 4.0*vcheck/3.0 !!$ the previous point was r=M in isotropic coords = r=9/4M in schwarzschild; this one is r=1.01M in isotropic @@ -176,8 +177,9 @@ j = floor((log10(rnew) - logrmin) / dlogr + 1.0) !!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck2 = rho_bondi(j) - call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, & - rho_bondi(j),rho_bondi(j+1)) +!!$ call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, & +!!$ rho_bondi(j),rho_bondi(j+1)) + call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot) uisocheck2 = vcheck2 / (1.0 - 1.0/2.02) / (1.0+ 1.0/2.02) drhodr = 100.0*(rhocheck2-rhocheck)/M @@ -224,8 +226,9 @@ rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*& (rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb)) - call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, & - rho_bondi(jb),rho_bondi(jb+1)) +!!$ call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, & +!!$ rho_bondi(jb),rho_bondi(jb+1)) + call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) rho(i,j,k) = rhotmp uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso) else

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 177)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_Bondi_new.F90 Log: GRHydro_InitData: improve error messages for find_bondi_solution_bracket func. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Bondi_new.F90 Delta lines: +14 -8 =================================================================== --- trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:00 UTC (rev 176) +++ trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:02 UTC (rev 177) @@ -164,7 +164,8 @@ !!$ // find the derivative near r=M rnew = 2.25 * M - j = floor ((log10(rnew) - logrmin) / dlogr )+1.0 + j = floor ((log10(rnew) - logrmin) / dlogr + 1.0) +!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck = rho_bondi(j) call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, & rho_bondi(j),rho_bondi(j+1)) @@ -172,7 +173,8 @@ !!$ the previous point was r=M in isotropic coords = r=9/4M in schwarzschild; this one is r=1.01M in isotropic rnew = 0.25 * 3.02**2 * M/1.01 - j = floor((log10(rnew) - logrmin) / dlogr )+1.0 + j = floor((log10(rnew) - logrmin) / dlogr + 1.0) +!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck2 = rho_bondi(j) call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, & rho_bondi(j),rho_bondi(j+1)) @@ -215,11 +217,13 @@ if(roverm > ONEmTINY) then rsch = 0.25 * ( 2.*riso + M)**2 / riso !!$ jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) - jb = floor( (log10(rsch) - logrmin) / dlogr ) +1.0 + jb = floor( (log10(rsch) - logrmin) / dlogr + 1.0) +!!$ jb = NINT((log10(rsch) - logrmin) / dlogr + 1.0) - if(jb >= N_points)jb = N_points-1 + if(jb >= N_points) jb = N_points-1 - rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*(rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb)) + rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*& + (rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb)) call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, & rho_bondi(jb),rho_bondi(jb+1)) rho(i,j,k) = rhotmp @@ -424,9 +428,11 @@ CCTK_INT :: max_newt_iter_b, extra_newt_iter_b if(rho.gt.rho1.or.rho.lt.rho2) then - write(6,*)'Very Bad rho! ',rho1,rho2,rho + write(6,*)'find_bondi_solution_bracket: Very bad rho! (rholow,rhoup,rho) = ',rho1,rho2,rho stop endif + + write(6,*)'find_bondi_solution_bracket: (rholow,rhoup,rho) = ',rho1,rho2,rho max_newt_iter_b = 30 newt_tol_b = 1.0e-15 @@ -498,7 +504,7 @@ endif if(x.gt.rho1.or.x.lt.rho2) then - write(6,*)'Bad rho, bad! ',rho1,rho2,rho + write(6,*)'find_bondi_solution_bracket: Bad rho, bad! (rholow,rhoup,rho) = ',rho1,rho2,rho stop endif @@ -533,7 +539,7 @@ if( ((abs(errx) <= newt_tol_b).and.(doing_extra == 0)) .or. & (i_extra > extra_newt_iter_b) .or. (n_iter >= (max_newt_iter_b-1)) ) then keep_iterating = 0 - if(n_iter>=(max_newt_iter_b-1))write(6,*)'Extra iterations!' + if(n_iter>=(max_newt_iter_b-1))write(6,*)'find_bondi_solution_bracket: Extra iterations!' endif !!$ if(r>8.1043 .and. r<8.1044)write(*,*)'iter4:',errx,newt_tol_b,keep_iterating, &

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 176)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_Bondi_new.F90 Log: Add find_bondi_solution_bracket and velocity smoothing from Josh's patch. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Bondi_new.F90 Delta lines: +188 -12 =================================================================== --- trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:58:58 UTC (rev 175) +++ trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:00 UTC (rev 176) @@ -58,7 +58,7 @@ CCTK_REAL :: r_bondi(N_points), logr_bondi(N_points), rho_bondi(N_points), u_bondi(N_points), v_bondi(N_points) CCTK_REAL :: drhodr, det, rhocheck, rhocheck2, riso, rnew, rsch, ucheck CCTK_REAL :: uiso, uisocheck, vcheck, ucheck2, vcheck2, xhat,yhat, zhat, xp, yp, zp - CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm + CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm,dudr,uisocheck2,auiso,buiso character(400) :: debug_message @@ -112,7 +112,8 @@ endif enddo - rhotmp = -1. !!$ start with guess +!!$ rhotmp = -1. !!$ start with guess + rhotmp=rhos !!$ start with value at sonic point! do i=imin,N_points rspher = r_bondi(i) @@ -126,7 +127,8 @@ v_bondi(i) = vtmp end do - rhotmp = -1. +!!$ rhotmp = -1. + rhotmp=rhos !!$ start with value at sonic point! do i=imin-1,1,-1 rspher = r_bondi(i) @@ -162,23 +164,39 @@ !!$ // find the derivative near r=M rnew = 2.25 * M - j = floor ( 0.5 + (log10(rnew) - logrmin) / dlogr ) + j = floor ((log10(rnew) - logrmin) / dlogr )+1.0 rhocheck = rho_bondi(j) - call find_bondi_solution(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot ) + call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, & + rho_bondi(j),rho_bondi(j+1)) uisocheck = 4.0*vcheck/3.0 +!!$ the previous point was r=M in isotropic coords = r=9/4M in schwarzschild; this one is r=1.01M in isotropic rnew = 0.25 * 3.02**2 * M/1.01 - j = floor( 0.5 + (log10(rnew) - logrmin) / dlogr ) + j = floor((log10(rnew) - logrmin) / dlogr )+1.0 rhocheck2 = rho_bondi(j) - call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot ) + call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, & + rho_bondi(j),rho_bondi(j+1)) + uisocheck2 = vcheck2 / (1.0 - 1.0/2.02) / (1.0+ 1.0/2.02) + drhodr = 100.0*(rhocheck2-rhocheck)/M +!!$ Don't divide by M here, to simplify the math + dudr = 100.0*(uisocheck2-uisocheck) + write(debug_message,'(a,3f22.14)') 'Rhocheck:',rhocheck,rhocheck2,drhodr call CCTK_INFO(debug_message) + write(debug_message,'(a,3f22.14)') 'Ucheck:',uisocheck,uisocheck2,dudr + call CCTK_INFO(debug_message) nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) + + + write(debug_message,'(a,3f22.14)') 'Lower bound coordinates:',x(1,1,1),y(1,1,1),z(1,1,1) + call CCTK_INFO(debug_message) + write(debug_message,'(a,3f22.14)') 'Upper bound coordinates:',x(nx,ny,nz),y(nx,ny,nz),z(nx,ny,nz) + call CCTK_INFO(debug_message) do i=1,nx do j=1,ny @@ -196,12 +214,14 @@ if(roverm > ONEmTINY) then rsch = 0.25 * ( 2.*riso + M)**2 / riso - jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) +!!$ jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) + jb = floor( (log10(rsch) - logrmin) / dlogr ) +1.0 - if(jb > N_points)jb = N_points + if(jb >= N_points)jb = N_points-1 - rhotmp = rho_bondi(jb) - call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) + rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*(rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb)) + call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, & + rho_bondi(jb),rho_bondi(jb+1)) rho(i,j,k) = rhotmp uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso) else @@ -211,7 +231,14 @@ rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0 endif utmp = Kval * rho(i,j,k)**( gam ) / (gam - 1.) - uiso = uisocheck * riso / M + +!!$ match to uiso and dudr at roverm=1 +!!$ a R + b R^3 ---> a+b = uisocheck; a+3b = dudr +!!$ b = (dudr-uisocheck)/2; a=3*uisocheck-dudr)/2 + + auiso = 1.5*uisocheck-0.5*dudr + buiso = 0.5*dudr-0.5*uisocheck + uiso = roverm*(auiso+buiso*roverm**2) endif eps(i,j,k) = utmp/rhotmp @@ -382,3 +409,152 @@ end subroutine find_bondi_solution + subroutine find_bondi_solution_bracket(r, rho, u, v, rs, rhos, M, Mdot, Kval, gam, Qdot,rho1,rho2 ) + implicit none + + CCTK_REAL rho1,rho2 + CCTK_REAL :: r, rho, u, v, rs, rhos, M, Mdot, Kval, gam, Qdot + CCTK_REAL :: ur,r_sol, rho_old + CCTK_REAL :: f, df, dx, x_old, resid, jac + CCTK_REAL :: errx, x_orig, x + CCTK_INT :: n_iter, i_extra, doing_extra, keep_iterating, i_increase + CCTK_REAL :: vp, h, hp, term + + CCTK_REAL :: newt_tol_b,small_bondi + CCTK_INT :: max_newt_iter_b, extra_newt_iter_b + + if(rho.gt.rho1.or.rho.lt.rho2) then + write(6,*)'Very Bad rho! ',rho1,rho2,rho + stop + endif + + max_newt_iter_b = 30 + newt_tol_b = 1.0e-15 + extra_newt_iter_b = 2 + small_bondi = 1.0e-20 + +!!$ if(r>8.1043 .and. r<8.1044)write(*,*)'init guess:',r,rho +!!$ write(*,*)'init guess:',r,rho + +!!$ if (rho < 0.) then +!!$ if( r > 0.9*rs .and. r < 1.1*rs ) then +!!$ rho = rhos +!!$ else +!!$ if(r < rs) then +!!$ ur = r**(-0.5) +!!$ else +!!$ ur = 0.5*r**(-1.5) +!!$ endif +!!$ rho = Mdot / (4.*M_PI * r * r * ur) +!!$ endif +!!$ endif + +!!$ if(r>8.1043 .and. r<8.1044) +!!$ write(*,*)'init guess:',r,ur,rho,rs,rhos,Mdot +!!$ if(r<1.0001e-10)write(*,*)'init guess:',r,ur,rho,rs,rhos,Mdot +!!$ write(*,*)'init guess:',r,ur,rho,rs,rhos,Mdot + +!!$ set global variables needed by residual function: + + r_sol = r + + +!!$ Use Newton's method to find rho: +!!$ gnr_bondi( rho, NEWT_DIM_B, bondi_resid) + errx = 1.0 + df=1.0 + f=1.0 + doing_extra = 0 + + rho_old=rho + x=rho + + n_iter = 0 + +!!$ Start the Newton-Raphson iterations : + keep_iterating = 1 + do while( keep_iterating == 1 ) + + hp = Kval * gam * x**(gam - 2.) !!$ // dh/drho + h = 1. + hp * x / ( gam - 1. ) + v = Mdot / ( 4. * M_PI * r_sol * r_sol * x ) + vp = -v / x !!$ // dv/drho + term = 1. - 2.*M/r_sol + v*v + resid = -Qdot + h * h * term + jac = 2. * h *( hp*term + h*v*vp ) + dx = -resid / jac + f = 0.5*resid**2 + df = -2. * f + + errx = 0. + x_old = x + + x=x+dx + + if(x.gt.rho1.or.x.lt.rho2) then +!!$ write(6,*)'Bad rho! ',rho1,rho2,rho + if(x.gt.rho1)x=0.5*(x_old+rho1) + if(x.lt.rho2)x=0.5*(x_old+rho2) + endif + + if(x.gt.rho1.or.x.lt.rho2) then + write(6,*)'Bad rho, bad! ',rho1,rho2,rho + stop + endif + +!!$ if(r>8.1043 .and. r<8.1044)write(*,*)'iter:',x,dx,resid,jac + +!!$ /****************************************/ +!!$ /* Calculate the convergence criterion */ +!!$ /****************************************/ + +!!$ /* For the new criterion, always look at relative error in indep. variable: */ +!!$ // METHOD specific: + if(x==0) then + errx = abs(dx) + x = small_bondi + else + errx = abs(dx/x) + endif + +!!$ /*****************************************************************************/ +!!$ /* If we've reached the tolerance level, then just do a few extra iterations */ +!!$ /* before stopping */ +!!$ /*****************************************************************************/ + +!!$ if(r>8.1043 .and. r<8.1044)write(*,*)'iter3:',errx,newt_tol_b,keep_iterating, & +!!$ doing_extra,i_extra + + if((abs(errx)<=newt_tol_b) .and. (doing_extra == 0) .and. (extra_newt_iter_b > 0)) & + doing_extra=1 + + if( doing_extra == 1 ) i_extra=i_extra+1 + + if( ((abs(errx) <= newt_tol_b).and.(doing_extra == 0)) .or. & + (i_extra > extra_newt_iter_b) .or. (n_iter >= (max_newt_iter_b-1)) ) then + keep_iterating = 0 + if(n_iter>=(max_newt_iter_b-1))write(6,*)'Extra iterations!' + endif + +!!$ if(r>8.1043 .and. r<8.1044)write(*,*)'iter4:',errx,newt_tol_b,keep_iterating, & +!!$ doing_extra,i_extra + + + n_iter=n_iter+1 + + end do + + rho=x + +!!$ Calculate other quantities: + u = Kval * rho**( gam ) / (gam - 1.) + v = Mdot / ( 4. * M_PI * r * r * rho ) + +!!$ if(r>8.1043 .and. r<8.1044) +!!$ write(*,*)'final:',r,rho,u,v + + + return + + end subroutine find_bondi_solution_bracket +

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/par/ (Rev. 175)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:58 PM Added: /trunk/par/ bondiM.par Log: GRHydro_InitData: magnetized Bondi par file. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/par/ ====================== File [added]: bondiM.par Delta lines: +728 -0 =================================================================== --- trunk/par/bondiM.par (rev 0) +++ trunk/par/bondiM.par 2012-11-09 01:58:58 UTC (rev 175) @@ -0,0 +1,728 @@ +#============================================================================== +# Magnetized bondi accretion into a BH represented by a smoothed puncture: +# McLachlan GRHydro_InitData SmoothedPuncture +#============================================================================== + +#------------------------------------------------------------------------------ +# Cactus parameters: +#------------------------------------------------------------------------------ +Cactus::cctk_run_title = "Magnetized bondi on smoothed puncture data" +Cactus::cctk_full_warnings = "yes" + +#Cactus::terminate = "never" +#Cactus::terminate = "time" +#Cactus::cctk_final_time = 1500.0 +#Cactus::terminate = "iteration" +#Cactus::cctk_itlast = 1 +#Cactus::cctk_itlast = 128 +#Cactus::cctk_itlast = 2560 +Cactus::terminate="runtime" +#Cactus::max_runtime=28800 # 46 days +#Cactus::max_runtime=1410 # 23.5 hours +#Cactus::max_runtime=1380 # 23 hours +#Cactus::max_runtime=690 # 11.5 hours +#Cactus::max_runtime=450 # 7.5 hours +Cactus::max_runtime=210 # 3.5 hours + +#------------------------------------------------------------------------------ +# Activate all necessary thorns: +#------------------------------------------------------------------------------ + +ActiveThorns = "Boundary CartGrid3D CoordBase Fortran InitBase IOUtil LocalReduce SymBase Time" +ActiveThorns = "AEILocalInterp LocalInterp" +#ActiveThorns = "MoL Slab SpaceMask SphericalSurface" +#ActiveThorns = "MoL ReflectionSymmetry Slab SpaceMask SphericalSurface" +#ActiveThorns = "MoL ReflectionSymmetry RotatingSymmetry180 Slab SpaceMask SphericalSurface" +ActiveThorns = "MoL ReflectionSymmetry RotatingSymmetry90 Slab SpaceMask SphericalSurface" +#ActiveThorns = "MoL ReflectionSymmetry Slab SpaceMask SphericalSurface" +ActiveThorns = "Carpet CarpetInterp CarpetIOASCII CarpetIOScalar CarpetLib CarpetIOBasic CarpetReduce CarpetRegrid2 CarpetSlab CarpetTracker CarpetMask LoopControl CarpetIOHDF5 CarpetEvolutionMask" +#ActiveThorns = "Formaline NaNChecker TerminationTrigger TimerReport" +ActiveThorns = "NaNChecker TerminationTrigger TimerReport" +ActiveThorns = "ADMbase ADMcoupling ADMmacros CoordGauge StaticConformal" +#ActiveThorns = "PunctureTracker" +#ActiveThorns = "ADManalysis Constants TmunuBase HydroBase " +ActiveThorns = "Constants TmunuBase HydroBase " +#ActiveThorns = "Constants TmunuBase " +ActiveThorns = "QuasiLocalMeasures" +#ActiveThorns = "ADMConstraints" +#ActiveThorns = "EOS_Base EOS_IdealFluid EOS_Polytrope" +ActiveThorns = "EOS_Omni" +ActiveThorns = "GRHydro GRHydro_InitData" +ActiveThorns = "SmoothedPuncture" +ActiveThorns = "BLAS LAPACK GSL" +#ActiveThorns = "TwoPunctures" +ActiveThorns = "SummationByParts" +#ActiveThorns = "TOVSolver" +ActiveThorns = "GenericFD NewRad" +ActiveThorns = "ML_BSSN ML_BSSN_Helper ML_ADMConstraints" +#ActiveThorns = "Hydro_Analysis" +ActiveThorns = "Dissipation" +ActiveThorns = "AHFinderDirect" +#ActiveThorns = "WeylScal4 Multipole" +#ActiveThorns = "SetMask_SphericalSurface" +#ActiveThorns = "SystemStatistics" + +#------------------------------------------------------------------------------ +# Diagnostic parameters: +#------------------------------------------------------------------------------ +AHFinderDirect::verbose_level = "physics details" + +#Carpet::output_timers_every = 1024 +#Carpet::output_timer_tree_every = 256 +#Carpet::output_initialise_timer_tree = "yes" +Carpet::storage_verbose = "no" +#Carpet::timers_verbose = "no" +#Carpet::verbose = "yes" # only for Initial Data +Carpet::verbose = "no" +Carpet::veryverbose = "no" +Carpet::grid_structure_filename = "carpet-grid-structure.asc" +Carpet::grid_coordinates_filename = "carpet-grid-coordinates.asc" + +CarpetLib::output_bboxes = "no" +CarpetLib::print_memstats_every = 1024 +CarpetLib::memstat_file = "carpetlib-memory-statistics.asc" +#CarpetLib::print_timestats_every = 1024 +#CarpetLib::timestat_file = "carpetlib-timing-statistics.asc" + +#CarpetEvolutionMask::verbose = "yes" + +CarpetMask::verbose = "yes" +CarpetReduce::verbose = "no" +#CarpetRegrid2::verbose = "yes" +#CarpetRegrid2::veryverbose = "yes" + +#CarpetTracker::verbose = "yes" + +Dissipation::verbose = "no" + +#NaNChecker::verbose = "all" + +#PunctureTracker::verbose = "yes" + +TimerReport::out_every = 1024 +TimerReport::out_filename = "TimerReport" +#TimerReport::output_all_timers = "yes" +#TimerReport::output_all_timers_together = "yes" +#TimerReport::output_all_timers_readable = "yes" +#TimerReport::before_checkpoint = "yes" +TimerReport::n_top_timers = 40 + +#TwoPunctures::verbose = "yes" + +#QuasiLocalMeasures::verbose = "yes" +#SphericalSurface::verbose = "yes" + +#------------------------------------------------------------------------------ +# Utility parameters: +#------------------------------------------------------------------------------ + +#NaNChecker::check_every = 16 +NaNChecker::check_every = 1 +NaNChecker::check_vars = " + ADMBase::curv + ADMBase::metric + ADMBase::lapse + ADMBase::shift + HydroBase::rho + HydroBase::eps + HydroBase::press + HydroBase::vel + ML_BSSN::ML_log_confac + ML_BSSN::ML_metric + ML_BSSN::ML_trace_curv + ML_BSSN::ML_curv + ML_BSSN::ML_Gamma + ML_BSSN::ML_lapse + ML_BSSN::ML_shift + ML_BSSN::ML_dtlapse + ML_BSSN::ML_dtshift +" +NaNChecker::action_if_found = "terminate" +#NaNChecker::action_if_found = "abort" +#NaNChecker::action_if_found = "just warn" #"terminate", "just warn", "abort" +#NaNChecker::check_after=0 + +#TerminationTrigger::max_walltime = 0.0 # hours +#TerminationTrigger::max_walltime = 24.0 # hours +#TerminationTrigger::on_remaining_walltime = 10.0 # minutes +TerminationTrigger::termination_from_file = "yes" +TerminationTrigger::create_termination_file = "yes" +TerminationTrigger::termination_file = "../TERMINATE" + + +#------------------------------------------------------------------------------ +# Run parameters: +#------------------------------------------------------------------------------ + +#------ +# Grid: +#------ + +Time::dtfac = 0.25 +MoL::ODE_Method = "rk4" +MoL::MoL_Intermediate_Steps = 4 +MoL::MoL_Num_Scratch_Levels = 1 + +CartGrid3D::type = "coordbase" +CartGrid3D::domain = "full" +CartGrid3D::avoid_origin = "no" + +CoordBase::domainsize = "minmax" +CoordBase::spacing = "gridspacing" # "gridspacing" or "numcells" + +CoordBase::xmin = 0.00 +CoordBase::ymin = 0.00 +CoordBase::zmin = 0.00 +CoordBase::xmax = 61.44 +CoordBase::ymax = 61.44 +CoordBase::zmax = 61.44 +CoordBase::dx = 1.28 # dx or ncells_x +CoordBase::dy = 1.28 # dy or ncells_y +CoordBase::dz = 1.28 # dz or ncells_z + +CoordBase::boundary_size_x_lower = 3 +CoordBase::boundary_size_y_lower = 3 +CoordBase::boundary_size_z_lower = 3 +CoordBase::boundary_size_x_upper = 3 +CoordBase::boundary_size_y_upper = 3 +CoordBase::boundary_size_z_upper = 3 +CoordBase::boundary_shiftout_x_lower = 1 +CoordBase::boundary_shiftout_y_lower = 1 +CoordBase::boundary_shiftout_z_lower = 1 +CoordBase::boundary_shiftout_x_upper = 0 +CoordBase::boundary_shiftout_y_upper = 0 +CoordBase::boundary_shiftout_z_upper = 0 + +Driver::ghost_size = 3 +GRHydro::GRHydro_stencil = 3 + +ReflectionSymmetry::reflection_x = "no" +ReflectionSymmetry::reflection_y = "no" +ReflectionSymmetry::reflection_z = "yes" +ReflectionSymmetry::avoid_origin_x = "no" +ReflectionSymmetry::avoid_origin_y = "no" +ReflectionSymmetry::avoid_origin_z = "no" + + +# General Carpet parameters: +Carpet::domain_from_coordbase = "yes" +Carpet::enable_all_storage = "no" +Carpet::use_buffer_zones = "yes" +Carpet::schedule_barriers = "no" + +Carpet::poison_new_timelevels = "yes" +Carpet::check_for_poison = "no" +#Carpet::poison_value = 113 # git version only! + +Carpet::init_3_timelevels = "no" +Carpet::init_fill_timelevels = "yes" + +CarpetEvolutionMask::enlarge_evolved_region_by = 0 +CarpetEvolutionMask::writeNaNs = "yes" +CarpetEvolutionMask::enforce_mask = "yes" +CarpetEvolutionMask::enforce_vars = " + ML_BSSN::ML_lapse + ML_BSSN::ML_shift + GRHydro::dens + GRHydro::scon + GRHydro::tau + GRHydro::Bcons +" + +CarpetLib::poison_new_memory = "yes" +CarpetLib::poison_value = 114 +CarpetLib::check_bboxes = "no" +CarpetLib::interleave_communications = "yes" +#CarpetLib::vary_tags = "yes" # git version only! +CarpetLib::combine_sends = "yes" +#CarpetLib::reduce_mpi_waitall = "yes" # git version only! + +CarpetInterp::tree_search = "yes" +CarpetInterp::check_tree_search = "no" + +CarpetRegrid2::freeze_unaligned_levels = "yes" +#CarpetRegrid2::snap_to_coarse = "yes" +CarpetRegrid2::symmetry_rotating90 = "yes" + +# System specific Carpet parameters: +Carpet::max_refinement_levels = 7 +Carpet::prolongation_order_space = 5 +Carpet::prolongation_order_time = 2 + +CarpetRegrid2::regrid_every = 0 +CarpetRegrid2::num_centres = 1 + +CarpetRegrid2::num_levels_1 = 7 +CarpetRegrid2::position_x_1 = 0.0 +CarpetRegrid2::position_y_1 = 0.0 +CarpetRegrid2::position_z_1 = 0.0 +CarpetRegrid2::radius_1[1] = 30.72 # 4 +CarpetRegrid2::radius_1[2] = 15.36 # 2 +CarpetRegrid2::radius_1[3] = 7.68 # 1 +CarpetRegrid2::radius_1[4] = 3.84 # 0.5 +CarpetRegrid2::radius_1[5] = 1.92 # 0.25 +CarpetRegrid2::radius_1[6] = 0.96 # 0.125 + +#CarpetRegrid2::radius_1[1] = 30.72 # 4 +#CarpetRegrid2::radius_1[2] = 15.36 # 2 +#CarpetRegrid2::radius_1[3] = 7.68 # 1 +#CarpetRegrid2::radius_1[4] = 3.84 # 0.5 +#CarpetRegrid2::radius_1[5] = 1.92 # 0.25 +#CarpetRegrid2::radius_1[6] = 0.96 # 0.125 + +#CarpetMask::excluded_surface [0] = 0 +#CarpetMask::excluded_surface_factor[0] = 1.0 + +#CarpetTracker::surface[0] = 0 + +#------ +# MODEL: +#------ + +HydroBase::initial_hydro = "magnetized_bondi_solution_iso" + +#GRHydro_InitData::mdot_sonicpt_bondi = 1.e-6 +GRHydro_InitData::mdot_sonicpt_bondi = 3.e-3 +#GRHydro_InitData::mdot_sonicpt_bondi = 1.0 +GRHydro_InitData::r_sonicpt_bondi = 8.1 +GRHydro_InitData::bondi_coordinates = "Isotropic" +GRHydro_InitData::bondi_rmax[0] = 100. +GRHydro_InitData::bondi_rmin[0] = 0.1 +GRHydro_InitData::bondi_bmag = 0.01 +#GRHydro_InitData::bondi_bmag = 1.0e-9 +#GRHydro_InitData::bondi_bmag = 0.0 +#GRHydro_InitData::bondi_bmag = 1.0e-3 + +EOS_Omni::gl_gamma = 1.666666666666667 +GRHydro::initial_Gamma = 1.666666666666667 + +ADMBase::initial_data = "smoothedpuncture" +ADMBase::initial_lapse = "smoothedpuncture" + +ADMBase::metric_type = "physical" +#ADMBase::initial_data = "twopunctures" +#ADMBase::initial_lapse = "twopunctures-averaged" +ADMBase::initial_shift = "zero" +ADMBase::initial_dtlapse = "zero" +ADMBase::initial_dtshift = "zero" +# needed for AHFinderDirect: +ADMBase::metric_timelevels = 3 + +#TwoPunctures::schedule_in_ADMBase_InitialData = "no" +#TwoPunctures::keep_u_around = "yes" + +#TwoPunctures::npoints_A = 16 +#TwoPunctures::npoints_B = 16 +#TwoPunctures::npoints_phi = 14 +#TwoPunctures::Newton_tol = 1.e-10 +#TwoPunctures::grid_setup_method = "evaluation" +##TwoPunctures::grid_setup_method = "Taylor expansion" +# +#TwoPunctures::par_b = 1.0 +# +#TwoPunctures::par_m_plus = 1.0 +#TwoPunctures::par_m_minus = 0.0 +# +#TwoPunctures::Newton_maxit = 10 +#TwoPunctures::do_residuum_debug_output = "no" +#TwoPunctures::do_initial_debug_output = "no" +##TwoPunctures::TP_Tiny = 1.0e-10 +#TwoPunctures::TP_Tiny = 1.0e-2 +#TwoPunctures::center_offset[0] = -1.0 +#TwoPunctures::center_offset[1] = 0 +#TwoPunctures::center_offset[2] = 0 +# +#---------- +# Numerics: +#---------- + +InitBase::initial_data_setup_method = "init_all_levels" + +TmunuBase::stress_energy_storage = "no" +TmunuBase::stress_energy_at_RHS = "no" +TmunuBase::timelevels = 1 +TmunuBase::prolongation_type = "none" +TmunuBase::support_old_CalcTmunu_mechanism = "no" + +HydroBase::timelevels = 3 + +ADMMacros::spatial_order = 4 +SummationByParts::order = 4 + +SpaceMask::use_mask = "yes" + +#Dissipation::ah_slope = 0.2 +#Dissipation::epsdis = 0.01 +#Dissipation::extra_dissipation_in_horizons = "yes" +#Dissipation::ah_max_epsdis=0.1 +#Dissipation::update_ah_every = 16 +#Dissipation::surface_number[0] = 0 +#Dissipation::horizon_number[0] = 1 +Dissipation::order = 5 +Dissipation::vars = " + ML_BSSN::ML_log_confac + ML_BSSN::ML_metric + ML_BSSN::ML_curv + ML_BSSN::ML_trace_curv + ML_BSSN::ML_Gamma + ML_BSSN::ML_lapse + ML_BSSN::ML_shift + ML_BSSN::ML_dtlapse + ML_BSSN::ML_dtshift +" + +SphericalSurface::nsurfaces = 1 +SphericalSurface::maxntheta = 39 +SphericalSurface::maxnphi = 76 + +SphericalSurface::ntheta [0] = 39 +SphericalSurface::nphi [0] = 76 +SphericalSurface::nghoststheta[0] = 2 +SphericalSurface::nghostsphi [0] = 2 +#SphericalSurface::radius [0] = 0.5 + +#SetMask_SphericalSurface::SetMask_SurfaceIndex[0] = 0 +#SetMask_SphericalSurface::SetMask_RadiusFactor[0] = 0.75 + + +#----------- +# Evolution: +#----------- + +HydroBase::evolution_method = "GRHydro" +HydroBase::Bvec_evolution_method = "GRHydro" + +GRHydro::method_type = "RSA FV" +GRHydro::recon_method = "tvd" +GRHydro::tvd_limiter = "vanleerMC2" +GRHydro::recon_vars = "primitive" +GRHydro::riemann_solver = "HLLE" +GRHydro::decouple_normal_Bfield = "yes" +GRHydro::bound = "none" +GRHydro::GRHydro_eos_type = "General" +GRHydro::GRHydro_eos_table = "Ideal_Fluid" + +GRHydro::rho_abs_min = 1.0e-13 +Grhydro::initial_rho_abs_min = 1.0e-13 +Grhydro::grhydro_atmo_tolerance = 0.0001 + +GRHydro::clean_divergence = "yes" +GRHydro::psidcspeed = "light speed" +GRHydro::kap_dc = 10.0 +GRHydro::track_divB = "yes" +GRHydro::calculate_bcom = "yes" + + +GRHydro::c2p_reset_pressure = "yes" +#GRHydro::GRHydro_enable_internal_excision = "false" +# +GRHydro::Grhydro_MaxNumConstrainedVars = 23 +GRHydro::GRHydro_MaxNumEvolvedVars = 10 + +GRHydro::GRHydro_MaxNumSandRVars = 16 +#GRHydro::GRHydro_MaxNumConstrainedVars = 12 +#GRHydro::GRHydro_MaxNumEvolvedVars = 0 +#GRHydro::GRHydro_MaxNumSandRVars = 0 +#GRHydro::GRHydro_MaxNumEvolvedVars = 6 # When 1 tracer is evolved! + +GRHydro::use_evolution_mask = "yes" + +ADMBase::evolution_method = "ML_BSSN" +ADMBase::lapse_evolution_method = "ML_BSSN" +ADMBase::shift_evolution_method = "ML_BSSN" +ADMBase::dtlapse_evolution_method= "ML_BSSN" +ADMBase::dtshift_evolution_method= "ML_BSSN" + +ML_BSSN::timelevels = 3 + +#ML_BSSN::conformalMethod = 1 # 1 for W, 0 for phi + +ML_BSSN::harmonicN = 1 # 1+log +ML_BSSN::harmonicF = 2.0 # 1+log +ML_BSSN::ShiftGammaCoeff = 0.75 +ML_BSSN::BetaDriver = 1.0 +ML_BSSN::LapseAdvectionCoeff = 1.0 +ML_BSSN::ShiftAdvectionCoeff = 1.0 + +ML_BSSN::MinimumLapse = 1.0e-8 + +ML_BSSN::my_initial_boundary_condition = "extrapolate-gammas" +ML_BSSN::my_rhs_boundary_condition = "NewRad" +Boundary::radpower = 2 + +ML_BSSN::ML_log_confac_bound = "none" +ML_BSSN::ML_metric_bound = "none" +ML_BSSN::ML_Gamma_bound = "none" +ML_BSSN::ML_trace_curv_bound = "none" +ML_BSSN::ML_curv_bound = "none" +ML_BSSN::ML_lapse_bound = "none" +ML_BSSN::ML_dtlapse_bound = "none" +ML_BSSN::ML_shift_bound = "none" +ML_BSSN::ML_dtshift_bound = "none" + +#------------------------------------------------------------------------------ +# Analysis: +#------------------------------------------------------------------------------ +AHFinderDirect::find_every = 512 + +#AHFinderDirect::run_at_CCTK_ANALYSIS = "yes" +#AHFinderDirect::run_at_CCTK_POSTSTEP = "no" +#AHFinderDirect::run_at_CCTK_POST_RECOVER_VARIABLES = "no" + +#AHFinderDirect::move_origins = "yes" +#AHFinderDirect::reshape_while_moving = "yes" +#AHFinderDirect::predict_origin_movement = "yes" + +# Hermite to order 3 to avoid discontinuities in the metric spatial derivatives: +AHFinderDirect::geometry_interpolator_name = "Hermite polynomial interpolation" +AHFinderDirect::geometry_interpolator_pars = "order=3" +AHFinderDirect::surface_interpolator_name = "Hermite polynomial interpolation" +AHFinderDirect::surface_interpolator_pars = "order=3" +#AHFinderDirect::geometry_interpolator_name = "Lagrange polynomial interpolation" +#AHFinderDirect::geometry_interpolator_pars = "order=4" +#AHFinderDirect::surface_interpolator_name = "Lagrange polynomial interpolation" +#AHFinderDirect::surface_interpolator_pars = "order=4" + +AHFinderDirect::output_h_every = 0 + +AHFinderDirect::N_horizons = 1 + +AHFinderDirect::reset_horizon_after_not_finding [1] = "no" +AHFinderDirect::initial_guess__coord_sphere__radius [1] = 0.45 +AHFinderDirect::origin_x [1] = 0.0 +AHFinderDirect::origin_y [1] = 0.0 +AHFinderDirect::origin_z [1] = 0.0 +AHFinderDirect::initial_guess__coord_sphere__x_center[1] = 0.0 +AHFinderDirect::initial_guess__coord_sphere__y_center[1] = 0.0 +AHFinderDirect::initial_guess__coord_sphere__z_center[1] = 0.0 +#AHFinderDirect::track_origin_from_grid_scalar [1] = "yes" +#AHFinderDirect::track_origin_source_x [1] = "PunctureTracker::pt_loc_x[0]" +#AHFinderDirect::track_origin_source_y [1] = "PunctureTracker::pt_loc_y[0]" +#AHFinderDirect::track_origin_source_z [1] = "PunctureTracker::pt_loc_z[0]" +AHFinderDirect::which_surface_to_store_info [1] = 0 + + +#Hydro_Analysis::Hydro_Analysis_comp_rho_max = "true" + +#PunctureTracker::track [0] = "yes" +#PunctureTracker::initial_x [0] = 10.64502257 +#PunctureTracker::initial_y [0] = 0 +#PunctureTracker::initial_z [0] = 0 +#PunctureTracker::which_surface_to_store_info[0] = 0 + +#QuasiLocalMeasures::num_surfaces = 1 +#QuasiLocalMeasures::spatial_order = 4 +#QuasiLocalMeasures::interpolator = "Lagrange polynomial interpolation" +#QuasiLocalMeasures::interpolator_options = "order=4" +#QuasiLocalMeasures::surface_index [0] = 0 + +#Multipole::nradii = 4 +#Multipole::radius[0] = 80 +#Multipole::radius[1] = 100 +#Multipole::radius[2] = 120 +#Multipole::radius[3] = 140 +#Multipole::ntheta = 120 +#Multipole::nphi = 240 +#Multipole::variables = "WeylScal4::Psi4r{sw=-2 cmplx='WeylScal4::Psi4i' name='psi4'}" +#Multipole::out_every = 4 +#Multipole::l_max = 4 +#Multipole::integration_method = "Simpson" + +#WeylScal4::fd_order = "4th" + +#------------------------------------------------------------------------------ +# Output: +#------------------------------------------------------------------------------ + +IO::out_dir = $parfile +IO::out_fileinfo = "none" + +IOBasic::outInfo_every = 1 +IOBasic::outInfo_reductions = "maximum" +IOBasic::outInfo_vars = " + Carpet::physical_time_per_hour + GRHydro::bcom_sq + HydroBase::press + HydroBase::rho + HydroBase::eps + ML_ADMConstraints::ML_Ham +" + +IOScalar::outScalar_every = 128 +#IOScalar::outScalar_every = 1 +IOScalar::one_file_per_group = "yes" +IOScalar::outScalar_reductions = "minimum maximum average norm1 norm2" +IOScalar::outScalar_vars = " + ADMBase::curv + ADMBase::lapse + ADMBase::metric + ADMBase::shift + GRHydro::bcom + GRHydro::bcom0 + GRHydro::bcom_sq + GRHydro::Bcons + GRHydro::dens + GrHydro::divB + GRHydro::scon + GRHydro::tau + HydroBase::Bvec + HydroBase::eps + HydroBase::press + HydroBase::rho + HydroBase::vel + HydroBase::w_lorentz + ML_ADMConstraints::ML_Ham + ML_ADMConstraints::ML_mom + ML_BSSN::ML_Ham + ML_BSSN::ML_mom +# SystemStatistics::process_memory_mb +" + +IOASCII::one_file_per_group = "yes" +#IOASCII::output_symmetry_points = "no" +#IOASCII::output_ghost_points = "no" #hg version only! +#IOASCII::output_boundary_points = "no" #hg version only! +#IOASCII::out3D_ghosts = "no" #hg (DEPRECATED) +#IOASCII::out3D_outer_ghosts = "no" #hg (DEPRECATED) + +IOASCII::out0D_every = 256 +IOASCII::out0D_vars = " + Carpet::timing +" + +#Set these IOASCII options for initial data only: +IOASCII::out1D_every = 640 +#IOASCII::out1D_every = 1 +IOASCII::out1D_d = "no" +IOASCII::out1D_vars = " + ADMBase::curv + ADMBase::lapse + ADMBase::metric + ADMBase::shift + GRHydro::bcom + GRHydro::bcom0 + GRHydro::bcom_sq + GRHydro::Bcons + GRHydro::dens + GrHydro::divB + GRHydro::scon + GRHydro::tau + HydroBase::Bvec + HydroBase::eps + HydroBase::press + HydroBase::rho + HydroBase::vel + HydroBase::w_lorentz + ML_BSSN::ML_Ham + ML_BSSN::ML_mom + ML_ADMConstraints::ML_Ham + ML_ADMConstraints::ML_mom +" + +#CarpetIOASCII::one_file_per_group = "yes" +#CarpetIOASCII::open_one_input_file_at_a_time = "yes" +#CarpetIOASCII::out2D_every = 256 +#CarpetIOASCII::out2D_xy = "yes" +#CarpetIOASCII::out2D_xz = "no" +#CarpetIOASCII::out2D_yz = "no" +#CarpetIOASCII::out2D_xyplane_z = 0.0 +#CarpetIOASCII::out2D_vars = " +# CarpetReduce::weight +# GRHydro::GRHydro_tracers +# Grid::coordinates +# HydroBase::eps +# HydroBase::press +# HydroBase::rho +# HydroBase::vel +# HydroBase::w_lorentz +# ADMBase::curv +# ADMBase::lapse +# ADMBase::metric +# ADMBase::shift +# ML_ADMConstraints::ML_Ham +# ML_ADMConstraints::ML_Mom +# WeylScal4::Psi4i +# WeylScal4::Psi4r +# " + +IOHDF5::one_file_per_group = "yes" +IOHDF5::open_one_input_file_at_a_time = "yes" +IOHDF5::output_boundary_points = "yes" +IOHDF5::output_buffer_points = "yes" +IOHDF5::output_ghost_points = "yes" +IOHDF5::output_symmetry_points = "yes" + +IOHDF5::out2D_every = 128 +#IOHDF5::out2D_every = 1 +IOHDF5::out2D_vars = " + HydroBase::Bvec + HydroBase::eps + HydroBase::press + HydroBase::rho + HydroBase::vel + HydroBase::w_lorentz + GRHydro::bcom + GRHydro::bcom0 + GRHydro::bcom_sq + GRHydro::Bcons + GRHydro::dens + GRHydro::divB + GrHydro::psidc + GRHydro::scon + GRHydro::tau + ADMBase::curv + ADMBase::lapse + ADMBase::metric + ADMBase::shift + ML_ADMConstraints::ML_Ham + ML_ADMConstraints::ML_mom +" + +IOHDF5::out3D_every = 128 +#IOHDF5::out3D_every = 1 +IOHDF5::out3D_vars = " + ADMBase::curv + ADMBase::lapse + ADMBase::metric + ADMBase::shift + CarpetReduce::weight + GRHydro::bcom + GRHydro::bcom0 + GRHydro::bcom_sq + GRHydro::Bcons + GRHydro::dens + GRHydro::divB + GRHydro::GRHydro_atmosphere_mask + GRHydro::GRHydro_atmosphere_mask_real + GrHydro::psidc + GRHydro::scon + GRHydro::tau + Grid::coordinates + HydroBase::Bvec + HydroBase::eps + HydroBase::press + HydroBase::rho + HydroBase::vel + HydroBase::w_lorentz + ML_ADMConstraints::ML_Ham + ML_ADMConstraints::ML_Mom +# TwoPunctures::puncture_u +" + +#------------------------------------------------------------------------------ +# Checkpoint/Recovery: +#------------------------------------------------------------------------------ +IOHDF5::checkpoint = "yes" +#IO::checkpoint_dir = $parfile +IO::checkpoint_dir = "./bondiM_checkpoints" +IO::checkpoint_ID = "yes" +IO::checkpoint_every = 6144 +#IO::checkpoint_every_walltime_hours = 6.0 +IO::checkpoint_keep=3 +IO::checkpoint_on_terminate = "yes" + +IO::recover = "autoprobe" +#IO::recover_dir = $parfile +IO::recover_dir = "./bondiM_checkpoints" + +

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 439)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:53 PM Modified: /trunk/ param.ccl /trunk/src/ GRHydro_Con2PrimM.F90 Log: GRHydro: Introduce inequalities cons variables needs to satisfy when B/=0. Introduce temporary prim variables as well. * Please refer to Illinois' paper (arXiv:1112.0568) appendices for details. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +198 -35 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2012-11-09 01:53:21 UTC (rev 438) +++ trunk/src/GRHydro_Con2PrimM.F90 2012-11-09 01:53:23 UTC (rev 439) @@ -15,6 +15,9 @@ #include "GRHydro_InterfacesM.h" #include "GRHydro_Macros.h" +#define ITER_TOL (1.0e-8) +#define MAXITER (50) + /*@@ @routine Conservative2PrimitiveM @date Sep 3, 2010 @@ -58,8 +61,11 @@ DECLARE_CCTK_FUNCTIONS integer :: i, j, k, itracer, nx, ny, nz - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin(1), epsmin(1) - CCTK_REAL :: b2 + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, sdet, pmin(1), epsmin(1) + CCTK_REAL :: oob, b2, d2, s2, bscon, bxhat, byhat, bzhat, bhatscon + CCTK_REAL :: Wm, Wm0, Wm_resid, Wmold + CCTK_REAL :: s2m, s2m0, s2m_resid, s2mold, s2max, taum + CCTK_INT :: niter CCTK_INT :: epsnegative character(len=100) warnline @@ -81,6 +87,17 @@ CCTK_REAL :: g11c, g12c, g13c, g22c, g23c, g33c CCTK_REAL :: tmp1 + ! Save the primitive variables to temporary functions before calling the + ! con2prim pointwise routines: + CCTK_REAL :: rho_tmp, press_tmp, eps_tmp + CCTK_REAL :: velx_tmp, vely_tmp, velz_tmp, w_lorentz_tmp + CCTK_REAL :: Bvecx_tmp, Bvecy_tmp, Bvecz_tmp + + ! Assume 3-metric is positive definite. Check deep inside the horizon + ! if this is actually satisfied and if it is not then cast the metric + !as conformally flat only for con2prim inversion purposes. + posdef = .true. + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then g11 => gaa g12 => gab @@ -145,7 +162,11 @@ !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr,keytemp, & - !$OMP local_perc_ptol,posdef,g11c,g12c,g13c,g22c,g23c,g33c,tmp1) + !$OMP local_perc_ptol,posdef,g11c,g12c,g13c,g22c,g23c,g33c,tmp1, & + !$OMP sdet,d2,s2,oob,bscon,bxhat,byhat,bzhat, & + !$OMP bhatscon,Wm,Wm0,Wm_resid,Wmold,s2m,s2m0,s2m_resid,s2mold,s2max, & + !$OMP taum,niter,rho_tmp,press_tmp,eps_tmp,velx_tmp,vely_tmp,velz_tmp, & + !$OMP w_lorentz_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp) do k = 1, nz do j = 1, ny do i = 1, nx @@ -157,6 +178,7 @@ epsnegative = 0 det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k)) + sdet = sqrt(det) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),& @@ -182,8 +204,14 @@ Y_e(i,j,k) = max(min(Y_e_con(i,j,k) / dens(i,j,k),GRHydro_Y_e_max),& GRHydro_Y_e_min) endif + + + b2=g11(i,j,k)*Bprim(i,j,k,1)**2+g22(i,j,k)*Bprim(i,j,k,2)**2+g33(i,j,k)*Bprim(i,j,k,3)**2+ & + 2.0*(g12(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,2)+g13(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,3)+ & + g23(i,j,k)*Bprim(i,j,k,2)*Bprim(i,j,k,3)) + - if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + if ( dens(i,j,k) .le. sdet*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then !call CCTK_WARN(1,"Con2Prim: Resetting to atmosphere") !write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k) @@ -192,11 +220,8 @@ ! temperature(i,j,k),y_e(i,j,k) !call CCTK_WARN(1,warnline) - b2=g11(i,j,k)*Bprim(i,j,k,1)**2+g22(i,j,k)*Bprim(i,j,k,2)**2+g33(i,j,k)*Bprim(i,j,k,3)**2+ & - 2.0*(g12(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,2)+g13(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,3)+ & - g23(i,j,k)*Bprim(i,j,k,2)*Bprim(i,j,k,3)) - dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + dens(i,j,k) = sdet*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 vup(i,j,k,:) = 0.d0 @@ -229,29 +254,119 @@ !!$ tau does need to take into account the existing B-field !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens [Press drops out] - tau(i,j,k) = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) + tau(i,j,k) = sdet * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) + if(tau(i,j,k).le.sdet*b2*0.5d0)then + tau(i,j,k) = GRHydro_tau_min + sdet*b2*0.5d0 + endif + cycle end if + if(evolve_temper.eq.0) then - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + + + if(sdet.ge.sqrtdet_thr) then + d2 = dens(i,j,k)**2 + s2 = uxx*scon(i,j,k,1)**2 + uyy*scon(i,j,k,2)**2 & + + uzz*scon(i,j,k,3)**2 & + + 2.0d0*uxy*scon(i,j,k,1)*scon(i,j,k,2) & + + 2.0d0*uxz*scon(i,j,k,1)*scon(i,j,k,3) & + + 2.0d0*uyz*scon(i,j,k,2)*scon(i,j,k,3) + oob = 1.0d0/sqrt(b2) + bxhat = oob*Bprim(i,j,k,1) + byhat = oob*Bprim(i,j,k,2) + bzhat = oob*Bprim(i,j,k,3) + bhatscon = bxhat*scon(i,j,k,1)+byhat*scon(i,j,k,2) & + +bzhat*scon(i,j,k,3) + bscon = Bprim(i,j,k,1)*scon(i,j,k,1) & + + Bprim(i,j,k,2)*scon(i,j,k,2) & + + Bprim(i,j,k,3)*scon(i,j,k,3) + ! Initial guesses for iterative procedure to find Wm: + Wm0 = sdet*sqrt(bhatscon**2+d2) + s2m0 = (Wm0**2*s2+bhatscon**2*(b2+2.0d0*Wm0)) & + / (Wm0+b2)**2 + Wm = sdet*sqrt(s2m0+d2) + s2m = (Wm**2*s2+bscon**2*(b2+2.0d0*Wm)) & + / (Wm+b2)**2 + s2m_resid = 1.0d60 + Wm_resid = 1.0d60 + niter = 0 + do while((s2m_resid.ge.ITER_TOL.and.Wm_resid.ge.ITER_TOL).and.& + niter.le.MAXITER) + Wmold = Wm + s2mold = s2m + Wm = sdet*sqrt(s2m+d2) + s2m = (Wm**2*s2+bscon**2*(b2+2.0d0*Wm)) & + / (Wm+b2)**2 + Wm_resid = abs(Wmold-Wm) + s2m_resid = abs(s2mold-s2m) + niter = niter + 1 + end do + !TODO: abort execution if niter .eq. MAXITER and warn user + taum = tau(i,j,k) - 0.5d0*sdet*b2 -0.5d0*(b2*s2-bscon**2)/ & + (sdet*(Wm+b2)**2) + s2max = taum*(taum+2.0d0*dens(i,j,k)) + if(taum.lt.GRHydro_tau_min)then + tau(i,j,k) = GRHydro_tau_min + 0.5d0*sdet*b2 + 0.5d0* & + (b2*s2-bscon**2)/(sdet*(Wm+b2)**2) + end if + if(s2.gt.s2max) then + scon(i,j,k,1) = scon(i,j,k,1)*sqrt(s2max/s2) + scon(i,j,k,2) = scon(i,j,k,2)*sqrt(s2max/s2) + scon(i,j,k,3) = scon(i,j,k,3)*sqrt(s2max/s2) + end if + endif + + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + + keytemp = 0 + + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & + GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & - Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), xye(1), xtemp(1), rho(i,j,k),& - vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),press(i,j,k),& - Bprim(i,j,k,1), Bprim(i,j,k,2), Bprim(i,j,k,3),b2, w_lorentz(i,j,k),& - g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),xye(1), & + xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& + eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& + w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k)) - else + + else ! if(evolve_temper.eq.0) then + + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + keytemp = 0 - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & + GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & - Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), Y_e(i,j,k), temperature(i,j,k), rho(i,j,k),& - vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),press(i,j,k),& - Bprim(i,j,k,1), Bprim(i,j,k,2), Bprim(i,j,k,3), b2, w_lorentz(i,j,k),& - g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), Y_e(i,j,k), & + temperature(i,j,k),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& + eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& + w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k)) if(GRHydro_C2P_failed(i,j,k).ne.0) then @@ -260,12 +375,25 @@ ! accuracy requirement; if it fails again, we abort GRHydro_C2P_failed(i,j,k) = 0 local_perc_ptol = GRHydro_eos_rf_prec*100.0d0 - call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, local_perc_ptol, local_gam(1), dens(i,j,k), & + ! Use the previous primitive values as initial guesses + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & + local_perc_ptol, local_gam(1), dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & - Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), Y_e(i,j,k), temperature(i,j,k), rho(i,j,k),& - vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),press(i,j,k),& - Bprim(i,j,k,1), Bprim(i,j,k,2), Bprim(i,j,k,3), b2, w_lorentz(i,j,k),& - g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), Y_e(i,j,k), & + temperature(i,j,k),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& + eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& + w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k)) if(GRHydro_C2P_failed(i,j,k).ne.0) then @@ -284,7 +412,7 @@ !$OMP END CRITICAL endif endif - endif + endif ! if(evolve_temper.eq.0) then if (epsnegative .ne. 0) then @@ -350,19 +478,50 @@ local_pgam=log(xpress(1)/local_K)/log(xrho(1)) sc = local_K*dens(i,j,k) - call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), & + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & - Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),& - vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),press(i,j,k),& - Bprim(i,j,k,1), Bprim(i,j,k,2), Bprim(i,j,k,3),b2, w_lorentz(i,j,k),& + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,& + velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& + Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k)) - + + +! call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), & +! scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & +! Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),& +! vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),press(i,j,k),& +! Bprim(i,j,k,1), Bprim(i,j,k,2), Bprim(i,j,k,3),b2, w_lorentz(i,j,k),& +! g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & +! uxx,uxy,uxz,uyy,uyz,uzz,det, & + #endif - end if - + end if ! if (epsnegative .ne. 0) then + + rho(i,j,k) = rho_tmp + press(i,j,k) = press_tmp + eps(i,j,k) = eps_tmp + vup(i,j,k,1) = velx_tmp + vup(i,j,k,2) = vely_tmp + vup(i,j,k,3) = velz_tmp + w_lorentz(i,j,k) = w_lorentz_tmp + Bprim(i,j,k,1) = Bvecx_tmp + Bprim(i,j,k,2) = Bvecy_tmp + Bprim(i,j,k,3) = Bvecz_tmp + if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance)) then ! if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) .or. GRHydro_C2P_failed(i,j,k) .ge. 1) then @@ -370,7 +529,7 @@ 2.0*(g12(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,2)+g13(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,3)+ & g23(i,j,k)*Bprim(i,j,k,2)*Bprim(i,j,k,3)) - dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + dens(i,j,k) = sdet*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 vup(i,j,k,:) = 0.d0 @@ -397,8 +556,12 @@ !!$ tau does need to take into account the existing B-field !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens [Press drops out] - tau(i,j,k) = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) + tau(i,j,k) = sdet * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) + if(tau(i,j,k).le.sdet*b2*0.5d0)then + tau(i,j,k) = GRHydro_tau_min + sdet*b2*0.5d0 + endif + end if ! ! Again, reset 3-metric only for con2prim inversion. Restoring Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +5 -0 =================================================================== --- trunk/param.ccl 2012-11-09 01:53:21 UTC (rev 438) +++ trunk/param.ccl 2012-11-09 01:53:23 UTC (rev 439) @@ -351,7 +351,12 @@ 0:1 :: "Either 0 or 1" } 0 +REAL sqrtdet_thr "Threshold to apply cons rescalings deep inside the horizon" STEERABLE=ALWAYS +{ + 1.0: :: "Larger values guarantees this sort of rescaling only deep inside the horizon" +} 30.0 + boolean c2p_reset_pressure "If the pressure guess is unphysical should we recompute it?" { } "no"

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 438)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:53 PM Modified: /trunk/src/ GRHydro_Con2PrimM.F90 Log: GRHydro: Check if 3-metric is definite positive before con2prim inversion. Only do this deep inside the horizon. Reset it to conformally flat if it isn't. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +21 -2 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2012-11-09 01:53:19 UTC (rev 437) +++ trunk/src/GRHydro_Con2PrimM.F90 2012-11-09 01:53:21 UTC (rev 438) @@ -76,6 +76,11 @@ CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim + logical :: posdef + + CCTK_REAL :: g11c, g12c, g13c, g22c, g23c, g33c + CCTK_REAL :: tmp1 + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then g11 => gaa g12 => gab @@ -139,7 +144,8 @@ !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & - !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr,keytemp,local_perc_ptol) + !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr,keytemp, & + !$OMP local_perc_ptol,posdef,g11c,g12c,g13c,g22c,g23c,g33c,tmp1) do k = 1, nz do j = 1, ny do i = 1, nx @@ -151,6 +157,7 @@ epsnegative = 0 det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k)) + call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),& g23(i,j,k),g33(i,j,k)) @@ -393,7 +400,19 @@ tau(i,j,k) = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) end if - + +! ! Again, reset 3-metric only for con2prim inversion. Restoring +! ! the current values for the 3-metric: +! if(.not.posdef)then +! g11(i,j,k) = g11c +! g12(i,j,k) = g12c +! g13(i,j,k) = g13c +! g22(i,j,k) = g22c +! g23(i,j,k) = g23c +! g33(i,j,k) = g33c +! posdef = .true. +! endif +! end do end do end do

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 437)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:53 PM Modified: /trunk/ interface.ccl, param.ccl /trunk/src/ GRHydro_Minima.F90, GRHydro_Prim2ConM.F90 Log: GRHydro: Introduce GRHydro_tau_min = 10^-10 max(tau(t=0)) to be used in tau rescaling in con2prim failures. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Minima.F90 Delta lines: +2 -0 =================================================================== --- trunk/src/GRHydro_Minima.F90 2012-11-09 01:53:16 UTC (rev 436) +++ trunk/src/GRHydro_Minima.F90 2012-11-09 01:53:19 UTC (rev 437) @@ -46,6 +46,8 @@ if (initial_atmosphere_factor > 0.0) GRHydro_rho_min = GRHydro_rho_min * initial_atmosphere_factor + GRHydro_tau_min = tau_rel_min + return end subroutine GRHydro_Rho_Minima_Setup File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +19 -1 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2012-11-09 01:53:16 UTC (rev 436) +++ trunk/src/GRHydro_Prim2ConM.F90 2012-11-09 01:53:19 UTC (rev 437) @@ -453,10 +453,13 @@ CCTK_REAL :: xtemp(1) CCTK_INT :: i, j, k CCTK_REAL :: det + CCTK_REAL :: maxtau0 + + maxtau0 = -1.0d60 if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(k,j,i,det) + !$OMP PARALLEL DO PRIVATE(k,j,i,det), REDUCTION(MAX:maxtau0) 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 @@ -473,10 +476,25 @@ eps(i,j,k),press(i,j,k),Bvecx(i,j,k), & Bvecy(i,j,k), Bvecz(i,j,k), w_lorentz(i,j,k)) + maxtau0 = max(maxtau0,tau(i,j,k)) + end do end do end do !$OMP END PARALLEL DO + + ! TODO: to actually reduce GRHydro_tau_min over all Carpet components + ! we need to modify Carpet looping to allow this function to be called + ! on all AMR levels before calling any other function. The best would be + ! to set a special bin where functions would be called on all levels first + ! instead of calling all functions per level. The workaround for this problem + ! is to set GRHydro_tau_min to a user specified value as it was set in + ! GRHydro_Minima.F90. Once this issue is solved, then uncomment the line + ! below and create two other routines to be run in global mode so that + ! GRHydro_tau_min can be properly initialized and reduced. + + !GRHydro_tau_min = GRHydro_tau_min * maxtau0 + else !$OMP PARALLEL DO PRIVATE(k,j,i,det) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +1 -0 =================================================================== --- trunk/interface.ccl 2012-11-09 01:53:16 UTC (rev 436) +++ trunk/interface.ccl 2012-11-09 01:53:19 UTC (rev 437) @@ -359,6 +359,7 @@ CCTK_REAL GRHydro_minima type = SCALAR { GRHydro_rho_min + GRHydro_tau_min # GRHydro_dens_min } "Atmosphere values" File [modified]: param.ccl Delta lines: +5 -0 =================================================================== --- trunk/param.ccl 2012-11-09 01:53:16 UTC (rev 436) +++ trunk/param.ccl 2012-11-09 01:53:19 UTC (rev 437) @@ -404,6 +404,11 @@ ##Parameters for the atmosphere. +REAL tau_rel_min "A minimum relative tau (taumin = maxtau(t=0) * tau_rel_min) below which tau is reschaled" +{ + 0: :: "" +} 1.e-10 + REAL rho_abs_min "A minimum rho below which evolution is turned off (atmosphere). If negative, the relative minimum will be used instead." { -1.0: :: ""

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 436)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:53 PM Added: /trunk/src/ GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed, GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c Log: This routine hopefully fixes the polytrope calls to magnetized Con2Prim while leaving the hybrid EOS behavior unchanged From: Josh Faber <jfaberrit(a)jfaber3.local> File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed Delta lines: <Unavailable> File [added]: GRHydro_Con2PrimM_pt_EOSOmni.c.hopefullyfixed.c Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 435)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:53 PM Modified: /trunk/ interface.ccl Log: ng the interface for magentized con2prim for consistency From: Josh Faber <jfaberrit(a)jfaber3.local> File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +2 -2 =================================================================== --- trunk/interface.ccl 2012-11-09 01:53:11 UTC (rev 434) +++ trunk/interface.ccl 2012-11-09 01:53:13 UTC (rev 435) @@ -68,10 +68,10 @@ CCTK_INT INOUT GRHydro_reflevel, \ CCTK_REAL OUT retval) -void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ +void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ - CCTK_REAL INOUT rho, \ + CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \

1 0

← Newer
1
2
3
4
5
6
7
8
9
Older →

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits November 2012