[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 319)

schnetter at cct.lsu.edu schnetter at cct.lsu.edu
Fri May 11 17:46:37 CDT 2012


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                                                    &&\



More information about the Commits mailing list