Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] TwoPunctures/trunk/src/ (Rev. 118)
by knarf＠cct.lsu.edu 10 Jan '11

10 Jan '11

User: knarf Date: 2011/01/10 01:20 PM Modified: /trunk/src/ TwoPunctures.c Log: patch from Bernard Kelly: Output the individual puncture ADM mass post-solve in TwoPunctures File Changes: Directory: /trunk/src/ ====================== File [modified]: TwoPunctures.c Delta lines: +12 -1 =================================================================== --- trunk/src/TwoPunctures.c 2010-11-21 03:29:55 UTC (rev 117) +++ trunk/src/TwoPunctures.c 2011-01-10 19:20:14 UTC (rev 118) @@ -219,6 +219,7 @@ CCTK_REAL admMass; if (! F) { + CCTK_REAL Mp_adm, Mm_adm, up, um; /* Solve only when called for the first time */ F = dvector (0, ntotal - 1); allocate_derivs (&u, ntotal); @@ -255,7 +256,7 @@ target ADM masses target_M_plus and target_M_minus and with initial guesses given by par_m_plus and par_m_minus. */ if(!(give_bare_mass)) { - CCTK_REAL tmp, Mp_adm, Mm_adm, Mp_adm_err, Mm_adm_err, up, um; + CCTK_REAL tmp, Mp_adm_err, Mm_adm_err; char valbuf[100]; CCTK_REAL M_p = target_M_plus; @@ -316,6 +317,16 @@ "The two puncture masses are mp=%.17g and mm=%.17g", (double) *mp, (double) *mm); + up = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v, par_b, 0., 0.); + um = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v,-par_b, 0., 0.); + + /* Calculate the ADM masses from the current bare mass guess */ + Mp_adm = (1 + up) * *mp + *mp * *mm / (4. * par_b); + Mm_adm = (1 + um) * *mm + *mp * *mm / (4. * par_b); + + CCTK_VInfo (CCTK_THORNSTRING, "Puncture 1 ADM mass is %g", (double)Mp_adm); + CCTK_VInfo (CCTK_THORNSTRING, "Puncture 2 ADM mass is %g", (double)Mm_adm); + /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */ admMass = (*mp + *mm - 4*par_b*PunctEvalAtArbitPosition(v.d0, 0, 1, 0, 0, nvar, n1, n2, n3));

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 210)
by knarf＠cct.lsu.edu 07 Jan '11

07 Jan '11

User: knarf Date: 2011/01/07 09:59 PM Modified: /trunk/src/ GRHydro_Prim2Con.F90 Log: correct variable declaration File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2011-01-08 01:37:53 UTC (rev 209) +++ trunk/src/GRHydro_Prim2Con.F90 2011-01-08 03:59:58 UTC (rev 210) @@ -46,7 +46,7 @@ integer :: i, j, k CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr - CCTK_REAL :: xtemp + CCTK_REAL :: xtemp(1) if(evolve_temper.ne.1) then !$OMP PARALLEL DO PRIVATE(i, j, avg_detl, avg_detr,&

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 209)
by cott＠tapir.caltech.edu 07 Jan '11

07 Jan '11

User: cott Date: 2011/01/07 07:37 PM Modified: /trunk/src/ GRHydro_Reconstruct.F90 Log: * fix a very nasty bug in PPM reconstruction that sneeked in with the very recent PPM changes. Old version of the code were not affected. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Reconstruct.F90 Delta lines: +3 -3 =================================================================== --- trunk/src/GRHydro_Reconstruct.F90 2011-01-06 01:13:02 UTC (rev 208) +++ trunk/src/GRHydro_Reconstruct.F90 2011-01-08 01:37:53 UTC (rev 209) @@ -361,7 +361,7 @@ if(clean_divergence.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),psidc(j,:,k),eps(j,:,k),press(j,:,k),& @@ -389,7 +389,7 @@ else !clean_divergence !$OMP PARALLEL DO PRIVATE(i, j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),dum(j,:,k),eps(j,:,k),press(j,:,k),& @@ -418,7 +418,7 @@ else !evolve_mhd !$OMP PARALLEL DO PRIVATE(i, j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 call SimplePPM_1d(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),&

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 208)
by bcmsma＠astro.rit.edu 05 Jan '11

05 Jan '11

User: bmundim Date: 2011/01/05 07:13 PM Modified: /trunk/src/ GRHydro_Con2PrimM.F90, GRHydro_Con2PrimM_polytype_pt.c, GRHydro_Con2PrimM_pt.c Log: RIT MHD development: Fix race conditions arising in the Con2Prim routines. Now both GNU and Intel compilers agree on the results for the magnetic tests. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +7 -5 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2011-01-04 14:40:06 UTC (rev 207) +++ trunk/src/GRHydro_Con2PrimM.F90 2011-01-06 01:13:02 UTC (rev 208) @@ -94,7 +94,7 @@ !$OMP PARALLEL DO PRIVATE(i,j,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & - !$OMP b2,local_gam,xrho,xpress,local_K,local_pgam,sc) + !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr) do k = 1, nz do j = 1, ny do i = 1, nx @@ -142,6 +142,7 @@ scon(i,j,k,:) = 0.d0 vel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 + xtemp = 0.0d0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) @@ -225,13 +226,14 @@ call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype.') !$OMP END CRITICAL - xrho=10.0d0 + xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& - 1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr) - local_K = xpress + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress; + xrho=10.0d0; xeps=1.0d0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& - xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr) + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) local_pgam=log(xpress/local_K)/log(xrho) sc = local_K*dens(i,j,k) File [modified]: GRHydro_Con2PrimM_polytype_pt.c Delta lines: +83 -166 =================================================================== --- trunk/src/GRHydro_Con2PrimM_polytype_pt.c 2011-01-04 14:40:06 UTC (rev 207) +++ trunk/src/GRHydro_Con2PrimM_polytype_pt.c 2011-01-06 01:13:02 UTC (rev 208) @@ -71,43 +71,33 @@ ***************************************************/ /* Local Globals */ -static CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq,rho_gm1; -static CCTK_REAL W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old, drho_dW; -static CCTK_REAL t2,t4,t7,t24,two_Bsq,t300,t400,s200,s100; +struct LocGlob { + CCTK_REAL Bsq,QdotBsq,Qtsq,D; + CCTK_REAL W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old; + CCTK_REAL t4,t7,t24,two_Bsq,t300,t400,s200,s100; +} ; // Declarations: -static CCTK_REAL vsq_calc(CCTK_REAL W); - -static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, +static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][1], - CCTK_REAL *, CCTK_REAL *, CCTK_INT, CCTK_REAL) ); + CCTK_REAL *, CCTK_REAL *, CCTK_INT, CCTK_REAL, + struct LocGlob *lgp) ); static void func_W( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][1], - CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos); -static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) ; + CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob * lgp); -static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, +static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT, CCTK_REAL) ); + CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT, CCTK_REAL, + struct LocGlob *lgp) ); + static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos); + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob *lgp); -static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos); -/* pressure as a function of rho0 and u */ -static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u, CCTK_REAL gammaeos) -{ - return((gammaeos - 1.)*u) ; -} - -/* Pressure as a function of rho0 and w = rho0 + u + p */ -static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w, CCTK_REAL gammaeos) -{ - return((gammaeos-1.)*(w - rho0)/gammaeos) ; -} - - void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( CCTK_INT *handle, CCTK_REAL *gamma_eos, CCTK_REAL *dens_in, @@ -220,8 +210,12 @@ CCTK_REAL detg = (*det); CCTK_REAL sqrt_detg = sqrt(detg); CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; + CCTK_REAL t2, rho_gm1; CCTK_INT i,j, i_increase,ntries ; + struct LocGlob lg; + + /* Assume ok initially: */ *retval = 0.; *epsnegative = 0; @@ -275,7 +269,7 @@ // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: - Bsq = + lg.Bsq = (*gxx) * (*Bx) * (*Bx) + (*gyy) * (*By) * (*By) + (*gzz) * (*Bz) * (*Bz) + @@ -285,25 +279,23 @@ (*gyz) * (*By) * (*Bz) ); QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; - QdotBsq = QdotB*QdotB ; + lg.QdotBsq = QdotB*QdotB ; - Qtsq = (usx * sx + usy * sy + usz * sz) ; + lg.Qtsq = (usx * sx + usy * sy + usz * sz) ; - D = dens; + lg.D = dens; - half_Bsq = 0.5*Bsq; + t2 = lg.D*lg.D; + lg.t4 = lg.QdotBsq*t2; + lg.t7 = lg.Bsq*lg.Bsq; + lg.t24 = 1/t2; + lg.two_Bsq = lg.Bsq + lg.Bsq; + lg.t300 = lg.QdotBsq*lg.Bsq*t2; + lg.t400 = lg.Qtsq*t2; - t2 = D*D; - t4 = QdotBsq*t2; - t7 = Bsq*Bsq; - t24 = 1/t2; - two_Bsq = Bsq + Bsq; - t300 = QdotBsq*Bsq*t2; - t400 = Qtsq*t2; - - s200 = D*gammaeos*(*sc_in); + lg.s200 = lg.D*gammaeos*(*sc_in); CCTK_REAL g_o_gm1 = (gammaeos/(gammaeos-1.)); - s100 = g_o_gm1*(*sc_in); + lg.s100 = g_o_gm1*(*sc_in); /* calculate W from last timestep and use for guess */ vsq = @@ -328,7 +320,7 @@ // Always calculate rho from D and gamma so that using D in EOS remains consistent // i.e. you don't get positive values for dP/d(vsq) . - rho0 = D / gamma ; + rho0 = lg.D / gamma ; CCTK_REAL gm1 = gammaeos-1.; rho_gm1 = pow(rho0,gm1); p = (*sc_in) * rho_gm1 / gamma; @@ -340,20 +332,20 @@ // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; - while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) - - QdotBsq*(2.*W_last + Bsq) ) <= W_last*W_last*(Qtsq-Bsq*Bsq)) + while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) + - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) && (i_increase < 10) ) { W_last *= 10.; i_increase++; } - W_for_gnr2 = W_for_gnr2_old = W_last; - rho_for_gnr2 = rho_for_gnr2_old = rho0; + lg.W_for_gnr2 = lg.W_for_gnr2_old = W_last; + lg.rho_for_gnr2 = lg.rho_for_gnr2_old = rho0; // Calculate W: x_1d[0] = W_last; - *retval = 1.0*oned_newton_raphson( x_1d, 1, gammaeos, func_W ) ; + *retval = 1.0*oned_newton_raphson( x_1d, 1, gammaeos, &lg, func_W ) ; W = x_1d[0]; @@ -369,15 +361,15 @@ } } - rho_g = x_rho[0] = rho_for_gnr2; + rho_g = x_rho[0] = lg.rho_for_gnr2; ntries = 0; - while ( (*retval = gnr2( x_rho, 1, gammaeos, func_rho)) && ( ntries++ < 10 ) ) { + while ( (*retval = gnr2( x_rho, 1, gammaeos, &lg, func_rho)) && ( ntries++ < 10 ) ) { rho_g *= 10.; x_rho[0] = rho_g; } - rho_for_gnr2 = x_rho[0]; + lg.rho_for_gnr2 = x_rho[0]; if( (*retval != 0) ) { *retval = 10; @@ -385,10 +377,10 @@ } // Calculate v^2 : - rho0 = rho_for_gnr2; + rho0 = lg.rho_for_gnr2; rho_gm1 = pow(rho0,gm1); - utsq = (D-rho0)*(D+rho0)/(rho0*rho0); + utsq = (lg.D-rho0)*(lg.D+rho0)/(rho0*rho0); gamma_sq = 1.+utsq; gamma = sqrt(gamma_sq); @@ -421,9 +413,9 @@ *w_lorentz = gamma; *pressure = p ; - g_o_WBsq = 1./(W+Bsq); + g_o_WBsq = 1./(W+lg.Bsq); QdB_o_W = QdotB / W; - *bsq = Bsq * (1.-vsq) + QdB_o_W*QdB_o_W; + *bsq = lg.Bsq * (1.-vsq) + QdB_o_W*QdB_o_W; *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ; *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; @@ -449,49 +441,8 @@ } -/**********************************************************************/ -/**************************************************************************** - vsq_calc(): - - -- evaluate v^2 (spatial, normalized velocity) from - W = \gamma^2 w - -****************************************************************************/ -static CCTK_REAL vsq_calc(CCTK_REAL W) -{ - CCTK_REAL Wsq,Xsq,Bsq_W; - - Wsq = W*W ; - Bsq_W = (Bsq + W); - Xsq = Bsq_W * Bsq_W; - - return( ( Wsq * Qtsq + QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); -} - - /******************************************************************** - x1_of_x0(): - - -- calculates v^2 from W with some physical bounds checking; - -- asumes x0 is already physical - -- makes v^2 physical if not; - -*********************************************************************/ - -static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) -{ - CCTK_REAL x1,vsq; - CCTK_REAL dv = 1.e-15; - - vsq = fabs(vsq_calc(x0)) ; // guaranteed to be positive - - return( ( vsq > 1. ) ? (1.0 - dv) : vsq ); - -} - -/******************************************************************** - validate_x(): -- makes sure that x[0,1] have physical values, based upon @@ -524,10 +475,10 @@ -- inspired in part by Num. Rec.'s routine newt(); *****************************************************************/ -static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, +static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][1], CCTK_REAL *, - CCTK_REAL *, CCTK_INT, CCTK_REAL) ) + CCTK_REAL *, CCTK_INT, CCTK_REAL, struct LocGlob *) ) { CCTK_REAL f, df, dx[1], x_old[1], resid[1], jac[1][1]; @@ -554,7 +505,7 @@ keep_iterating = 1; while( keep_iterating ) { - (*funcd) (x, dx, resid, jac, &f, &df, n, gammaeos); /* returns with new dx, f, df */ + (*funcd) (x, dx, resid, jac, &f, &df, n, gammaeos, lgp); /* returns with new dx, f, df */ /* Save old values before calculating the new: */ errx = 0.; @@ -568,8 +519,8 @@ // //METHOD specific: // i_increase = 0; -// while( (( x[0]*x[0]*x[0] * ( x[0] + 2.*Bsq ) - -// QdotBsq*(2.*x[0] + Bsq) ) <= x[0]*x[0]*(Qtsq-Bsq*Bsq)) +// while( (( x[0]*x[0]*x[0] * ( x[0] + 2.*lgp->Bsq ) - +// lgp->QdotBsq*(2.*x[0] + lgp->Bsq) ) <= x[0]*x[0]*(lgp->Qtsq-lgp->Bsq*lgp->Bsq)) // && (i_increase < 10) ) { // x[0] -= (1.*i_increase) * dx[0] / 10. ; // i_increase++; @@ -615,7 +566,7 @@ if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { #if(DEBUG_CON2PRIMM) fprintf(stderr,"\ngnr not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n", - f,df,x[0],x_old[0],W_for_gnr2_old,W_for_gnr2,rho_for_gnr2_old,rho_for_gnr2); fflush(stderr); + f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); #endif return(2); } @@ -652,52 +603,51 @@ n = dimension of x[]; *********************************************************************************/ static void func_W(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos) + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob *lgp) { CCTK_INT retval, ntries; CCTK_REAL W, x_rho[1], rho, rho_g ; - CCTK_REAL t15, t200,t1000 ; + CCTK_REAL t15, t200,t1000, drho_dW, rho_gm1 ; W = x[0]; - W_for_gnr2_old = W_for_gnr2; - W_for_gnr2 = W; + lgp->W_for_gnr2_old = lgp->W_for_gnr2; + lgp->W_for_gnr2 = W; // get rho from NR: - rho_g = x_rho[0] = rho_for_gnr2; + rho_g = x_rho[0] = lgp->rho_for_gnr2; ntries = 0; - while ( (retval = gnr2( x_rho, 1, gammaeos, func_rho)) && ( ntries++ < 10 ) ) { + while ( (retval = gnr2( x_rho, 1, gammaeos, lgp, func_rho)) && ( ntries++ < 10 ) ) { rho_g *= 10.; x_rho[0] = rho_g; } #if(DEBUG_CON2PRIMM) if( x_rho[0] <= 0. ) { - fprintf(stderr,"gnr2 neg rho = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], rho_for_gnr2, rho_for_gnr2_old, x[0], W_for_gnr2_old); + fprintf(stderr,"gnr2 neg rho = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], lgp->rho_for_gnr2, lgp->rho_for_gnr2_old, x[0], lgp->W_for_gnr2_old); fflush(stderr); } if( retval ) { - fprintf(stderr,"gnr2 retval = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], rho_for_gnr2, rho_for_gnr2_old, x[0], W_for_gnr2_old); + fprintf(stderr,"gnr2 retval = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], lgp->rho_for_gnr2, lgp->rho_for_gnr2_old, x[0], lgp->W_for_gnr2_old); fflush(stderr); } #endif - rho_for_gnr2_old = rho_for_gnr2; - rho = rho_for_gnr2 = x_rho[0]; + lgp->rho_for_gnr2_old = lgp->rho_for_gnr2; + rho = lgp->rho_for_gnr2 = x_rho[0]; CCTK_REAL gm1 = gammaeos-1.; rho_gm1 = pow(rho,gm1); - drho_dW = -rho*rho/( -rho_gm1*s200 + W*rho); + drho_dW = -rho*rho/( -rho_gm1*lgp->s200 + W*rho); - // t6 = rho*rho; - // t100 = (D-rho)*(D+rho); // t2 - t6; - t15 = -(D-rho)*(D+rho); // t6-t2 - t200 = W + two_Bsq; + t15 = -(lgp->D-rho)*(lgp->D+rho); // t6-t2 + t200 = W + lgp->two_Bsq; t1000 = rho*drho_dW; - resid[0] = (t300+(t4+t4+(t400+t15*(t7+(t200)*W))*W)*W)*t24; - jac[0][0] = 2*(t4+(t400+t15*t7+(3.0*t15*Bsq+t7*t1000+(t15+t15+t1000*(t200))*W)*W)*W)*t24; + resid[0] = (lgp->t300+(lgp->t4+lgp->t4+(lgp->t400+t15*(lgp->t7+(t200)*W))*W)*W)*lgp->t24; + jac[0][0] = 2*(lgp->t4+(lgp->t400+t15*lgp->t7+(3.0*t15*lgp->Bsq+lgp->t7*t1000+(t15+t15+t1000*(t200))*W)*W)*W)*lgp->t24; dx[0] = -resid[0]/jac[0][0]; @@ -705,11 +655,11 @@ *f = -0.5*(*df); -// fprintf(stdout,"QdotBsq = %28.18e ; \n",QdotBsq ); +// fprintf(stdout,"QdotBsq = %28.18e ; \n",lgp->QdotBsq ); // fprintf(stdout,"Sc = %28.18e ; \n",Sc ); -// fprintf(stdout,"Bsq = %28.18e ; \n",Bsq ); -// fprintf(stdout,"Qtsq = %28.18e ; \n",Qtsq ); -// fprintf(stdout,"Dc = %28.18e ; \n",D ); +// fprintf(stdout,"Bsq = %28.18e ; \n",lgp->Bsq ); +// fprintf(stdout,"Qtsq = %28.18e ; \n",lgp->Qtsq ); +// fprintf(stdout,"Dc = %28.18e ; \n",lgp->D ); // fprintf(stdout,"drhodW = %28.18e ; \n",drho_dW ); // fprintf(stdout,"W = %28.18e ; \n",W ); // fprintf(stdout,"rho = %28.18e ; \n",rho ); @@ -730,9 +680,9 @@ -- used to calculate rho from W *****************************************************************/ -static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, +static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT,CCTK_REAL) ) + CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT,CCTK_REAL, struct LocGlob *lgp) ) { CCTK_REAL f, df, dx[1], x_old[1], resid[1], jac[1][1]; @@ -757,7 +707,7 @@ keep_iterating = 1; while( keep_iterating ) { - (*funcd) (x, dx, resid, jac, &f, &df, n,gammaeos); /* returns with new dx, f, df */ + (*funcd) (x, dx, resid, jac, &f, &df, n,gammaeos, lgp); /* returns with new dx, f, df */ /* Save old values before calculating the new: */ //-fast errx = 0.; @@ -801,7 +751,7 @@ if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { #if(DEBUG_CON2PRIMM) fprintf(stderr,"\ngnr2 not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n", - f,df,x[0],x_old[0],W_for_gnr2_old,W_for_gnr2,rho_for_gnr2_old,rho_for_gnr2); fflush(stderr); + f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); #endif return(2); } @@ -839,7 +789,8 @@ *********************************************************************************/ // for the isentropic version: eq. (27) static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos) + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob *lgp) { CCTK_REAL A, B, C, rho, W, B0; @@ -848,13 +799,13 @@ CCTK_REAL gm1 = gammaeos-1.; rho = x[0]; - W = W_for_gnr2; + W = lgp->W_for_gnr2; - rho_gm1 = t40 = pow(rho,gm1); + t40 = pow(rho,gm1); - resid[0] = (rho*W+(-t40*s100-D)*D); + resid[0] = (rho*W+(-t40*lgp->s100-lgp->D)*lgp->D); t14 = t40/rho; // rho^(g-2) - jac[0][0] = -t14*s200 + W; + jac[0][0] = -t14*lgp->s200 + W; // drho_dW = -rho/jac[0][0]; dx[0] = -resid[0]/jac[0][0]; @@ -862,7 +813,7 @@ *f = -0.5*(*df); // fprintf(stdout,"deriv3 %g %g %g %g %g \n",rho,W,resid[0],jac[0][0],dx[0]); -// fprintf(stdout,"Dc := %28.18e ; \n",D); +// fprintf(stdout,"Dc := %28.18e ; \n",lgp->D); // fprintf(stdout,"Sc := %28.18e ; \n",Sc); // fprintf(stdout,"rho := %28.18e ; \n",rho); // fprintf(stdout,"W := %28.18e ; \n",W); @@ -874,40 +825,6 @@ } -/********************************************************************** - ********************************************************************** - - The following routines specify the equation of state. All routines - above here should be indpendent of EOS. If the user wishes - to use another equation of state, the below functions must be replaced - by equivalent routines based upon the new EOS. - - ********************************************************************** -**********************************************************************/ - -/**********************************************************************/ -/********************************************************************** - eos_info(): - - -- returns with all the EOS-related values needed; - **********************************************************************/ -static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos) -{ - register double ftmp,gtmp; - - ftmp = 1. - vsq; - gtmp = sqrt(ftmp); - - CCTK_REAL gam_m1_o_gam = (gammaeos-1.)/gammaeos; - - *dpdw = gam_m1_o_gam * ftmp ; - *dpdvsq = gam_m1_o_gam * ( 0.5 * D/gtmp - W ) ; - - return( gam_m1_o_gam * ( W * ftmp - D * gtmp ) ); // p - -} - - /****************************************************************************** END ******************************************************************************/ File [modified]: GRHydro_Con2PrimM_pt.c Delta lines: +45 -40 =================================================================== --- trunk/src/GRHydro_Con2PrimM_pt.c 2011-01-04 14:40:06 UTC (rev 207) +++ trunk/src/GRHydro_Con2PrimM_pt.c 2011-01-06 01:13:02 UTC (rev 208) @@ -71,23 +71,26 @@ ***************************************************/ /* Local Globals */ -static CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ; +struct LocGlob { + CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ; +} ; // Declarations: -static CCTK_REAL vsq_calc(CCTK_REAL W); +static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp); -static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, +static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *, CCTK_REAL *, CCTK_REAL) ); + CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); static void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos); -static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) ; + CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); +static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) ; + // EOS STUFF: -static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos); +static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos, struct LocGlob *lgp); /* pressure as a function of rho0 and u */ static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u, CCTK_REAL gammaeos) { @@ -213,6 +216,8 @@ CCTK_REAL sqrt_detg = sqrt(detg); CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_INT i,j, i_increase ; + + struct LocGlob lg; gammaeos = *gamma_eos; @@ -268,7 +273,7 @@ // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: - Bsq = + lg.Bsq = (*gxx) * (*Bx) * (*Bx) + (*gyy) * (*By) * (*By) + (*gzz) * (*Bz) * (*Bz) + @@ -278,15 +283,15 @@ (*gyz) * (*By) * (*Bz) ); QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; - QdotBsq = QdotB*QdotB ; + lg.QdotBsq = QdotB*QdotB ; - Qdotn = -(tau + dens) ; + lg.Qdotn = -(tau + dens) ; - Qtsq = (usx * sx + usy * sy + usz * sz) ; + lg.Qtsq = (usx * sx + usy * sy + usz * sz) ; - D = dens; + lg.D = dens; - half_Bsq = 0.5*Bsq; + lg.half_Bsq = 0.5*lg.Bsq; /* calculate W from last timestep and use for guess */ vsq = @@ -311,7 +316,7 @@ // Always calculate rho from D and gamma so that using D in EOS remains consistent // i.e. you don't get positive values for dP/d(vsq) . - rho0 = D / gamma ; + rho0 = lg.D / gamma ; u = (*epsilon) * rho0; p = pressure_rho0_u(rho0,u,gammaeos) ; // EOS w = rho0 + u + p ; @@ -321,8 +326,8 @@ // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; - while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) - - QdotBsq*(2.*W_last + Bsq) ) <= W_last*W_last*(Qtsq-Bsq*Bsq)) + while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) + - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) && (i_increase < 10) ) { W_last *= 10.; i_increase++; @@ -330,8 +335,8 @@ // Calculate W and vsq: x_2d[0] = fabs( W_last ); - x_2d[1] = x1_of_x0( W_last ) ; - *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, func_vsq ) ; + x_2d[1] = x1_of_x0( W_last, &lg ) ; + *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, &lg, func_vsq ) ; W = x_2d[0]; vsq = x_2d[1]; @@ -357,7 +362,7 @@ // Recover the primitive variables from the scalars and conserved variables: gtmp = sqrt(1. - vsq); gamma = 1./gtmp ; - rho0 = D * gtmp; + rho0 = lg.D * gtmp; w = W * (1. - vsq) ; p = pressure_rho0_w(rho0,w,gammaeos) ; // EOS @@ -375,9 +380,9 @@ *w_lorentz = gamma; *pressure = p ; - g_o_WBsq = 1./(W+Bsq); + g_o_WBsq = 1./(W+lg.Bsq); QdB_o_W = QdotB / W; - *bsq = Bsq * (1.-vsq) + QdB_o_W*QdB_o_W; + *bsq = lg.Bsq * (1.-vsq) + QdB_o_W*QdB_o_W; *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ; *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; @@ -411,15 +416,15 @@ W = \gamma^2 w ****************************************************************************/ -static CCTK_REAL vsq_calc(CCTK_REAL W) +static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp) { CCTK_REAL Wsq,Xsq,Bsq_W; Wsq = W*W ; - Bsq_W = (Bsq + W); + Bsq_W = (lgp->Bsq + W); Xsq = Bsq_W * Bsq_W; - return( ( Wsq * Qtsq + QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); + return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); } @@ -433,12 +438,12 @@ *********************************************************************/ -static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) +static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) { CCTK_REAL x1,vsq; CCTK_REAL dv = 1.e-15; - vsq = fabs(vsq_calc(x0)) ; // guaranteed to be positive + vsq = fabs(vsq_calc(x0,lgp)) ; // guaranteed to be positive return( ( vsq > 1. ) ? (1.0 - dv) : vsq ); @@ -478,10 +483,10 @@ -- inspired in part by Num. Rec.'s routine newt(); *****************************************************************/ -static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, +static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], CCTK_REAL *, - CCTK_REAL *, CCTK_REAL) ) + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) { CCTK_REAL f, df, dx[2], x_old[2]; CCTK_REAL resid[2], jac[2][2]; @@ -507,7 +512,7 @@ keep_iterating = 1; while( keep_iterating ) { - (*funcd) (x, dx, resid, jac, &f, &df, gammaeos); /* returns with new dx, f, df */ + (*funcd) (x, dx, resid, jac, &f, &df, gammaeos, lgp); /* returns with new dx, f, df */ /* Save old values before calculating the new: */ @@ -599,7 +604,7 @@ *********************************************************************************/ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos) + CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) { @@ -612,22 +617,22 @@ Wsq = W*W; - p_tmp = eos_info(W, vsq, &dPdW, &dPdvsq, gammaeos); + p_tmp = eos_info(W, vsq, &dPdW, &dPdvsq, gammaeos, lgp); // These expressions were calculated using Mathematica, but made into efficient // code using Maple. Since we know the analytic form of the equations, we can // explicitly calculate the Newton-Raphson step: - t2 = -half_Bsq+dPdvsq; - t3 = Bsq+W; + t2 = -lgp->half_Bsq+dPdvsq; + t3 = lgp->Bsq+W; t4 = t3*t3; t23 = 1/W; - t16 = QdotBsq*t23*t23; - t11 = Qtsq-vsq*t4+t16*(Bsq+W+W); - t18 = -Qdotn-half_Bsq*(1.0+vsq)+0.5*t16-W+p_tmp; + t16 = lgp->QdotBsq*t23*t23; + t11 = lgp->Qtsq-vsq*t4+t16*(lgp->Bsq+W+W); + t18 = -lgp->Qdotn-lgp->half_Bsq*(1.0+vsq)+0.5*t16-W+p_tmp; t24 = t16*t23; t25 = -1.0+dPdW-t24; - t35 = t25*t3+(Bsq-2.0*dPdvsq)*(t16+vsq*W)*t23; + t35 = t25*t3+(lgp->Bsq-2.0*dPdvsq)*(t16+vsq*W)*t23; // t21 = 1/t3; // t36 = 1/t35; t21 = 1/(t3*t35); @@ -666,7 +671,7 @@ -- returns with all the EOS-related values needed; **********************************************************************/ -static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos) +static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos, struct LocGlob *lgp) { register double ftmp,gtmp; @@ -676,9 +681,9 @@ double gam_m1_o_gam = ((gammaeos-1.)/gammaeos); *dpdw = gam_m1_o_gam * ftmp ; - *dpdvsq = gam_m1_o_gam * ( 0.5 * D/gtmp - W ) ; + *dpdvsq = gam_m1_o_gam * ( 0.5 * lgp->D/gtmp - W ) ; - return( gam_m1_o_gam * ( W * ftmp - D * gtmp ) ); // p + return( gam_m1_o_gam * ( W * ftmp - lgp->D * gtmp ) ); // p }

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 207)
by cott＠tapir.caltech.edu 04 Jan '11

04 Jan '11

User: cott Date: 2011/01/04 08:40 AM Modified: /trunk/src/ GRHydro_Prim2Con.F90 Log: * Update error handling in _hot routine of Prim2Con. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +25 -3 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2011-01-02 15:54:48 UTC (rev 206) +++ trunk/src/GRHydro_Prim2Con.F90 2011-01-04 14:40:06 UTC (rev 207) @@ -231,9 +231,11 @@ ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: temp0(1) n = 1;keytemp = 0;anyerr = 0;keyerr(1) = 0 ! end EOS Omni vars + temp0 = temp w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz & *dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz& *dvely*dvelz)) @@ -262,19 +264,39 @@ call CCTK_WARN(1,warnline) endif else - call CCTK_WARN(1,"EOS error in prim2con_hot") + ! This is a way of recovering even on finer refinement levels: + ! Use the average temperature at the interface instead of the + ! reconstructed specific internal energy. + call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z call CCTK_WARN(1,warnline) - write(warnline,"(1P10E15.6)") drho,deps,temp,ye + write(warnline,"(1P10E15.6)") drho,deps,temp0,ye call CCTK_WARN(1,warnline) write(warnline,"(A7,i8)") "code: ",keyerr(1) call CCTK_WARN(1,warnline) write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(1,warnline) - call CCTK_WARN(0,"Aborting!!!") + keytemp=1 + temp = temp0 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + keytemp=0 + if(anyerr.ne.0) then + call CCTK_WARN(1,"EOS error in prim2con_hot") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif endif endif + vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz

1 0

[svn:einsteintoolkit] Multipole/trunk/ (Rev. 74)
by ian.hinder＠aei.mpg.de 03 Jan '11

03 Jan '11

User: hinder Date: 2011/01/03 12:50 PM Added: /trunk/ configuration.ccl Modified: /trunk/ interface.ccl, param.ccl /trunk/src/ multipole.cc, utils.cc, utils.hh Log: Add HDF5 output support One HDF5 file per variable, and one extensible dataset per radius per mode. HDF5 is required only optionally, so this thorn can be compiled without it. File Changes: Directory: /trunk/src/ ====================== File [modified]: multipole.cc Delta lines: +21 -8 =================================================================== --- trunk/src/multipole.cc 2011-01-03 18:50:31 UTC (rev 73) +++ trunk/src/multipole.cc 2011-01-03 18:50:53 UTC (rev 74) @@ -97,16 +97,29 @@ static void output_mode(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, CCTK_REAL real_lm, CCTK_REAL imag_lm, int l, int m) { - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; - if (CCTK_MyProc(cctkGH) == 0) + if (output_ascii) { - ostringstream name; - name << "mp_" << v->name << "_l" << l << "_m" << m << - "_r" << setiosflags(ios::fixed) << setprecision(2) << rad << ".asc"; - Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name.str(), real_lm, imag_lm); + if (CCTK_MyProc(cctkGH) == 0) + { + ostringstream name; + name << "mp_" << v->name << "_l" << l << "_m" << m << + "_r" << setiosflags(ios::fixed) << setprecision(2) << rad << ".asc"; + Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name.str(), real_lm, imag_lm); + } } + if (output_hdf5) + { + if (CCTK_MyProc(cctkGH) == 0) + { + ostringstream h5datasetname; + h5datasetname << "l" << l << "_m" << m << "_r" << setiosflags(ios::fixed) << setprecision(2) << rad; + Multipole_OutputComplexToH5File(CCTK_PASS_CTOC, "mp_" + v->name + ".h5", h5datasetname.str(), + real_lm, imag_lm); + } + } } static void output_1D(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, @@ -116,7 +129,7 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - if (CCTK_MyProc(cctkGH) == 0) + if (CCTK_MyProc(cctkGH) == 0 && output_ascii) { if (out_1d_every != 0 && cctk_iteration % out_1d_every == 0) { File [modified]: utils.cc Delta lines: +120 -0 =================================================================== --- trunk/src/utils.cc 2011-01-03 18:50:31 UTC (rev 73) +++ trunk/src/utils.cc 2011-01-03 18:50:53 UTC (rev 74) @@ -2,12 +2,25 @@ #include <string.h> #include <math.h> #include <string> +#include <sys/stat.h> +#include <errno.h> #include <iostream> #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" +#ifdef HAVE_CAPABILITY_HDF5 +// We currently support the HDF5 1.6 API (and when using 1.8 the +// compatibility mode introduced by H5_USE_16_API). Several machines +// in SimFactory use HDF5 1.6, so we cannot drop support for it. It +// seems it is hard to support both the 1.6 and 1.8 API +// simultaneously; for example H5Fopen takes a different number of +// arguments in the two versions. +#define H5_USE_16_API +#include <hdf5.h> +#endif + #include "utils.hh" #include "integrate.hh" @@ -133,7 +146,114 @@ } } +//////////////////////////////////////////////////////////////////////// +// HDF5 complex output +//////////////////////////////////////////////////////////////////////// +#ifdef HAVE_CAPABILITY_HDF5 + +bool dataset_exists(hid_t file, const string &dataset_name) +{ + // To test whether a dataset exists, the recommended way in API 1.6 + // is to use H5Gget_objinfo, but this prints an error to stderr if + // the dataset does not exist. We explicitly avoid this by wrapping + // the call in H5E_BEGIN_TRY/H5E_END_TRY statements. In 1.8, + // H5Gget_objinfo is deprecated, and H5Lexists does the job. See + // http://www.mail-archive.com/hdf-forum@hdfgroup.org/msg00125.html + + #if 1 + bool exists; + H5E_BEGIN_TRY + { + exists = H5Gget_objinfo(file, dataset_name.c_str(), 1, NULL) >= 0; + } + H5E_END_TRY; + return exists; + #else + return H5Lexists(file, dataset_name.c_str(), H5P_DEFAULT); + #endif +} + +void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &basename, const string &datasetname, + CCTK_REAL redata, CCTK_REAL imdata) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + string output_name = out_dir + string("/") + basename; + int status = 0; + static int first_time = true; + + hid_t file = + (first_time && IO_TruncateOutputFiles(cctkGH)) ? + H5Fcreate(output_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) : + H5Fopen(output_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + + first_time = false; + + hid_t dataset = -1; + + if (dataset_exists(file, datasetname)) + { + dataset = H5Dopen(file, datasetname.c_str()); + } + else + { + hsize_t dims[2] = {0,3}; + hsize_t maxdims[2] = {H5S_UNLIMITED, 3}; + hid_t dataspace = H5Screate_simple(2, dims, maxdims); + + hid_t cparms = -1; + hsize_t chunk_dims[2] = {hdf5_chunk_size,3}; + cparms = H5Pcreate (H5P_DATASET_CREATE); + status = H5Pset_chunk(cparms, 2, chunk_dims); + + dataset = H5Dcreate(file, datasetname.c_str(), H5T_NATIVE_DOUBLE, dataspace, cparms); + H5Pclose(cparms); + } + + hid_t filespace = H5Dget_space (dataset); + hsize_t filedims[2]; + hsize_t maxdims[2]; + status = H5Sget_simple_extent_dims(filespace, filedims, maxdims); + + filedims[0] += 1; + hsize_t size[2] = {filedims[0], filedims[1]}; + status = H5Dextend (dataset, size); + H5Sclose(filespace); + + /* Select a hyperslab */ + hsize_t offset[2] = {filedims[0]-1, 0}; + hsize_t dims2[2] = {1,3}; + filespace = H5Dget_space (dataset); + status = H5Sselect_hyperslab (filespace, H5S_SELECT_SET, offset, NULL, + dims2, NULL); + + CCTK_REAL data[] = {cctk_time, redata, imdata}; + + hid_t memdataspace = H5Screate_simple(2, dims2, NULL); + + /* Write the data to the hyperslab */ + status = H5Dwrite (dataset, H5T_NATIVE_DOUBLE, memdataspace, filespace, + H5P_DEFAULT, data); + + H5Dclose(dataset); + H5Sclose(filespace); + H5Sclose(memdataspace); + + status = H5Fclose(file); +} + +#else + +void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &basename, const string &datasetname, + CCTK_REAL redata, CCTK_REAL imdata) +{ + CCTK_WARN(0,"HDF5 output has been requested but Cactus has been compiled without HDF5 support"); +} + +#endif + //////////////////////////////////////////////////////////////////////// // Coordinates //////////////////////////////////////////////////////////////////////// File [modified]: utils.hh Delta lines: +3 -0 =================================================================== --- trunk/src/utils.hh 2011-01-03 18:50:31 UTC (rev 73) +++ trunk/src/utils.hh 2011-01-03 18:50:53 UTC (rev 74) @@ -20,6 +20,9 @@ void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL redata, CCTK_REAL imdata); +void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &filename, const string &dataset, + CCTK_REAL redata, CCTK_REAL imdata); + void Multipole_CoordSetup(int ntheta, int nphi, CCTK_REAL xhat[], CCTK_REAL yhat[], CCTK_REAL zhat[], CCTK_REAL th[], Directory: /trunk/ ================== File [added]: configuration.ccl Delta lines: +3 -0 =================================================================== --- trunk/configuration.ccl (rev 0) +++ trunk/configuration.ccl 2011-01-03 18:50:53 UTC (rev 74) @@ -0,0 +1,3 @@ +OPTIONAL HDF5 +{ +} File [modified]: interface.ccl Delta lines: +5 -0 =================================================================== --- trunk/interface.ccl 2011-01-03 18:50:31 UTC (rev 73) +++ trunk/interface.ccl 2011-01-03 18:50:53 UTC (rev 74) @@ -2,6 +2,11 @@ implements: multipole inherits: Grid +CCTK_INT FUNCTION IO_TruncateOutputFiles \ + (CCTK_POINTER_TO_CONST IN GH) + +REQUIRES FUNCTION IO_TruncateOutputFiles + public: CCTK_REAL harmonics type=GF timelevels=1 File [modified]: param.ccl Delta lines: +13 -0 =================================================================== --- trunk/param.ccl 2011-01-03 18:50:31 UTC (rev 73) +++ trunk/param.ccl 2011-01-03 18:50:53 UTC (rev 74) @@ -137,3 +137,16 @@ { -100:* :: "Deprecated" } -100 + +CCTK_BOOLEAN output_ascii "Output a simple ASCII file for each mode at each radius" +{ +} "yes" + +CCTK_BOOLEAN output_hdf5 "Output an HDF5 file for each variable containing one dataset per mode at each radius" +{ +} "no" + +CCTK_INT hdf5_chunk_size "How many iterations to preallocate in extensible HDF5 datasets" +{ + 1:* :: "Any integer" +} 200

1 0

[svn:einsteintoolkit] Multipole/trunk/ (Rev. 73)
by ian.hinder＠aei.mpg.de 03 Jan '11

03 Jan '11

User: hinder Date: 2011/01/03 12:50 PM Modified: /trunk/ param.ccl /trunk/src/ utils.cc Log: Use IO::out_dir parameter by sharing from implementation rather than via CCTK_ParameterGet File Changes: Directory: /trunk/src/ ====================== File [modified]: utils.cc Delta lines: +1 -8 =================================================================== --- trunk/src/utils.cc 2011-01-03 18:50:17 UTC (rev 72) +++ trunk/src/utils.cc 2011-01-03 18:50:31 UTC (rev 73) @@ -24,16 +24,9 @@ bool first_time = cctk_iteration == 0; const char *mode = first_time ? "w" : "a"; - CCTK_STRING *out_dir = (CCTK_STRING *) CCTK_ParameterGet("out_dir", "IOUtil", NULL); - if (*out_dir == 0) - { - CCTK_WARN(1,"Parameter IOUtil::out_dir not found"); - return 0; - } + string output_name(string(out_dir) + "/" + name); - string output_name(string(*out_dir) + "/" + name); - FILE *fp = fopen(output_name.c_str(), mode); if (fp == 0) Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +4 -0 =================================================================== --- trunk/param.ccl 2011-01-03 18:50:17 UTC (rev 72) +++ trunk/param.ccl 2011-01-03 18:50:31 UTC (rev 73) @@ -1,5 +1,9 @@ #param.ccl for thorn Multipole +shares: IO + +USES STRING out_dir + restricted: #Interpolator params: CCTK_STRING interpolator_name "Which interpolator should I use"

1 0

[svn:einsteintoolkit] Multipole/trunk/src/ (Rev. 72)
by ian.hinder＠aei.mpg.de 03 Jan '11

03 Jan '11

User: hinder Date: 2011/01/03 12:50 PM Modified: /trunk/src/ multipole.cc, utils.cc, utils.hh Log: Convert to C++ strings File Changes: Directory: /trunk/src/ ====================== File [modified]: multipole.cc Delta lines: +26 -26 =================================================================== --- trunk/src/multipole.cc 2010-11-21 03:29:03 UTC (rev 71) +++ trunk/src/multipole.cc 2011-01-03 18:50:17 UTC (rev 72) @@ -1,6 +1,8 @@ #include <stdio.h> #include <string.h> #include <assert.h> +#include <string> +#include <iomanip> #include "cctk.h" #include "cctk_Arguments.h" @@ -12,6 +14,8 @@ #include "utils.hh" #include "sphericalharmonic.hh" +using namespace std; + static const int var_name_length = 30; static const int max_vars = 10; @@ -24,7 +28,7 @@ int index; int imag_index; int spin_weight; - char name[var_name_length]; + string name; } variable_desc; @@ -51,7 +55,7 @@ // not present v->imag_index = -1; v->spin_weight = 0; - strcpy(v->name, CCTK_VarName(v->index)); + v->name = string(CCTK_VarName(v->index)); if (optstring != 0) { @@ -60,30 +64,31 @@ if (table >= 0) { const int buffer_length = 256; - char imag_name[buffer_length]; + char buffer[buffer_length]; Util_TableGetInt(table, &v->spin_weight , "sw"); ///////////////////////////////////////////////////////////// CCTK_VInfo(CCTK_THORNSTRING,"spinweight %d", v->spin_weight); ///////////////////////////////////////////////////////////// - if (Util_TableGetString(table, buffer_length, imag_name , "cmplx") >= 0) + if (Util_TableGetString(table, buffer_length, buffer , "cmplx") >= 0) { - v->imag_index = CCTK_VarIndex(imag_name); + v->imag_index = CCTK_VarIndex(buffer); } - Util_TableGetString(table, var_name_length, v->name , "name"); + Util_TableGetString(table, buffer_length, buffer , "name"); + v->name = string(buffer); } } vs->n_vars++; } -static void parse_variables_string(const char *var_string, variable_desc v[max_vars], int *n_variables) +static void parse_variables_string(const string &var_string, variable_desc v[max_vars], int *n_variables) { variables_desc vars; vars.n_vars = 0; vars.vars = v; - int ierr = CCTK_TraverseString(var_string, fill_variable, &vars, CCTK_GROUP_OR_VAR); + int ierr = CCTK_TraverseString(var_string.c_str(), fill_variable, &vars, CCTK_GROUP_OR_VAR); assert(ierr > 0); *n_variables = vars.n_vars; @@ -95,12 +100,12 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - char name_tmp[1000]; - if (CCTK_MyProc(cctkGH) == 0) { - sprintf(name_tmp, "mp_%s_l%d_m%d_r%1.2f.asc", v->name, l, m, rad); - Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name_tmp, real_lm, imag_lm); + ostringstream name; + name << "mp_" << v->name << "_l" << l << "_m" << m << + "_r" << setiosflags(ios::fixed) << setprecision(2) << rad << ".asc"; + Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name.str(), real_lm, imag_lm); } } @@ -111,26 +116,21 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - char name_tmp[1000]; - if (CCTK_MyProc(cctkGH) == 0) { if (out_1d_every != 0 && cctk_iteration % out_1d_every == 0) { - const char *real_name = CCTK_VarName(v->index); - sprintf(name_tmp, "mp_%s_r%1.2f.th.asc", real_name, rad); - Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_theta, real); - sprintf(name_tmp, "mp_%s_r%1.2f.ph.asc", real_name, rad); - Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_phi, real); + ostringstream real_base; + real_base << "mp_" << string(CCTK_VarName(v->index)) << "_r" << setiosflags(ios::fixed) << setprecision(2) << rad; + Multipole_Output1D(CCTK_PASS_CTOC, real_base.str()+string(".th.asc"), array_size, th, ph, mp_theta, real); + Multipole_Output1D(CCTK_PASS_CTOC, real_base.str()+string(".ph.asc"), array_size, th, ph, mp_phi, real); if (v->imag_index != -1) { - const char *imag_name = CCTK_VarName(v->imag_index); - sprintf(name_tmp, "mp_%s_r%1.2f.th.asc", imag_name, rad); - Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_theta, imag); - - sprintf(name_tmp, "mp_%s_r%1.2f.ph.asc", imag_name, rad); - Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_phi, imag); + ostringstream imag_base; + imag_base << "mp_" << string(CCTK_VarName(v->imag_index)) << "_r" << setiosflags(ios::fixed) << setprecision(2) << rad; + Multipole_Output1D(CCTK_PASS_CTOC, imag_base.str()+string(".th.asc"), array_size, th, ph, mp_theta, imag); + Multipole_Output1D(CCTK_PASS_CTOC, imag_base.str()+string(".ph.asc"), array_size, th, ph, mp_phi, imag); } } } @@ -285,7 +285,7 @@ yhat = new CCTK_REAL[array_size]; zhat = new CCTK_REAL[array_size]; - parse_variables_string(variables, vars, &n_variables); + parse_variables_string(string(variables), vars, &n_variables); get_spin_weights(vars, n_variables, spin_weights, &n_spin_weights); Multipole_CoordSetup(ntheta, nphi, xhat, yhat, zhat, th, ph); setup_harmonics(spin_weights, n_spin_weights, lmax, th, ph, File [modified]: utils.cc Delta lines: +11 -15 =================================================================== --- trunk/src/utils.cc 2010-11-21 03:29:03 UTC (rev 71) +++ trunk/src/utils.cc 2011-01-03 18:50:17 UTC (rev 72) @@ -1,6 +1,8 @@ #include <stdio.h> #include <string.h> #include <math.h> +#include <string> +#include <iostream> #include "cctk.h" #include "cctk_Arguments.h" @@ -9,19 +11,19 @@ #include "utils.hh" #include "integrate.hh" +using namespace std; + //////////////////////////////////////////////////////////////////////// // File manipulation //////////////////////////////////////////////////////////////////////// -FILE *Multipole_OpenOutputFile(CCTK_ARGUMENTS, char const *name) +FILE *Multipole_OpenOutputFile(CCTK_ARGUMENTS, const string &name) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - const int buf_size = 1024; bool first_time = cctk_iteration == 0; const char *mode = first_time ? "w" : "a"; - char output_name[buf_size]; CCTK_STRING *out_dir = (CCTK_STRING *) CCTK_ParameterGet("out_dir", "IOUtil", NULL); if (*out_dir == 0) @@ -30,19 +32,13 @@ return 0; } - if ((int)strlen(*out_dir)+1 > buf_size) - { - CCTK_WARN(1,"Parameter IOUtil::out_dir is too long"); - return 0; - } + string output_name(string(*out_dir) + "/" + name); - sprintf(output_name, "%s/%s", *out_dir, name); + FILE *fp = fopen(output_name.c_str(), mode); - FILE *fp = fopen(output_name, mode); - if (fp == 0) { - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "Could not open output file %s", output_name); + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, (string("Could not open output file ") + output_name).c_str()); } return fp; @@ -72,7 +68,7 @@ } -void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, char *name, int array_size, +void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, const string &name, int array_size, CCTK_REAL const th[], CCTK_REAL const ph[], CCTK_REAL const xs[], CCTK_REAL const ys[], CCTK_REAL const zs[], CCTK_REAL const data[]) @@ -90,7 +86,7 @@ // Misc //////////////////////////////////////////////////////////////////////// -void Multipole_Output1D(CCTK_ARGUMENTS, char const *name, int array_size, +void Multipole_Output1D(CCTK_ARGUMENTS, const string &name, int array_size, CCTK_REAL const th[], CCTK_REAL const ph[], mp_coord coord, CCTK_REAL const data[]) { @@ -133,7 +129,7 @@ fprintf(fp, "%f %.19g %.19g\n", cctk_time, redata, imdata); } -void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, char const *name, CCTK_REAL redata, CCTK_REAL imdata) +void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL redata, CCTK_REAL imdata) { DECLARE_CCTK_ARGUMENTS; File [modified]: utils.hh Delta lines: +6 -3 =================================================================== --- trunk/src/utils.hh 2010-11-21 03:29:03 UTC (rev 71) +++ trunk/src/utils.hh 2011-01-03 18:50:17 UTC (rev 72) @@ -3,19 +3,22 @@ #define __utils_h #include "cctk.h" +#include <string> +using namespace std; + enum mp_coord {mp_theta, mp_phi}; -void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, char const *name, int array_size, +void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, const string &name, int array_size, CCTK_REAL const th[], CCTK_REAL const ph[], CCTK_REAL const x[], CCTK_REAL const y[], CCTK_REAL const z[], CCTK_REAL const data[]); -void Multipole_Output1D(CCTK_ARGUMENTS, char const *name, int array_size, +void Multipole_Output1D(CCTK_ARGUMENTS, const string &name, int array_size, CCTK_REAL const th[], CCTK_REAL const ph[], mp_coord coord, CCTK_REAL const data[]); -void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, char const *name, CCTK_REAL redata, CCTK_REAL imdata); +void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL redata, CCTK_REAL imdata); void Multipole_CoordSetup(int ntheta, int nphi, CCTK_REAL xhat[], CCTK_REAL yhat[],

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 206)
by cott＠tapir.caltech.edu 02 Jan '11

02 Jan '11

User: cott Date: 2011/01/02 09:54 AM Modified: /trunk/ schedule.ccl Log: * change scheduling of GRHydro_EOSHandle from local to global. File Changes: Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +2 -0 =================================================================== --- trunk/schedule.ccl 2011-01-01 00:04:02 UTC (rev 205) +++ trunk/schedule.ccl 2011-01-02 15:54:48 UTC (rev 206) @@ -235,11 +235,13 @@ schedule GRHydro_EOSHandle AT CCTK_INITIAL before HydroBase_Initial { LANG: C + OPTIONS: global } "Set the EOS number" schedule GRHydro_EOSHandle AT CCTK_POST_RECOVER_VARIABLES { LANG: C + OPTIONS: global } "Set the EOS number" ####################################

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/par/ (Rev. 120)
by bcmsma＠astro.rit.edu 31 Dec '10

31 Dec '10

User: bmundim Date: 2010/12/31 06:10 PM Added: /trunk/par/ balsara1a.par, c2pM_omni.par, p2cM_omni.par, p2cM_polytype.par Log: Add a few parameter files for MHD testing purposes. File Changes: Directory: /trunk/par/ ====================== File [added]: balsara1a.par Delta lines: +61 -0 =================================================================== --- trunk/par/balsara1a.par (rev 0) +++ trunk/par/balsara1a.par 2011-01-01 00:10:16 UTC (rev 120) @@ -0,0 +1,61 @@ +ActiveThorns = "time MoL cartgrid3d carpetioascii ioutil Fortran boundary hydrobase grhydro grhydro_initdata nanchecker eos_omni CoordBase LoopControl CarpetLib Carpet CarpetIOHDF5 CarpetReduce spacemask ADMCoupling ADMMacros ADMBASE STATICCONFORMAL TMUNUBASE SymBase localreduce CarpetIOBasic initbase GRHydro_InitData" + +#This is a hack. ghost_size should be the same as grhydro::grhydro_stencil. + +driver::ghost_size=3 +grhydro::grhydro_stencil=3 + +time::dtfac = 0.25 + +methodoflines::ODE_Method = "rk2" +methodoflines::MoL_Intermediate_Steps=2 + +#admbase::initial_data="prim2con2prim_test" +#GRHydro_InitData::Bx_init=0.1 +#GRHydro_InitData::By_init=0.12 +#GRHydro_InitData::Bz_init=0.14 +HydroBase::initial_hydro="shocktube" +GRHydro_InitData::shocktube_type="xshock" +GRHydro_InitData::shock_case="Balsara1" + +HydroBase::timelevels=2 +HydroBase::evolution_method="grhydro" +HydroBase::Bvec_evolution_method="grhydro" + +SpaceMask::use_mask="yes" +EOS_Omni::gl_gamma=2.0 +grhydro::recon_method="tvd" +grhydro::method_type="RSA FV" +grhydro::recon_vars="primitive" +grhydro::riemann_solver="HLLE" +grhydro::tvd_limiter="vanleerMC" +grhydro::bound="flat" +grhydro::grhydro_eos_type = "General" +grhydro::grhydro_eos_table="Ideal_Fluid" +grhydro::Grhydro_MaxNumConstrainedVars=13 +grhydro::GRHydro_MaxNumEvolvedVars=10 + +grid::type = "BySpacing" +grid::domain = "full" +grid::dxyz = 0.01 + +driver::global_nx = 150 +driver::global_ny = 10 +driver::global_nz = 10 + +#cactus::cctk_itlast = 40 +cactus::cctk_itlast = 200 + +IO::out_dir = "balsara1a" +#IOBasic::outInfo_every = 1 +#IOBasic::outInfo_vars = "HydroBase::rho" +CarpetIOBasic::outInfo_vars="hydrobase::rho" +CarpetIOBasic::outInfo_every=1 +CarpetIOASCII::out1D_criterion = "time" +CarpetIOASCII::out1D_dt = 0.01 +CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec" + +CarpetIOHDF5::out_every = 1 +CarpetIOHDF5::out_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec" + + File [added]: c2pM_omni.par Delta lines: +54 -0 =================================================================== --- trunk/par/c2pM_omni.par (rev 0) +++ trunk/par/c2pM_omni.par 2011-01-01 00:10:16 UTC (rev 120) @@ -0,0 +1,54 @@ +ActiveThorns = "time MoL cartgrid3d carpetioascii ioutil Fortran boundary hydrobase grhydro grhydro_initdata nanchecker eos_omni CoordBase LoopControl CarpetLib Carpet CarpetReduce spacemask ADMCoupling ADMMacros ADMBASE STATICCONFORMAL TMUNUBASE SymBase localreduce CarpetIOBasic initbase GRHydro_InitData" + +#This is a hack. ghost_size should be the same as grhydro::grhydro_stencil. + +driver::ghost_size=3 +grhydro::grhydro_stencil=3 + +time::dtfac = 0.25 + +methodoflines::ODE_Method = "rk2" +methodoflines::MoL_Intermediate_Steps=2 + +admbase::initial_data="con2prim2con_test" +GRHydro_InitData::Bx_init=0.1 +GRHydro_InitData::By_init=0.12 +GRHydro_InitData::Bz_init=0.14 +HydroBase::initial_hydro="shocktube" +GRHydro_InitData::shocktube_type="xshock" +GRHydro_InitData::shock_case="Simple" + +HydroBase::timelevels=2 +HydroBase::evolution_method="grhydro" +HydroBase::Bvec_evolution_method="grhydro" + +SpaceMask::use_mask="yes" +EOS_Omni::gl_gamma=2.5 +#EOS_Omni::gl_gamma=1.66666666666666 +grhydro::recon_method="ppm" +grhydro::method_type="RSA FV" +grhydro::recon_vars="primitive" +grhydro::riemann_solver="HLLE" +grhydro::tvd_limiter="vanleerMC" +grhydro::bound="flat" +grhydro::grhydro_eos_table="Ideal_Fluid" +grhydro::Grhydro_MaxNumConstrainedVars=13 +grhydro::GRHydro_MaxNumEvolvedVars=9 + +grid::type = "BySpacing" +grid::domain = "full" +grid::dxyz = 0.01 + +driver::global_nx = 200 +driver::global_ny = 7 +driver::global_nz = 7 + +cactus::cctk_itlast = 200 + +IO::out_dir = "grhydro_2c_omni" +#IOBasic::outInfo_every = 1 +#IOBasic::outInfo_vars = "HydroBase::rho" +CarpetIOBasic::outInfo_vars="hydrobase::rho" +CarpetIOBasic::outInfo_every=1 +CarpetIOASCII::out1D_every = 1 +CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon" File [added]: p2cM_omni.par Delta lines: +54 -0 =================================================================== --- trunk/par/p2cM_omni.par (rev 0) +++ trunk/par/p2cM_omni.par 2011-01-01 00:10:16 UTC (rev 120) @@ -0,0 +1,54 @@ +ActiveThorns = "time MoL cartgrid3d carpetioascii ioutil Fortran boundary hydrobase grhydro grhydro_initdata nanchecker eos_omni CoordBase LoopControl CarpetLib Carpet CarpetReduce spacemask ADMCoupling ADMMacros ADMBASE STATICCONFORMAL TMUNUBASE SymBase localreduce CarpetIOBasic initbase GRHydro_InitData" + +#This is a hack. ghost_size should be the same as grhydro::grhydro_stencil. + +driver::ghost_size=3 +grhydro::grhydro_stencil=3 + +time::dtfac = 0.25 + +methodoflines::ODE_Method = "rk2" +methodoflines::MoL_Intermediate_Steps=2 + +admbase::initial_data="prim2con2prim_test" +GRHydro_InitData::Bx_init=0.1 +GRHydro_InitData::By_init=0.12 +GRHydro_InitData::Bz_init=0.14 +HydroBase::initial_hydro="shocktube" +GRHydro_InitData::shocktube_type="xshock" +GRHydro_InitData::shock_case="Simple" + +HydroBase::timelevels=2 +HydroBase::evolution_method="grhydro" +HydroBase::Bvec_evolution_method="grhydro" + +SpaceMask::use_mask="yes" +EOS_Omni::gl_gamma=2.5 +#EOS_Omni::gl_gamma=1.66666666666666 +grhydro::recon_method="ppm" +grhydro::method_type="RSA FV" +grhydro::recon_vars="primitive" +grhydro::riemann_solver="HLLE" +grhydro::tvd_limiter="vanleerMC" +grhydro::bound="flat" +grhydro::grhydro_eos_table="Ideal_Fluid" +grhydro::Grhydro_MaxNumConstrainedVars=13 +grhydro::GRHydro_MaxNumEvolvedVars=9 + +grid::type = "BySpacing" +grid::domain = "full" +grid::dxyz = 0.01 + +driver::global_nx = 200 +driver::global_ny = 7 +driver::global_nz = 7 + +cactus::cctk_itlast = 200 + +IO::out_dir = "grhydro_2c_omni" +#IOBasic::outInfo_every = 1 +#IOBasic::outInfo_vars = "HydroBase::rho" +CarpetIOBasic::outInfo_vars="hydrobase::rho" +CarpetIOBasic::outInfo_every=1 +CarpetIOASCII::out1D_every = 1 +CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon" File [added]: p2cM_polytype.par Delta lines: +56 -0 =================================================================== --- trunk/par/p2cM_polytype.par (rev 0) +++ trunk/par/p2cM_polytype.par 2011-01-01 00:10:16 UTC (rev 120) @@ -0,0 +1,56 @@ +ActiveThorns = "time MoL cartgrid3d carpetioascii ioutil Fortran boundary hydrobase grhydro grhydro_initdata nanchecker eos_omni CoordBase LoopControl CarpetLib Carpet CarpetReduce spacemask ADMCoupling ADMMacros ADMBASE STATICCONFORMAL TMUNUBASE SymBase localreduce CarpetIOBasic initbase GRHydro_InitData" + +#This is a hack. ghost_size should be the same as grhydro::grhydro_stencil. + +driver::ghost_size=3 +grhydro::grhydro_stencil=3 + +time::dtfac = 0.25 + +methodoflines::ODE_Method = "rk2" +methodoflines::MoL_Intermediate_Steps=2 + +admbase::initial_data="prim2con2prim_polytype_test" +GRHydro_InitData::Bx_init=0.1 +GRHydro_InitData::By_init=0.12 +GRHydro_InitData::Bz_init=0.14 +HydroBase::initial_hydro="shocktube" +GRHydro_InitData::shocktube_type="xshock" +GRHydro_InitData::shock_case="Simple" + +HydroBase::timelevels=2 +HydroBase::evolution_method="grhydro" +HydroBase::Bvec_evolution_method="grhydro" + +SpaceMask::use_mask="yes" +EOS_Omni::poly_gamma=2.5 +EOS_Omni::poly_gamma_ini=2.5 +EOS_Omni::poly_k = 4.2 +#EOS_Omni::gl_gamma=1.66666666666666 +grhydro::recon_method="ppm" +grhydro::method_type="RSA FV" +grhydro::recon_vars="primitive" +grhydro::riemann_solver="HLLE" +grhydro::tvd_limiter="vanleerMC" +grhydro::bound="flat" +grhydro::grhydro_eos_table="Ideal_Fluid" +grhydro::Grhydro_MaxNumConstrainedVars=13 +grhydro::GRHydro_MaxNumEvolvedVars=9 + +grid::type = "BySpacing" +grid::domain = "full" +grid::dxyz = 0.01 + +driver::global_nx = 200 +driver::global_ny = 7 +driver::global_nz = 7 + +cactus::cctk_itlast = 200 + +IO::out_dir = "grhydro_2c_omni" +#IOBasic::outInfo_every = 1 +#IOBasic::outInfo_vars = "HydroBase::rho" +CarpetIOBasic::outInfo_vars="hydrobase::rho" +CarpetIOBasic::outInfo_every=1 +CarpetIOASCII::out1D_every = 1 +CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon"

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits