NYU-CAL · ajdittmann · Jun 8, 2023 · Jul 12, 2023 · Jul 12, 2023 · Aug 27, 2023
diff --git a/Hydro/euler.c b/Hydro/euler.c
@@ -146,7 +146,7 @@ void cons2prim( const double * cons , double * prim , const double * x ,
    }
 
    if( isothermal ){
-      double cs2 = get_cs2( x );
+      double cs2 = get_cs2( x, rho );
       Pp = cs2*rho/gamma_law;
    }
 

diff --git a/Hydro/euler_altvisc.c b/Hydro/euler_altvisc.c
@@ -117,7 +117,7 @@ void cons2prim( const double * cons , double * prim , const double * x , double
 
    if( Pp  < PRE_FLOOR*rho ) Pp = PRE_FLOOR*rho;
    if( isothermal ){
-      double cs2 = get_cs2( x );
+      double cs2 = get_cs2( x, rho );
       Pp = cs2*rho/gamma_law;
    }
 

diff --git a/Hydro/euler_cart.c b/Hydro/euler_cart.c
@@ -110,7 +110,7 @@ void cons2prim( const double * cons , double * prim , const double * x , double
 
    if( Pp  < PRE_FLOOR*rho ) Pp = PRE_FLOOR*rho;
    if( isothermal ){
-      double cs2 = get_cs2( x );
+      double cs2 = get_cs2( x, rho );
       Pp = cs2*rho/gamma_law;
    }
 

diff --git a/Hydro/euler_cartInterp.c b/Hydro/euler_cartInterp.c
@@ -168,7 +168,7 @@ void cons2prim( const double * cons , double * prim , const double * x ,
    }
 
    if( isothermal ){
-      double cs2 = get_cs2( x );
+      double cs2 = get_cs2( x, rho );
       Pp = cs2*rho/gamma_law;
    }
 

diff --git a/Hydro/euler_cart_cons.c b/Hydro/euler_cart_cons.c
@@ -122,7 +122,7 @@ void cons2prim( const double * cons , double * prim , const double * x , double
 
    if( Pp  < PRE_FLOOR*rho ) Pp = PRE_FLOOR*rho;
    if( isothermal ){
-      double cs2 = get_cs2( x );
+      double cs2 = get_cs2( x, rho );
       Pp = cs2*rho/gamma_law;
    }
 

diff --git a/Hydro/euler_pro.c b/Hydro/euler_pro.c
@@ -137,7 +137,7 @@ void cons2prim(const double *cons, double *prim, const double *x, double dV,
         Pp = PRE_FLOOR*rho;
 
     if( isothermal ){
-        double cs2 = get_cs2( x );
+        double cs2 = get_cs2( x, rho );
         Pp = cs2*rho/gamma_law;
     }
 

diff --git a/Hydro/euler_sph.c b/Hydro/euler_sph.c
@@ -124,7 +124,7 @@ void cons2prim( const double * cons , double * prim , const double * x , double
 
    if( Pp  < PRE_FLOOR*rho ) Pp = PRE_FLOOR*rho;
    if( isothermal ){
-      double cs2 = get_cs2( x );
+      double cs2 = get_cs2( x, rho );
       Pp = cs2*rho/gamma_law;
    }
 

diff --git a/Hydro/greuler.c b/Hydro/greuler.c
@@ -652,7 +652,7 @@ void cons2prim_solve_isothermal(const double *cons, double *prim, const double *
     double s2 = 0.0;
     double Us = 0.0;
 
-    double cs2N = get_cs2(x);
+    double cs2N = get_cs2(x, 0.0);
     double P_o_rho = cs2N / gamma_law;
     double h = 1.0 + gamma_law * P_o_rho / (gamma_law-1.0);
     //double P_o_rhoh = cs2N / gamma_law;

diff --git a/Hydro/mhd.c b/Hydro/mhd.c
@@ -169,7 +169,7 @@ void cons2prim( const double * cons , double * prim , const double * x , double
 
    if( Pp  < PRE_FLOOR*rho ) Pp = PRE_FLOOR*rho;
    if( isothermal ){
-      double cs2 = get_cs2( x );
+      double cs2 = get_cs2( x, rho );
       Pp = cs2*rho/gamma_law;
    }
 

diff --git a/Hydro/mhd_cart.c b/Hydro/mhd_cart.c
@@ -151,7 +151,7 @@ void cons2prim( const double * cons , double * prim , const double * x , double
 
    if( Pp  < PRE_FLOOR*rho ) Pp = PRE_FLOOR*rho;
    if( isothermal ){
-      double cs2 = get_cs2( x );
+      double cs2 = get_cs2( x, rho );
       Pp = cs2*rho/gamma_law;
    }
 

diff --git a/Hydro/mhd_sph.c b/Hydro/mhd_sph.c
@@ -166,7 +166,7 @@ void cons2prim( const double * cons , double * prim , const double * x , double
 
    if( Pp  < PRE_FLOOR*rho ) Pp = PRE_FLOOR*rho;
    if( isothermal ){
-      double cs2 = get_cs2( x );
+      double cs2 = get_cs2( x, rho );
       Pp = cs2*rho/gamma_law;
    }
 

diff --git a/Initial/binary.c b/Initial/binary.c
@@ -1,4 +1,4 @@
-
+#include "../omega.h"
 #include "../paul.h"
 
 static double gam  = 0.0;
@@ -11,15 +11,13 @@ void setICparams( struct domain * theDomain ){
    Mach = theDomain->theParList.Disk_Mach;
 }
 
-double get_cs2( double );
-
 void initial( double * prim , double * x ){
 
    double r = x[0];
    double z = x[2];
 
    double sint = z/sqrt(r*r+z*z);
-   double cs2 = get_cs2( r );
+   double cs2 = get_cs2( x );
 
    double rho = 1.0*exp(-sint*sint*Mach*Mach);
    double Pp  = rho*cs2/gam;

diff --git a/Initial/bl.c b/Initial/bl.c
@@ -1,4 +1,5 @@
 #include "../paul.h"
+#include "../omega.h"
 
 static double M = 0.0;
 static double as = 0.0;
@@ -10,8 +11,6 @@ static double sig_atm = 0.0;
 static double mach = 0.0;
 static int isothermal = 0;
 
-double get_cs2(double r);
-
 double x_func_sin2(double x, void *args);
 double r_func_schw_sc(double r, void *args);
 double int_sim_ad(double a, double b, double tol, 
@@ -45,7 +44,7 @@ void initial( double * prim , double * x ){
 
    double cs20;
    if(isothermal)
-      cs20 = get_cs2(r);
+      cs20 = get_cs2(x);
    else
       cs20 = M/(r0 * mach*mach);
 

diff --git a/Initial/blBRS.c b/Initial/blBRS.c
@@ -23,7 +23,6 @@ void initial( double * prim , double * x ){
    get_rpz(x, rpz);
    double r = rpz[0];
 
-   double cs2 = get_cs2(x);
    double omega = 1.0/sqrt(r*r*r);
    double rho = 1.0;
    double A, B, C;
@@ -45,6 +44,7 @@ void initial( double * prim , double * x ){
      rho = exp(gam*Mach*Mach*(1./r + om0*om0*0.5*r*r + off2 - off1 - off3));
    }
 
+   double cs2 = get_cs2(x, rho);
    double Pp = rho*cs2/gam;
 
    double X = 0.0;

diff --git a/Initial/blBaro.c b/Initial/blBaro.c
@@ -0,0 +1,71 @@
+#include "../paul.h"
+#include "../geometry.h"
+#include "../omega.h"
+
+
+static double gam  = 0.0;
+static double Mach = 0.0;
+
+static double rbl = 0.0;
+static double dbl = 0.0;
+static double om0 = 0.0;
+static double n = 0.0;
+
+void setICparams( struct domain * theDomain ){
+   gam  = theDomain->theParList.Adiabatic_Index;
+   Mach = theDomain->theParList.Disk_Mach;
+   rbl = theDomain->theParList.initPar1;
+   dbl = theDomain->theParList.initPar2;
+   om0 = theDomain->theParList.initPar3;
+   n = theDomain->theParList.Cs2_Par;
+}
+
+void initial( double * prim , double * x ){
+   double rpz[3];
+   get_rpz(x, rpz);
+   double r = rpz[0];
+
+   double omega = 1.0/sqrt(r*r*r);
+   double rho = 1.0;
+   double A, B, C;
+   double off1, off2, off3;
+
+   A = rbl - dbl*0.5;
+   B = (1.0/pow(rbl + dbl*0.5, 1.5) - om0)/dbl;
+   C = om0 - B*A;
+   off1 = 1./(rbl+dbl*0.5) + C*C*0.5*pow(rbl+dbl*0.5, 2.0) + 2.0*C*B*pow(rbl+dbl*0.5, 3.0)/3.0+ B*B*0.25*pow(rbl+dbl*0.5, 4.0);
+   off2 = 1./(rbl-dbl*0.5) + C*C*0.5*pow(rbl-dbl*0.5, 2.0) + 2.0*C*B*pow(rbl-dbl*0.5, 3.0)/3.0+ B*B*0.25*pow(rbl-dbl*0.5, 4.0);
+   off3 = 1./(rbl - dbl*0.5) + om0*om0*0.5*pow(rbl-dbl*0.5, 2.0);
+
+   if (r <= rbl + 0.5*dbl) {
+     omega = C + r*B;
+     //rho = exp(gam*Mach*Mach*(1./r + C*C*0.5*r*r + 2.0*C*B*r*r*r/3.0 + B*B*0.25*r*r*r*r - off1));
+     rho = 1.0 + pow(gam*Mach*Mach*(1-1/n)*(1./r + C*C*0.5*r*r + 2.0*C*B*r*r*r/3.0 + B*B*0.25*r*r*r*r - off1), 1/(n-1) );
+   }
+   if (r < rbl - 0.5*dbl) {
+     omega = om0;
+     rho = 1.0 + pow(gam*Mach*Mach*(1-1/n)*(1./r + om0*om0*0.5*r*r + off2 - off1 - off3), 1/(n-1));
+     //rho = exp(gam*Mach*Mach*(1./r + om0*om0*0.5*r*r + off2 - off1 - off3));
+   }
+
+   double cs2 = get_cs2(x, rho);
+   double Pp = rho*cs2/gam;
+
+   double X = 0.0;
+   if( r > rbl - 0.5*dbl ) X = (r-(rbl-dbl*0.5))/dbl;
+   if( r > rbl + 0.5*dbl ) X = 1.0;
+
+   double vr = 0.0;
+   double Vrpz[3] = {vr, r*omega, 0.0};
+   double V[3];
+   get_vec_from_rpz(x, Vrpz, V);
+   get_vec_contravariant(x, V, V);
+
+   prim[RHO] = rho;
+   prim[PPP] = Pp;
+   prim[URR] = V[0];
+   prim[UPP] = V[1];
+   prim[UZZ] = V[2];
+   if( NUM_N>0 ) prim[NUM_C] = X;
+
+}
diff --git a/Initial/blob.c b/Initial/blob.c
@@ -71,7 +71,7 @@ void initial( double * prim , double * x )
 
     if(isothermal_flag)
     {
-        double cs2 = get_cs2(x);
+        double cs2 = get_cs2(x, rho);
         P = rho * cs2 / gam;
     }
 

diff --git a/Initial/cb.c b/Initial/cb.c
@@ -7,9 +7,6 @@ static double Mdot = 0.0;
 static double R0 = 0.0;
 static double rho0 = 0.0;
 
-double get_cs2(double *);
-double get_nu(double *, double *);
-
 void setICparams( struct domain * theDomain )
 {
     gam = theDomain->theParList.Adiabatic_Index;
@@ -24,7 +21,6 @@ void initial(double *prim, double *x)
     double r = x[0];
     double phi = x[1];
 
-    double cs2 = get_cs2(x);
     double om;
     if(r>1.0)
         om = pow(r, -1.5);
@@ -34,6 +30,7 @@ void initial(double *prim, double *x)
     double nu = get_nu(x, prim);
 
     double rho = Mdot/(3*M_PI*nu);
+    double cs2 = get_cs2(x, rho);
     if(nu == 0.0)
         rho = Mdot/(3*M_PI * 0.1*cs2/om);
     double v = -1.5*nu/r;

diff --git a/Initial/cbJdot.c b/Initial/cbJdot.c
@@ -35,8 +35,6 @@ void initial(double *prim, double *x)
     double R = r + 0.001;
     double phi = x[1];
 
-    double cs2 = get_cs2(x);
-
     double rho, efact, fth, dfth;
 
     //double om = 1.0;

diff --git a/Initial/cbMMLth.c b/Initial/cbMMLth.c
@@ -36,7 +36,6 @@ void initial(double *prim, double *x)
     double R = r + 0.05;
     double phi = x[1];
 
-    double cs2 = get_cs2(x);
 
     double rho, efact;
 
@@ -61,6 +60,8 @@ void initial(double *prim, double *x)
     rho = sig0*efact + epsfl;
     double drho = sig0*efact*xi*pow((R/redge),-xi)/R;
 
+    double cs2 = get_cs2(x, rho);
+
     double v = -1.5*nu/(R);
     double P = rho*cs2/gam;
 

diff --git a/Initial/cb_cc.c b/Initial/cb_cc.c
@@ -38,7 +38,7 @@ void initial(double *prim, double *x)
     double r = rpz[0];
     double phi = rpz[1];
 
-    double cs2 = get_cs2(x);
+    double cs2 = get_cs2(x, 1.0);
 
     double f = 1.0;
     if(Rout > 0.0)

diff --git a/Initial/cb_duff.c b/Initial/cb_duff.c
@@ -12,9 +12,6 @@ static struct planet *thePlanets = NULL;
 static double Mach = 0.0;
 static int Npl = 0;
 
-double get_cs2(double *);
-double get_nu(double *, double *);
-
 void setICparams( struct domain * theDomain )
 {
     gam = theDomain->theParList.Adiabatic_Index;
@@ -35,7 +32,6 @@ void initial(double *prim, double *x)
     double R = r/R0;
     double phi = x[1];
 
-    double cs2 = get_cs2(x);
 
     double rho, fact;
 
@@ -47,6 +43,7 @@ void initial(double *prim, double *x)
       om = pow(r,-1.5);
     }
 
+    double cs2 = get_cs2(x, rho);
     double v = -1.5*nu/(r+0.0001);
     double P = rho*cs2/gam;
 

diff --git a/Initial/cylinder.c b/Initial/cylinder.c
@@ -45,7 +45,7 @@ void initial(double *prim, double *x)
     double rho = rho_ref;
     double P;
     if(isothermal_flag)
-        P = get_cs2(x) * rho / gam;
+        P = get_cs2(x, rho) * rho / gam;
     else
         P = cs_ref*cs_ref * rho / gam;
 

diff --git a/Initial/gapGlobalDisk.c b/Initial/gapGlobalDisk.c
@@ -17,8 +17,6 @@ static int cs2choice = 0;
 static double a = 1.0;
 static double om_bary = 0.0;
 
-double get_cs2(double *);
-double get_nu(double *, double *);
 double phigrav( double , double , double , int); //int here is type
 double fgrav( double , double , double , int);
 
@@ -65,7 +63,7 @@ void initial(double *prim, double *x)
    //homtot +=  fgrav( thePlanets[1].M, d1p , thePlanets[1].eps, thePlanets[1].type)/fmax(d1p, thePlanets[1].eps*0.01);
    //homtot +=  fgrav( thePlanets[2].M, d2p , thePlanets[2].eps, thePlanets[2].type)/fmax(d2p, thePlanets[2].eps*0.01);
 
-   double cs2 = get_cs2(x);
+   double cs2 = get_cs2(x, 1.0);
 
    double alpha = get_nu(x, prim);
    alpha = alpha*sqrt(homtot)/cs2;

diff --git a/Initial/gapLocalDisk.c b/Initial/gapLocalDisk.c
@@ -17,8 +17,6 @@ static int cs2choice = 0;
 static double a = 1.0;
 static double om_bary = 0.0;
 
-double get_cs2(double *);
-double get_nu(double *, double *);
 double phigrav( double , double , double , int); //int here is type
 double fgrav( double , double , double , int);
 
@@ -64,10 +62,6 @@ void initial(double *prim, double *x)
 
    double homtot = 0.0;
    homtot +=  fgrav( thePlanets[0].M, R , thePlanets[0].eps, thePlanets[0].type)/R;
-   double cs2 = get_cs2(x);
-
-   double alpha = get_nu(x,prim);
-   alpha = alpha*sqrt(homtot)/cs2;
 
    double dphi = 0.0;
    if (cs2choice == 5){
@@ -81,8 +75,9 @@ void initial(double *prim, double *x)
    double dSigma = -xi*exp(-pow(r/rin,xi))*pow(r/rin, xi-1.0)/rin;
 
    double rho = Sigma;
-   double P = rho*cs2/gam;
+   double cs2 = get_cs2(x, rho);
 
+   double P = rho*cs2/gam;
    double vr	= -1.5*nu/r;
    double vp = sqrt(fabs((massq/(1.0+massq))/(r*r*r) + (cs2*dSigma/gam + Sigma*dphi/gam)/(r*Sigma)))+om_bary;