diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index b5aa9ac8..e2fd16dc 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -8,6 +8,7 @@ #include // For 3d rectangular grids #include +#include #include #include #include @@ -301,7 +302,7 @@ namespace BasisClasses { auto [ret, worst, lworst] = orthoCompute(orthoCheck(knots)); // For the CTest log - std::cout << "Spherical::orthoTest: worst=" << worst << std::endl; + std::cout << "---- Spherical::orthoTest: worst=" << worst << std::endl; return ret; } @@ -536,7 +537,7 @@ namespace BasisClasses { auto [ret, worst, lworst] = orthoCompute(orthoCheck()); // For the CTest log - std::cout << "FlatDisk::orthoTest: worst=" << worst << std::endl; + std::cout << "---- FlatDisk::orthoTest: worst=" << worst << std::endl; return ret; } @@ -741,18 +742,20 @@ namespace BasisClasses //@{ //! Basis construction parameters - double aratio, hratio, dweight, rwidth, ashift, rfactor, rtrunc, ppow; + double aratio, hratio, dweight, rwidth, ashift, rfactor, rtrunc, ppow, Mfac, HERNA; bool Ignore, deproject; + std::string pyname; //! DiskType support // - enum DiskType { constant, gaussian, mn, exponential, doubleexpon }; + enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; DiskType DTYPE; std::string mtype, dtype, dmodel; static const std::map dtlookup; double DiskDens(double R, double z, double phi); double dcond(double R, double z, double phi, int M); + std::shared_ptr pyDens; protected: @@ -837,7 +840,7 @@ namespace BasisClasses { auto [ret, worst, lworst] = orthoCompute(sl->orthoCheck()); // For the CTest log - std::cout << "Cylindrical::orthoTest: worst=" << worst << std::endl; + std::cout << "---- Cylindrical::orthoTest: worst=" << worst << std::endl; return ret; } }; @@ -981,7 +984,7 @@ namespace BasisClasses } } - std::cout << "Slab::orthoTest: worst=" << worst << std::endl; + std::cout << "---- Slab::orthoTest: worst=" << worst << std::endl; if (worst > __EXP__::orthoTol) return false; return true; } @@ -1103,7 +1106,7 @@ namespace BasisClasses } } - std::cout << "Cube::orthoTest: worst=" << worst << std::endl; + std::cout << "---- Cube::orthoTest: worst=" << worst << std::endl; if (worst > __EXP__::orthoTol) return false; return true; } diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index cd38730c..ecc320a9 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -60,6 +60,7 @@ namespace BasisClasses "self_consistent", "cachename", "modelname", + "pyname", "rnum" }; @@ -860,7 +861,9 @@ namespace BasisClasses {"gaussian", DiskType::gaussian}, {"mn", DiskType::mn}, {"exponential", DiskType::exponential}, - {"doubleexpon", DiskType::doubleexpon} + {"doubleexpon", DiskType::doubleexpon}, + {"diskbulge", DiskType::diskbulge}, + {"python", DiskType::python} }; const std::set @@ -911,6 +914,8 @@ namespace BasisClasses "aratio", "hratio", "dweight", + "Mfac", + "HERNA", "rwidth", "rfactor", "rtrunc", @@ -921,7 +926,8 @@ namespace BasisClasses "self_consistent", "playback", "coefCompute", - "coefMaster" + "coefMaster", + "pyname" }; Cylindrical::Cylindrical(const YAML::Node& CONF) : @@ -980,6 +986,22 @@ namespace BasisClasses w2*exp(-R/a2)/(4.0*M_PI*a2*a2*h2*f2*f2) ; } break; + + case DiskType::diskbulge: + { + double f = cosh(z/hcyl); + double rr = pow(pow(R, 2) + pow(z,2), 0.5); + double w1 = Mfac; + double w2 = 1.0 - Mfac; + double as = HERNA; + + ans = w1*exp(-R/acyl)/(4.0*M_PI*acyl*acyl*hcyl*f*f) + + w2*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0) ; + } + break; + case DiskType::python: + ans = (*pyDens)(R, z, phi); + break; case DiskType::exponential: default: { @@ -1068,6 +1090,12 @@ namespace BasisClasses // Ratio of second disk relative to the first disk for disk basis construction with double-exponential dweight = 1.0; + // mass fraction for disk for diskbulge + Mfac = 1.0; + + // Hernquist scale a disk basis construction with diskbulge + HERNA = 0.10; + // Width for erf truncation for EOF conditioning density (ignored if zero) rwidth = 0.0; @@ -1124,16 +1152,19 @@ namespace BasisClasses if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); if (conf["aratio" ]) aratio = conf["aratio" ].as(); - if (conf["hratio" ]) aratio = conf["hratio" ].as(); - if (conf["dweight" ]) aratio = conf["dweight" ].as(); - if (conf["rwidth" ]) aratio = conf["rwidth" ].as(); - if (conf["ashift" ]) aratio = conf["ashift" ].as(); + if (conf["hratio" ]) hratio = conf["hratio" ].as(); + if (conf["dweight" ]) dweight = conf["dweight" ].as(); + if (conf["Mfac" ]) Mfac = conf["Mfac" ].as(); + if (conf["HERNA" ]) HERNA = conf["HERNA" ].as(); + if (conf["rwidth" ]) rwidth = conf["rwidth" ].as(); + if (conf["ashift" ]) ashift = conf["ashift" ].as(); if (conf["rfactor" ]) rfactor = conf["rfactor" ].as(); if (conf["rtrunc" ]) rtrunc = conf["rtrunc" ].as(); if (conf["pow" ]) ppow = conf["ppow" ].as(); if (conf["mtype" ]) mtype = conf["mtype" ].as(); if (conf["dtype" ]) dtype = conf["dtype" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); + if (conf["pyname" ]) pyname = conf["pyname" ].as(); // Deprecation warning if (conf["density" ]) { @@ -1261,7 +1292,7 @@ namespace BasisClasses DTYPE = dtlookup.at(dtype); // key is not in the map. if (myid==0) // Report DiskType - std::cout << "DiskType is <" << dtype << ">" << std::endl; + std::cout << "---- DiskType is <" << dtype << ">" << std::endl; } catch (const std::out_of_range& err) { if (myid==0) { @@ -1273,6 +1304,12 @@ namespace BasisClasses throw std::runtime_error("Cylindrical::initialize: invalid DiskType"); } + // Check for and initialize the Python density type + // + if (DTYPE == DiskType::python) { + pyDens = std::make_shared(pyname); + } + // Use these user models to deproject for the EOF spherical basis // if (deproject) { diff --git a/exputil/CMakeLists.txt b/exputil/CMakeLists.txt index c0722c39..94c947ef 100644 --- a/exputil/CMakeLists.txt +++ b/exputil/CMakeLists.txt @@ -35,13 +35,14 @@ set(QPDISTF_SRC QPDistF.cc qld.c) set(SLEDGE_SRC sledge.f) set(PARTICLE_SRC Particle.cc ParticleReader.cc header.cc) set(CUDA_SRC cudaParticle.cu cudaSLGridMP2.cu) +set(PYWRAP_SRC DiskDensityFunc.cc) set(exputil_SOURCES ${ODE_SRC} ${ROOT_SRC} ${QUAD_SRC} ${RANDOM_SRC} ${UTIL_SRC} ${SPECFUNC_SRC} ${PHASE_SRC} ${SYMP_SRC} ${INTERP_SRC} ${MASSMODEL_SRC} ${ORBIT_SRC} ${BIORTH_SRC} ${POLY_SRC} ${GAUSS_SRC} ${QPDISTF_SRC} ${BESSEL_SRC} ${OPTIMIZATION_SRC} - ${SLEDGE_SRC} ${PARTICLE_SRC} ${CUDA_SRC}) + ${SLEDGE_SRC} ${PARTICLE_SRC} ${CUDA_SRC} ${PYWRAP_SRC}) set(common_INCLUDE_DIRS $ $ @@ -53,7 +54,7 @@ set(common_INCLUDE_DIRS $ set(common_LINKLIB ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp ${VTK_LIBRARIES} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES} - ${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB}) + ${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB} pybind11::embed) if(ENABLE_CUDA) list(APPEND common_LINKLIB CUDA::toolkit CUDA::cudart) diff --git a/exputil/DiskDensityFunc.cc b/exputil/DiskDensityFunc.cc new file mode 100644 index 00000000..2c6f00bc --- /dev/null +++ b/exputil/DiskDensityFunc.cc @@ -0,0 +1,31 @@ +#include + +namespace py = pybind11; + +DiskDensityFunc::DiskDensityFunc(const std::string& modulename, + const std::string& funcname) + : funcname(funcname) +{ + // Check if a Python interpreter exists + if (Py_IsInitialized() == 0) { + py::initialize_interpreter(); + // Mark the interpreter as started by this instance + started = true; + } + + // Bind the disk_density function from Python + disk_density = + py::reinterpret_borrow + ( py::module::import(modulename.c_str()).attr(funcname.c_str()) ); +} + +DiskDensityFunc::~DiskDensityFunc() +{ + // Only end the interpreter if it was started by this instance + if (started) py::finalize_interpreter(); +} + +double DiskDensityFunc::operator() (double R, double z, double phi) +{ + return disk_density(R, z, phi).cast(); +} diff --git a/include/DiskDensityFunc.H b/include/DiskDensityFunc.H new file mode 100644 index 00000000..0469214c --- /dev/null +++ b/include/DiskDensityFunc.H @@ -0,0 +1,61 @@ +#ifndef _DiskDensityFunc_H_ +#define _DiskDensityFunc_H_ + +#include + +#include +#include + +/** Python disk-density function wrapper + + The wrapper class will take a user supplied Python module that + must supply the 'disk-density' function + + An example disk-density Python module: + --------------------------cut here-------------------------------- + import math + + def disk_density(R, z, phi): + """Computes the disk density at a given radius, height, and + azimuth. This would be the user's target density for basis + creation. + + """ + a = 1.0 # Scale radius + h = 0.2 # Scale height + f = math.exp(-math.fabs(z)/h) # Prevent overflows + sech = 2.0*f / (1.0 + f*f) # + return math.exp(-R/a)*sech*sech/(4*math.pi*h*a*a) + + --------------------------cut here-------------------------------- + + */ +class __attribute__((visibility("default"))) +DiskDensityFunc +{ +private: + + //! The disk-density function + pybind11::function disk_density; + + //! Interpreter started? + bool started = false; + + //! Name of density function + std::string funcname; + +public: + + //! Constructor + DiskDensityFunc(const std::string& modulename, + const std::string& funcname="disk_density"); + + //! Destructor + ~DiskDensityFunc(); + + //! Evaluate the density + double operator() (double R, double z, double phi); + +}; + +#endif diff --git a/include/DiskModels.H b/include/DiskModels.H index 3440794f..2b4c75b0 100644 --- a/include/DiskModels.H +++ b/include/DiskModels.H @@ -2,7 +2,7 @@ #define _DISK_MODELS_H #include - +#include //! A Ferrers Ellipsoid + Evacuated Exponential Disc (semi-realistic //! bar+disk model) @@ -223,5 +223,58 @@ public: }; -#endif +//! The usual exponential disk + a Hernquist bulge +class diskbulge : public EmpCylSL::AxiDisk +{ +private: + double a, h, as, Mfac, d1, d2; + +public: + + diskbulge(double a, double h, double as, double Mfac, double M=1) : + a(a), h(h), as(as), Mfac(Mfac), EmpCylSL::AxiDisk(M, "diskbulge") + { + params.push_back(a); + params.push_back(h); + params.push_back(as); + params.push_back(Mfac); + } + + double operator()(double R, double z, double phi=0.) + { + double sh = 1.0/cosh(z/h); + d1 = 0.25*M*Mfac/(M_PI*a*a*h) * exp(-R/a) * sh * sh; + + double rr = pow(pow(R, 2) + pow(z,2), 0.5); + d2 = M*(1-Mfac)*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0); + + return d1 + d2; + + } + +}; + + +//! A user-defined Python model +class PyModel : public EmpCylSL::AxiDisk +{ +private: + + std::shared_ptr pyDens; + +public: + + PyModel(std::string& pyname) + { + pyDens = std::make_shared(pyname); + } + + double operator()(double R, double z, double phi=0.) + { + return (*pyDens)(R, z, phi); + } + +}; + +#endif diff --git a/src/Cylinder.H b/src/Cylinder.H index 08e012c5..0faf7dda 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -103,6 +103,8 @@ class MixtureBasis; @param playback file reads a coefficient file and uses it to compute the basis function output for resimiulation + @param python is the file name of Python module which supplies the 'disk_density' function for conditioning the cylindrical basis. A non-null string triggers the use of the Python interpreter to evaluate the target density function. + */ class Cylinder : public Basis { @@ -130,7 +132,7 @@ private: int ncylnx, ncylny, ncylr; double hcyl, hexp, snr, rem; int nmax, ncylodd, ncylrecomp, npca, npca0, nvtk, cmapR, cmapZ; - std::string cachename; + std::string cachename, pyname; bool self_consistent, logarithmic, pcavar, pcainit, pcavtk, pcadiag, pcaeof; bool try_cache, firstime, dump_basis, compute, firstime_coef; diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 80eecc83..af8103e0 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -10,52 +10,17 @@ #include #include #include +#include #include #include #include -Timer timer_debug; - -double EXPSCALE=1.0, HSCALE=1.0, ASHIFT=0.25; - - //@{ //! These are for testing exclusively (should be set false for production) static bool cudaAccumOverride = false; static bool cudaAccelOverride = false; //@} -double DiskDens(double R, double z, double phi) -{ - double f = cosh(z/HSCALE); - return exp(-R/EXPSCALE)/(4.0*M_PI*EXPSCALE*EXPSCALE*HSCALE*f*f); -} - - -double dcond(double R, double z, double phi, int M) -{ - // - // No shift for M==0 - // - if (M==0) return DiskDens(R, z, phi); - - // - // Fold into [-PI/M, PI/M] for M>=1 - // - double dmult = M_PI/M, phiS; - if (phi>M_PI) - phiS = phi + dmult*(int)((2.0*M_PI - phi)/dmult); - else - phiS = phi - dmult*(int)(phi/dmult); - - // - // Apply a shift along the x-axis - // - double x = R*cos(phiS) - ASHIFT*EXPSCALE; - double y = R*sin(phiS); - return DiskDens(sqrt(x*x + y*y), z, atan2(y, x)); -} - const std::set Cylinder::valid_keys = { @@ -107,6 +72,7 @@ Cylinder::valid_keys = { "playback", "coefCompute", "coefMaster", + "pyname", "dumpbasis" }; @@ -238,9 +204,6 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : // Set parameters for external dcond function // - EXPSCALE = acyl; - HSCALE = hcyl; - ASHIFT = ashift; eof = 0; bool cache_ok = false; @@ -300,6 +263,45 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : << "this step will take many minutes." << std::endl; + std::shared_ptr dfunc; + if (pyname.size()) dfunc = std::make_shared(pyname); + + // Use either the user-supplied Python target or the default + // exponential disk model + auto DiskDens = [&](double R, double z, double phi) + { + if (dfunc) return (*dfunc)(R, z, phi); + double f = exp(-fabs(z)/hcyl); // Overflow prevention + double s = 2.0*f/(1.0 + f*f); // in sech computation + return exp(-R/acyl)*s*s/(4.0*M_PI*acyl*acyl*hcyl); + }; + + // The conditioning function for the EOF with an optional shift + // for M>0 + auto dcond = [&](double R, double z, double phi, int M) + { + // + // No shift for M==0 + // + if (M==0) return DiskDens(R, z, phi); + + // + // Fold into [-PI/M, PI/M] for M>=1 + // + double dmult = M_PI/M, phiS; + if (phi>M_PI) + phiS = phi + dmult*(int)((2.0*M_PI - phi)/dmult); + else + phiS = phi - dmult*(int)(phi/dmult); + + // + // Apply a shift along the x-axis + // + double x = R*cos(phiS) - ashift*acyl; + double y = R*sin(phiS); + return DiskDens(sqrt(x*x + y*y), z, atan2(y, x)); + }; + ortho->generate_eof(rnum, pnum, tnum, dcond); } diff --git a/utils/ICs/EllipsoidForce.H b/utils/ICs/EllipsoidForce.H index 5ffcb68c..6bc56a69 100644 --- a/utils/ICs/EllipsoidForce.H +++ b/utils/ICs/EllipsoidForce.H @@ -7,6 +7,8 @@ #include #include +#include + #include #include @@ -66,7 +68,7 @@ private: vector dX, table; double rtmin, rtmax; bool quadr, tmade, tlog; - bool nindx(vector& x, vector& n); + bool nindx(Eigen::Vector3d& x, Eigen::Vector3i& n); void write_cache(); bool read_cache(); diff --git a/utils/ICs/EllipsoidForce.cc b/utils/ICs/EllipsoidForce.cc index c3f58bbc..3d918d65 100644 --- a/utils/ICs/EllipsoidForce.cc +++ b/utils/ICs/EllipsoidForce.cc @@ -54,7 +54,7 @@ EllipsoidForce::EllipsoidForce(double a0, double a1, double a2, double mass, if (quadr) { - vector x(3, 0); + Eigen::Vector3d x; x.setZero(); lrmin = log(rmin); lrmax = log(rmax); ldr = (lrmax - lrmin)/n; @@ -494,7 +494,7 @@ void EllipsoidForce::MakeTable(int n1, int n2, int n3) tmade = true; } -bool EllipsoidForce::nindx(vector& x, vector& n) +bool EllipsoidForce::nindx(Eigen::Vector3d& x, Eigen::Vector3i& n) { for (int i=0; i<3; i++) { if (tlog) { @@ -536,13 +536,14 @@ void EllipsoidForce::TableEval(vector x, vector& force) double f; int indx; - vector x1(x), a(3); - vector n(3), n1(3); + Eigen::Vector3d x1, a; + Eigen::Vector3i n, n1; // // Put lower limit on the grid interpolation // for (int i=0; i<3; i++) { + x1[i] = x[i]; if (tlog) x1[i] = max(fabs(x1[i]), exp(rtmin)); else diff --git a/utils/ICs/check_coefs.cc b/utils/ICs/check_coefs.cc index a9faf3af..02d75a71 100644 --- a/utils/ICs/check_coefs.cc +++ b/utils/ICs/check_coefs.cc @@ -515,7 +515,7 @@ main(int ac, char **av) // Report dtype // if (myid==0) - std::cout << "DiskType is <" << disktype << ">" << std::endl; + std::cout << "---- DiskType is <" << disktype << ">" << std::endl; //==================== // OpenMP control diff --git a/utils/ICs/cylcache.cc b/utils/ICs/cylcache.cc index fd78cdbf..a4da8071 100644 --- a/utils/ICs/cylcache.cc +++ b/utils/ICs/cylcache.cc @@ -165,20 +165,25 @@ void set_fpu_gdb_handler(void) // Global variables // -enum DiskType { constant, gaussian, mn, exponential, doubleexpon }; +enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; std::map dtlookup = { {"constant", DiskType::constant}, {"gaussian", DiskType::gaussian}, {"mn", DiskType::mn}, {"exponential", DiskType::exponential}, - {"doubleexpon", DiskType::doubleexpon} + {"doubleexpon", DiskType::doubleexpon}, + {"diskbulge", DiskType::diskbulge}, + {"python", DiskType::python} }; DiskType DTYPE; +std::string pyname; double ASCALE; double ASHIFT; double HSCALE; +double HERNA; +double Mfac; double RTRUNC = 1.0; double RWIDTH = 0.0; double ARATIO = 1.0; @@ -187,6 +192,8 @@ double DWEIGHT = 1.0; #include +static std::shared_ptr pydens; + double DiskDens(double R, double z, double phi) { double ans = 0.0; @@ -230,6 +237,24 @@ double DiskDens(double R, double z, double phi) w2*exp(-R/a2)/(4.0*M_PI*a2*a2*h2*f2*f2) ; } break; + case DiskType::diskbulge: + { + double f = cosh(z/HSCALE); + double rr = pow(pow(R, 2) + pow(z,2), 0.5); + double w1 = Mfac; + double w2 = (1-Mfac); + double as = HERNA; + + ans = w1*exp(-R/ASCALE)/(4.0*M_PI*ASCALE*ASCALE*HSCALE*f*f) + + w2*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0) ; + } + break; + case DiskType::python: + { + if (not pydens) pydens = std::make_shared(pyname); + ans = (*pydens)(R, z, phi); + } + break; case DiskType::exponential: default: { @@ -327,7 +352,7 @@ main(int ac, char **av) ("ctype", "DiskHalo radial coordinate scaling type (one of: Linear, Log, Rat)", cxxopts::value(ctype)->default_value("Log")) ("LMAXFID", "Maximum angular order for spherical basis in adaptive construction of the cylindrical basis", - cxxopts::value(LMAXFID)->default_value("72")) + cxxopts::value(LMAXFID)->default_value("32")) ("NMAXFID", "Maximum radial order for the spherical basis in adapative construction of the cylindrical basis", cxxopts::value(NMAXFID)->default_value("32")) ("MMAX", "Maximum azimuthal order for the cylindrical basis", @@ -341,7 +366,7 @@ main(int ac, char **av) ("NCYLODD", "Number of vertically odd basis functions per harmonic order", cxxopts::value(NCYLODD)->default_value("6")) ("NMAX", "Total number of basis functions per harmonic order", - cxxopts::value(NMAX)->default_value("18")) + cxxopts::value(NMAX)->default_value("20")) ("VFLAG", "Diagnostic flag for EmpCylSL", cxxopts::value(VFLAG)->default_value("31")) ("expcond", "Use analytic target density rather than particle distribution", @@ -368,6 +393,10 @@ main(int ac, char **av) cxxopts::value(PPower)->default_value("4.0")) ("DWEIGHT", "Ratio of second disk relative to the first disk for disk basis construction with double-exponential", cxxopts::value(DWEIGHT)->default_value("1.0")) + ("Mfac", "Mass fraction of the disk for diskbulge model, sets disk/bulge mass (e.g. 0.75 -> 75% of mass in disk, 25% of mass in bulge)", + cxxopts::value(Mfac)->default_value("1.0")) + ("HERNA", "Hernquist scale a parameter in diskbulge model", + cxxopts::value(HERNA)->default_value("0.10")) ("RTRUNC", "Maximum disk radius for erf truncation of EOF conditioning density", cxxopts::value(RTRUNC)->default_value("0.1")) ("RWIDTH", "Width for erf truncation for EOF conditioning density (ignored if zero)", @@ -384,6 +413,8 @@ main(int ac, char **av) cxxopts::value(CMAPR)->default_value("1")) ("CMAPZ", "Vertical coordinate mapping type for cylindrical grid (0=none, 1=rational fct)", cxxopts::value(CMAPZ)->default_value("1")) + ("pyname", "The name of the Python module supplying the disk_density function", + cxxopts::value(pyname)) ; cxxopts::ParseResult vm; @@ -490,7 +521,7 @@ main(int ac, char **av) DTYPE = dtlookup.at(dtype); // key is not in the map. if (myid==0) // Report DiskType - std::cout << "DiskType is <" << dtype << ">" << std::endl; + std::cout << "---- DiskType is <" << dtype << ">" << std::endl; } catch (const std::out_of_range& err) { if (myid==0) { @@ -556,6 +587,8 @@ main(int ac, char **av) << " rmax=" << EmpCylSL::RMAX << " a=" << ASCALE << " h=" << HSCALE + << " as=" << HERNA + << " Mfac=" << Mfac << " nmax2=" << NMAXFID << " lmax2=" << LMAXFID << " mmax=" << MMAX diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index 5a8825aa..f0890ed2 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -92,6 +92,7 @@ #include #include #include +#include #include #include #include @@ -239,17 +240,20 @@ const double f_H = 0.76; // Global variables // -enum DiskType { constant, gaussian, mn, exponential, doubleexpon }; +enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; std::map dtlookup = { {"constant", DiskType::constant}, {"gaussian", DiskType::gaussian}, {"mn", DiskType::mn}, {"exponential", DiskType::exponential}, - {"doubleexpon", DiskType::doubleexpon} + {"doubleexpon", DiskType::doubleexpon}, + {"diskbulge", DiskType::diskbulge}, + {"python", DiskType::python} }; DiskType DTYPE; +std::string pyname; double ASCALE; double ASHIFT; double HSCALE; @@ -258,12 +262,15 @@ double RWIDTH = 0.0; double ARATIO = 1.0; double HRATIO = 1.0; double DWEIGHT = 1.0; +double Mfac = 1.0; +double HERNA = 0.10; #include double DiskDens(double R, double z, double phi) { double ans = 0.0; + static std::shared_ptr dfunc; switch (DTYPE) { @@ -304,6 +311,27 @@ double DiskDens(double R, double z, double phi) w2*exp(-R/a2)/(4.0*M_PI*a2*a2*h2*f2*f2) ; } break; + + case DiskType::diskbulge: + { + double acyl = ASCALE; + double hcyl = HSCALE; + double f = cosh(z/HSCALE); + double rr = pow(pow(R, 2) + pow(z,2), 0.5); + double w1 = Mfac; + double w2 = (1-Mfac); + double as = HERNA; + + ans = w1*exp(-R/acyl)/(4.0*M_PI*acyl*acyl*hcyl*f*f) + + w2*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0); + } + break; + case DiskType::python: + { + if (!dfunc) dfunc = std::make_shared(pyname); + ans = (*dfunc)(R, z, phi); + } + break; case DiskType::exponential: default: { @@ -529,6 +557,10 @@ main(int ac, char **av) cxxopts::value(ARATIO)->default_value("1.0")) ("HRATIO", "Vertical scale height ratio for disk basis construction with doubleexpon", cxxopts::value(HRATIO)->default_value("1.0")) + ("Mfac", "Mass fraction of the disk for diskbulge model, sets disk/bulge mass (e.g. 0.75 -> 75% of mass in disk, 25% of mass in bulge)", + cxxopts::value(Mfac)->default_value("1.0")) + ("HERNA", "Hernquist scale a parameter in diskbulge model", + cxxopts::value(HERNA)->default_value("0.10")) ("DWEIGHT", "Ratio of second disk relative to the first disk for disk basis construction with double-exponential", cxxopts::value(DWEIGHT)->default_value("1.0")) ("RTRUNC", "Maximum disk radius for erf truncation of EOF conditioning density", @@ -576,6 +608,8 @@ main(int ac, char **av) ("allow", "Allow multimass algorithm to generature negative masses for testing") ("nomono", "Allow non-monotonic mass interpolation") ("diskmodel", "Table describing the model for the disk plane") + ("pyname", "Name of module with the user-specified target disk density", + cxxopts::value(pyname)) ; cxxopts::ParseResult vm; @@ -685,7 +719,7 @@ main(int ac, char **av) DTYPE = dtlookup.at(dtype); // key is not in the map. if (myid==0) // Report DiskType - std::cout << "DiskType is <" << dtype << ">" << std::endl; + std::cout << "---- DiskType is <" << dtype << ">" << std::endl; } catch (const std::out_of_range& err) { if (myid==0) {