From 5cd05f78f9a47a1452d4eafcdf44f047e55f8889 Mon Sep 17 00:00:00 2001 From: CarrieFilion <15cfilion@gmail.com> Date: Mon, 23 Dec 2024 16:01:30 -0500 Subject: [PATCH 01/14] adding an exponential disk + Hernquist bulge model --- expui/BiorthBasis.H | 4 ++-- expui/BiorthBasis.cc | 35 ++++++++++++++++++++++++++++++----- include/DiskModels.H | 33 ++++++++++++++++++++++++++++++++- utils/ICs/cylcache.cc | 33 +++++++++++++++++++++++++++++---- utils/ICs/initial.cc | 26 ++++++++++++++++++++++++-- 5 files changed, 117 insertions(+), 14 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index a4af27c1a..6d3b41500 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -590,12 +590,12 @@ 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; //! DiskType support // - enum DiskType { constant, gaussian, mn, exponential, doubleexpon }; + enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge }; DiskType DTYPE; std::string mtype, dtype, dmodel; diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index e39bfe825..62522170c 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -860,7 +860,8 @@ namespace BasisClasses {"gaussian", DiskType::gaussian}, {"mn", DiskType::mn}, {"exponential", DiskType::exponential}, - {"doubleexpon", DiskType::doubleexpon} + {"doubleexpon", DiskType::doubleexpon}, + {"diskbulge", DiskType::diskbulge} }; const std::set @@ -911,6 +912,8 @@ namespace BasisClasses "aratio", "hratio", "dweight", + "Mfac", + "HERNA", "rwidth", "rfactor", "rtrunc", @@ -980,6 +983,19 @@ 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-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::exponential: default: { @@ -1062,12 +1078,19 @@ namespace BasisClasses // Radial scale length ratio for disk basis construction with doubleexpon aratio = 1.0; + // Vertical scale height ratio for disk basis construction with doubleexpon hratio = 1.0; // 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,10 +1147,12 @@ 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(); diff --git a/include/DiskModels.H b/include/DiskModels.H index 3440794fe..6c6ed3c3d 100644 --- a/include/DiskModels.H +++ b/include/DiskModels.H @@ -223,5 +223,36 @@ 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; + + } + +}; + +#endif \ No newline at end of file diff --git a/utils/ICs/cylcache.cc b/utils/ICs/cylcache.cc index fd78cdbf8..20671a4c9 100644 --- a/utils/ICs/cylcache.cc +++ b/utils/ICs/cylcache.cc @@ -165,20 +165,23 @@ void set_fpu_gdb_handler(void) // Global variables // -enum DiskType { constant, gaussian, mn, exponential, doubleexpon }; +enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge }; std::map dtlookup = { {"constant", DiskType::constant}, {"gaussian", DiskType::gaussian}, {"mn", DiskType::mn}, {"exponential", DiskType::exponential}, - {"doubleexpon", DiskType::doubleexpon} + {"doubleexpon", DiskType::doubleexpon}, + {"diskbulge", DiskType::diskbulge} }; DiskType DTYPE; double ASCALE; double ASHIFT; double HSCALE; +double HERNA; +double Mfac; double RTRUNC = 1.0; double RWIDTH = 0.0; double ARATIO = 1.0; @@ -230,6 +233,18 @@ 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::exponential: default: { @@ -327,7 +342,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 +356,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 +383,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 +403,10 @@ 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")) + ("HERNA", "The Hernquist 'a' scale radius, diskbulge model only", + cxxopts::value(dmodel)->default_value("0.1")) + ("Mfac", "Fraction of mass in the disk, remaining mass will be in the bulge (i.e. Mfac = .75 would have 75% of the mass be disk, 25% bulge). diskbulge model only", + cxxopts::value(dmodel)->default_value("0.1")) ; cxxopts::ParseResult vm; @@ -556,6 +579,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 5a8825aa6..7cbc9f5b8 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -239,14 +239,15 @@ const double f_H = 0.76; // Global variables // -enum DiskType { constant, gaussian, mn, exponential, doubleexpon }; +enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge }; std::map dtlookup = { {"constant", DiskType::constant}, {"gaussian", DiskType::gaussian}, {"mn", DiskType::mn}, {"exponential", DiskType::exponential}, - {"doubleexpon", DiskType::doubleexpon} + {"doubleexpon", DiskType::doubleexpon}, + {"diskbulge", DiskType::diskbulge} }; DiskType DTYPE; @@ -258,6 +259,8 @@ 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 @@ -304,6 +307,21 @@ 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::exponential: default: { @@ -529,6 +547,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", From 4fb989465a3cf3a2f9d577f407bed30dc864afae Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Dec 2024 17:18:57 -0500 Subject: [PATCH 02/14] Implementation of a Python functor for the target density function in 'pyEXP' and 'exp' --- expui/BiorthBasis.H | 5 ++- expui/BiorthBasis.cc | 34 +++++++++++------ exputil/CMakeLists.txt | 5 ++- include/DiskDensityFunc.H | 33 +++++++++++++++++ include/DiskModels.H | 26 ++++++++++++- src/Cylinder.H | 4 +- src/Cylinder.cc | 77 +++++++++++++++++++-------------------- utils/ICs/cylcache.cc | 15 +++++++- 8 files changed, 140 insertions(+), 59 deletions(-) create mode 100644 include/DiskDensityFunc.H diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 6d3b41500..41d32d278 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -8,6 +8,7 @@ #include // For 3d rectangular grids #include +#include #include #include #include @@ -592,16 +593,18 @@ namespace BasisClasses //! Basis construction parameters 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, diskbulge }; + 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: diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 62522170c..e2bba7ed1 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -60,6 +60,7 @@ namespace BasisClasses "self_consistent", "cachename", "modelname", + "pyname", "rnum" }; @@ -861,7 +862,8 @@ namespace BasisClasses {"mn", DiskType::mn}, {"exponential", DiskType::exponential}, {"doubleexpon", DiskType::doubleexpon}, - {"diskbulge", DiskType::diskbulge} + {"diskbulge", DiskType::diskbulge}, + {"python", DiskType::python} }; const std::set @@ -924,7 +926,8 @@ namespace BasisClasses "self_consistent", "playback", "coefCompute", - "coefMaster" + "coefMaster", + "pyname" }; Cylindrical::Cylindrical(const YAML::Node& CONF) : @@ -986,16 +989,19 @@ namespace BasisClasses 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-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) ; + 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: { @@ -1078,7 +1084,6 @@ namespace BasisClasses // Radial scale length ratio for disk basis construction with doubleexpon aratio = 1.0; - // Vertical scale height ratio for disk basis construction with doubleexpon hratio = 1.0; @@ -1159,6 +1164,7 @@ namespace BasisClasses 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" ]) { @@ -1298,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 c785e5817..e1e5db1c7 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/include/DiskDensityFunc.H b/include/DiskDensityFunc.H new file mode 100644 index 000000000..e9b0f1cb3 --- /dev/null +++ b/include/DiskDensityFunc.H @@ -0,0 +1,33 @@ +#ifndef _DiskDensityFunc_H_ +#define _DiskDensityFunc_H_ + +#include + +#include +#include + +// Convenience +namespace py = pybind11; + +//! Python disk-density function wrapper +class DiskDensityFunc +{ +private: + + //! The disk-density function + py::function disk_density; + +public: + + //! Constructor + DiskDensityFunc(const std::string& modulename); + + //! 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 6c6ed3c3d..2b4c75b08 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) @@ -255,4 +255,26 @@ public: }; -#endif \ No newline at end of file + +//! 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 08e012c58..0faf7ddaa 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 e355bcdf9..ad369eb60 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 = { @@ -106,7 +71,8 @@ Cylinder::valid_keys = { "self_consistent", "playback", "coefCompute", - "coefMaster" + "coefMaster", + "pyname" }; Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : @@ -237,9 +203,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; @@ -299,6 +262,40 @@ 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); + + auto DiskDens = [&](double R, double z, double phi) + { + if (dfunc) return (*dfunc)(R, z, phi); + double f = cosh(z/hcyl); + return exp(-R/acyl)/(4.0*M_PI*acyl*acyl*hcyl*f*f); + }; + + 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/cylcache.cc b/utils/ICs/cylcache.cc index 20671a4c9..95daf0089 100644 --- a/utils/ICs/cylcache.cc +++ b/utils/ICs/cylcache.cc @@ -165,7 +165,7 @@ void set_fpu_gdb_handler(void) // Global variables // -enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge }; +enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; std::map dtlookup = { {"constant", DiskType::constant}, @@ -173,10 +173,12 @@ std::map dtlookup = {"mn", DiskType::mn}, {"exponential", DiskType::exponential}, {"doubleexpon", DiskType::doubleexpon}, - {"diskbulge", DiskType::diskbulge} + {"diskbulge", DiskType::diskbulge}, + {"python", DiskType::python} }; DiskType DTYPE; +std::string pyname; double ASCALE; double ASHIFT; double HSCALE; @@ -190,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; @@ -245,6 +249,11 @@ double DiskDens(double R, double z, double phi) 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); + } case DiskType::exponential: default: { @@ -407,6 +416,8 @@ main(int ac, char **av) cxxopts::value(dmodel)->default_value("0.1")) ("Mfac", "Fraction of mass in the disk, remaining mass will be in the bulge (i.e. Mfac = .75 would have 75% of the mass be disk, 25% bulge). diskbulge model only", cxxopts::value(dmodel)->default_value("0.1")) + ("pyname", "The name of the python module supplying the disk_density", + cxxopts::value(pyname)) ; cxxopts::ParseResult vm; From 4ec255e2117666e83debcb9ed4efe4afcf932cd6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Dec 2024 23:05:35 -0500 Subject: [PATCH 03/14] Missing commit of this source file, ooops. --- exputil/DiskDensityFunc.cc | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 exputil/DiskDensityFunc.cc diff --git a/exputil/DiskDensityFunc.cc b/exputil/DiskDensityFunc.cc new file mode 100644 index 000000000..02b439df1 --- /dev/null +++ b/exputil/DiskDensityFunc.cc @@ -0,0 +1,23 @@ +#include + +namespace py = pybind11; // Convenience + + +DiskDensityFunc::DiskDensityFunc(const std::string& modulename) +{ + py::initialize_interpreter(); + + disk_density = + py::reinterpret_borrow + ( py::module::import(modulename.c_str()).attr("disk_density") ); +} + +DiskDensityFunc::~DiskDensityFunc() +{ + py::finalize_interpreter(); +} + +double DiskDensityFunc::operator() (double R, double z, double phi) +{ + return disk_density(R, z, phi).cast(); +} From dc566edf82bdcd5ea58526bce37176a37d06d787 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 25 Dec 2024 10:50:19 -0500 Subject: [PATCH 04/14] Missing parameter in getHeader_hdf5 --- exputil/EmpCylSL.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 01b941501..17c9b51e3 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -1253,6 +1253,7 @@ YAML::Node EmpCylSL::getHeader_hdf5(const std::string& cachefile) node["mmax"] = getInt("mmax"); node["numx"] = getInt("numx"); node["numy"] = getInt("numy"); + node["lmaxfid"] = getInt("lmaxfid"); node["nmaxfid"] = getInt("nmaxfid"); node["nmax"] = getInt("nmax"); node["neven"] = getInt("neven"); From 92c31dd81cc8777987db76dad96f524fc79da21c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 25 Dec 2024 10:51:42 -0500 Subject: [PATCH 05/14] Added documenation and an example Python module --- include/DiskDensityFunc.H | 33 +++++++++++++++++++++++++++------ 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/include/DiskDensityFunc.H b/include/DiskDensityFunc.H index e9b0f1cb3..e8bb1eb3a 100644 --- a/include/DiskDensityFunc.H +++ b/include/DiskDensityFunc.H @@ -6,16 +6,37 @@ #include #include -// Convenience -namespace py = pybind11; - -//! Python disk-density function wrapper -class DiskDensityFunc +/** 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 - py::function disk_density; + pybind11::function disk_density; public: From 6c2fbcdce37e594b71adfbe14380825e7334c274 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 25 Dec 2024 10:52:20 -0500 Subject: [PATCH 06/14] Duplicated Hernquist parameters for diskbulge --- utils/ICs/cylcache.cc | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/utils/ICs/cylcache.cc b/utils/ICs/cylcache.cc index 95daf0089..79959422f 100644 --- a/utils/ICs/cylcache.cc +++ b/utils/ICs/cylcache.cc @@ -254,6 +254,7 @@ double DiskDens(double R, double z, double phi) if (not pydens) pydens = std::make_shared(pyname); ans = (*pydens)(R, z, phi); } + break; case DiskType::exponential: default: { @@ -412,11 +413,7 @@ 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")) - ("HERNA", "The Hernquist 'a' scale radius, diskbulge model only", - cxxopts::value(dmodel)->default_value("0.1")) - ("Mfac", "Fraction of mass in the disk, remaining mass will be in the bulge (i.e. Mfac = .75 would have 75% of the mass be disk, 25% bulge). diskbulge model only", - cxxopts::value(dmodel)->default_value("0.1")) - ("pyname", "The name of the python module supplying the disk_density", + ("pyname", "The name of the Python module supplying the disk_density function", cxxopts::value(pyname)) ; From ce62bb47f25d46abe86413bbb6edc12ea4078220 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 25 Dec 2024 13:15:43 -0500 Subject: [PATCH 07/14] Check for existence of interpreter and only start a new one if needed --- exputil/DiskDensityFunc.cc | 10 ++++++---- include/DiskDensityFunc.H | 3 +++ 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/exputil/DiskDensityFunc.cc b/exputil/DiskDensityFunc.cc index 02b439df1..61247a74e 100644 --- a/exputil/DiskDensityFunc.cc +++ b/exputil/DiskDensityFunc.cc @@ -1,11 +1,13 @@ #include -namespace py = pybind11; // Convenience - +namespace py = pybind11; DiskDensityFunc::DiskDensityFunc(const std::string& modulename) { - py::initialize_interpreter(); + if (Py_IsInitialized() == 0) { + py::initialize_interpreter(); + started = true; + } disk_density = py::reinterpret_borrow @@ -14,7 +16,7 @@ DiskDensityFunc::DiskDensityFunc(const std::string& modulename) DiskDensityFunc::~DiskDensityFunc() { - py::finalize_interpreter(); + if (started) py::finalize_interpreter(); } double DiskDensityFunc::operator() (double R, double z, double phi) diff --git a/include/DiskDensityFunc.H b/include/DiskDensityFunc.H index e8bb1eb3a..1ad58620a 100644 --- a/include/DiskDensityFunc.H +++ b/include/DiskDensityFunc.H @@ -38,6 +38,9 @@ private: //! The disk-density function pybind11::function disk_density; + //! Interpreter started? + bool started = false; + public: //! Constructor From d593f2d71359739e8332db3a7b026bf205529979 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 25 Dec 2024 13:17:07 -0500 Subject: [PATCH 08/14] For consistency in stdout info --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index e2bba7ed1..80d6f3bdc 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1292,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) { From e3e78f4b09c07eb8d935fe98063fd5f40f8d8272 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 25 Dec 2024 16:10:57 -0500 Subject: [PATCH 09/14] Consistency in stdout info --- expui/BiorthBasis.H | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 41d32d278..ee31673d1 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -302,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; } @@ -537,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; } @@ -689,7 +689,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; } }; @@ -833,7 +833,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; } @@ -955,7 +955,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; } From 86c8c4e74e925089e9066c6e721f9e5fb581a0c6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 25 Dec 2024 16:11:54 -0500 Subject: [PATCH 10/14] Allow the function name to be user settable, although we are not currently using this feature --- exputil/DiskDensityFunc.cc | 10 ++++++++-- include/DiskDensityFunc.H | 6 +++++- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/exputil/DiskDensityFunc.cc b/exputil/DiskDensityFunc.cc index 61247a74e..2c6f00bc5 100644 --- a/exputil/DiskDensityFunc.cc +++ b/exputil/DiskDensityFunc.cc @@ -2,20 +2,26 @@ namespace py = pybind11; -DiskDensityFunc::DiskDensityFunc(const std::string& modulename) +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("disk_density") ); + ( 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(); } diff --git a/include/DiskDensityFunc.H b/include/DiskDensityFunc.H index 1ad58620a..0469214cb 100644 --- a/include/DiskDensityFunc.H +++ b/include/DiskDensityFunc.H @@ -41,10 +41,14 @@ private: //! Interpreter started? bool started = false; + //! Name of density function + std::string funcname; + public: //! Constructor - DiskDensityFunc(const std::string& modulename); + DiskDensityFunc(const std::string& modulename, + const std::string& funcname="disk_density"); //! Destructor ~DiskDensityFunc(); From b5789498b77ff998821e5954c165b069a55e806e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 26 Dec 2024 10:53:54 -0500 Subject: [PATCH 11/14] Missing instantiation in 'initial.cc' --- src/Cylinder.cc | 9 +++++++-- utils/ICs/initial.cc | 16 ++++++++++++++-- 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/Cylinder.cc b/src/Cylinder.cc index ad369eb60..98773a411 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -265,13 +265,18 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : 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 = cosh(z/hcyl); - return exp(-R/acyl)/(4.0*M_PI*acyl*acyl*hcyl*f*f); + 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) { // diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index 7cbc9f5b8..f0ad4b6d2 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -92,6 +92,7 @@ #include #include #include +#include #include #include #include @@ -239,7 +240,7 @@ const double f_H = 0.76; // Global variables // -enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge }; +enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; std::map dtlookup = { {"constant", DiskType::constant}, @@ -247,10 +248,12 @@ std::map dtlookup = {"mn", DiskType::mn}, {"exponential", DiskType::exponential}, {"doubleexpon", DiskType::doubleexpon}, - {"diskbulge", DiskType::diskbulge} + {"diskbulge", DiskType::diskbulge}, + {"python", DiskType::python} }; DiskType DTYPE; +std::string pyname; double ASCALE; double ASHIFT; double HSCALE; @@ -267,6 +270,7 @@ double HERNA = 0.10; double DiskDens(double R, double z, double phi) { double ans = 0.0; + static std::shared_ptr dfunc; switch (DTYPE) { @@ -322,6 +326,12 @@ double DiskDens(double R, double z, double phi) 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: { @@ -598,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; From 14e5da81791181f6cc08ec429aec43ea1d006b76 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 26 Dec 2024 11:58:18 -0500 Subject: [PATCH 12/14] Eliminate std::vector assignment warnings --- utils/ICs/EllipsoidForce.H | 4 +++- utils/ICs/EllipsoidForce.cc | 9 +++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/utils/ICs/EllipsoidForce.H b/utils/ICs/EllipsoidForce.H index 5ffcb68c4..6bc56a695 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 c3f58bbcc..3d918d650 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 From 3b84992d7fe34fe17ddfe41d1b8b70bd21aa0bde Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 26 Dec 2024 15:14:57 -0500 Subject: [PATCH 13/14] stdout log formating only --- utils/ICs/check_coefs.cc | 2 +- utils/ICs/cylcache.cc | 2 +- utils/ICs/initial.cc | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/utils/ICs/check_coefs.cc b/utils/ICs/check_coefs.cc index a9faf3af5..02d75a714 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 79959422f..a4da8071b 100644 --- a/utils/ICs/cylcache.cc +++ b/utils/ICs/cylcache.cc @@ -521,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) { diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index f0ad4b6d2..f0890ed2b 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -719,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) { From 94375ff95bb13aa72af436a2afe9c3f76281f3b2 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 28 Dec 2024 09:45:43 -0500 Subject: [PATCH 14/14] Fix typo in merge edit --- src/Cylinder.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cylinder.cc b/src/Cylinder.cc index a5a11ac56..af8103e0d 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -1,4 +1,4 @@ -not#include +#include #include #include #include