From 0de2558ec21b2f723acdefc1e3b11251a1a7849e Mon Sep 17 00:00:00 2001 From: ajdittmann Date: Wed, 12 Jul 2023 12:40:32 -0400 Subject: [PATCH 1/9] tentative lin/log gridding updates --- gridsetup.c | 23 ++++++++++++----------- paul.h | 1 - 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/gridsetup.c b/gridsetup.c index 8f84101..134db95 100644 --- a/gridsetup.c +++ b/gridsetup.c @@ -103,9 +103,20 @@ void setupGrid( struct domain * theDomain ){ theDomain->r_jph[j] = Rmin + x*(Rmax-Rmin); }else if( LogZoning == 1 ){ theDomain->r_jph[j] = Rmin*pow(Rmax/Rmin,x); - }else{ + }else if( LogZoning == 2 ){ theDomain->r_jph[j] = R0*(pow(Rmax/R0,x)-1) + Rmin + (R0-Rmin)*x; } + }else{ + double x1 = 1.0/(1.0 - R0*log(R0/Rmax)/(R0-Rmin)); + double b = exp((1-Rmin/R0)/x1); + double a = log(b)*Rmax/b; + if( x>x1 ){ + theDomain->r_jph[j] = a*pow(b, x)/log(b); + } + else{ + theDomain->r_jph[j] = a*pow(b, x1)*x + Rmin; + } + } } double dz = (Zmax-Zmin)/(double)Num_Z; double z0 = Zmin + (double)N0z*dz; @@ -113,16 +124,6 @@ void setupGrid( struct domain * theDomain ){ theDomain->z_kph[k] = z0 + ((double)k+1.)*dz; } - double dr0; - if(LogZoning == 0) - dr0 = (Rmax-Rmin) / (double) Num_R; - else if(LogZoning == 1) - dr0 = Rmin * (pow(Rmax/Rmin, 1.0/Num_R) - 1.0); - else - dr0 = R0*(pow(Rmax/R0,1.0/Num_R)-1) + (R0-Rmin)/Num_R; - theDomain->dr0 = dr0; - - theDomain->phi_max = theDomain->theParList.phimax; setGeometryParams( theDomain ); diff --git a/paul.h b/paul.h index 04a62f8..be832f1 100644 --- a/paul.h +++ b/paul.h @@ -227,7 +227,6 @@ struct domain{ double phi_max; int * fIndex_r; int * fIndex_z; - double dr0; time_t Wallt_init; int rank,size; From 4d0d7bbeec6d97ca9e8d6817ebbcfd3a7d5117c7 Mon Sep 17 00:00:00 2001 From: ajdittmann Date: Wed, 12 Jul 2023 15:41:33 -0400 Subject: [PATCH 2/9] typo fix --- gridsetup.c | 1 - 1 file changed, 1 deletion(-) diff --git a/gridsetup.c b/gridsetup.c index 134db95..7389b78 100644 --- a/gridsetup.c +++ b/gridsetup.c @@ -105,7 +105,6 @@ void setupGrid( struct domain * theDomain ){ theDomain->r_jph[j] = Rmin*pow(Rmax/Rmin,x); }else if( LogZoning == 2 ){ theDomain->r_jph[j] = R0*(pow(Rmax/R0,x)-1) + Rmin + (R0-Rmin)*x; - } }else{ double x1 = 1.0/(1.0 - R0*log(R0/Rmax)/(R0-Rmin)); double b = exp((1-Rmin/R0)/x1); From a9063b2554ad68064992ab6ec8298ee93dd1c1d1 Mon Sep 17 00:00:00 2001 From: Alexander Dittmann Date: Wed, 30 Aug 2023 23:49:16 -0400 Subject: [PATCH 3/9] tweaks for mixed softening setup --- Initial/keplerIso.c | 12 ++++--- Planets/cb_mixed.c | 57 +++++++++++++++++++++++++++++++++ Planets/cb_varyq_mixed.c | 68 ++++++++++++++++++++++++++++++++++++++++ Planets/cmA.c | 29 +++++++++-------- 4 files changed, 146 insertions(+), 20 deletions(-) create mode 100644 Planets/cb_mixed.c create mode 100644 Planets/cb_varyq_mixed.c diff --git a/Initial/keplerIso.c b/Initial/keplerIso.c index 40bc514..d9704c3 100644 --- a/Initial/keplerIso.c +++ b/Initial/keplerIso.c @@ -3,6 +3,7 @@ #include "../geometry.h" static double gam = 0.0; +static double nu = 1.0; static double beta = 0.0; static double Mach = 0.0; static double eps = 0.0; @@ -19,25 +20,26 @@ void setICparams( struct domain * theDomain ){ viscChoice = theDomain->theParList.visc_profile; viscPar = theDomain->theParList.visc_par; - + nu = theDomain->theParList.viscosity; } void initial( double * prim , double * x ){ + double rs = sqrt(x[0]*x[0] + eps*eps); double r = x[0]; - double nu = get_nu(x, prim); - double rho = 1.0; + //double nu = get_nu(x, prim); + double rho = 1.0; ///nu; double Pp = rho*get_cs2(x)/gam; double dlogrho = 0.0; if (viscChoice == 1) dlogrho = -0.5; if (viscChoice == 2) dlogrho = -1*viscPar; - double omega2 = (1.0/(r*r*r))*(1.0 - 1.0/(Mach*Mach) + dlogrho/(Mach*Mach)); + double omega2 = (1.0/(rs*rs*r))*(1.0 - 1.0/(Mach*Mach) + dlogrho/(Mach*Mach)); double omega = sqrt(omega2); - double Vrpz[3] = {-1.5*nu/r, r*omega, 0.0}; + double Vrpz[3] = {-1.5*nu/rs, r*omega, 0.0}; double V[3]; get_vec_from_rpz(x, Vrpz, V); get_vec_contravariant(x, V, V); diff --git a/Planets/cb_mixed.c b/Planets/cb_mixed.c new file mode 100644 index 0000000..db1bed9 --- /dev/null +++ b/Planets/cb_mixed.c @@ -0,0 +1,57 @@ + +#include "../paul.h" + +static double q_planet = 1.0; +static double Mach = 1.0; +static double eps = 0.0; +static double nu = 1.0; +static double tau = 2.0; + +void setPlanetParams( struct domain * theDomain ){ + + theDomain->Npl = 2; + q_planet = theDomain->theParList.Mass_Ratio; + Mach = theDomain->theParList.Disk_Mach; + eps = theDomain->theParList.grav_eps; + nu = theDomain->theParList.viscosity; +} + +int planet_motion_analytic( void ){ + return(1); +} + +void initializePlanets( struct planet * thePlanets ){ + + double a = 1.0; + + double q = q_planet; + double mu = q/(1.+q); + + double om = pow( a , -1.5 ); + + thePlanets[0].M = 1.-mu; + thePlanets[0].vr = 0.0; + thePlanets[0].omega = om; + thePlanets[0].vz = 0.0; + thePlanets[0].r = a*mu; + thePlanets[0].phi = M_PI; + thePlanets[0].z = 0.0; + thePlanets[0].eps = eps; + thePlanets[0].type = PLSPLINE; + + thePlanets[1].M = mu; + thePlanets[1].vr = 0.0; + thePlanets[1].omega = om; + thePlanets[1].vz = 0.0; + thePlanets[1].r = a*(1.-mu); + thePlanets[1].phi = 0.0; + thePlanets[1].z = 0.0; + thePlanets[1].eps = eps; + thePlanets[1].type = PLPOINTMASS; +} + +void movePlanets( struct planet * thePlanets , double t , double dt ){ + thePlanets[0].phi += thePlanets[0].omega*dt; + thePlanets[1].phi += thePlanets[1].omega*dt; +} + diff --git a/Planets/cb_varyq_mixed.c b/Planets/cb_varyq_mixed.c new file mode 100644 index 0000000..4f33416 --- /dev/null +++ b/Planets/cb_varyq_mixed.c @@ -0,0 +1,68 @@ + +#include "../paul.h" + +static double q_planet = 1.0; +static double Mach = 1.0; +static double eps = 0.0; +static double nu = 1.0; +static double tau = 2.0; + +void setPlanetParams( struct domain * theDomain ){ + + theDomain->Npl = 2; + q_planet = theDomain->theParList.Mass_Ratio; + Mach = theDomain->theParList.Disk_Mach; + eps = theDomain->theParList.grav_eps; + nu = theDomain->theParList.viscosity; +} + +int planet_motion_analytic( void ){ + return(1); +} + +void initializePlanets( struct planet * thePlanets ){ + + double a = 1.0; + + double q = q_planet; + double mu = q/(1.+q); + + double om = pow( a , -1.5 ); + + thePlanets[0].M = 1.-mu; + thePlanets[0].vr = 0.0; + thePlanets[0].omega = om; + thePlanets[0].vz = 0.0; + thePlanets[0].r = a*mu; + thePlanets[0].phi = M_PI; + thePlanets[0].z = 0.0; + thePlanets[0].eps = eps; + thePlanets[0].type = PLSPLINE; + + thePlanets[1].M = mu; + thePlanets[1].vr = 0.0; + thePlanets[1].omega = om; + thePlanets[1].vz = 0.0; + thePlanets[1].r = a*(1.-mu); + thePlanets[1].phi = 0.0; + thePlanets[1].z = 0.0; + thePlanets[1].eps = eps; + thePlanets[1].type = PLPOINTMASS; +} + +void movePlanets( struct planet * thePlanets , double t , double dt ){ + thePlanets[0].phi += thePlanets[0].omega*dt; + thePlanets[1].phi += thePlanets[1].omega*dt; + + double T = tau*2.0*M_PI/nu; + double q = q_planet*(1.+pow(t/T, 2.)); + double mu = q/(1.+q); + + thePlanets[0].M = 1.-mu; + thePlanets[1].M = mu; + + thePlanets[0].r = mu; + thePlanets[1].r = 1.-mu; + +} + diff --git a/Planets/cmA.c b/Planets/cmA.c index 218585f..a1b5878 100644 --- a/Planets/cmA.c +++ b/Planets/cmA.c @@ -18,30 +18,30 @@ int planet_motion_analytic( void ){ void initializePlanets( struct planet * thePlanets ){ double a = 1.0; - + double q = q_planet; double mu = q/(1.+q); double om = pow( a , -1.5 ); - thePlanets[0].M = 1.-mu; - thePlanets[0].vr = 0.0; - thePlanets[0].omega = om; - thePlanets[0].vz = 0.0; - thePlanets[0].r = a*mu; - thePlanets[0].phi = M_PI; - thePlanets[0].z = 0.0; + thePlanets[0].M = 1.-mu; + thePlanets[0].vr = 0.0; + thePlanets[0].omega = om; + thePlanets[0].vz = 0.0; + thePlanets[0].r = a*mu; + thePlanets[0].phi = M_PI; + thePlanets[0].z = 0.0; thePlanets[0].eps = eps; thePlanets[0].type = PLPOINTMASS; thePlanets[1].M = mu; - thePlanets[1].vr = 0.0; - thePlanets[1].omega = om; - thePlanets[1].vz = 0.0; - thePlanets[1].r = a*(1.-mu); - thePlanets[1].phi = 0.0; - thePlanets[1].z = 0.0; + thePlanets[1].vr = 0.0; + thePlanets[1].omega = om; + thePlanets[1].vz = 0.0; + thePlanets[1].r = a*(1.-mu); + thePlanets[1].phi = 0.0; + thePlanets[1].z = 0.0; thePlanets[1].eps = eps; thePlanets[1].type = PLPOINTMASS; } @@ -51,4 +51,3 @@ void movePlanets( struct planet * thePlanets , double t , double dt ){ thePlanets[0].phi += thePlanets[0].omega*dt; thePlanets[1].phi += thePlanets[1].omega*dt; } - From d02a5a61e17a1b910633d57b12d331e7f13cbbd8 Mon Sep 17 00:00:00 2001 From: ajdittmann Date: Wed, 2 Oct 2024 17:37:46 -0400 Subject: [PATCH 4/9] adding energy diagnostic to reports --- Reports/rep_cb.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Reports/rep_cb.c b/Reports/rep_cb.c index 18ba530..8b0b186 100644 --- a/Reports/rep_cb.c +++ b/Reports/rep_cb.c @@ -2,7 +2,7 @@ #include "../geometry.h" #include "../planet.h" -#define N_AUX_PER_PLANET 14 +#define N_AUX_PER_PLANET 15 static double gamma_law = 0.0; static int Npl = 0; @@ -58,7 +58,7 @@ void get_distributed_aux_reports(double *Q, struct domain *theDomain) PL_SNK_M, PL_GRV_JZ, PL_SNK_JZ, PL_GRV_PX, PL_GRV_PY, PL_SNK_PX, PL_SNK_PY, PL_GRV_K, PL_SNK_K, PL_SNK_MX, PL_SNK_MY, PL_SNK_SZ, - PL_GRV_U, PL_SNK_U}; + PL_GRV_U, PL_SNK_U, PL_SNK_UGAS}; int q; for(q=0; q Date: Thu, 7 Nov 2024 16:55:26 -0500 Subject: [PATCH 5/9] eccentric setup --- Initial/keplerQP.c | 44 ++++++++++++ Planets/planet_ecc.c | 4 ++ Planets/planet_excised.c | 4 ++ Templates/kep_ecc.in | 30 +++++++++ Templates/kep_ecc.par | 140 +++++++++++++++++++++++++++++++++++++++ omega.c | 10 ++- 6 files changed, 230 insertions(+), 2 deletions(-) create mode 100644 Initial/keplerQP.c create mode 100644 Templates/kep_ecc.in create mode 100644 Templates/kep_ecc.par diff --git a/Initial/keplerQP.c b/Initial/keplerQP.c new file mode 100644 index 0000000..12d6408 --- /dev/null +++ b/Initial/keplerQP.c @@ -0,0 +1,44 @@ +#include "../paul.h" +#include "../omega.h" +#include "../geometry.h" + +static double gam = 0.0; +static double nu = 1.0; +static double Mach = 0.0; +static double eps = 0.0; +static double p = 0.0; +static double q = 0.0; + +void setICparams( struct domain * theDomain ){ + gam = theDomain->theParList.Adiabatic_Index; + Mach = theDomain->theParList.Disk_Mach; + eps = theDomain->theParList.grav_eps; + + q = theDomain->theParList.Cs2_Par; //q + p = theDomain->theParList.visc_par; //p + nu = theDomain->theParList.viscosity; +} + +void initial( double * prim , double * x ){ + + double r = x[0]; + + nu *= pow(fmax(x[0],1e-10), p); + + double rho = 1.0*pow(x[0], -p); + double Pp = rho*get_cs2(x)/gam; + + double omega2 = (1.0/(r*r*r))*(1.0 - 1.0/(Mach*Mach)*(p+q)*pow(r, 1-q) ); + double omega = sqrt(omega2); + + double Vrpz[3] = {-1.5*nu/r, 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]; +} diff --git a/Planets/planet_ecc.c b/Planets/planet_ecc.c index 8a9f781..6058020 100644 --- a/Planets/planet_ecc.c +++ b/Planets/planet_ecc.c @@ -38,16 +38,20 @@ void initializePlanets( struct planet * thePlanets ){ thePlanets[0].M = 1.0 - mu; thePlanets[0].vr = 0.0; thePlanets[0].omega = 0.0; + thePlanets[0].vz = 0.0; thePlanets[0].r = R*mu; thePlanets[0].phi = M_PI; + thePlanets[0].z = 0.0; thePlanets[0].eps = 0.0; thePlanets[0].type = PLPOINTMASS; thePlanets[1].M = mu; thePlanets[1].vr = 0.0; thePlanets[1].omega = om; + thePlanets[1].vz = 0.0; thePlanets[1].r = R*(1.0-mu); thePlanets[1].phi = 0.0; + thePlanets[1].z = 0.0; thePlanets[1].eps = eps; thePlanets[1].type = PLPOINTMASS; diff --git a/Planets/planet_excised.c b/Planets/planet_excised.c index ac3e928..ac5776f 100644 --- a/Planets/planet_excised.c +++ b/Planets/planet_excised.c @@ -27,8 +27,10 @@ void initializePlanets( struct planet * thePlanets ){ thePlanets[0].M = 1.-mu; thePlanets[0].vr = 0.0; thePlanets[0].omega = om; + thePlanets[0].vz = 0.0; thePlanets[0].r = a*mu; thePlanets[0].phi = M_PI; + thePlanets[0].z = 0.0; thePlanets[0].eps = 0.0; thePlanets[0].type = PLPOINTMASS; @@ -36,8 +38,10 @@ void initializePlanets( struct planet * thePlanets ){ thePlanets[1].M = mu; thePlanets[1].vr = 0.0; thePlanets[1].omega = om; + thePlanets[1].vz = 0.0; thePlanets[1].r = a*(1.-mu); thePlanets[1].phi = 0.0; + thePlanets[1].z = 0.0; thePlanets[1].eps = eps; thePlanets[1].type = PLPOINTMASS; } diff --git a/Templates/kep_ecc.in b/Templates/kep_ecc.in new file mode 100644 index 0000000..53d8548 --- /dev/null +++ b/Templates/kep_ecc.in @@ -0,0 +1,30 @@ +### Makefile_opt.in for Disco ### + +# This file contains user-specified build options. + +# Build options +INITIAL = keplerQP +HYDRO = euler +GEOMETRY = cylindrical +BOUNDARY = fixed3d +OUTPUT = h5out +RESTART = h5in +PLANETS = planet_ecc +HLLD = hlld +ANALYSIS = diag_cb +REPORT = rep_cb + +METRIC = flat_rpz +FRAME = normal + +NUM_C = 5 # Number of conserved variables (Hydro: 5, MHD: 8) +NUM_N = 0 # Number of passive variables (>=0) +CT_MODE = 0 # Constained Transport mode. None: 0, 2D: 1, 3D: 2. + +ENABLE_CART_INTERP = 0 # 0: Default. Cartesian Interpolation routines will + # NOT be compiled into Disco. The Cartesian_Interp + # parameters will do nothing, but the code may run + # faster. + # 1: Cartesian Interpolation routines will be compiled + # into Disco. The Cartesian_Interp parameters will + # be active. diff --git a/Templates/kep_ecc.par b/Templates/kep_ecc.par new file mode 100644 index 0000000..00b170e --- /dev/null +++ b/Templates/kep_ecc.par @@ -0,0 +1,140 @@ + +//Disco Parameter File// +Restart 0 //Restart from input.h5 + +//Time Range// +T_Start 0.0 //1.0 +T_End 100.0 //0.12 //100000.0 +T_Times_2pi 1 + +//Output Frequency// +Num_Reports 1e4 +Num_Snapshots 0 +Num_Checkpoints 100 +Use_Logtime 0 + +//Grid Dimensions// +Num_R 1152 +Num_Z 1 +aspect 1.0 //Cell Aspect Ratio, for calculating N_phi +Max_Aspect_Short 1.5 +Max_Aspect_Long 1.5 //Aspect Ratio AMR Refinement Criteria + +//Domain Parameters// +R_Min 0.2 +R_Max 5.0 +Z_Min -0.5 //-0.0625 +Z_Max 0.5 //0.0625 +R_Periodic 0 +Z_Periodic 0 +Phi_Max 1.0 +P_Times_2pi 1 +Z_Times_pi 0 +Log_Zoning 1 //0=Uniform, 1=Log, 2=Hybrid +Log_Radius 0.25 //Only for option 2 of Log_Zoning +NoBC_Rmin 0 //0 = Regular boundary. Ghost zones added and +NoBC_Rmax 0 // boundary condition applied. +NoBC_Zmin 0 //1 = Coordinate boundary (like r=0 origin). +NoBC_Zmax 0 // No ghost zones or explicit BC applied. + +//Numerical Parameters +CFL 0.5 +PLM 1.5 +Max_DT 0.0 //If > 0.0, the maximum allowed time step. +Timestep 0 // 0=RK2, 1=FE +Riemann_Solver 1 //0=HLL, 1=HLLC, 2=HLLD +Mesh_Motion 3 //0=Fixed, 1=By Cell, 2=Riemann Solver, 3=Avg +Exact_Mesh_Omega 0 //0=Fixed, 1=Rigid, 2=Kep, 3=Rig2Kep +Exact_Mesh_Omega_Par 0.0 //unused so far +Absorbing_BC 0 +Initial_Regrid 0 +Density_Floor 1e-5 +Pressure_Floor 1e-5 +Constrained_Transport 0 +Cartesian_Interp 0 // Fixes the pole at the cost of disks 0 off, 1 on +Cartesian_Interp_R0 0.1 // Radius within which Cartesian_Interp acts. + +//Hydro Parameters +Adiabatic_Index 1.00001 //1.66666666667 +Isothermal 1 +Cs2_Profile 6 //1=Flat,4=ConstMach,5=Binary,6=PowerLaw +Cs2_Par 1.0 //PowerLaw for profile 6 +Energy_Omega 0 //0=Fixed, 1=Rigid, 2=Kep, 3=Rig2Kep, 4=Gaussian +Energy_Omega_Par 0.0 //unused so far +Use_Viscosity 1 +Viscosity 0.0 +Viscosity_Profile 2 //0=Flat, 1=alpha, 2=powerlaw, 3=power law (binary) +Viscosity_Par 0.0 //PowerLaw for profiles 2,3 + + +//Planet Parameters +Mass_Ratio 1e-5 +Eccentricity 0.01 +Drift_Rate 0.0 //-2.5e-4 +Drift_Exp 0.66666667 +Grav2D 1 +Softening 0.015 +//Rotating Frame +RotFrame 0 //0=off, 1=on +RotOmega 0.0 //Angular speed of frame +RotD 0.0 //Distance from rotational axis to grid axis. + +//Disk Parameters +Mach_Number 20.0 +Include_Atmos 1 +//Initial Condition Parameters +Init_Par0 0 // int +Init_Par1 0.0 // doub +Init_Par2 0.0 // doub +Init_Par3 0.0 // doub +Init_Par4 0.0 // doub +Init_Par5 0.0 // doub +Init_Par6 0.0 // doub +Init_Par7 0.0 // doub +Init_Par8 0.0 // doub + +//Noise Parameters +Noise_Type 0 // 0: No Noise, 1: URR and UPP +Noise_Abs 0.0 // Amplitude of Absolute Noise +Noise_Rel 0.0 // Amplitude of Relative Noise +//Sink & Nozzle Parameters +Sink_Type 0 //0: No Sink +Sink_Number 0 //0: all planets have sinks. else: first N planets have sinks +Sink_Par1 0.0 //(double) +Sink_Par2 0.0 //(double) +Sink_Par3 0.0 //(double) +Sink_Par4 0.0 //(double) +Sink_Par5 0.0 //(double) +Nozzle_Type 0 //0: No Nozzle, 1: Circular Nozzle +Nozzle_Par1 0.0 //(double) +Nozzle_Par2 0.0 //(double) +Nozzle_Par3 0.0 //(double) +Nozzle_Par4 0.0 //(double) + +//Damping parameters +DampInner_Type 0 //0: No damping, 1: "time" is actual time, 2: "time" is fraction of orbital time +DampInner_Time 1.e-2 //(double) +DampInner_Len 0.01 //(double) +DampOuter_Type 0 //0: No damping, 1: "time" is actual time, 2: "time" is fraction of orbital time +DampOuter_Time 1.e-2 //(double) +DampOuter_Len 0.167 //(double) +DampLower_Type 0 //0: No damping, 1: "time" is actual time, 2: "time" is fraction of orbital time +DampLower_Time 1.e-2 //(double) +DampLower_Len 1.0 //(double) +DampUpper_Type 0 //0: No damping, 1: "time" is actual time, 2: "time" is fraction of orbital time +DampUpper_Time 1.e-2 //(double) +DampUpper_Len 1.0 //(double) + +//Cooling parameters +Cool_Type 0 //0: No cooling, 1: beta cooling, 2: beta relaxation +Cool_Par1 1.e-2 //(double) +Cool_Par2 0.0 //(double) +Cool_Par3 0.0 //(double) +Cool_Par4 0.0 //(double) +//Metric Parameters +Metric_Par0 0 //(int) +Metric_Par1 0.0 //Om_rot (double) +Metric_Par2 1.0 //M (double) +Metric_Par3 0.0 //a (double) +Metric_Par4 0.0 + diff --git a/omega.c b/omega.c index c26f8ba..b11f442 100644 --- a/omega.c +++ b/omega.c @@ -193,12 +193,18 @@ double get_cs2( const double *x ){ int pi; double phip = 0.0; - + for(pi = 0; pi < Npl; pi++) phip -= planetaryPotential(thePlanets+pi, xyz); cs2 = fabs(phip)/(Mach*Mach); } + else if(cs2Choice == 6) + { + double r = x[0]; + double scaling = pow(r, -cs2Par); // + cs2 = scaling/(Mach*Mach); + } else cs2 = 1.0; @@ -212,7 +218,7 @@ double get_nu(const double x[], const double prim[]){ double c2 = gamma_law*prim[PPP]/prim[RHO]; nu *= c2/get_height_om(x); } - //generic power law w.r.t. r=0 + //generic power law w.r.t. r=1 if (viscChoice == 2){ nu *= pow(fmax(x[0],1e-10), viscPar); } From 9ba1de62807722ec9698166e20f20b71263210c4 Mon Sep 17 00:00:00 2001 From: ajdittmann Date: Thu, 5 Dec 2024 21:09:27 -0500 Subject: [PATCH 6/9] astrocentric planet setup for comparison with simple linear theory --- Planets/fakeAstrocentric.c | 116 +++++++++++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) create mode 100644 Planets/fakeAstrocentric.c diff --git a/Planets/fakeAstrocentric.c b/Planets/fakeAstrocentric.c new file mode 100644 index 0000000..3ba19d4 --- /dev/null +++ b/Planets/fakeAstrocentric.c @@ -0,0 +1,116 @@ + +#include "../paul.h" + +static double q_planet = 1.0; +static double e_planet = 0.0; +static double eps = 0.05; + +void setPlanetParams( struct domain * theDomain ){ + + theDomain->Npl = 2; + q_planet = theDomain->theParList.Mass_Ratio; + e_planet = theDomain->theParList.Eccentricity; + eps = theDomain->theParList.grav_eps; +} + +double root0( double E , double e ){ + return( E - e*sin(E) ); +} + +double root1( double E , double e ){ + return( 1. - e*cos(E) ); +} + +int planet_motion_analytic( void ){ + return(1); +} + +void initializePlanets( struct planet * thePlanets ){ + + double a = 1.0; + double e = e_planet; + double R = a*(1.-e); + double om = pow( a , -1.5 )*sqrt(1.-e*e)/(1.-e)/(1.-e); + + double q = q_planet; + double mu = q/(1.+q); + + thePlanets[0].M = 1.0 - mu; + thePlanets[0].vr = 0.0; + thePlanets[0].omega = 0.0; + thePlanets[0].vz = 0.0; + thePlanets[0].r = 0.0; + thePlanets[0].phi = M_PI; + thePlanets[0].z = 0.0; + thePlanets[0].eps = 0.0; + thePlanets[0].type = PLPOINTMASS; + + thePlanets[1].M = mu; + thePlanets[1].vr = 0.0; + thePlanets[1].omega = om; + thePlanets[1].vz = 0.0; + thePlanets[1].r = R; + thePlanets[1].phi = 0.0; + thePlanets[1].z = 0.0; + thePlanets[1].eps = eps; + thePlanets[1].type = PLPOINTMASS; + + +} + +void movePlanets( struct planet * thePlanets , double t , double dt ){ + + double TOL = 1e-8; + + double r0 = thePlanets[0].r + thePlanets[1].r; + double phi0 = thePlanets[1].phi; + + double vr = thePlanets[0].vr + thePlanets[1].vr; + double omega = thePlanets[1].omega; + + double l = r0*r0*omega; + double en = 0.5*vr*vr - 1./r0 + 0.5*l*l/r0/r0; + + double a = 1./2./fabs(en); + double b = l/sqrt(2.*fabs(en)); + double f = sqrt(fabs(a*a-b*b)); + double e = f/a; + + double x0 = r0*cos(phi0); + double y0 = r0*sin(phi0); + + double E0 = atan2( y0/b , (x0+f)/a ); + double M0 = E0 - e*sin(E0); + double M = M0 + l*dt/a/b; + +//Newton-Rapheson to solve M = E - e*sin(E) + double E = M; //Guess value for E is M. + double ff = root0( E , e ) - M; + while( fabs(ff) > TOL ){ + double dfdE = root1( E , e ); + double dE = -ff/dfdE; + E += dE; + ff = root0( E , e ) - M; + } + double x = a*cos(E)-f; + double y = b*sin(E); + double R = sqrt(x*x+y*y); + double phi = atan2(y,x); + + vr = sqrt( fabs( 2.*en + 2./R - l*l/R/R ) ); + if( y<0.0 ) vr *= -1.; + + double mu = q_planet/(1.+q_planet); + + thePlanets[1].r = R; + thePlanets[1].phi = phi; + thePlanets[1].omega = l/R/R; + thePlanets[1].vr = vr; + + thePlanets[0].r = 0.0; + thePlanets[0].phi = 0.0; + thePlanets[0].omega = 0.0; + thePlanets[0].vr = 0.0; + +} + From d2e172a645a01c29e617d8ac91122a8d4b735084 Mon Sep 17 00:00:00 2001 From: ajdittmann Date: Thu, 5 Dec 2024 21:09:42 -0500 Subject: [PATCH 7/9] minor text alignment --- Planets/planet_ecc.c | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/Planets/planet_ecc.c b/Planets/planet_ecc.c index 6058020..130bd3f 100644 --- a/Planets/planet_ecc.c +++ b/Planets/planet_ecc.c @@ -35,23 +35,23 @@ void initializePlanets( struct planet * thePlanets ){ double q = q_planet; double mu = q/(1.+q); - thePlanets[0].M = 1.0 - mu; - thePlanets[0].vr = 0.0; - thePlanets[0].omega = 0.0; - thePlanets[0].vz = 0.0; - thePlanets[0].r = R*mu; - thePlanets[0].phi = M_PI; - thePlanets[0].z = 0.0; + thePlanets[0].M = 1.0 - mu; + thePlanets[0].vr = 0.0; + thePlanets[0].omega = om; + thePlanets[0].vz = 0.0; + thePlanets[0].r = R*mu; + thePlanets[0].phi = M_PI; + thePlanets[0].z = 0.0; thePlanets[0].eps = 0.0; thePlanets[0].type = PLPOINTMASS; - thePlanets[1].M = mu; - thePlanets[1].vr = 0.0; - thePlanets[1].omega = om; - thePlanets[1].vz = 0.0; - thePlanets[1].r = R*(1.0-mu); - thePlanets[1].phi = 0.0; - thePlanets[1].z = 0.0; + thePlanets[1].M = mu; + thePlanets[1].vr = 0.0; + thePlanets[1].omega = om; + thePlanets[1].vz = 0.0; + thePlanets[1].r = R*(1.0-mu); + thePlanets[1].phi = 0.0; + thePlanets[1].z = 0.0; thePlanets[1].eps = eps; thePlanets[1].type = PLPOINTMASS; From c24ddbac7456fb51b55c34e35209c9efcaba63a9 Mon Sep 17 00:00:00 2001 From: ajdittmann Date: Thu, 5 Dec 2024 21:09:59 -0500 Subject: [PATCH 8/9] new damping options --- paul.h | 1 + readpar.c | 1 + sink.c | 57 +++++++++++++++++++++++++++++++------------------------ 3 files changed, 34 insertions(+), 25 deletions(-) diff --git a/paul.h b/paul.h index be832f1..9ff6aa9 100644 --- a/paul.h +++ b/paul.h @@ -156,6 +156,7 @@ struct param_list{ double coolPar3; double coolPar4; + int dampType; int dampInnerType; int dampOuterType; int dampUpperType; diff --git a/readpar.c b/readpar.c index 5653ac7..628797b 100644 --- a/readpar.c +++ b/readpar.c @@ -167,6 +167,7 @@ int read_par_file( struct domain * theDomain ){ err += readvar(pfile, "Cool_Par2", VAR_DOUB, &(theList->coolPar2)); err += readvar(pfile, "Cool_Par3", VAR_DOUB, &(theList->coolPar3)); err += readvar(pfile, "Cool_Par4", VAR_DOUB, &(theList->coolPar4)); + err += readvar(pfile, "Damp_Type", VAR_INT , &(theList->dampType)); err += readvar(pfile, "DampInner_Type", VAR_INT , &(theList->dampInnerType)); err += readvar(pfile, "DampInner_Time", VAR_DOUB, &(theList->dampTimeInner)); err += readvar(pfile, "DampInner_Len", VAR_DOUB, &(theList->dampLenInner)); diff --git a/sink.c b/sink.c index 8068a43..d4c480c 100644 --- a/sink.c +++ b/sink.c @@ -37,6 +37,7 @@ static double rmin; static double zmax; static double zmin; +static int dampType = 0; static int DAMP_OUTER = 0; static int DAMP_INNER = 0; static int DAMP_LOWER = 0; @@ -87,6 +88,7 @@ void setSinkParams(struct domain *theDomain) thePlanets = theDomain->thePlanets; Npl = theDomain->Npl; + dampType = theDomain->theParList.dampType; DAMP_INNER = theDomain->theParList.dampInnerType; DAMP_OUTER = theDomain->theParList.dampOuterType; DAMP_UPPER = theDomain->theParList.dampUpperType; @@ -396,26 +398,26 @@ void damping(const double *prim, double *cons, const double *xp, double omtot = 1.0; if (DAMP_INNER==2 || DAMP_OUTER==2 || DAMP_LOWER==2 || DAMP_UPPER==2){ - double gx = xyz[0]; - double gy = xyz[1]; - - int pi; - omtot = 0.0; - double px, py, dx, dy, mag; - double Fxyz[3]; - for (pi=0; pi Date: Fri, 17 Jan 2025 13:30:36 -0500 Subject: [PATCH 9/9] updated get_cs2 calls, *untested* --- Hydro/euler.c | 2 +- Hydro/euler_altvisc.c | 2 +- Hydro/euler_cart.c | 2 +- Hydro/euler_cartInterp.c | 2 +- Hydro/euler_cart_cons.c | 2 +- Hydro/euler_pro.c | 2 +- Hydro/euler_sph.c | 2 +- Hydro/greuler.c | 2 +- Hydro/mhd.c | 2 +- Hydro/mhd_cart.c | 2 +- Hydro/mhd_sph.c | 2 +- Initial/binary.c | 6 ++-- Initial/bl.c | 5 ++- Initial/blBRS.c | 2 +- Initial/blBaro.c | 71 +++++++++++++++++++++++++++++++++++++++ Initial/blob.c | 2 +- Initial/cb.c | 5 +-- Initial/cbJdot.c | 2 -- Initial/cbMMLth.c | 3 +- Initial/cb_cc.c | 2 +- Initial/cb_duff.c | 5 +-- Initial/cylinder.c | 2 +- Initial/gapGlobalDisk.c | 4 +-- Initial/gapLocalDisk.c | 9 ++--- Initial/ggt.c | 6 ++-- Initial/kep_half.c | 6 ++-- Initial/kepler.c | 2 +- Initial/keplerCool.c | 2 +- Initial/keplerIso.c | 2 +- Initial/keplerQP.c | 2 +- Initial/kepler_cav.c | 2 +- Initial/kepler_gap.c | 5 ++- Initial/kepler_gap_2.c | 7 ++-- Initial/kepler_mhd.c | 7 ++-- Initial/kepler_softened.c | 2 +- Initial/rigid_perturb.c | 2 +- Initial/shear_cart.c | 5 ++- Initial/singleSinkTest.c | 2 +- Initial/warp.c | 2 +- omega.c | 6 +++- omega.h | 2 +- sink.c | 2 +- 42 files changed, 127 insertions(+), 77 deletions(-) create mode 100644 Initial/blBaro.c diff --git a/Hydro/euler.c b/Hydro/euler.c index 61943b8..03e0dab 100644 --- 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 index e2aa31a..53268e1 100644 --- 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 index 53693cd..f8636d0 100644 --- 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 index 552ddc0..bef37e9 100644 --- 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 index 98f7006..40a7243 100644 --- 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 index bdbf12a..52839ee 100644 --- 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 index efdef2a..8669a3e 100644 --- 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 index 305d10f..eb6369c 100644 --- 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 index 39f0380..deb31d1 100644 --- 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 index 7a85d98..5096e0a 100644 --- 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 index 399c12d..8b2291f 100644 --- 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 index 3dd2a79..a2b1c41 100644 --- 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 index 940c799..2ed1378 100644 --- 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 index 8a3afde..a3f538a 100644 --- 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 new file mode 100644 index 0000000..1b92ff3 --- /dev/null +++ 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 index ab97b10..770c69e 100644 --- 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 index 2ed031c..079d135 100644 --- 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 index cba39ed..038eb99 100644 --- 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 index b9d4ae3..c16dedb 100644 --- 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 index 70d7f86..3a0aeda 100644 --- 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 index 337a533..9e6d0af 100644 --- 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 index c252463..6bb3e76 100644 --- 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 index ac2673e..8b0fe66 100644 --- 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 index e0bcb7e..419841a 100644 --- 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; diff --git a/Initial/ggt.c b/Initial/ggt.c index 60ec821..9223eb6 100644 --- a/Initial/ggt.c +++ b/Initial/ggt.c @@ -1,4 +1,4 @@ - +#include "../omega.h" #include "../paul.h" static double gam = 0.0; @@ -13,20 +13,18 @@ void setICparams( struct domain * theDomain ){ rho_floor = theDomain->theParList.Density_Floor; } -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 rho = 1.0*exp(-sint*sint*Mach*Mach); double rho = pow(r,-1.5); // rho *= exp(-pow(r/10.,8.)); if( rho < rho_floor ) rho = rho_floor; + double cs2 = get_cs2( x, rho ); double Pp = rho*cs2/gam; //double n = 1.5; //double omega = ( pow( r , n-1.5 ) + 1. )/( pow( r , n ) + 1. ); diff --git a/Initial/kep_half.c b/Initial/kep_half.c index a057a4b..2265d0e 100644 --- a/Initial/kep_half.c +++ b/Initial/kep_half.c @@ -1,4 +1,4 @@ - +#include "../omega.h" #include "../paul.h" static double gam = 0.0; @@ -11,8 +11,6 @@ void setICparams( struct domain * theDomain ){ Mach = theDomain->theParList.Disk_Mach; } -double get_cs2( double ); - void initial( double * prim , double * x ){ double r = x[0]; @@ -24,7 +22,7 @@ void initial( double * prim , double * x ){ double k = 0.5; double rho = 1.0/pow(r,k) ;//+ 1.*exp(-20.*sc_r2); - double cs2 = get_cs2( r ); + double cs2 = get_cs2( x, rho ); double Pp = rho*cs2/gam; double omega02 = 0.5/pow(r,3.) ;//+.5-.5/pow(1.-r,3.); double omegaP2 = omega02*( (1.+k)/Mach/Mach ); diff --git a/Initial/kepler.c b/Initial/kepler.c index c69aa7b..1386dc1 100644 --- a/Initial/kepler.c +++ b/Initial/kepler.c @@ -33,7 +33,7 @@ void initial( double * prim , double * x ){ double rho = 1.0; double Pp = rho/(gam * Mach * Mach); if(isothermal_flag) - Pp = rho * get_cs2(x); + Pp = rho * get_cs2(x, rho); double visc = 0.0; //get_nu(x, prim); //if (nu > 0.0) rho = rho/nu; diff --git a/Initial/keplerCool.c b/Initial/keplerCool.c index 54e4cb8..06074de 100644 --- a/Initial/keplerCool.c +++ b/Initial/keplerCool.c @@ -29,7 +29,7 @@ void initial( double * prim , double * x ){ double nu = get_nu(x, prim); double rho = 1.0; - double Pp = rho*(get_cs2(x)/gam + (gam-1)*9*nu*beta*omega/4); + double Pp = rho*(get_cs2(x, rho)/gam + (gam-1)*9*nu*beta*omega/4); double Vrpz[3] = {-1.5*nu/r, r*omega, 0.0}; diff --git a/Initial/keplerIso.c b/Initial/keplerIso.c index d9704c3..3eca041 100644 --- a/Initial/keplerIso.c +++ b/Initial/keplerIso.c @@ -30,7 +30,7 @@ void initial( double * prim , double * x ){ //double nu = get_nu(x, prim); double rho = 1.0; ///nu; - double Pp = rho*get_cs2(x)/gam; + double Pp = rho*get_cs2(x,rho)/gam; double dlogrho = 0.0; if (viscChoice == 1) dlogrho = -0.5; diff --git a/Initial/keplerQP.c b/Initial/keplerQP.c index 12d6408..f13b191 100644 --- a/Initial/keplerQP.c +++ b/Initial/keplerQP.c @@ -26,7 +26,7 @@ void initial( double * prim , double * x ){ nu *= pow(fmax(x[0],1e-10), p); double rho = 1.0*pow(x[0], -p); - double Pp = rho*get_cs2(x)/gam; + double Pp = rho*get_cs2(x, rho)/gam; double omega2 = (1.0/(r*r*r))*(1.0 - 1.0/(Mach*Mach)*(p+q)*pow(r, 1-q) ); double omega = sqrt(omega2); diff --git a/Initial/kepler_cav.c b/Initial/kepler_cav.c index e11f8fd..454a108 100644 --- a/Initial/kepler_cav.c +++ b/Initial/kepler_cav.c @@ -43,7 +43,7 @@ void initial( double * prim , double * x ){ double rho = efact*(1-epsfl) + epsfl; double drho = (1-epsfl)*efact*xi*pow((R/redge),-xi)/R; - double P = rho*get_cs2(x)/gam; + double P = rho*get_cs2(x, rho)/gam; double addom = rho*dphitot + phitot*drho; addom *= 1.0/(Mach*Mach*r*rho); diff --git a/Initial/kepler_gap.c b/Initial/kepler_gap.c index e963ddb..1c49bae 100644 --- a/Initial/kepler_gap.c +++ b/Initial/kepler_gap.c @@ -1,4 +1,4 @@ - +#include "../omega.h" #include "../paul.h" static double gam = 0.0; @@ -22,7 +22,6 @@ void setICparams( struct domain * theDomain ){ alpha_csd = theDomain->theParList.initPar2; } -double get_cs2(double *); void initial( double * prim , double * x ){ @@ -70,7 +69,7 @@ void initial( double * prim , double * x ){ //double Pp = rho/(gam*Mach*Mach); //double Pp = rho/(gam*mach_csd*mach_csd); - double Pp = rho * get_cs2(x) / gam; + double Pp = rho * get_cs2(x, rho) / gam; double X = 0.0; if( r*cos(x[1]) > 0.0 ) X = 1.0; diff --git a/Initial/kepler_gap_2.c b/Initial/kepler_gap_2.c index e6fb1bd..43ea1eb 100644 --- a/Initial/kepler_gap_2.c +++ b/Initial/kepler_gap_2.c @@ -1,4 +1,4 @@ - +#include "../omega.h" #include "../paul.h" static double gam = 0.0; @@ -27,7 +27,6 @@ void setICparams( struct domain * theDomain ){ struct Gap_Vals gap_density(double r, double d, double M, double alpha, double q, double rho_0); //double gap_density(double r, double d, double M, double alpha, double q, double rho_0); //double d_gap_density(double r, double d, double M, double alpha, double q, double rho_0); -double get_cs2(double *); /* double gap_density(double r, double a, double Mach, double alpha, double q, double rho_0){ @@ -126,9 +125,9 @@ void initial( double * prim , double * x ){ double d_rho = vals.d_r_rho; //double rho = gap_density(R, a, mach_csd, alpha_csd, q_planet, rho_0); //double d_rho = d_gap_density(R, a, mach_csd, alpha_csd, q_planet, rho_0); - double Pp = rho * get_cs2(x) / gam; + double Pp = rho * get_cs2(x, rho) / gam; - double omega2 = (1./(R*R*R)) - (1/(R*gam*rho))*get_cs2(x)*d_rho; + double omega2 = (1./(R*R*R)) - (1/(R*gam*rho))*get_cs2(x, rho)*d_rho; double X = 0.0; if( r*cos(x[1]) > 0.0 ) X = 1.0; diff --git a/Initial/kepler_mhd.c b/Initial/kepler_mhd.c index 55e08e5..109193c 100644 --- a/Initial/kepler_mhd.c +++ b/Initial/kepler_mhd.c @@ -1,4 +1,4 @@ - +#include "../omega.h" #include "../paul.h" /* @@ -33,7 +33,6 @@ void get_rpz(double *x, double *rpz); void get_coords_from_rpz(double *rpz, double *x); void get_vec_from_rpz(double *x, double *Vrpz, double *V); void get_vec_contravariant(double *x, double *V, double *Vc); -double get_cs2(double *x); void initial( double * prim , double * X ) { @@ -125,8 +124,8 @@ void initial( double * prim , double * X ) double x0[3]; get_coords_from_rpz(rpz0, x0); - cs20 = get_cs2(x0); - cs2 = get_cs2(X); + cs20 = get_cs2(x0, 1.0); + cs2 = get_cs2(X, 1.0); rho = rho_atm0 * exp(h/cs2); rho0 = rho_atm0 * exp(h0/cs20); diff --git a/Initial/kepler_softened.c b/Initial/kepler_softened.c index 0cb1720..afb68e5 100644 --- a/Initial/kepler_softened.c +++ b/Initial/kepler_softened.c @@ -24,7 +24,7 @@ void initial( double * prim , double * x ){ double Pp; if(isothermal_flag) { - double cs2 = get_cs2(x); + double cs2 = get_cs2(x, rho); Pp = rho*cs2/gam; } else diff --git a/Initial/rigid_perturb.c b/Initial/rigid_perturb.c index 11e6959..37008a9 100644 --- a/Initial/rigid_perturb.c +++ b/Initial/rigid_perturb.c @@ -21,7 +21,7 @@ void initial( double * prim , double * x ){ int m = 4; rho += 0.1 * sin(m*(phi-omega*(*t))) * exp(-(r-0.5)*(r-0.5)/(2*0.1*0.1)); - double Pp = 0.01; //get_cs2(x) * rho; // 0.01; + double Pp = 0.01; double X = 0.0; diff --git a/Initial/shear_cart.c b/Initial/shear_cart.c index c3e927c..b69b114 100644 --- a/Initial/shear_cart.c +++ b/Initial/shear_cart.c @@ -1,4 +1,4 @@ - +#include "../omega.h" #include "../paul.h" static double nu = 0.0; @@ -20,7 +20,6 @@ void get_xyz(double *, double *); void get_vec_from_xyz(double *, double *, double *); void get_vec_contravariant(double *, double *, double *); -double get_cs2(double *x); void setICparams( struct domain * theDomain ) { @@ -55,7 +54,7 @@ void initial(double *prim, double *x) rho = exp(-rho_kappa * X); if(isothermal_flag) - P = get_cs2(x) * rho / gam; + P = get_cs2(x, rho) * rho / gam; else P = cs0*cs0 * rho / gam; diff --git a/Initial/singleSinkTest.c b/Initial/singleSinkTest.c index 4c497a5..757e6f6 100644 --- a/Initial/singleSinkTest.c +++ b/Initial/singleSinkTest.c @@ -115,7 +115,7 @@ void initial( double * prim , double * x ) double gr = grpz[0]; - double cs2 = get_cs2(x); + double cs2 = get_cs2(x, sig); double Pp = sig * cs2 / gam; double vr = 0.0; diff --git a/Initial/warp.c b/Initial/warp.c index 417c84d..df913cf 100644 --- a/Initial/warp.c +++ b/Initial/warp.c @@ -47,9 +47,9 @@ void initial( double * prim , double * x ){ double dz = z0-z; double xi = sqrt(ds*ds+dz*dz); - double cs2 = get_cs2(x); double visc = get_nu(x, prim); double sigma = 1.0/visc; + double cs2 = get_cs2(x, sigma); double h = sqrt(2.0*cs2)/omega double rho = sigma*exp(-xi*xi/(h*h)); diff --git a/omega.c b/omega.c index b11f442..d73065b 100644 --- a/omega.c +++ b/omega.c @@ -165,7 +165,7 @@ double get_height_om( const double *x){ return sqrt(omtot2); } -double get_cs2( const double *x ){ +double get_cs2( const double *x, double rho ){ double cs2; if(cs2Choice == 1) @@ -205,6 +205,10 @@ double get_cs2( const double *x ){ double scaling = pow(r, -cs2Par); // cs2 = scaling/(Mach*Mach); } + else if(cs2Choice == 7) + { + cs2 = pow(rho, cs2Par - 1.0)/(Mach*Mach); + } else cs2 = 1.0; diff --git a/omega.h b/omega.h index e6f0cef..a064319 100644 --- a/omega.h +++ b/omega.h @@ -12,7 +12,7 @@ double get_om1( const double *x); double get_om2( const double *x); double get_height_om( const double *x); -double get_cs2( const double *x ); +double get_cs2( const double *x, double rho ); double get_nu( const double *x , const double *prim); #endif diff --git a/sink.c b/sink.c index d4c480c..ad851d8 100644 --- a/sink.c +++ b/sink.c @@ -352,7 +352,7 @@ void cooling(const double *prim, double *cons, double beta = coolPar1; double x[3]; get_centroid_arr(xp, xm, x); - double enTarget = get_cs2(x)/(gamma_law*gm1); + double enTarget = get_cs2(x, rho)/(gamma_law*gm1); double enCurrent = press/(rho*gm1); double omtot = 0.0;