Commits May 2012

commits@lists.einsteintoolkit.org

11 participants
183 discussions

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 325)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:12 PM Modified: /trunk/src/ GRHydro_Boundaries.F90 Log: GRHydro: Fixed minor bugs in registration in Boundaries Some variables were stil set to "Flat" boundaries even though "None" is supposed to be set. Patch by Christian Reisswig. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Boundaries.F90 Delta lines: +9 -3 =================================================================== --- trunk/src/GRHydro_Boundaries.F90 2012-05-13 19:12:00 UTC (rev 324) +++ trunk/src/GRHydro_Boundaries.F90 2012-05-13 19:12:48 UTC (rev 325) @@ -260,8 +260,14 @@ "HydroBase::press", "None") ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::eps", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "None") + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lvel", "None") + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::vel", "None") + endif + if(evolve_mhd.ne.0) then if (general_coordinates .ne. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & @@ -299,7 +305,7 @@ "HydroBase::entropy", "None") endif - end if + endif if (CCTK_EQUALS(bound,"scalar")) then call CCTK_WARN(0, "Until somebody uses this I see no reason to support it")

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 324)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:12 PM Modified: /trunk/src/ GRHydro_Con2PrimM_pt_EOSOmni.c Log: Adding Josh's latest changes to the general EOS Con2PrimM routine. patch by Philipp Moesta File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM_pt_EOSOmni.c Delta lines: +5 -5 =================================================================== --- trunk/src/GRHydro_Con2PrimM_pt_EOSOmni.c 2012-05-13 19:10:12 UTC (rev 323) +++ trunk/src/GRHydro_Con2PrimM_pt_EOSOmni.c 2012-05-13 19:12:00 UTC (rev 324) @@ -89,7 +89,7 @@ CCTK_REAL [], CCTK_REAL [][2], CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); -static CCTK_INT threed_newton_raphson( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, +static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][3], CCTK_REAL *, CCTK_REAL *, @@ -386,7 +386,7 @@ } else { //USE 3d NR for non-polytropes! x_3d[2] = u; - *retval = 1.0*threed_newton_raphson( x_3d, &eosvars, &lg, func_vsq_eosomni ) ; + *retval = 1.0*threed_newton_raphson_omni( x_3d, &eosvars, &lg, func_vsq_eosomni ) ; } W = x_3d[0]; @@ -646,14 +646,14 @@ /************************************************************ - threed_newton_raphson(): + threed_newton_raphson_omni(): -- performs Newton-Rapshon method on an 2d system for polytropes. -- inspired in part by Num. Rec.'s routine newt(); *****************************************************************/ -static CCTK_INT threed_newton_raphson( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, +static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *eosvars, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][3], CCTK_REAL *, CCTK_REAL *, struct eosomnivars *, struct LocGlob *) ) @@ -900,7 +900,7 @@ j01 = -1.0*B2plusW_sq; - //Eq. 5: Qdotn + B^2(1+vsq)/2 - QdotBsq/2/W^2 + vsq W + rho_0(vsq) + u = 0 + //Eq. 5: -Qdotn - B^2(1+vsq)/2 + QdotBsq/2/W^2 - vsq W - rho_0(vsq) - u = 0 //rho0 = D * sqrt(1-vsq) resid[1] = -lgp->Qdotn - lgp->half_Bsq*(1.0+vsq) + 0.5*QB2Winv2 - vsq*W - rho0 - u;

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 323)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:10 PM Modified: /trunk/src/ GRHydro_Con2Prim.F90 Log: * improve error handling/reporting in hot C2P * set up directory for performance workshop Patch by Christian Ott. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +86 -79 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2012-05-13 19:08:58 UTC (rev 322) +++ trunk/src/GRHydro_Con2Prim.F90 2012-05-13 19:10:12 UTC (rev 323) @@ -194,7 +194,7 @@ z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) else - call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),& + call Con2Prim_pt_hot(cctk_iteration,i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),& scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),vup(i,j,k,1),& vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& press(i,j,k),w_lorentz(i,j,k), & @@ -207,7 +207,7 @@ ! accuracy requirement; if it fails again, we abort GRHydro_C2P_failed(i,j,k) = 0 local_perc_ptol = GRHydro_perc_ptol*100.0d0 - call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),& + call Con2Prim_pt_hot(cctk_iteration,i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),& scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),& vup(i,j,k,1),vup(i,j,k,2), vup(i,j,k,3),eps(i,j,k),& temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), & @@ -563,7 +563,8 @@ end subroutine Con2Prim_pt -subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, ye_con, rho, velx, vely, & +subroutine Con2Prim_pt_hot(cctk_iteration, ii,jj,kk,handle, dens, & + sx, sy, sz, tau, ye_con, rho, velx, vely, & velz, epsilon, temp, ye, press, w_lorentz, uxx, uxy, uxz, uyy, & uyz, uzz, det, x, y, z, r, epsnegative, GRHydro_rho_min, pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed, local_perc_ptol) @@ -577,7 +578,7 @@ y, z, r, GRHydro_rho_min CCTK_REAL temp, ye CCTK_REAL s2, c0, c1, c2, c3, c4, f, df, ftol, v2, w, vlowx, vlowy, vlowz - CCTK_INT ii,jj,kk,count, i, handle, GRHydro_reflevel + CCTK_INT cctk_iteration, ii,jj,kk,count, i, handle, GRHydro_reflevel CCTK_REAL GRHydro_C2P_failed CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, & w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin @@ -649,32 +650,34 @@ enddo write(warnline,"(A30,i5)") "Iterations of heat injection: ",nf call CCTK_WARN(failinfomode,warnline) - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + write(warnline,"(5i5,1P10E15.6)") cctk_iteration,GRHydro_Reflevel,ii,jj,kk,x,y,z call CCTK_WARN(failinfomode,warnline) - if(nf.gt.nfudgemax.and.GRHydro_c2p_reset_eps_tau_hot_eos.ne.1) then - call CCTK_WARN(failinfomode,"EOS error in c2p 1: injected heat too many times") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) - call CCTK_WARN(failwarnmode,"Aborting!!!") - else - call CCTK_WARN(failinfomode,"EOS error in c2p 1: LAST RESORT -- reset eps and tau based on old temp") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - keytemp = 1 - call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& - rho,epsilon,temp,ye,xpress,keyerr,anyerr) - tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - keytemp = 0 + if(nf.gt.nfudgemax) then + if(GRHydro_c2p_reset_eps_tau_hot_eos.ne.1) then + call CCTK_WARN(failinfomode,"EOS error in c2p 0: injected heat too many times") + write(warnline,"(5i5,1P10E15.6)") cctk_iteration,GRHydro_reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel, iteration: ", GRHydro_reflevel, cctk_iteration + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + else + call CCTK_WARN(failinfomode,"EOS error in c2p 0: LAST RESORT -- reset eps and tau based on old temp") + write(warnline,"(5i5,1P10E15.6)") cctk_iteration,GRHydro_reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + keytemp = 1 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,epsilon,temp,ye,xpress,keyerr,anyerr) + tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + keytemp = 0 + endif endif else call CCTK_WARN(failinfomode,"EOS error in c2p 0") @@ -740,32 +743,34 @@ enddo write(warnline,"(A30,i5)") "Iterations of heat injection: ",nf call CCTK_WARN(failinfomode,warnline) - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + write(warnline,"(5i5,1P10E15.6)") cctk_iteration,GRHydro_Reflevel,ii,jj,kk,x,y,z call CCTK_WARN(failinfomode,warnline) - if(nf.gt.nfudgemax.and.GRHydro_c2p_reset_eps_tau_hot_eos.ne.1) then - call CCTK_WARN(failinfomode,"EOS error in c2p 1: injected heat too many times") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) - call CCTK_WARN(failwarnmode,"Aborting!!!") - else - call CCTK_WARN(failinfomode,"EOS error in c2p 1: LAST RESORT -- reset eps and tau based on old temp") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - keytemp = 1 - call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& - rho,epsilon,temp,ye,xpress,keyerr,anyerr) - tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - keytemp = 0 + if(nf.gt.nfudgemax) then + if(GRHydro_c2p_reset_eps_tau_hot_eos.ne.1) then + call CCTK_WARN(failinfomode,"EOS error in c2p 1: injected heat too many times") + write(warnline,"(5i5,1P10E15.6)") cctk_iteration,GRHydro_reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel, iteration: ", GRHydro_reflevel, cctk_iteration + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + else + call CCTK_WARN(failinfomode,"EOS error in c2p 1: LAST RESORT -- reset eps and tau based on old temp") + write(warnline,"(5i5,1P10E15.6)") cctk_iteration,GRHydro_reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + keytemp = 1 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,epsilon,temp,ye,xpress,keyerr,anyerr) + tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + keytemp = 0 + endif endif else call CCTK_WARN(failinfomode,"EOS error in c2p 1") @@ -849,32 +854,34 @@ enddo write(warnline,"(A30,i5)") "Iterations of heat injection: ",nf call CCTK_WARN(failinfomode,warnline) - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + write(warnline,"(5i5,1P10E15.6)") cctk_iteration,GRHydro_Reflevel,ii,jj,kk,x,y,z call CCTK_WARN(failinfomode,warnline) - if(nf.gt.nfudgemax.and.GRHydro_c2p_reset_eps_tau_hot_eos.ne.1) then - call CCTK_WARN(failinfomode,"EOS error in c2p 1: injected heat too many times") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) - call CCTK_WARN(failwarnmode,"Aborting!!!") - else - call CCTK_WARN(failinfomode,"EOS error in c2p 1: LAST RESORT -- reset eps and tau based on old temp") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - keytemp = 1 - call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& - rho,epsilon,temp,ye,xpress,keyerr,anyerr) - tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - keytemp = 0 + if(nf.gt.nfudgemax) then + if(GRHydro_c2p_reset_eps_tau_hot_eos.ne.1) then + call CCTK_WARN(failinfomode,"EOS error in c2p 2: injected heat too many times") + write(warnline,"(5i5,1P10E15.6)") cctk_iteration,GRHydro_reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel, iteration: ", GRHydro_reflevel, cctk_iteration + call CCTK_WARN(failinfomode,warnline) + call CCTK_WARN(failwarnmode,"Aborting!!!") + else + call CCTK_WARN(failinfomode,"EOS error in c2p 2: LAST RESORT -- reset eps and tau based on old temp") + write(warnline,"(5i5,1P10E15.6)") cctk_iteration,GRHydro_reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + keytemp = 1 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + rho,epsilon,temp,ye,xpress,keyerr,anyerr) + tau = ( (rho + rho*epsilon) + xpress ) * w_lorentz**2 - xpress - rho*w_lorentz + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + keytemp = 0 + endif endif else call CCTK_WARN(failinfomode,"EOS error in c2p 2")

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 322)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:08 PM Modified: /trunk/src/ GRHydro_Con2Prim.F90 Log: Error check for C2P_failed mask greater or equal 1, do not abort code if in the atmosphere though. patch by Philipp Moesta File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +4 -1 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2012-05-13 19:03:31 UTC (rev 321) +++ trunk/src/GRHydro_Con2Prim.F90 2012-05-13 19:08:58 UTC (rev 322) @@ -1980,7 +1980,7 @@ do j = 1, ny do i = 1, nx - if (GRHydro_C2P_failed(i,j,k) == 1) then + if (GRHydro_C2P_failed(i,j,k) .ge. 1) then if(con2prim_oct_hack.ne.0.and.& (x(i,j,k) .lt. -1.0d-12 .or.& @@ -1989,6 +1989,9 @@ cycle endif + if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + cycle + end if !$OMP CRITICAL call CCTK_WARN(1,'Con2Prim failed; stopping the code.')

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 321)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:03 PM Modified: /trunk/src/ GRHydro_Con2PrimM_pt.c Log: Consistent atmosphere treatment. patch by Philipp Moesta File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM_pt.c Delta lines: +13 -2 =================================================================== --- trunk/src/GRHydro_Con2PrimM_pt.c 2012-05-13 19:01:37 UTC (rev 320) +++ trunk/src/GRHydro_Con2PrimM_pt.c 2012-05-13 19:03:31 UTC (rev 321) @@ -45,6 +45,7 @@ #include <complex.h> #include "cctk.h" +#include "cctk_Parameters.h" /* Set this to be 1 if you want debug output */ #define DEBUG_CON2PRIMM (0) @@ -215,8 +216,10 @@ CCTK_REAL detg = (*det); CCTK_REAL sqrt_detg = sqrt(detg); CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; - CCTK_INT i,j, i_increase ; - + CCTK_INT i,j, i_increase; + + DECLARE_CCTK_PARAMETERS; + struct LocGlob lg; gammaeos = *gamma_eos; @@ -395,7 +398,15 @@ *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; *velz = g_o_WBsq * ( usz + QdB_o_W*(*Bz) ) ; + if (*rho <= rho_abs_min*(1.0+GRHydro_atmo_tolerance) ) { + *rho = rho_abs_min; + *velx = 0.0; + *vely = 0.0; + *velz = 0.0; + *w_lorentz = 1.0; + } + #if(DEBUG_CON2PRIMM) fprintf(stdout,"rho = %26.16e \n",*rho ); fprintf(stdout,"epsilon = %26.16e \n",*epsilon );

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 320)
by roland.haas＠physics.gatech.edu 13 May '12

13 May '12

User: rhaas Date: 2012/05/13 02:01 PM Modified: /trunk/src/ GRHydro_Con2PrimM.F90 Log: Error checking for C2P_failed > 1, not only equals 1 and check for atmosphere using density after Con2PrimM_pt was performed. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +28 -1 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2012-05-11 22:46:36 UTC (rev 319) +++ trunk/src/GRHydro_Con2PrimM.F90 2012-05-13 19:01:37 UTC (rev 320) @@ -243,6 +243,34 @@ #endif end if + + if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) .or. GRHydro_C2P_failed(i,j,k) .ge. 1) then + + b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ & + 2.0*(gxy(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,2)+gxz(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,3)+ & + gyz(i,j,k)*Bvec(i,j,k,2)*Bvec(i,j,k,3)) + + dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + rho(i,j,k) = GRHydro_rho_min + 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) + + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + + ! w_lorentz=1, so the expression for tau reduces to: + + !!$ tau does need to take into account the existing B-field + !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens [Press drops out] + tau(i,j,k) = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) + + end if + end do end do end do @@ -487,7 +515,6 @@ !$OMP END CRITICAL endif endif - end do end do end do

1 0

[svn:einsteintoolkit] AHFinderDirect/trunk/src/include/ (Rev. 1564)
by schnetter＠cct.lsu.edu 11 May '12

11 May '12

User: eschnett Date: 2012/05/11 07:40 PM Modified: /trunk/src/include/ config.h Log: Support CCTK_REAL4 File Changes: Directory: /trunk/src/include/ ============================== File [modified]: config.h Delta lines: +7 -0 =================================================================== --- trunk/src/include/config.h 2012-04-28 19:12:20 UTC (rev 1563) +++ trunk/src/include/config.h 2012-05-12 00:40:26 UTC (rev 1564) @@ -54,8 +54,15 @@ * FIXME: this assumes fp == C 'double' * CCTK_REAL_PRECISION_{4,8,16} are helpful, but not quite enough) */ +#include <cctk.h> #undef FP_IS_FLOAT +#undef FP_IS_DOUBLE +#ifdef CCTK_REAL_PRECISION_4 +#define FP_IS_FLOAT +#endif +#ifdef CCTK_REAL_PRECISION_8 #define FP_IS_DOUBLE +#endif #if defined(FP_IS_FLOAT) #define FP_SCANF_FORMAT "%f"

1 0

New changesets in Carpet Mercurial
by Carpet mercurial repositories 11 May '12

11 May '12

changeset: 3549:e64eee13a84d tag: tip user: Roland Haas <rhaas(a)caltech.edu> date: Wed May 09 16:56:06 2012 -0700 summary: CarpetLib: disable debug output

1 0

[svn:einsteintoolkit] ADMBase/trunk/ (Rev. 66)
by schnetter＠cct.lsu.edu 11 May '12

11 May '12

User: eschnett Date: 2012/05/11 06:15 PM Modified: /trunk/ param.ccl /trunk/src/ InitSymBound.c Log: Allow different boundary conditions in ADMBase Introduce a parameter ADMBase::admbase_boundary_condition selection which boundary condition ADMBase applies. File Changes: Directory: /trunk/src/ ====================== File [modified]: InitSymBound.c Delta lines: +22 -10 =================================================================== --- trunk/src/InitSymBound.c 2011-10-25 18:44:57 UTC (rev 65) +++ trunk/src/InitSymBound.c 2012-05-11 23:15:39 UTC (rev 66) @@ -123,9 +123,22 @@ } /* A macro for selecting boundary conditions and checking for errors */ -#define ADMBASE_BC(NAME, STENCIL) \ - {if (Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, STENCIL, -1, NAME, "flat") < 0) \ - CCTK_WARN(0, "Failed to register BC for "NAME"!");} +static void select_bc(cGH const* const cctkGH, char const* const groupname) +{ + DECLARE_CCTK_PARAMETERS; + + int const ierr = + Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, + groupname, admbase_boundary_condition); + if (ierr < 0) + { + CCTK_VWarn (CCTK_WARN_ABORT, + __LINE__, __FILE__, CCTK_THORNSTRING, + "Failed to select boundary conditions for group %s", + groupname); + } +} + /* Select boundary conditions on ADMBase variables */ void ADMBase_Boundaries(CCTK_ARGUMENTS) { @@ -153,23 +166,22 @@ if (CCTK_EQUALS(evolution_method, "none" ) || CCTK_EQUALS(evolution_method, "static")) { - ADMBASE_BC("ADMBase::metric", stencil); - ADMBASE_BC("ADMBase::curv", stencil); + select_bc(cctkGH, "ADMBase::metric"); + select_bc(cctkGH, "ADMBase::curv"); } if (CCTK_EQUALS(lapse_evolution_method, "static")) - ADMBASE_BC("ADMBase::lapse", stencil); + select_bc(cctkGH, "ADMBase::lapse"); if (!CCTK_EQUALS(initial_dtlapse, "none") && CCTK_EQUALS(dtlapse_evolution_method, "static")) - ADMBASE_BC("ADMBase::dtlapse", stencil); + select_bc(cctkGH, "ADMBase::dtlapse"); if (!CCTK_EQUALS(initial_shift, "none") && CCTK_EQUALS(shift_evolution_method, "static")) - ADMBASE_BC("ADMBase::shift", stencil); + select_bc(cctkGH, "ADMBase::shift"); if (!CCTK_EQUALS(initial_dtshift, "none") && CCTK_EQUALS(dtshift_evolution_method, "static")) - ADMBASE_BC("ADMBase::dtshift", stencil); + select_bc(cctkGH, "ADMBase::dtshift"); } -#undef ADMBase_BC Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +7 -0 =================================================================== --- trunk/param.ccl 2011-10-25 18:44:57 UTC (rev 65) +++ trunk/param.ccl 2012-05-11 23:15:39 UTC (rev 66) @@ -61,6 +61,13 @@ +STRING admbase_boundary_condition "Boundary condition for ADMBase variables" +{ + "" :: "must be a registered boundary condition" +} "flat" + + + KEYWORD metric_type "The semantics of the metric variables (physical, static conformal, etc)" { "physical" :: "metric and extrinsic curvature are the physical ones"

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 319)
by schnetter＠cct.lsu.edu 11 May '12

11 May '12

User: eschnett Date: 2012/05/11 05:46 PM Modified: /trunk/src/ GRHydro_Con2Prim.F90, GRHydro_Eigenproblem_Marquina.F90, GRHydro_FluxSplit.F90, GRHydro_PPM.F90 Log: Support real*16 (and real*4) in GRHydro Support using precisions other than real*8 in GRHydro by removing all explicit references to double precision functions and constants, and using type-generic functions and constants instead. In particular, use "one" and "half" as constants in some places, and use "abs", "max" etc. instead of "dabs", "dmax" etc. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: +14 -12 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2012-05-11 22:17:18 UTC (rev 318) +++ trunk/src/GRHydro_Con2Prim.F90 2012-05-11 22:46:36 UTC (rev 319) @@ -1026,6 +1026,8 @@ DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS + CCTK_REAL, parameter :: half = 0.5d0 + integer :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: uxxl, uxyl, uxzl, uyyl, uyzl, uzzl,& uxxr, uxyr, uxzr, uyyr, uyzr, uzzr, pmin, epsmin @@ -1135,9 +1137,9 @@ velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),& pressminus(i,j,k),w_lorentzminus(i,j,k),& uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& - x(i,j,k)-0.5d0*CCTK_DELTA_SPACE(1),& - y(i,j,k)-0.5d0*CCTK_DELTA_SPACE(2), & - z(i,j,k)-0.5d0*CCTK_DELTA_SPACE(3),r(i,j,k),& + x(i,j,k)-half*CCTK_DELTA_SPACE(1),& + y(i,j,k)-half*CCTK_DELTA_SPACE(2), & + z(i,j,k)-half*CCTK_DELTA_SPACE(3),r(i,j,k),& epsnegative,GRHydro_rho_min,pmin, epsmin, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) if (epsnegative) then !$OMP CRITICAL @@ -1149,9 +1151,9 @@ velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),& pressminus(i,j,k),w_lorentzminus(i,j,k),& uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& - x(i,j,k)-0.5d0*CCTK_DELTA_SPACE(1),& - y(i,j,k)-0.5d0*CCTK_DELTA_SPACE(2), & - z(i,j,k)-0.5d0*CCTK_DELTA_SPACE(3),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) + x(i,j,k)-half*CCTK_DELTA_SPACE(1),& + y(i,j,k)-half*CCTK_DELTA_SPACE(2), & + z(i,j,k)-half*CCTK_DELTA_SPACE(3),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) end if if (epsminus(i,j,k) .lt. 0.0d0) then @@ -1182,9 +1184,9 @@ velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),& pressplus(i,j,k),w_lorentzplus(i,j,k),& uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& - x(i,j,k)+0.5d0*CCTK_DELTA_SPACE(1),& - y(i,j,k)+0.5d0*CCTK_DELTA_SPACE(2), & - z(i,j,k)+0.5d0*CCTK_DELTA_SPACE(3),r(i,j,k),& + x(i,j,k)+half*CCTK_DELTA_SPACE(1),& + y(i,j,k)+half*CCTK_DELTA_SPACE(2), & + z(i,j,k)+half*CCTK_DELTA_SPACE(3),r(i,j,k),& epsnegative,GRHydro_rho_min,pmin, epsmin, GRHydro_reflevel,GRHydro_C2P_failed(i,j,k)) if (epsnegative) then !$OMP CRITICAL @@ -1196,9 +1198,9 @@ velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),& pressplus(i,j,k),w_lorentzplus(i,j,k),& uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& - x(i,j,k)+0.5d0*CCTK_DELTA_SPACE(1),& - y(i,j,k)+0.5d0*CCTK_DELTA_SPACE(2), & - z(i,j,k)+0.5d0*CCTK_DELTA_SPACE(3),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) + x(i,j,k)+half*CCTK_DELTA_SPACE(1),& + y(i,j,k)+half*CCTK_DELTA_SPACE(2), & + z(i,j,k)+half*CCTK_DELTA_SPACE(3),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) end if if (epsplus(i,j,k) .lt. 0.0d0) then File [modified]: GRHydro_Eigenproblem_Marquina.F90 Delta lines: +9 -9 =================================================================== --- trunk/src/GRHydro_Eigenproblem_Marquina.F90 2012-05-11 22:17:18 UTC (rev 318) +++ trunk/src/GRHydro_Eigenproblem_Marquina.F90 2012-05-11 22:46:36 UTC (rev 319) @@ -200,11 +200,11 @@ !!$ FINAL - lam1 = dmax1(dabs(lam1l),dabs(lam1r)) + lam1 = max(abs(lam1l),abs(lam1r)) lam2 = lam1 lam3 = lam1 - lamp = dmax1(dabs(lampl),dabs(lampr)) - lamm = dmax1(dabs(lamml),dabs(lammr)) + lamp = max(abs(lampl),abs(lampr)) + lamm = max(abs(lamml),abs(lammr)) !!$ lam(1) = lamm !!$ lam(2) = lam1 @@ -918,11 +918,11 @@ !!$ FINAL - lam1 = dmax1(dabs(lam1l),dabs(lam1r)) + lam1 = max(abs(lam1l),abs(lam1r)) lam2 = lam1 lam3 = lam1 - lamp = dmax1(dabs(lampl),dabs(lampr)) - lamm = dmax1(dabs(lamml),dabs(lammr)) + lamp = max(abs(lampl),abs(lampr)) + lamm = max(abs(lamml),abs(lammr)) !!$ lam(1) = lamm !!$ lam(2) = lam1 @@ -1610,11 +1610,11 @@ !!$ FINAL - lam1 = dmax1(dabs(lam1l),dabs(lam1r)) + lam1 = max(abs(lam1l),abs(lam1r)) lam2 = lam1 lam3 = lam1 - lamp = dmax1(dabs(lampl),dabs(lampr)) - lamm = dmax1(dabs(lamml),dabs(lammr)) + lamp = max(abs(lampl),abs(lampr)) + lamm = max(abs(lamml),abs(lammr)) !!$ lam(1) = lamm !!$ lam(2) = lam1 File [modified]: GRHydro_FluxSplit.F90 Delta lines: +17 -15 =================================================================== --- trunk/src/GRHydro_FluxSplit.F90 2012-05-11 22:17:18 UTC (rev 318) +++ trunk/src/GRHydro_FluxSplit.F90 2012-05-11 22:46:36 UTC (rev 319) @@ -437,6 +437,8 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS + CCTK_REAL, parameter :: half = 0.5d0 + CCTK_INT :: i, nx, handle, ll CCTK_REAL, dimension(nx) :: gxx, gxy, gxz, gyy, gyz, gzz, & u, det, alp, beta, & @@ -493,21 +495,21 @@ !!$ Calculate the maximum eigenvalue and put it here call eigenproblem_leftright(handle, & - 0.5d0 * (rho (i) + rho (i+1)), & - 0.5d0 * (velx1 (i) + velx1 (i+1)), & - 0.5d0 * (vely1 (i) + vely1 (i+1)), & - 0.5d0 * (velz1 (i) + velz1 (i+1)), & - 0.5d0 * (eps (i) + eps (i+1)), & - 0.5d0 * (w_lorentz(i) + w_lorentz(i+1)), & - 0.5d0 * (gxx (i) + gxx (i+1)), & - 0.5d0 * (gxy (i) + gxy (i+1)), & - 0.5d0 * (gxz (i) + gxz (i+1)), & - 0.5d0 * (gyy (i) + gyy (i+1)), & - 0.5d0 * (gyz (i) + gyz (i+1)), & - 0.5d0 * (gzz (i) + gzz (i+1)), & - 0.5d0 * (u (i) + u (i+1)), & - 0.5d0 * (alp (i) + alp (i+1)), & - 0.5d0 * (beta (i) + beta (i+1)), & + half * (rho (i) + rho (i+1)), & + half * (velx1 (i) + velx1 (i+1)), & + half * (vely1 (i) + vely1 (i+1)), & + half * (velz1 (i) + velz1 (i+1)), & + half * (eps (i) + eps (i+1)), & + half * (w_lorentz(i) + w_lorentz(i+1)), & + half * (gxx (i) + gxx (i+1)), & + half * (gxy (i) + gxy (i+1)), & + half * (gxz (i) + gxz (i+1)), & + half * (gyy (i) + gyy (i+1)), & + half * (gyz (i) + gyz (i+1)), & + half * (gzz (i) + gzz (i+1)), & + half * (u (i) + u (i+1)), & + half * (alp (i) + alp (i+1)), & + half * (beta (i) + beta (i+1)), & lambda,& levec,& revec) File [modified]: GRHydro_PPM.F90 Delta lines: +26 -20 =================================================================== --- trunk/src/GRHydro_PPM.F90 2012-05-11 22:17:18 UTC (rev 318) +++ trunk/src/GRHydro_PPM.F90 2012-05-11 22:46:36 UTC (rev 319) @@ -66,6 +66,8 @@ DECLARE_CCTK_PARAMETERS + CCTK_REAL, parameter :: one = 1 + CCTK_INT :: handle,poly,nx CCTK_REAL :: dx CCTK_REAL, dimension(nx) :: rho,velx,vely,velz,eps @@ -112,10 +114,10 @@ !use_std_ppm = 1 #define STEEP(x,dx,dmx) \ - if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0.d0 ) then &&\ - dmx(i) = sign(1.d0, dx(i)) * \ - min(abs(dx(i)), 2.d0 * abs(x(i) - x(i-1)), \ - 2.d0 * abs(x(i+1) - x(i))) &&\ + if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0 ) then &&\ + dmx(i) = sign(one, dx(i)) * \ + min(abs(dx(i)), 2 * abs(x(i) - x(i-1)), \ + 2 * abs(x(i+1) - x(i))) &&\ else &&\ dmx(i) = 0.d0 &&\ end if @@ -205,7 +207,7 @@ D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ - D2aLim = sign(1.0d0, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ + D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ else &&\ @@ -478,7 +480,7 @@ end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(1.d0, -dpress(i)) + s=sign(one, -dpress(i)) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) if (abs(1.d0 - flatten) > 0.d0) then trivial_rp(i-1) = .false. @@ -543,8 +545,8 @@ D3aL = D2aC - D2aL &&\ D3aR = D2aRR - D2aR &&\ D3aLL = D2aL - D2aLL &&\ - if (sign(1.0d0, D2a) .eq. sign(1.0d0, D2aC) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aL) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aR)) then &&\ - D2aLim = sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ else &&\ D2aLim = 0 &&\ end if &&\ @@ -608,9 +610,9 @@ ! D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ ! D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ ! D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ -! if (D2a .ne. 0 .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aC) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aL) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aR)) then &&\ -! aplus(i) = a(i) + (aplus(i) -a(i)) * sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) / D2a &&\ -! aminus(i) = a(i) + (aminus(i)-a(i)) * sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) / D2a &&\ +! if (D2a .ne. 0 .and. sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ +! aplus(i) = a(i) + (aplus(i) -a(i)) * sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) / D2a &&\ +! aminus(i) = a(i) + (aminus(i)-a(i)) * sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) / D2a &&\ ! else &&\ ! aplus(i) = a(i) &&\ ! aminus(i) = a(i) &&\ @@ -637,8 +639,8 @@ D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - if (sign(1.0d0, D2a) .eq. sign(1.0d0, D2aC) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aL) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aR)) then &&\ - D2aLim = sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ else &&\ D2aLim = 0 &&\ end if &&\ @@ -898,6 +900,8 @@ DECLARE_CCTK_PARAMETERS + CCTK_REAL, parameter :: one = 1 + CCTK_INT :: nx CCTK_REAL :: dx CCTK_REAL, dimension(nx) :: rho,velx,vely,velz @@ -934,7 +938,7 @@ do i = 2, nx - 1 if( (tracer(i+1,itracer) - tracer(i,itracer)) * & (tracer(i,itracer) - tracer(i-1,itracer)) > 0.0d0 ) then - dmtracer(i,itracer) = sign(1.0d0,dtracer(i,itracer)) * & + dmtracer(i,itracer) = sign(one,dtracer(i,itracer)) * & min(abs(dtracer(i,itracer)), 2.0d0 * & abs(tracer(i,itracer) - tracer(i-1,itracer)), & 2.0d0 * abs(tracer(i+1,itracer) - tracer(i,itracer))) @@ -1111,6 +1115,8 @@ DECLARE_CCTK_PARAMETERS + CCTK_REAL, parameter :: one = 1 + CCTK_INT :: nx CCTK_REAL :: dx CCTK_REAL, dimension(nx) :: rho,velx,vely,velz @@ -1151,7 +1157,7 @@ do i = 2, nx - 1 if( (Y_e(i+1) - Y_e(i)) * & (Y_e(i) - Y_e(i-1)) > 0.0d0 ) then - dmY_e(i) = sign(1.0d0,dY_e(i)) * & + dmY_e(i) = sign(one,dY_e(i)) * & min(abs(dY_e(i)), 2.0d0 * & abs(Y_e(i) - Y_e(i-1)), & 2.0d0 * abs(Y_e(i+1) - Y_e(i))) @@ -1194,7 +1200,7 @@ D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ - alim = 0.5d0*(a(i)+a(i+1)) - sign(1.0d0, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ + alim = 0.5d0*(a(i)+a(i+1)) - sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ else &&\ alim = 0.5d0*(a(i)+a(i+1)) &&\ end if &&\ @@ -1292,8 +1298,8 @@ D3aL = D2aC - D2aL &&\ D3aR = D2aRR - D2aR &&\ D3aLL = D2aL - D2aLL &&\ - if (sign(1.0d0, D2a) .eq. sign(1.0d0, D2aC) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aL) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aR)) then &&\ - D2aLim = sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ else &&\ D2aLim = 0 &&\ end if &&\ @@ -1358,8 +1364,8 @@ D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ - if (sign(1.0d0, D2a) .eq. sign(1.0d0, D2aC) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aL) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aR)) then &&\ - D2aLim = sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ else &&\ D2aLim = 0 &&\ end if &&\

1 0

← Newer
1
...
9
10
11
12
13
14
15
...
19
Older →

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits May 2012