Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 207)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:49 PM Modified: /trunk/ param.ccl Log: GRHydro_InitData: fix some text only typos From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +2 -2 =================================================================== --- trunk/param.ccl 2013-03-28 01:49:08 UTC (rev 206) +++ trunk/param.ccl 2013-03-28 01:49:10 UTC (rev 207) @@ -18,7 +18,7 @@ "simple_wave" :: "Set initial data from Anile Miller Motta, Phys.Fluids. 26, 1450 (1983)" "monopole" :: "Monopole at the center" "cylexp" :: "Cylindrical Explosion" - "rotor" :: "Magnetic Rotor test from DelZanna,Bucciantini, and Londrillo A&A 400, 397â413 (2003)" + "rotor" :: "Magnetic Rotor test from DelZanna,Bucciantini, and Londrillo A&A 400, 397-413 (2003)" "advectedloop":: "Magnetic advected loop test" "alfvenwave" :: "Circularly polarized Alfven wave" "hydro_bondi_solution" :: "Spherical single black hole Bondi solution" @@ -275,7 +275,7 @@ (0:* :: "Any positive number" } 0.8 -CCTK_REAL cyl_r_outer "Inner Radius" +CCTK_REAL cyl_r_outer "Outer Radius" { (0:* :: "Any positive number" } 1.0

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 206)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:49 PM Modified: /trunk/src/ GRHydro_BondiM.c Log: GRHydro_Init_Data: remove sincos from Bondi ID routine From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +5 -11 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:49:05 UTC (rev 205) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:49:08 UTC (rev 206) @@ -67,13 +67,12 @@ #define DLOOP1 for(i=0 ;i<NDIM ;i++) #define DLOOP2 for(i=0 ;i<NDIM ;i++) for(j=0 ;j<NDIM ;j++) +static inline void sincos(double theta, CCTK_REAL *sth, CCTK_REAL *cth) +{ + *sth = sin(theta); + *cth = cos(theta); +} -#if !defined(__INTEL_COMPILER) -# define LOCAL_SINCOS -# define sincos( theta_ , sth_ , cth_ ) { *(sth_) = sin((theta_)) ; *(cth_) = cos((theta_)) ; } -#endif - - //Newton-Raphson parameters: #define NEWT_DIM_B (1 ) #define MAX_NEWT_ITER_B (30 ) /* Max. # of Newton-Raphson iterations for find_root_2D(); */ @@ -1349,8 +1348,3 @@ #undef COORD_ISOTROPIC #undef DLOOP1 #undef DLOOP2 - -#ifdef LOCAL_SINCOS -# undef sincos -# undef LOCAL_SINCOS -#endif

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 205)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:49 PM Modified: /trunk/ param.ccl, schedule.ccl /trunk/src/ GRHydro_BondiM.c Log: GRHydro_Init_Data: add routine to force boundary data to Bondi From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +46 -8 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:49:03 UTC (rev 204) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:49:05 UTC (rev 205) @@ -890,22 +890,41 @@ -- driver routine for the Bondi solution; ***********************************************************************************/ -static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max); +static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, + CCTK_REAL range_min, CCTK_REAL range_max, + const int range_imin[], const int range_imax[]); void GRHydro_BondiM(CCTK_ARGUMENTS) { - GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100); + DECLARE_CCTK_ARGUMENTS; + const int null[3] = {0,0,0}; + GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100, cctk_lsh, null); } void GRHydro_BondiM_Range(CCTK_ARGUMENTS) { DECLARE_CCTK_PARAMETERS; + DECLARE_CCTK_ARGUMENTS; + const int null[3] = {0,0,0}; GRHydro_BondiM_Internal(CCTK_PASS_CTOC, bondi_freeze_inner_radius, - bondi_freeze_outer_radius); + bondi_freeze_outer_radius, + cctk_lsh, null); } +void GRHydro_BondiM_Boundary(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_PARAMETERS; + DECLARE_CCTK_ARGUMENTS; + const int imax[3] = {cctk_lsh[0]-cctk_nghostzones[0], + cctk_lsh[1]-cctk_nghostzones[1], + cctk_lsh[2]-cctk_nghostzones[2]}; + GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100, + cctk_nghostzones, imax); +} -static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max) +static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, + CCTK_REAL range_min, CCTK_REAL range_max, + const int range_imin[], const int range_imax[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; @@ -920,7 +939,7 @@ int nbondi; - const int size = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]; + const int size = cctk_ash[0] * cctk_ash[1] * cctk_ash[2]; int i, imin, j; # define velx(i_) vel[i_ + 0*size] @@ -952,7 +971,7 @@ } else if( CCTK_EQUALS(bondi_Bvec_method, "transform")) { bvec_method = BVEC_TRANSFORM; } else { - CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, "Unknown bvec setup method '%s'", bondi_Bvec_method); } @@ -1109,8 +1128,12 @@ *********************************************************************************/ xpos[TT] = 0.; - for(i=0; i < size; i++) { + for(int kk = 0 ; kk < cctk_lsh[2] ; kk++) { + for(int jj = 0 ; jj < cctk_lsh[1] ; jj++) { + for(int ii = 0 ; ii < cctk_lsh[0] ; ii++) { + i = CCTK_GFINDEX3D(cctkGH, ii,jj,kk); + xpos[XX] = x[i] ; xpos[YY] = y[i] ; xpos[ZZ] = z[i] ; @@ -1152,8 +1175,21 @@ rspher = x_spher[RR]; + // the conditions here are a bit strange. The logic goes like this: + // range_min,range_max (and range_imin, range_imax) are radii (and indices) + // bounding the volume where we want to evolve the fields. So during + // evolution we only reset to ID data if we are NOT inside this volume. For + // ID we pass special values for min/max such that these conditions never + // trigger. + // TODO: do this more nicely. The fact that a comment is needed should tell + // me that this is getting too complex. if(rspher > range_min && rspher < range_max) continue; + if(ii >= range_imin[0] && ii < range_imax[0] && + jj >= range_imin[1] && jj < range_imax[1] && + kk >= range_imin[2] && kk < range_imax[2]) { + continue; + } /* Find nearest point in the Bondi solution : */ j = (int) ( 0.5 + (log10(rspher) - logrmin) / dlogr ) ; @@ -1266,7 +1302,9 @@ if(psidc) psidc[i] = 0.; - } + } // ii + } // jj + } // kk free(r_bondi); free(logr_bondi); free(rho_bondi); free(u_bondi); free(v_bondi); Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +4 -0 =================================================================== --- trunk/param.ccl 2013-03-28 01:49:03 UTC (rev 204) +++ trunk/param.ccl 2013-03-28 01:49:05 UTC (rev 205) @@ -420,6 +420,10 @@ *:* :: "any value" } 1e300 +CCTK_BOOLEAN bondi_overwrite_boundary "reset data to initial data in outer boundary in boundary condition" +{ +} "no" + # For Poloidal Magnetic field test: CCTK_REAL poloidal_A_b "Vector potential strength" File [modified]: schedule.ccl Delta lines: +10 -0 =================================================================== --- trunk/schedule.ccl 2013-03-28 01:49:03 UTC (rev 204) +++ trunk/schedule.ccl 2013-03-28 01:49:05 UTC (rev 205) @@ -257,6 +257,16 @@ } "force analytic solution outside anulus" } +if(bondi_overwrite_boundary ) { +# this is a terrible HACK. We really should at least schedule ourselves in +# Boundaries and look at cctk_bbox and the symmetry flags +# by now Cactus is convenient CCTK_LOOP macros for this... +SCHEDULE GRHydro_BondiM_Boundary IN HydroBase_Boundaries BEFORE HydroBase_Select_Boundaries +{ + LANG: C +} "force analytic solution in boundaries" +} + if (CCTK_EQUALS(initial_Bvec, "poloidalmagfield") || CCTK_EQUALS(initial_Avec, "poloidalmagfield")) { # SCHEDULE GRHydro_PoloidalMagFieldM AT CCTK_INITIAL AFTER IN HydroBase_Initial AFTER rnsid_init AFTER TOV_Initial_Data after CCCC_StarMapper_InitialData

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 204)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:49 PM Modified: /trunk/ interface.ccl, param.ccl /trunk/src/ GRHydro_C2P2CM.F90 Log: Finish the implementation of con2prim using the entropy equation. * Tested successfully pointwisely. It still needs to be tested on evolution. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_C2P2CM.F90 Delta lines: +87 -2 =================================================================== --- trunk/src/GRHydro_C2P2CM.F90 2013-03-28 01:49:02 UTC (rev 203) +++ trunk/src/GRHydro_C2P2CM.F90 2013-03-28 01:49:03 UTC (rev 204) @@ -44,10 +44,12 @@ CCTK_REAL :: gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send CCTK_REAL :: dens_send,sx_send,sy_send,sz_send,tau_send CCTK_REAL :: bconsx_send, bconsy_send, bconsz_send + CCTK_REAL :: entropy_send, entropycons_send; CCTK_REAL :: rho_send,velx_send,vely_send,velz_send,eps_send CCTK_REAL :: press_send,w_lorentz_send CCTK_REAL :: bvcx_send, bvcy_send, bvcz_send, b2_send CCTK_REAL :: C2P_failed + CCTK_REAL :: BdotB, Bdotv, b2, beta_mag CCTK_INT :: epsnegative ! begin EOS Omni vars @@ -119,6 +121,7 @@ ) w_lorentz_send = 1.0d0/w_lorentz_send + epsnegative = 0 xrho = 1.0d-10 @@ -136,8 +139,46 @@ xrho,epsval,xtemp,xye,local_gam,keyerr,anyerr) local_gam = local_gam + 1.0 + if(use_c2p_with_entropy_eqn.eq.1)then + entropy_send = ( local_gam(1) - 1.0d0 ) * eps_send * & + rho_send**(2.0d0 - local_gam(1)) + entropycons_send = sdet * w_lorentz_send * entropy_send + endif + C2P_failed = 0.0d0 + write(*,*) 'C2P2CM test: metric values.' + write(*,*) ' gxx: ', gxx_send + write(*,*) ' gxy: ', gxy_send + write(*,*) ' gxz: ', gxz_send + write(*,*) ' gyy: ', gyy_send + write(*,*) ' gyz: ', gyz_send + write(*,*) ' gzz: ', gzz_send + write(*,*) ' uxx: ', uxx + write(*,*) ' uxy: ', uxy + write(*,*) ' uxz: ', uxz + write(*,*) ' uyy: ', uyy + write(*,*) ' uyz: ', uyz + write(*,*) ' uzz: ', uzz + write(*,*) ' det: ', det + + write(*,*) 'C2P2CM test: initial primitive guess values.' + write(*,*) ' primitive variables: ' + write(*,*) ' rho : ',rho_send + write(*,*) ' velx : ',velx_send + write(*,*) ' vely : ',vely_send + write(*,*) ' velz : ',velz_send + write(*,*) ' press : ',press_send + write(*,*) ' eps : ',eps_send + if(use_c2p_with_entropy_eqn.eq.1)then + write(*,*) ' entropy : ',entropy_send + endif + write(*,*) ' W : ',w_lorentz_send + write(*,*) ' Bvecx : ',bvcx_send + write(*,*) ' Bvecy : ',bvcy_send + write(*,*) ' Bvecz : ',bvcz_send + write(*,*) ' C2P_failed : ',C2P_failed + write(*,*) 'C2P2CM test: initial values.' write(*,*) ' conservative variables: ' write(*,*) ' dens: ',dens_send @@ -148,19 +189,57 @@ write(*,*) ' Bconsx : ',bconsx_send write(*,*) ' Bconsy : ',bconsy_send write(*,*) ' Bconsz : ',bconsz_send + if(use_c2p_with_entropy_eqn.eq.1)then + write(*,*) ' entropycons : ',entropycons_send + endif write(*,*) ' eps : ',eps_send write(*,*) ' W : ',w_lorentz_send write(*,*) ' Bvecx : ',bvcx_send write(*,*) ' Bvecy : ',bvcy_send write(*,*) ' Bvecz : ',bvcz_send + +! Is this point magnetically dominated? + BdotB = gxx_send*bvcx_send**2+gyy_send*bvcy_send**2+gzz_send*bvcz_send**2& + +2.0*(gxy_send*bvcx_send*bvcy_send+gxz_send*bvcx_send*bvcy_send+ & + gyz_send*bvcy_send*bvcz_send) + Bdotv = gxx_send*bvcx_send*velx_send+gyy_send*bvcy_send*vely_send & + +gzz_send*bvcz_send*velz_send & + +gxy_send*(bvcx_send*vely_send+bvcy_send*velx_send) & + +gxz_send*(bvcx_send*velz_send+bvcz_send*velx_send) & + +gyz_send*(bvcy_send*velz_send+bvcz_send*vely_send) + + b2 = BdotB/w_lorentz_send**2 + Bdotv**2 + + beta_mag = 2.0*press_send/b2; + + write(*,*) ' BdotB : ',BdotB + write(*,*) ' Bdotv : ',Bdotv + write(*,*) ' b2 : ',b2 + write(*,*) 'beta_mag : ',beta_mag - write(*,*) 'C2P2CM test: getting the associated primitive variables.' - call Con2PrimGenM(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,local_gam(1),dens_send,sx_send,sy_send,sz_send, & + if(use_c2p_with_entropy_eqn.eq.0)then + write(*,*) 'C2P2CM test: getting the associated primitive variables.' + call Con2PrimGenM(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,local_gam(1),dens_send,sx_send,sy_send,sz_send, & tau_send,bconsx_send,bconsy_send,bconsz_send,xtemp(1),xye(1),rho_send,velx_send,vely_send,velz_send, & eps_send,press_send,bvcx_send,bvcy_send,bvcz_send,b2_send,w_lorentz_send, & gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,& uxx,uxy,uxz,uyy,uyz,uzz,det,& epsnegative,C2P_failed) + else + write(*,*) 'C2P2CMee test: getting the associated primitive variables.' + call Con2PrimGenMee(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,local_gam(1),dens_send,sx_send,sy_send,sz_send, & + tau_send,bconsx_send,bconsy_send,bconsz_send,entropycons_send,xtemp(1),xye(1),rho_send,velx_send,vely_send,velz_send, & + eps_send,press_send,bvcx_send,bvcy_send,bvcz_send,b2_send,w_lorentz_send, & + gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,& + uxx,uxy,uxz,uyy,uyz,uzz,det,& + epsnegative,C2P_failed) + endif + + if(use_c2p_with_entropy_eqn.eq.1)then + entropy_send = ( local_gam(1) - 1.0d0 ) * eps_send * & + rho_send**(2.0d0 - local_gam(1)) + entropycons_send = sdet * w_lorentz_send * entropy_send + endif write(*,*) 'C2P2CM test: the primitive variables are' write(*,*) ' primitive variables: ' @@ -170,6 +249,9 @@ write(*,*) ' velz : ',velz_send write(*,*) ' press : ',press_send write(*,*) ' eps : ',eps_send + if(use_c2p_with_entropy_eqn.eq.1)then + write(*,*) ' entropy : ',entropy_send + endif write(*,*) ' W : ',w_lorentz_send write(*,*) ' Bvecx : ',bvcx_send write(*,*) ' Bvecy : ',bvcy_send @@ -191,6 +273,9 @@ write(*,*) ' Bconsx : ',bconsx_send write(*,*) ' Bconsy : ',bconsy_send write(*,*) ' Bconsz : ',bconsz_send + if(use_c2p_with_entropy_eqn.eq.1)then + write(*,*) ' entropycons : ',entropycons_send + endif write(*,*) ' eps : ',eps_send write(*,*) ' W : ',w_lorentz_send write(*,*) ' Bvecx : ',bvcx_send Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +22 -0 =================================================================== --- trunk/interface.ccl 2013-03-28 01:49:02 UTC (rev 203) +++ trunk/interface.ccl 2013-03-28 01:49:03 UTC (rev 204) @@ -131,6 +131,27 @@ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) +void FUNCTION Con2PrimGenMee(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 entropycons, \ + 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, \ + CCTK_REAL INOUT Bvecsq, \ + CCTK_REAL INOUT w_lorentz, \ + CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ + CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ + CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ + CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ + CCTK_REAL INOUT det, \ + CCTK_INT OUT epsnegative, \ + CCTK_REAL OUT retval) + void FUNCTION Con2PrimPolyM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT sc, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ @@ -157,6 +178,7 @@ USES FUNCTION Con2PrimPoly USES FUNCTION Con2PrimGen USES FUNCTION Con2PrimGenM +USES FUNCTION Con2PrimGenMee USES FUNCTION Con2PrimPolyM CCTK_INT FUNCTION EOS_Omni_GetHandle(CCTK_STRING IN name) File [modified]: param.ccl Delta lines: +5 -0 =================================================================== --- trunk/param.ccl 2013-03-28 01:49:02 UTC (rev 203) +++ trunk/param.ccl 2013-03-28 01:49:03 UTC (rev 204) @@ -188,6 +188,11 @@ } 6.666666666666667d-7 # Initial conservative values used by c2p tests: + +CCTK_BOOLEAN use_c2p_with_entropy_eqn "Use the con2prim routine that uses the entropy equation instead of the energy equation" +{ +} no + REAL dens_init "Initial conserved mass density" { (0:* :: "Anything positive."

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 203)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:49 PM Modified: /trunk/ param.ccl Log: GRHydro_InitData: Allow n_p=0 again for poloidal magnetic field From: Tanja Bode <tanja.bode(a)physics.gatech.edu> File Changes: Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +1 -1 =================================================================== --- trunk/param.ccl 2013-03-28 01:49:00 UTC (rev 202) +++ trunk/param.ccl 2013-03-28 01:49:02 UTC (rev 203) @@ -424,7 +424,7 @@ CCTK_INT poloidal_n_p "Vector potential strength" { - (0:* :: "Any positive integer." + 0:* :: "Any positive integer." } 3 CCTK_REAL poloidal_P_cut "Pressure used to confine the B field inside a star"

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 202)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:49 PM Modified: /trunk/src/ GRHydro_BondiM.c Log: GRHydro_Init_Data: damp hydro quantities to atmosphere inside of rKS < 1 M From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +13 -0 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:59 UTC (rev 201) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:49:00 UTC (rev 202) @@ -1238,6 +1238,19 @@ Bvecz(i) = bondi_bmag*SQR(M)*z[i]/sqrt(det)/CUBE(r[i]); } + // damp everything down to atmosphere inside of 1M + if(rspher < M) { + const CCTK_REAL smooth = 0.5*(1+tanh(tan(M_PI*(rspher/M-0.5)))); + rho[i] *= smooth; + velx(i) *= smooth; + vely(i) *= smooth; + velz(i) *= smooth; + eps[i] *= smooth; + Bvecx(i) *= smooth; + Bvecy(i) *= smooth; + Bvecz(i) *= smooth; + } + Prim2ConGenM(*GRHydro_eos_handle,gxx[i],gxy[i], gxz[i],gyy[i],gyz[i],gzz[i], det, &dens[i],&sx(i),&sy(i),&sz(i),

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 201)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/src/ GRHydro_BondiM.c Log: GRHydro_Init_Data: compute w_lorentz from uT and alp instead of vel later on prim2con will of course compute the grid function using vel From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +4 -12 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:57 UTC (rev 200) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:59 UTC (rev 201) @@ -1225,21 +1225,13 @@ vely(i) = (ucon[YY]/ucon[TT] + betay[i]) * inv_alp; velz(i) = (ucon[ZZ]/ucon[TT] + betaz[i]) * inv_alp; - w_lorentz[i] = 1./sqrt(1-(gxx[i]*SQR(velx(i))+gyy[i]*SQR(vely(i))+gzz[i]*SQR(velz(i))+ - 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i))) ); - if(isnan(w_lorentz[i])) { - CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING, - "nan in w_lorentz at (%g,%g,%g): v2 = %g", - x[i],y[i],z[i],(gxx[i]*SQR(velx(i))+gyy[i]*SQR(vely(i))+gzz[i]*SQR(velz(i))+ - 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i)))); - } - SpatialDet(gxx[i],gxy[i],gxz[i],gyy[i],gyz[i],gzz[i],&det); if(bvec_method == BVEC_TRANSFORM) { - Bvecx(i) = w_lorentz[i]*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; - Bvecy(i) = w_lorentz[i]*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; - Bvecz(i) = w_lorentz[i]*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; + const CCTK_REAL w = ucon[TT] * alp[i]; + Bvecx(i) = w*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; + Bvecy(i) = w*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; + Bvecz(i) = w*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; } else { Bvecx(i) = bondi_bmag*SQR(M)*x[i]/sqrt(det)/CUBE(r[i]); Bvecy(i) = bondi_bmag*SQR(M)*y[i]/sqrt(det)/CUBE(r[i]);

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 200)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/src/ GRHydro_BondiM.c Log: GRHydro_Init_Data: remove rho_guess from find_bondi_solution was used to call the newton solver twice with the very same arguments. Not sure what the idea was. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +1 -12 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:56 UTC (rev 199) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:57 UTC (rev 200) @@ -647,7 +647,7 @@ const int ntries = 10000; int itry; - CCTK_REAL rhotmp, rho_guess; + CCTK_REAL rhotmp; CCTK_REAL dr,ur; @@ -662,17 +662,12 @@ } else { // rhotmp = (sqrt(Qdot) - 1.) * (gamma_eos - 1.) / ( gamma_eos * K ); - // rho_guess = pow( rhotmp , (1./(gamma_eos - 1.)) ); if(r < rs) { ur = pow(r,-0.5) ; } else { ur = 0.5*pow(r,-1.5) ; } *rho = Mdot / (4.*M_PI * r * r * ur); } } - // safe guess value for multiple tries - rho_guess = *rho; - - // set global variables needed by residual function: r_sol = r ; @@ -681,12 +676,6 @@ retval = gnr_bondi( rho, NEWT_DIM_B, bondi_resid); - // first try guess if failure - if( retval ) { - *rho = rho_guess; - retval = gnr_bondi( rho, NEWT_DIM_B, bondi_resid); - } - // If we were unsure about the guess and solver fails, then creep from known solution to desired point: if( retval ) {

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 199)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/src/ GRHydro_BondiM.c Log: GRHydro_InitDat: clear psidc and c2p_failed in ID routine this is mostly important when we freeze the evolution afterwards From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +4 -0 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:55 UTC (rev 198) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:56 UTC (rev 199) @@ -1267,6 +1267,10 @@ eps[i],&press[i], Bvecx(i),Bvecy(i),Bvecz(i), &w_lorentz[i]); + + GRHydro_C2P_failed[i] = 0.; + if(psidc) + psidc[i] = 0.; }

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 198)
by rhaas＠tapir.caltech.edu 27 Mar '13

27 Mar '13

User: rhaas Date: 2013/03/27 08:48 PM Modified: /trunk/ param.ccl /trunk/src/ GRHydro_BondiM.c Log: GRHydro_InitData: set Bvec according to UIUC formula From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM.c Delta lines: +46 -26 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:53 UTC (rev 197) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:55 UTC (rev 198) @@ -901,6 +901,8 @@ -- driver routine for the Bondi solution; ***********************************************************************************/ +static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max); + void GRHydro_BondiM(CCTK_ARGUMENTS) { GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100); @@ -914,17 +916,18 @@ } -void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max) +static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - CCTK_REAL det; + CCTK_REAL det, inv_alp; CCTK_REAL rhotmp, utmp, vtmp, rspher, xpos[NDIM], x_spher[NDIM], ucon[NDIM], bcon[NDIM]; int retval; CCTK_REAL *r_bondi, *logr_bondi, *rho_bondi, *u_bondi, *v_bondi; CCTK_REAL dlogr,logrmin; CCTK_REAL rmin_bondi,rmax_bondi; -// CCTK_INT GRHydro_reflevel=0; + int bvec_method; + enum {BVEC_DIRECT, BVEC_TRANSFORM}; int nbondi; @@ -946,14 +949,26 @@ if( CCTK_EQUALS(bondi_coordinates,"Boyer-Lindquist") ) { coord_type = COORD_BOYERLINDQUIST ; - } - if( CCTK_EQUALS(bondi_coordinates,"Kerr-Schild") ) { + } else if( CCTK_EQUALS(bondi_coordinates,"Kerr-Schild") ) { coord_type = COORD_KERRSCHILD ; - } - if( CCTK_EQUALS(bondi_coordinates,"Isotropic") ) { + } else if( CCTK_EQUALS(bondi_coordinates,"Isotropic") ) { coord_type = COORD_ISOTROPIC ; + } else { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Unknown coordinate type '%s'", bondi_coordinates); } + if( CCTK_EQUALS(bondi_Bvec_method, "direct")) { + bvec_method = BVEC_DIRECT; + } else if( CCTK_EQUALS(bondi_Bvec_method, "transform")) { + bvec_method = BVEC_TRANSFORM; + } else { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Unknown bvec setup method '%s'", bondi_Bvec_method); + } + + + /* xyz location of the black hole : */ i = 0; x_cen = bh_bondi_pos_x[i] ; @@ -1215,11 +1230,11 @@ } } - det = 1./alp[i]; /* temp var */ + inv_alp = 1./alp[i]; - velx(i) = (ucon[XX]/ucon[TT] + betax[i]) * det; - vely(i) = (ucon[YY]/ucon[TT] + betay[i]) * det; - velz(i) = (ucon[ZZ]/ucon[TT] + betaz[i]) * det; + velx(i) = (ucon[XX]/ucon[TT] + betax[i]) * inv_alp; + vely(i) = (ucon[YY]/ucon[TT] + betay[i]) * inv_alp; + velz(i) = (ucon[ZZ]/ucon[TT] + betaz[i]) * inv_alp; w_lorentz[i] = 1./sqrt(1-(gxx[i]*SQR(velx(i))+gyy[i]*SQR(vely(i))+gzz[i]*SQR(velz(i))+ 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i))) ); @@ -1230,24 +1245,29 @@ 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i)))); } - - Bvecx(i) = w_lorentz[i]*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; - Bvecy(i) = w_lorentz[i]*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; - Bvecz(i) = w_lorentz[i]*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; - SpatialDet(gxx[i],gxy[i],gxz[i],gyy[i],gyz[i],gzz[i],&det); - Prim2ConGenM(*GRHydro_eos_handle,gxx[i],gxy[i], - gxz[i],gyy[i],gyz[i],gzz[i], - det, &dens[i],&sx(i),&sy(i),&sz(i), - &tau[i], - &Bconsx(i),&Bconsy(i),&Bconsz(i), - rho[i], - velx(i),vely(i),velz(i), - eps[i],&press[i], - Bvecx(i),Bvecy(i),Bvecz(i), - &w_lorentz[i]); + if(bvec_method == BVEC_TRANSFORM) { + Bvecx(i) = w_lorentz[i]*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; + Bvecy(i) = w_lorentz[i]*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; + Bvecz(i) = w_lorentz[i]*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; + } else { + Bvecx(i) = bondi_bmag*SQR(M)*x[i]/sqrt(det)/CUBE(r[i]); + Bvecy(i) = bondi_bmag*SQR(M)*y[i]/sqrt(det)/CUBE(r[i]); + Bvecz(i) = bondi_bmag*SQR(M)*z[i]/sqrt(det)/CUBE(r[i]); + } + Prim2ConGenM(*GRHydro_eos_handle,gxx[i],gxy[i], + gxz[i],gyy[i],gyz[i],gzz[i], + det, &dens[i],&sx(i),&sy(i),&sz(i), + &tau[i], + &Bconsx(i),&Bconsy(i),&Bconsz(i), + rho[i], + velx(i),vely(i),velz(i), + eps[i],&press[i], + Bvecx(i),Bvecy(i),Bvecz(i), + &w_lorentz[i]); + } free(r_bondi); free(logr_bondi); free(rho_bondi); free(u_bondi); free(v_bondi); Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +6 -0 =================================================================== --- trunk/param.ccl 2013-03-28 01:48:53 UTC (rev 197) +++ trunk/param.ccl 2013-03-28 01:48:55 UTC (rev 198) @@ -395,6 +395,12 @@ (0:* :: "positive" } 1.0 +CCTK_KEYWORD bondi_Bvec_method "how to compute the magnetic field vector" +{ + "direct" :: "directly from Cartesian metric" + "transform" :: "transform Schwarzschild solution to Kerr Schild" +} "direct" + CCTK_BOOLEAN bondi_evolve_only_annulus "reset to initial data outside of bondi_freeze_inner_radius and bondi_freeze_outer_radius" { } "no"

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits