Changes between Version 1 and Version 2 of u/bliu/processing

Timestamp:

06/28/18 12:35:57 (7 years ago)

Author:

Baowei Liu

Comment:

—

u/bliu/processing

@@ +1 @@
+== Compare SplinePotential(r,r_soft) with RelativePotential_1(mass,pos,r_soft) ==
+
+1. When r/r_soft > 0.5 SplinePotential
+
+{{{
+(-16/5+(y^2*(32/3+y*(48/5-32/15*y)))+1/15*y)/r_soft
+}}}
+Or
+
+{{{
+(-16/5+32/3*y^2 + 48/5*y^3 -32/15*y^4 +1/15*y)/r_soft
+}}}
+
+{{{
+FUNCTION SplinePotential(r,r_soft)
+      REAL(KIND=qPREC) :: r, r_soft,y
+      REAL(KIND=qPREC) :: SplinePotential
+      REAL(KIND=qPREC), PARAMETER :: &
+           GSC0=-14d0/5d0, &
+           GSC1=16d0/3d0, &
+           GSC2=-48d0/5d0, &
+           GSC3=32d0/5d0, &
+           GSC4=32d0/3d0, &
+           GSC5=-16d0, &
+           GSC6=48d0/5d0, &
+           GSC7=-32d0/15d0, &
+           GSC8=1d0/15d0, &
+           GSC9=-16d0/5d0
+      y=r/r_soft
+      IF (GravityDim == 3) THEN
+         IF (r >= r_soft) THEN
+            SplinePotential=-1d0/r
+         ELSE
+            y=r/r_soft
+            IF (r >= half*r_soft) THEN
+               SplinePotential=(GSC9+(y**2*(GSC4+y*(GSC5+y*(GSC6+y*GSC7)))+GSC8/y )) / r_soft
+            ELSE
+               SplinePotential=(GSC0 + (y**2*(GSC1+y**2*(GSC2+y*GSC3)))) / r_soft
+            END IF
+         END IF
+      ELSE
+         SplinePotential=PlummerPotential2D(r,r_soft)
+      END IF
+   END FUNCTION SplinePotential
+
+}}}
+
+2. When r/r_soft > 0.5 RelativePotential_1(mass,pos,r_soft) See below for the script from Luke
+
+{{{
+-1/15 - 32/3*y^3 + 16*y^4 - 48/5*y^5 + 32/15 *y^6 + 48/15*y)
+}}}
+
+
+{{{
+   FUNCTION RelativePotential_1(mass,pos,r_soft)
+      REAL(KIND=qPREC), DIMENSION(3) :: pos
+      REAL(KIND=qPREC) :: RelativePotential_1
+      REAL(KIND=qPREC) :: mass, r_soft, r, r2
+      REAL(KIND=qPREC) :: u_1, phi0_rel_1, spline_A_1, spline_B_1
+      REAL(KIND=qPREC) :: phi_rel_1_A, phi_rel_1_B, phi_rel_1_C
+
+      r2=sum(pos(1:nDim)**2)
+      r=sqrt(r2)
+
+      u_1=r/r_soft
+
+      phi0_rel_1=-1d0*mass*ScaleGrav/r
+      spline_A_1=-16d0/3d0*u_1**3 + 48d0/5d0*u_1**3 - 32d0/5d0*u_1**6 + 14d0/5d0*u_1
+      spline_B_1=-1d0/15d0 - 32d0/3d0*u_1**3 + 16d0*u_1**4 - 48d0/5d0*u_1**5 + 32d0/15d0*u_1**6 + 48d0/15d0*u_1
+
+      phi_rel_1_A=0d0
+      phi_rel_1_B=0d0
+      phi_rel_1_C=0d0
+      if(u_1<0.5) then
+         phi_rel_1_A=phi0_rel_1*spline_A_1
+      else if(u_1>1.0) then
+         phi_rel_1_C=phi0_rel_1
+      else
+         phi_rel_1_B=phi0_rel_1*spline_B_1
+      end if
+
+      RelativePotential_1=(phi_rel_1_A+phi_rel_1_B+phi_rel_1_C)
+
+   END FUNCTION RelativePotential_1
+
+}}}
+
+
+
+
+========
 1. Tests: with data Visit/2.13.0 failed opening data under /public/lchamand/Yisheng_chombo/out_143
   small memory test with Visit/2.12.3: succeeded with pseuducolor plot  rho, failed with px/rho