c Mahc_Evolve: Numerical Evolution of the General Relativistic Hydro Equations c Copyright (C) 2001 Mark Miller /*@@ @file Mahc_get_prim_var.F @date December 1999 @author Mark Miller @desc Calculate primative hydro variables. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_DefineThorn.h" /*@@ @routine get_prim_var @date December 1999 @author Mark Miller @desc Calculate primative hydro variables on all points of the local grid. @enddesc @calls @calledby @history @endhistory @@*/ subroutine get_prim_var(CCTK_FARGUMENTS) implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS integer CCTK_Equals c LOCAL VARS CCTK_REAL det,uxx,uxy,uxz,uyy,uyz,uzz integer i,j,k,nx,ny,nz logical some_hydro_dust,mahc_eos_tables,mahc_eos_polytrope logical mahc_eos_ideal some_hydro_dust = CCTK_EQUALS(some_hydro,'MahcHydroDust') mahc_eos_tables = CCTK_EQUALS(mahc_eos,'tables') mahc_eos_polytrope = CCTK_EQUALS(mahc_eos,'polytrope') mahc_eos_ideal = CCTK_EQUALS(mahc_eos,'ideal_fluid') nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) do k=1,nz do j=1,ny do i=1,nx c calculate u** at point i,j,k det = gxx(i,j,k)*gyy(i,j,k)*gzz(i,j,k) + . 2.0d0*gxy(i,j,k)*gxz(i,j,k)*gyz(i,j,k) - . gxx(i,j,k) * (gyz(i,j,k)**2) - . gyy(i,j,k) * (gxz(i,j,k)**2) - . gzz(i,j,k) * (gxy(i,j,k)**2) uxx=(-gyz(i,j,k)**2 + gyy(i,j,k)*gzz(i,j,k))/det uxy=(gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k))/det uyy=(-gxz(i,j,k)**2 + gxx(i,j,k)*gzz(i,j,k))/det uxz=(-gxz(i,j,k)*gyy(i,j,k) + gxy(i,j,k)*gyz(i,j,k))/det uyz=(gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k))/det uzz=(-gxy(i,j,k)**2 + gxx(i,j,k)*gyy(i,j,k))/det call get_prim_var_pt(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), . mahc_eos_gamma,uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), . z(i,j,k),r(i,j,k),i,j,k,mahc_eps_min,mahc_densmin, . mahc_eos_k,mahc_perc_ptol,mahc_del_ptol,some_hydro_dust, . mahc_eos_tables,mahc_eos_polytrope,mahc_eos_ideal, . mahc_eos_handle,mahc_primvar_solver,mahc_prim_energyadd) enddo enddo enddo return end /*@@ @routine get_prim_var_pt @date December 1999 @author Mark Miller @desc Calculate primative hydro variables on all points of the local grid. @enddesc @calls @calledby @history @endhistory @@*/ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccc GET_PRIM_VAR_PT cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccc Calculates the primative variables at a point. ccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c #define LORENTZ_SOLVE #define PRESSURE_SOLVE subroutine get_prim_var_pt(de,sx,sy,sz,en,rho, . velx,vely,velz,epsilon,press,w_lorentz,gamma, . uxx,uxy,uxz,uyy,uyz,uzz,det,x,y,z,r, . ival,jval,kval,mahc_eps_min,mahc_densmin, . mahc_eos_k,mahc_perc_ptol,mahc_del_ptol,some_hydro_dust, . mahc_eos_tables,mahc_eos_polytrope,mahc_eos_ideal, . mahc_eos_handle,mahc_primvar_solver,mahc_prim_energyadd) implicit none CCTK_REAL de,sx,sy,sz,en,rho,velx,vely,velz,epsilon,press, . uxx,uxy,uxz,uyy,uyz,uzz,gamma,det,w_lorentz,x,y,z,r CCTK_REAL mahc_eps_min,mahc_densmin,mahc_eos_k CCTK_REAL mahc_perc_ptol,mahc_del_ptol,mahc_prim_energyadd DECLARE_CCTK_PARAMETERS logical some_hydro_dust,mahc_eos_tables,mahc_eos_polytrope logical mahc_eos_ideal integer mahc_eos_handle,mahc_primvar_solver,ival,jval,kval c LOCAL VARS CCTK_REAL s2,c0,c1,c2,c3,c4,f,df,ftol,v2,w,vlowx,vlowy,vlowz CCTK_REAL de_prim,en_prim,sx_prim,sy_prim,sz_prim,test integer itnum,count,countmax,i,count_again,count_done CCTK_REAL ude,usx,usy,usz,uen,pold,pnew,epsold,epsnew,w2,w2mhalf CCTK_REAL s2byd2,s_reduct,unde,unsx,unsy,unsz,rhoold,rhonew,rhonew_gm1 CCTK_REAL eps_eos,press_eos,deps_by_drho,dpress_by_drho,unen CCTK_REAL f1,f2,df1bydrho,df1bydeps,df2bydrho,df2bydeps CCTK_REAL dpressbydrho,dpressbydeps,detmat,temp1 CCTK_REAL drhobydpress,depsbydpress CCTK_REAL sigmanew, sigmaold, rhoval, pressval, epsval,dtaudeps CCTK_REAL dpressdsig,drhodsig,dtau,sigma_c,w_max_temp CCTK_REAL fact,w_temp,tau_new CCTK_REAL unde_in,unen_in,s2_in,epsilon_in,press_in,rho_in #ifdef CACTUSEOS_EOS_BASE #include "EOS_Base.inc" #endif if(de .le. 0.0d0) then rho = 0.0d0 velx = 0.0d0 vely = 0.0d0 velz = 0.0d0 epsilon = 0.0d0 w_lorentz = 1.0d0 press = 0.0d0 return endif c PRIMATIVE VARIABLE SOLVE FOR DUST! if(some_hydro_dust) then if(de .le. 0.0d0) then rho = 0.0d0 velx = 0.0d0 vely = 0.0d0 velz = 0.0d0 epsilon = 0.0d0 press = 0.0d0 sx = 0.0d0 sy = 0.0d0 sz = 0.0d0 en = 0.0d0 else s2 = sx*sx*uxx + sy*sy*uyy + sz*sz*uzz + 2.0d0*sx*sy*uxy + . 2.0d0*sx*sz*uxz + 2.0d0*sy*sz*uyz s2byd2 = s2/(de*de) w_lorentz = sqrt(1.0d0 + s2byd2) press = 0.0d0 epsilon = 0.0d0 rho = de / (sqrt(det) * w_lorentz) vlowx = sx / (de * w_lorentz) vlowy = sy / (de * w_lorentz) vlowz = sz / (de * w_lorentz) en = de * (w_lorentz - 1.0d0) if(s2byd2 - 1.0d0 .le. 1.0d-9) then en = de * 0.5d0 * s2byd2 * (1.0d0 - 0.25d0 * s2byd2) endif velx = uxx * vlowx + uxy * vlowy + uxz * vlowz vely = uxy * vlowx + uyy * vlowy + uyz * vlowz velz = uxz * vlowx + uyz * vlowy + uzz * vlowz endif return endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccc primative variable solve from EOS tables cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(mahc_eos_tables) then #ifndef CACTUSEOS_EOS_BASE write(*,*) 'If you want to use tables, you have to compile' write(*,*) 'with EOS_Base!' call CCTK_WARN(0,"Mahc_get_prim_var") #else c IF(MAHC_PRIMVAR_SOLVER = SIGMA) if(mahc_primvar_solver .eq. 2) then unde = de/sqrt(det) unsx = sx/sqrt(det) unsy = sy/sqrt(det) unsz = sz/sqrt(det) unen = en/sqrt(det) s2 = unsx*unsx*uxx + unsy*unsy*uyy + unsz*unsz*uzz + . 2.0d0*unsx*unsy*uxy + . 2.0d0*unsx*unsz*uxz + 2.0d0*unsy*unsz*uyz c save some of these, so if things go wrong, we can c print them out to reinact the crime scene unde_in = unde unen_in = unen s2_in = s2 epsilon_in = epsilon press_in = press rho_in = rho sigmanew = epsilon + press/rho sigmaold = sigmanew count = 0 itnum = 0 count_again = 0 temp1 = sqrt(s2 + ((sigmanew + 1.0d0)*unde)**2) rhoval = (sigmanew + 1.0d0)*unde*unde/temp1 pressval = (sigmanew + 1.0d0)*unde*unde/rhoval - . (unde + unen) epsval = sigmanew - pressval/rhoval f = pressval - EOS_Pressure(mahc_eos_handle,rhoval,epsval) dpressbydrho = EOS_DPressByDRho(mahc_eos_handle,rhoval,epsval) dpressbydeps = EOS_DPressByDEps(mahc_eos_handle,rhoval,epsval) do while( (count .le. 3) ) itnum = itnum + 1 if(itnum .ge. 100) then write(*,*) 'sigma primvar solver: too many newton steps' call CCTK_WARN(0,'get_prim_var_pt') endif df = rhoval - ( dpressbydrho * s2 * unde*unde/(temp1**3) - . dpressbydeps*s2*( (unde+unen)/temp1 - 1.0d0)/ . (((sigmanew+1.0d0)*unde)**2) ) if( (abs((sigmaold - sigmanew)/(sigmanew+1.0d-30)) .le. . sig_rdiff_tol) .or. . ( (abs(sigmaold-sigmanew) .le. sig_adiff_tol) .or. . (abs(f) .lt. (unde+unen)*sig_f_tol) )) then count = count + 1 elseif ((abs(f/df*ude**2*s2/temp1**3) .lt. 1.0d-17) .and. . (abs(f/df*ude**2*s2/temp1**3*pressval/rhoval**2) .lt. . 1.0d-17)) then count = count + 1 endif sigmaold = sigmanew sigmanew = sigmaold - f/df if(itnum .eq. 99) then write(*,*) '************************************' write(*,*) 'itnum : ',itnum write(*,*) 'ijk : ',ival,jval,kval write(*,*) 'x : ',x write(*,*) 'y : ',y write(*,*) 'z : ',z write(*,*) 'r : ',r write(*,*) 'sigmaold: ',sigmaold write(*,*) 'sigmanew: ',sigmanew write(*,*) 'f : ',f write(*,*) 'df : ',df write(*,*) 'rhoold : ',rhoval write(*,*) 'epsold : ',epsval write(*,*) 'pressold: ',pressval write(*,*) 'temp1 : ',temp1 write(*,*) 'unde : ',unde write(*,*) 'unen : ',unen write(*,*) 's2 : ',s2 write(*,*) 'wlorentz: ',unde / rhoval write(*,*) 'rho,eps : ',rho,epsilon write(*,*) 'press in: ',press write(*,*) 'pressEOS: ',EOS_Pressure(mahc_eos_handle,rho,epsilon) write(*,*) 'input unde: ',unde_in write(*,*) 'input unen: ',unen_in write(*,*) 'input s2 : ',s2_in write(*,*) 'input eps : ',epsilon_in write(*,*) 'input pres: ',press_in write(*,*) 'input rho : ',rho_in itnum = 99 endif temp1 = sqrt(s2 + ((sigmanew + 1.0d0)*unde)**2) rhoval = (sigmanew + 1.0d0)*unde*unde/temp1 pressval = (sigmanew + 1.0d0)*unde*unde/rhoval - . (unde + unen) epsval = sigmanew - pressval/rhoval c OK, here are all of the nasty, dirty tricks that we will play. c These tricks are driven as indicated by the following c pseudo-code: c c if (dens < sigma_dens_trick) then c if(dens < sigma_dens_min) :set dens to be 1.05*sigma_dens_min c else if (W > sigma_w_max) : decrease sx, sy, and sz to make c W = 0.95 * sigma_w_max c else if (epsilon < sigma_eps_min) : increase tau to make c eps = 1.05 sigma_eps_min c endif c set aux variables c endif c First, make sure we are below some minimum density: if(unde .lt. sigma_dens_trick) then c define sigma_c used in some conditionals: sigma_c = sigmanew if(sigma_c .lt. 0.0d0) sigma_c = 0.0d0 c make sure density is at least some minimum value if(unde .lt. sigma_dens_min) then unde = 1.05d0 * sigma_dens_min mahc_prim_energyadd = mahc_prim_energyadd + . (unde*sqrt(det) - de) de = unde * sqrt(det) elseif(sqrt(1.0d0 + s2 / ( (sigma_c + 1.0d0)*unde)**2) .gt. . sigma_w_max) then c decrease s2 s.t. w_lorentz = 0.95*sigma_w_max w_max_temp = 0.95d0*sigma_w_max if(w_max_temp .lt. 1.0d0) w_max_temp = 1.0d0 fact = (((sigma_c + 1.0d0)*unde)**2)* . (w_max_temp**2 - 1.0d0) fact = sqrt(fact/s2) unsx = unsx * fact unsy = unsy * fact unsz = unsz * fact sx = unsx*sqrt(det) sy = unsy*sqrt(det) sz = unsz*sqrt(det) s2 = unsx*unsx*uxx + unsy*unsy*uyy + unsz*unsz*uzz + . 2.0d0*unsx*unsy*uxy + . 2.0d0*unsx*unsz*uxz + 2.0d0*unsy*unsz*uyz elseif(epsval .lt. sigma_eps_min) then w_temp = sqrt(1.0d0 + s2 / ( (sigma_c + 1.0d0)*unde)**2) tau_new = unde * (sigma_eps_min*1.05d0/w_temp - . sigma_c*(1.0d0 - w_temp**2)/w_temp + w_temp - 1.0d0) if(tau_new .gt. unen) then dtau = tau_new - unen mahc_prim_energyadd = mahc_prim_energyadd + . sqrt(det) * dtau unen = tau_new en = unen * sqrt(det) endif endif c now that we are done with tricks, we should reset the c auxilary variables temp1 = sqrt(s2 + ((sigmanew + 1.0d0)*unde)**2) rhoval = (sigmanew + 1.0d0)*unde*unde/temp1 pressval = (sigmanew + 1.0d0)*unde*unde/rhoval - . (unde + unen) epsval = sigmanew - pressval/rhoval endif f = pressval - EOS_Pressure(mahc_eos_handle,rhoval,epsval) dpressbydrho = EOS_DPressByDRho(mahc_eos_handle,rhoval,epsval) dpressbydeps = EOS_DPressByDEps(mahc_eos_handle,rhoval,epsval) enddo rho = rhoval epsilon = epsval press = pressval w_lorentz = unde/rho vlowx = unsx / ( (sigmanew + 1.0d0) * w_lorentz * unde) vlowy = unsy / ( (sigmanew + 1.0d0) * w_lorentz * unde) vlowz = unsz / ( (sigmanew + 1.0d0) * w_lorentz * unde) velx = uxx * vlowx + uxy * vlowy + uxz * vlowz vely = uxy * vlowx + uyy * vlowy + uyz * vlowz velz = uxz * vlowx + uyz * vlowy + uzz * vlowz return endif c REGULAR PRESSURE SOLVE c undensitize the evolved variables unde = de/sqrt(det) unsx = sx/sqrt(det) unsy = sy/sqrt(det) unsz = sz/sqrt(det) unen = en/sqrt(det) pnew = press pold = pnew s2 = unsx*unsx*uxx + unsy*unsy*uyy + unsz*unsz*uzz + . 2.0d0*unsx*unsy*uxy + . 2.0d0*unsx*unsz*uxz + 2.0d0*unsy*unsz*uyz temp1 = (unen+unde+pnew)**2 - s2 c make sure initial guess for pressure is sensible if (temp1 .le. 0.0d0) then pnew = sqrt(s2) - unde - unen + 1.0d-12*(unde+unen) temp1 = (unen+unde+pnew)**2 - s2 if(temp1 .le. 0.0d0) then write(*,*) 'get_prim_var_pt:alg error: EOS#1' call CCTK_WARN(0,'get_prim_var_pt') endif endif w_lorentz = (unde + unen + pnew) / sqrt(temp1) rho = unde / w_lorentz epsilon = (unen+unde+pnew)/(unde*w_lorentz) - pnew/rho - 1.0d0 if(epsilon .lt. 0.0d0) then write(*,*) 'get_prim_var_pt:alg error: EOS#2' write(*,*) ' initial epsilon is ',epsilon call CCTK_WARN(0,'get_prim_var_pt') endif f1 = pnew - EOS_Pressure(mahc_eos_handle,rho,epsilon) count = 0 count_done = 0 countmax = 100 do while ( count_done .le. 3) c stopping criterion if(abs(f1)/abs(pnew) .le. 1.0d-10) count_done = count_done+1 if ( ((abs(pnew - pold)/abs(pnew) .le. mahc_perc_ptol) .and. . (abs(pnew - pold) .le. mahc_del_ptol)) ) then count_done = count_done+1 endif count = count + 1 if (count .gt. countmax) then write(*,*) 'get_prim_var(press): error: did not converge in ' write(*,*) ' ',countmax,' steps' write(*,*) 'function is ',f1 write(*,*) 'pressure is ',pnew write(*,*) 'old pressure is ',pold write(*,*) 'mahc_perc_ptol is ',mahc_perc_ptol write(*,*) 'mahc_perc_ptol condition is ', . abs(pnew - pold)/abs(pnew) write(*,*) 'mahc_del_ptol is ',mahc_del_ptol write(*,*) 'mahc_del_ptol condition is ',abs(pnew - pold) call CCTK_WARN(0,'get_prim_var') endif dpressbydrho = EOS_DPressByDRho(mahc_eos_handle,rho,epsilon) dpressbydeps = EOS_DPressByDEps(mahc_eos_handle,rho,epsilon) drhobydpress = unde * s2 / (sqrt(temp1)*(unde+unen+pnew)**2) depsbydpress = pnew * s2 / (rho * (unde + unen + pnew) * temp1) df = 1.0d0 - dpressbydrho*drhobydpress - . dpressbydeps*depsbydpress pold = pnew pnew = pnew - f1/df if(pnew .le. 0.0d0) then write(*,*) 'get_prim_var_pt:alg error: EOS#3' call CCTK_WARN(0,'get_prim_var_pt') endif temp1 = (unen+unde+pnew)**2 - s2 if (temp1 .le. 0.0d0) then write(*,*) 'get_prim_var_pt:alg error: EOS#4' write(*,*) ' initial epsilon is ',epsilon call CCTK_WARN(0,'get_prim_var_pt') endif w_lorentz = (unde + unen + pnew) / sqrt(temp1) rho = unde / w_lorentz epsilon = (unen+unde+pnew)/(unde*w_lorentz) - pnew/rho - 1.0d0 if(epsilon .lt. 0.0d0) then write(*,*) 'get_prim_var_pt:alg error: EOS#5' write(*,*) ' initial epsilon is ',epsilon call CCTK_WARN(0,'get_prim_var_pt') endif f1 = pnew - EOS_Pressure(mahc_eos_handle,rho,epsilon) enddo c now, reconstruct the primative variables press = pnew vlowx = unsx / ( (rho + rho*epsilon + press)*w_lorentz**2) vlowy = unsy / ( (rho + rho*epsilon + press)*w_lorentz**2) vlowz = unsz / ( (rho + rho*epsilon + press)*w_lorentz**2) velx = uxx * vlowx + uxy * vlowy + uxz * vlowz vely = uxy * vlowx + uyy * vlowy + uyz * vlowz velz = uxz * vlowx + uyz * vlowy + uzz * vlowz return #endif endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccc primative variable solve using zero temperature polytropic EOS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(mahc_eos_polytrope) then c undensitize the evolved variables unde = de/sqrt(det) unsx = sx/sqrt(det) unsy = sy/sqrt(det) unsz = sz/sqrt(det) s2 = unsx*unsx*uxx + unsy*unsy*uyy + unsz*unsz*uzz + . 2.0d0*unsx*unsy*uxy + . 2.0d0*unsx*unsz*uxz + 2.0d0*unsy*unsz*uyz rhonew = rho rhoold = rho #define YES #ifdef NO STOP c0 = mahc_eos_k * gamma / (gamma - 1.0d0) c1 = 1.0d0 + c0*rhonew**(gamma-1.0d0) f = (c1)**2 * ((unde/rhonew)**2 - 1.0d0) - s2/(unde**2) df = 2.0d0 * c1 * mahc_eos_k * gamma * rhonew**(gamma-2.0d0) * . ((unde/rhonew)**2 - 1.0d0) - 2.0d0 * c1**2 * unde*unde / . (rhonew**3) count = 0 ftol = 1.0d-12 do while( (abs(f) .gt. ftol) .or. (count .lt. 5)) count = count + 1 rhoold = rhonew rhonew = rhonew - f/df if(rhonew .le. 0.0d0) then write(*,*) 'hyp_get_prim: rho is le zero!' endif c1 = 1.0d0 + c0*rhonew**(gamma-1.0d0) f = (c1)**2 * ((unde/rhonew)**2 - 1.0d0) - s2/(unde**2) df = 2.0d0 * c1 * mahc_eos_k * gamma * rhonew**(gamma-2.0d0) * . ((unde/rhonew)**2 - 1.0d0) - 2.0d0 * c1**2 * unde*unde / . (rhonew**3) if(count .ge. 100) then write(*,*) 'hyp_get_prim_var_pt: polytrope EOS:' write(*,*) 'rho not converging in 100 newton steps' write(*,*) 'unde,s2: ',unde,s2 write(*,*) 'rho: ',rhonew write(*,*) 'f, df: ',f,df write(*,*) 'ftol: ',ftol write(*,*) 'x,y,z,r: ',x,y,z,r call CCTK_WARN(0,"get_prim_var_pt, pt2") endif enddo #endif #ifdef YES c0 = mahc_eos_k * gamma / (gamma - 1.0d0) rhonew_gm1 = rhonew**(gamma-1.0d0) c1 = 1.0d0 + c0 * rhonew_gm1 w_lorentz = sqrt(1.0d0 + s2/( (unde*c1)**2)) f = rhonew * w_lorentz - unde count = 0 #ifndef SINGLE_PRECISION ftol = 1.0d-12 #else ftol = 1.0e-6 #endif do while( (abs(f) .gt. ftol*de) .or. (count .lt. 3)) count = count + 1 rhoold = rhonew df = w_lorentz - s2 * mahc_eos_k * gamma * rhonew_gm1 / . (w_lorentz * unde * unde * c1**3) rhonew = rhonew - f/df if(rhonew .le. 0.0d0) then rhonew = 1.0d-15 endif rhonew_gm1 = rhonew**(gamma-1.0d0) c1 = 1.0d0 + c0 * rhonew_gm1 w_lorentz = sqrt(1.0d0 + s2/( (unde*c1)**2)) f = rhonew * w_lorentz - unde if(count .ge. 100) then write(*,*) 'hyp_get_prim_var_pt: polytrope EOS:' write(*,*) 'rho not converging in 100 newton steps' write(*,*) 'unde,s2: ',unde,s2 write(*,*) 'rho: ',rhonew write(*,*) 'f, df: ',f,df write(*,*) 'ftol: ',ftol write(*,*) 'x,y,z,r: ',x,y,z,r call CCTK_WARN(0,"get_prim_var_pt, pt3") endif enddo if(rhonew .lt. mahc_densmin) then rhonew = mahc_densmin rhonew_gm1 = rhonew**(gamma-1.0d0) sx = 0.0d0 sy = 0.0d0 sz = 0.0d0 unsx = 0.0d0 unsy = 0.0d0 unsz = 0.0d0 endif #endif c now, reconstruct the primative variables, pressure, and tau rho = rhonew epsilon = mahc_eos_k * rhonew_gm1/(gamma - 1.0d0) press = mahc_eos_k * rho * rhonew_gm1 w_lorentz = unde/rho vlowx = unsx / ( (rho + rho*epsilon + press)*w_lorentz**2) vlowy = unsy / ( (rho + rho*epsilon + press)*w_lorentz**2) vlowz = unsz / ( (rho + rho*epsilon + press)*w_lorentz**2) velx = uxx * vlowx + uxy * vlowy + uxz * vlowz vely = uxy * vlowx + uyy * vlowy + uyz * vlowz velz = uxz * vlowx + uyz * vlowy + uzz * vlowz w_lorentz = 1.0d0/sqrt(1.0d0 - vlowx*velx - vlowy*vely - vlowz*velz) en = sqrt(det) * ( (rho + rho*epsilon + press)*w_lorentz**2 - . press) - de return endif if(.not. mahc_eos_ideal) then write(*,*) 'hyp_get_prim_var_pt: bad mahc_eos variable!' STOP endif #ifdef LORENTZ_SOLVE c STOPPING CONDITIONS HARD-CODED IN: countmax = 100 ftol = 1.0d-14 c set up coefficients for lorentz factor equation s2 = sx*sx*uxx + sy*sy*uyy + sz*sz*uzz + 2.*sx*sy*uxy + . 2.*sx*sz*uxz + 2.*sy*sz*uyz c0 = -((-1.d0+gamma)**2*(de**2+s2)) c1 = 2.d0*de*(de+en)*(-1.d0+gamma)*gamma c2 = de**2-2.d0*de**2*gamma- & 2.d0*de*en*gamma**2-en**2*gamma**2-2.d0*g & amma*s2+2.d0*gamma**2*s2 c3 = 2.d0*de*(de+en)*(1.d0-gamma)*gamma c4 = gamma**2*(de**2+2.d0*de*en+en**2-s2) c calculate function and derivative from initial guess w = 0.99 f = c4 * w**4 + c3 * w**3 + c2 * w**2 + c1 * w + c0 df = 4. * c4 * w**3 + 3. * c3 * w**2 + 2. * c2 * w + c1 c write(*,*) 'gpvp: initial f/de2 and df: ',f/(de * de),df if(df*f .gt. 0.0) then write(*,*) 'get_prim_var: error: must make algorithm smarter ...' write(*,*) ' initial guess not good enough !' write(*,*) 'de is ',de write(*,*) 'sx is ',sx write(*,*) 'sy is ',sy write(*,*) 'sz is ',sz write(*,*) 'en is ',en write(*,*) 'c0 is ',c0 write(*,*) 'c1 is ',c1 write(*,*) 'c2 is ',c2 write(*,*) 'c3 is ',c3 write(*,*) 'c4 is ',c4 write(*,*) 'w: is ',w call CCTK_WARN(0,"get_prim_var_pt, pt4") endif c start Newtons method: count = 0 do while( (abs(f/(de * de)) .gt. ftol) .or. (count .lt. 5)) count = count + 1 if (count .gt. countmax) then write(*,*) 'get_prim_var: error: did not converge in ' write(*,*) ' ',countmax,' steps' write(*,*) 'de is ',de write(*,*) 'sx is ',sx write(*,*) 'sy is ',sy write(*,*) 'sz is ',sz write(*,*) 'en is ',en call CCTK_WARN(0,"get_prim_var_pt, pt5") endif c calculate new w w = w - f/df c calculate function and derivative f = c4 * w**4 + c3 * w**3 + c2 * w**2 + c1 * w + c0 df = 4. * c4 * w**3 + 3. * c3 * w**2 + 2. * c2 * w + c1 c write(*,*) 'gpvp: w, f/de2 and df:',w,f/(de*de),df enddo c if close to 1, just make it 1. if(abs(1. - w) .lt. 1.0d-12) then w = 1.0d0 endif c write(*,*) 'get_prim_var_pt: w is ',w if(w .lt. 1.) then write(*,*) 'get_prim_var_pt: error: lorentz factor less than ' write(*,*) ' one: ',w write(*,*) 'gpvp: f and df:',f,df write(*,*) 'de is ',de write(*,*) 'sx is ',sx write(*,*) 'sy is ',sy write(*,*) 'sz is ',sz write(*,*) 'en is ',en write(*,*) 'c4 is ',c4 write(*,*) 'c3 is ',c3 write(*,*) 'c2 is ',c2 write(*,*) 'c1 is ',c1 write(*,*) 'c0 is ',c0 write(*,*) 'det is ',det call CCTK_WARN(0,"get_prim_var_pt, pt6") endif w_lorentz = w c now calculate all the primative variables: rho = de / (sqrt(det) * w) epsilon = w * (en/de + (1. - w)) / ( 1. + gamma * (w * w - 1.)) vlowx = sx / (de * w * (1.0 + gamma * epsilon)) vlowy = sy / (de * w * (1.0 + gamma * epsilon)) vlowz = sz / (de * w * (1.0 + gamma * epsilon)) velx = uxx * vlowx + uxy * vlowy + uxz * vlowz vely = uxy * vlowx + uyy * vlowy + uyz * vlowz velz = uxz * vlowx + uyz * vlowy + uzz * vlowz c write(*,*) 'get_prim_var_pt: rho is ',rho c write(*,*) 'get_prim_var_pt: velx is ',vx c write(*,*) 'get_prim_var_pt: vely is ',vy c write(*,*) 'get_prim_var_pt: velz is ',vz c write(*,*) 'get_prim_var_pt: epsilon is ',epsilon c as a consistency check, calculate residual and yell c if there is a problem de_prim = de - sqrt(det) * rho * w en_prim = en - sqrt(det) * (rho * (1.0d0 + epsilon * gamma) * w * w - . (gamma - 1.0d0)*epsilon*rho - rho * w) sx_prim = sx-sqrt(det)*rho*(1.0d0+epsilon*gamma)*w*w*vlowx sy_prim = sy-sqrt(det)*rho*(1.0d0+epsilon*gamma)*w*w*vlowy sz_prim = sz-sqrt(det)*rho*(1.0d0+epsilon*gamma)*w*w*vlowz press = (gamma - 1.0d0)*rho*epsilon test = max(abs(de_prim),abs(en_prim)) test = max(test,abs(sx_prim)) test = max(test,abs(sy_prim)) test = max(test,abs(sz_prim)) if(test .ge. 1.0d-10) then write(*,*) 'hyp_get_prim_var_pt: error: ' write(*,*) de_prim,en_prim,sx_prim,sy_prim,sz_prim call CCTK_WARN(0,"get_prim_var_pt, pt7") endif #endif #ifdef PRESSURE_SOLVE c SOLVE with P = (gamma-1)*eps*rho countmax = 100 c undensitize the variables ude = de /sqrt( det) usx = sx /sqrt( det) usy = sy /sqrt( det) usz = sz /sqrt( det) uen = en /sqrt( det) s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.*usx*usy*uxy + . 2.*usx*usz*uxz + 2.*usy*usz*uyz c set initial guess for pressure: pold = (gamma - 1.0d0) * rho * epsilon if( (uen + pold + ude)**2 - s2 .le. 0.0d0) then if(ude .lt. mahc_densmin) then s_reduct = ((uen + pold + ude)/sqrt(s2))/2.0d0 sx = sx * s_reduct sy = sy * s_reduct sz = sz * s_reduct usx = usx * s_reduct usy = usy * s_reduct usz = usz * s_reduct s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.0d0*usx*usy*uxy + . 2.0d0*usx*usz*uxz + 2.0d0*usy*usz*uyz endif endif if( (uen + pold + ude)**2 - s2 .le. 0.0d0) then write(*,*) 'hyp_get_prim_var: physical variables unphysical!' write(*,*) 'undensitized energy: ',uen write(*,*) 'undensitized density: ',ude write(*,*) 'undensitized s2: ',s2 write(*,*) 'pressure guess: ',pold write(*,*) 'mahc_densmin: ',mahc_densmin write(*,*) 'coordinates: x,y,z,r:',x,y,z,r call CCTK_WARN(0,"get_prim_var_pt, pt8") endif c calculate rho and epsilon as a function of evolved vars and press rho = ude * sqrt( (uen + pold + ude)**2 - s2)/(uen + pold + ude) w_lorentz = (uen + pold + ude) / sqrt( (uen + pold + ude)**2 - s2) epsilon = (sqrt( (uen + pold + ude)**2 - s2) - pold * w_lorentz - . ude)/ude c calculate the function f = pold - (gamma - 1.0d0) * rho * epsilon count = 0 pnew = pold do while ( ((abs(pnew - pold)/abs(pnew) .gt. mahc_perc_ptol) .and. . (abs(pnew - pold) .gt. mahc_del_ptol)) .or. . (count .lt. 5)) count = count + 1 if (count .gt. countmax) then write(*,*) 'get_prim_var(press): error: did not converge in ' write(*,*) ' ',countmax,' steps' write(*,*) 'de is ',de write(*,*) 'sx is ',sx write(*,*) 'sy is ',sy write(*,*) 'sz is ',sz write(*,*) 'en is ',en write(*,*) 'press func: ',f write(*,*) 'deriv of press func: ',df write(*,*) 'pnew: ',pnew write(*,*) 'pold: ',pold write(*,*) 'rho: ',rho write(*,*) 'epsilon: ',epsilon write(*,*) 'w_lorentz: ',w_lorentz write(*,*) 'xyz location: ',x,y,z write(*,*) 'radius: ',r call CCTK_WARN(0,"get_prim_var_pt, pt9") endif df = 1.0d0 - (gamma - 1.0d0) * ( . (uen + ude - rho * epsilon) / (uen + pnew + ude) - . ude / sqrt( (uen + pnew + ude)**2 - s2)) pold = pnew pnew = pold - f/df if( (uen + pnew + ude)**2 - s2 .le. 0.0d0) then if(ude .lt. mahc_densmin) then s_reduct = ((uen + pnew + ude)/sqrt(s2))/2.0d0 sx = sx * s_reduct sy = sy * s_reduct sz = sz * s_reduct usx = usx * s_reduct usy = usy * s_reduct usz = usz * s_reduct s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.0d0*usx*usy*uxy + . 2.0d0*usx*usz*uxz + 2.0d0*usy*usz*uyz endif endif c calculate primatives and function again rho = ude * sqrt( (uen + pnew + ude)**2 - s2)/(uen + pnew + ude) w_lorentz = (uen + pnew + ude) / sqrt( (uen + pnew + ude)**2 - s2) epsilon = (sqrt( (uen + pnew + ude)**2 - s2) - pnew * w_lorentz - . ude)/ude f = pnew - (gamma - 1.0d0) * rho * epsilon enddo c do a couple Newtons just for polishing do i=1,2 df = 1.0d0 - (gamma - 1.0d0) * ( . (uen + ude - rho * epsilon) / (uen + pnew + ude) - . ude / sqrt( (uen + pnew + ude)**2 - s2)) pold = pnew pnew = pold - f/df if( (uen + pnew + ude)**2 - s2 .le. 0.0d0) then if(ude .lt. mahc_densmin) then s_reduct = ((uen + pnew + ude)/sqrt(s2))/2.0d0 sx = sx * s_reduct sy = sy * s_reduct sz = sz * s_reduct usx = usx * s_reduct usy = usy * s_reduct usz = usz * s_reduct s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.0d0*usx*usy*uxy + . 2.0d0*usx*usz*uxz + 2.0d0*usy*usz*uyz endif endif c calculate primatives and function again rho = ude * sqrt( (uen + pnew + ude)**2 - s2)/(uen + pnew + ude) w_lorentz = (uen + pnew + ude) / sqrt( (uen + pnew + ude)**2 - s2) epsilon = (sqrt( (uen + pnew + ude)**2 - s2) - pnew * w_lorentz - . ude)/ude f = pnew - (gamma - 1.0d0) * rho * epsilon enddo c calculate primatives press = pnew vlowx = usx / ( (rho + rho*epsilon + press) * w_lorentz**2) vlowy = usy / ( (rho + rho*epsilon + press) * w_lorentz**2) vlowz = usz / ( (rho + rho*epsilon + press) * w_lorentz**2) velx = uxx * vlowx + uxy * vlowy + uxz * vlowz vely = uxy * vlowx + uyy * vlowy + uyz * vlowz velz = uxz * vlowx + uyz * vlowy + uzz * vlowz c Check and make sure we actually solved the equations. This c can someday be commented out. de_prim = de - sqrt(det) * rho * w_lorentz en_prim = en - sqrt(det) * ((rho * (1.0d0 + epsilon) + press) * . w_lorentz**2 - press - rho * w_lorentz) sx_prim = sx-sqrt(det)*(rho*(1.0d0+epsilon) + . press)*vlowx*w_lorentz**2 sy_prim = sy-sqrt(det)*(rho*(1.0d0+epsilon) + . press)*vlowy*w_lorentz**2 sz_prim = sz-sqrt(det)*(rho*(1.0d0+epsilon) + . press)*vlowz*w_lorentz**2 test = max(abs(de_prim),abs(en_prim)) test = max(test,abs(sx_prim)) test = max(test,abs(sy_prim)) test = max(test,abs(sz_prim)) #ifndef SINGLE_PRECISION if(test/rho .ge. 1.0d-10) then write(*,*) 'hyp_get_prim_var_pt: 1: error: ' write(*,*) de_prim,en_prim,sx_prim,sy_prim,sz_prim write(*,*) 'en : ',en write(*,*) 'det: ',det write(*,*) 'rho: ',rho write(*,*) 'epsilon: ',epsilon write(*,*) 'press: ',press write(*,*) 'w_lorentz: ',w_lorentz call CCTK_WARN(0,"get_prim_var_pt, pt10") endif #endif c Now, check for eps. If it is lower than some prescribed c threshold then do the thing all over again, except for c substituting the energy equation with P = eos_k * rho**gamma. if(epsilon .le. mahc_eps_min) then c set initial guess for eps (just use the eps from previous calc) epsold = epsilon c calculate w_lorentz and rho as a function of evolved vars and eps w2mhalf = sqrt(0.25d0 + s2 / ((ude * (1.0d0 + epsold*gamma))**2)) w2 = 0.5d0 + w2mhalf w_lorentz = sqrt(w2) rho = ude/w_lorentz c calculate the function f = (gamma - 1.0d0) * rho * epsold - mahc_eos_k * (rho ** gamma) count = 0 epsnew = epsold do while ( ((abs(epsnew - epsold)/abs(epsnew+1.0d-30) .gt. 1.0e-11) .and. . (abs(epsnew - epsold) .gt. 1.0d-22)) .or. . (count .lt. 5)) count = count + 1 if (count .gt. countmax) then write(*,*) 'get_prim_var (eps solve (stage 2)):' write(*,*) 'get_prim_var(eps): error: did not converge in ' write(*,*) ' ',countmax,' steps' write(*,*) 'de is ',de write(*,*) 'sx is ',sx write(*,*) 'sy is ',sy write(*,*) 'sz is ',sz write(*,*) 'en is ',en write(*,*) 'eps func: ',f write(*,*) 'deriv of eps func: ',df write(*,*) 'epsnew: ',epsnew write(*,*) 'epsold: ',epsold write(*,*) 'rho: ',rho write(*,*) 'w_lorentz: ',w_lorentz write(*,*) 'xyz location: ',x,y,z write(*,*) 'radius: ',r call CCTK_WARN(0,"get_prim_var_pt, pt11") endif df = (gamma - 1.0d0)*rho + gamma * s2 * . ( (gamma - 1.0d0)*epsnew - . mahc_eos_k * gamma * rho**(gamma-1.0d0))/ . (ude * w2 * w2mhalf * (1.0d0 + epsnew*gamma)**3) epsold = epsnew epsnew = epsold - f/df c calculate primatives and function again w2mhalf = sqrt(0.25d0 + s2 / ((ude * (1.0d0 + epsnew*gamma))**2)) w2 = 0.5d0 + w2mhalf w_lorentz = sqrt(w2) rho = ude/w_lorentz f = (gamma - 1.0d0) * rho * epsnew - mahc_eos_k * (rho ** gamma) enddo c do a couple Newtons just for polishing do i=1,2 df = (gamma - 1.0d0)*rho + gamma * s2 * . ( (gamma - 1.0d0)*epsnew - . mahc_eos_k * gamma * rho**(gamma-1.0d0))/ . (ude * w2 * w2mhalf * (1.0d0 + epsnew*gamma)**3) epsold = epsnew epsnew = epsold - f/df c calculate primatives and function again w2mhalf = sqrt(0.25d0 + s2 / ((ude * (1.0d0 + epsnew*gamma))**2)) w2 = 0.5d0 + w2mhalf w_lorentz = sqrt(w2) rho = ude/w_lorentz f = (gamma - 1.0d0) * rho * epsnew - mahc_eos_k * (rho ** gamma) enddo c calculate primatives epsilon = epsnew press = (gamma - 1.0d0) * rho * epsilon vlowx = usx / ( (rho + rho*epsilon + press) * w_lorentz**2) vlowy = usy / ( (rho + rho*epsilon + press) * w_lorentz**2) vlowz = usz / ( (rho + rho*epsilon + press) * w_lorentz**2) velx = uxx * vlowx + uxy * vlowy + uxz * vlowz vely = uxy * vlowx + uyy * vlowy + uyz * vlowz velz = uxz * vlowx + uyz * vlowy + uzz * vlowz c actually change the updated variable tau (here, called en) c to keep things consistent! en = sqrt(det) * ( (rho*(1.0d0+epsilon) + press) * w_lorentz**2 - . press - rho*w_lorentz) c as a consistency check, calculate residual and yell c if there is a problem de_prim = de - sqrt(det) * rho * w_lorentz en_prim = (gamma - 1.0d0) * rho * epsilon - . mahc_eos_k * (rho ** gamma) sx_prim = sx-sqrt(det)*(rho*(1.0d0+epsilon) + . press)*vlowx*w_lorentz**2 sy_prim = sy-sqrt(det)*(rho*(1.0d0+epsilon) + . press)*vlowy*w_lorentz**2 sz_prim = sz-sqrt(det)*(rho*(1.0d0+epsilon) + . press)*vlowz*w_lorentz**2 test = max(abs(de_prim),abs(en_prim)) test = max(test,abs(sx_prim)) test = max(test,abs(sy_prim)) test = max(test,abs(sz_prim)) if(test .ge. 1.0d-10) then write(*,*) 'hyp_get_prim_var_pt: error: ' write(*,*) ' stage 2: ' write(*,*) de_prim,en_prim,sx_prim,sy_prim,sz_prim call CCTK_WARN(0,"get_prim_var_pt, pt12") endif endif if(epsilon .lt. 0.0d0) then write(*,*) 'hyp_get_prim_var_pt: stopping the code.' write(*,*) ' specific internal energy just went below 0! ' write(*,*) 'xyz location: ',x,y,z write(*,*) 'radius: ',r write(*,*) 'velocities: ',velx,vely,velz call CCTK_WARN(0,"get_prim_var_pt, pt13") endif #endif return end