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*(-16+y*(48/5-32/15*y)))+1/15/y))/r_soft

Or

(-16/5+32/3*y^2 -16*y^3+48/5*y^4 -32/15*y^5 +1/15/y)/r_soft

Multiplying by y and changing the sign

(-1/15+16/5y-32/3*y^3 +16*y^4-48/5*y^5 +32/15*y^6)/r

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 + 48/15*y - 32/3*y^3 + 16*y^4 - 48/5*y^5 + 32/15 *y^6 )

   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 big memory test: succeeded with pseducolor plot px/rho

Calculation using Expressions with scalar/vector mesh variables uses a lot computing resources..

Visit 2.12.3 on Bluehive bhc0198 with ~500GB memory crashes when plot "En/E_pot_rel_1"

AddPlot("Pseudocolor","En/E_pot_rel_1")
DrawPlots()
Query("Weighted Variable Sum")
erg_Epot_gas_par1_box = GetQueryOutputValue()

Definition structure of variable En/E_pot_rel_1

Temp/x1sc: 0.0000
Temp/y1sc: 0.0000
Temp/z1sc: 0.0000

=>

Pos/rVec: coord(Mesh)
Part/rVec_part1: {<Temp/x1sc>,<Temp/y1sc>,<Temp/z1sc>}

=>

Pos/rVec_rel_1: <Pos/rVec> -<Part/rVec_part1>

=>

Pos/r_rel_1: magnitude(<Pos/rVec_rel_1>)
Param/r_soft: 167600000000.000000
Param/G_Newton: 6.674e-8
Param/M_1: 729799999999999961087848627568640.000000

=>

Temp/u_1: <Pos/r_rel_1>/<Param/r_soft>
Pot/Phi0_rel_1: -1.0*<Param/G_Newton>*<Param/M_1>/<Pos/r_rel_1>
Temp/spline_A_1: -16./3*<Temp/u_1>^3 +48./5*<Temp/u_1>^5 -32./5*<Temp/u_1>^6 +14./5*<Temp/u_1>
Temp/spline_B_1:-1./15 -32./3*<Temp/u_1>^3 +16.*<Temp/u_1>^4 -48./5*<Temp/u_1>^5 +32./15*<Temp/u_1>^6 +48./15*<Temp/u_1>

=>

Temp/Phi_rel_1_A: if(lt(<Temp/u_1>,0.5), 1.0, 0.0)*<Pot/Phi0_rel_1>*<Temp/spline_A_1>
Temp/Phi_rel_1_B: if(and(ge(<Temp/u_1>,0.5),lt(<Temp/u_1>,1.)), 1.0, 0.0)*<Pot/Phi0_rel_1>*<Temp/spline_B_1>
Temp/Phi_rel_1_C: if(ge(<Temp/u_1>,1.), 1.0, 0.0)*<Pot/Phi0_rel_1>

=>

Pot/Phi_rel_1: <Temp/Phi_rel_1_A> +<Temp/Phi_rel_1_B> +<Temp/Phi_rel_1_C>

=>

En/E_pot_rel_1: rho*<Pot/Phi_rel_1>