Skip to content

Commit a33c3a3

Browse files
The9CatThe9Cat
authored
Merge pull request #103 from CarrieFilion/diskbulge
2 parents 9db72e7 + 94375ff commit a33c3a3

13 files changed

+332
-72
lines changed

‎expui/BiorthBasis.H

+10-7
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include <unsupported/Eigen/CXX11/Tensor> // For 3d rectangular grids
99
#include <yaml-cpp/yaml.h>
1010

11+
#include <DiskDensityFunc.H>
1112
#include <ParticleReader.H>
1213
#include <Coefficients.H>
1314
#include <BiorthBess.H>
@@ -301,7 +302,7 @@ namespace BasisClasses
301302
{
302303
auto [ret, worst, lworst] = orthoCompute(orthoCheck(knots));
303304
// For the CTest log
304-
std::cout << "Spherical::orthoTest: worst=" << worst << std::endl;
305+
std::cout << "---- Spherical::orthoTest: worst=" << worst << std::endl;
305306
return ret;
306307
}
307308

@@ -536,7 +537,7 @@ namespace BasisClasses
536537
{
537538
auto [ret, worst, lworst] = orthoCompute(orthoCheck());
538539
// For the CTest log
539-
std::cout << "FlatDisk::orthoTest: worst=" << worst << std::endl;
540+
std::cout << "---- FlatDisk::orthoTest: worst=" << worst << std::endl;
540541
return ret;
541542
}
542543

@@ -741,18 +742,20 @@ namespace BasisClasses
741742

742743
//@{
743744
//! Basis construction parameters
744-
double aratio, hratio, dweight, rwidth, ashift, rfactor, rtrunc, ppow;
745+
double aratio, hratio, dweight, rwidth, ashift, rfactor, rtrunc, ppow, Mfac, HERNA;
745746
bool Ignore, deproject;
747+
std::string pyname;
746748

747749
//! DiskType support
748750
//
749-
enum DiskType { constant, gaussian, mn, exponential, doubleexpon };
751+
enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python };
750752
DiskType DTYPE;
751753
std::string mtype, dtype, dmodel;
752754

753755
static const std::map<std::string, DiskType> dtlookup;
754756
double DiskDens(double R, double z, double phi);
755757
double dcond(double R, double z, double phi, int M);
758+
std::shared_ptr<DiskDensityFunc> pyDens;
756759

757760
protected:
758761

@@ -837,7 +840,7 @@ namespace BasisClasses
837840
{
838841
auto [ret, worst, lworst] = orthoCompute(sl->orthoCheck());
839842
// For the CTest log
840-
std::cout << "Cylindrical::orthoTest: worst=" << worst << std::endl;
843+
std::cout << "---- Cylindrical::orthoTest: worst=" << worst << std::endl;
841844
return ret;
842845
}
843846
};
@@ -981,7 +984,7 @@ namespace BasisClasses
981984
}
982985
}
983986

984-
std::cout << "Slab::orthoTest: worst=" << worst << std::endl;
987+
std::cout << "---- Slab::orthoTest: worst=" << worst << std::endl;
985988
if (worst > __EXP__::orthoTol) return false;
986989
return true;
987990
}
@@ -1103,7 +1106,7 @@ namespace BasisClasses
11031106
}
11041107
}
11051108

1106-
std::cout << "Cube::orthoTest: worst=" << worst << std::endl;
1109+
std::cout << "---- Cube::orthoTest: worst=" << worst << std::endl;
11071110
if (worst > __EXP__::orthoTol) return false;
11081111
return true;
11091112
}

‎expui/BiorthBasis.cc

+44-7
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,7 @@ namespace BasisClasses
6060
"self_consistent",
6161
"cachename",
6262
"modelname",
63+
"pyname",
6364
"rnum"
6465
};
6566

@@ -860,7 +861,9 @@ namespace BasisClasses
860861
{"gaussian", DiskType::gaussian},
861862
{"mn", DiskType::mn},
862863
{"exponential", DiskType::exponential},
863-
{"doubleexpon", DiskType::doubleexpon}
864+
{"doubleexpon", DiskType::doubleexpon},
865+
{"diskbulge", DiskType::diskbulge},
866+
{"python", DiskType::python}
864867
};
865868

866869
const std::set<std::string>
@@ -911,6 +914,8 @@ namespace BasisClasses
911914
"aratio",
912915
"hratio",
913916
"dweight",
917+
"Mfac",
918+
"HERNA",
914919
"rwidth",
915920
"rfactor",
916921
"rtrunc",
@@ -921,7 +926,8 @@ namespace BasisClasses
921926
"self_consistent",
922927
"playback",
923928
"coefCompute",
924-
"coefMaster"
929+
"coefMaster",
930+
"pyname"
925931
};
926932

927933
Cylindrical::Cylindrical(const YAML::Node& CONF) :
@@ -980,6 +986,22 @@ namespace BasisClasses
980986
w2*exp(-R/a2)/(4.0*M_PI*a2*a2*h2*f2*f2) ;
981987
}
982988
break;
989+
990+
case DiskType::diskbulge:
991+
{
992+
double f = cosh(z/hcyl);
993+
double rr = pow(pow(R, 2) + pow(z,2), 0.5);
994+
double w1 = Mfac;
995+
double w2 = 1.0 - Mfac;
996+
double as = HERNA;
997+
998+
ans = w1*exp(-R/acyl)/(4.0*M_PI*acyl*acyl*hcyl*f*f) +
999+
w2*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0) ;
1000+
}
1001+
break;
1002+
case DiskType::python:
1003+
ans = (*pyDens)(R, z, phi);
1004+
break;
9831005
case DiskType::exponential:
9841006
default:
9851007
{
@@ -1068,6 +1090,12 @@ namespace BasisClasses
10681090
// Ratio of second disk relative to the first disk for disk basis construction with double-exponential
10691091
dweight = 1.0;
10701092

1093+
// mass fraction for disk for diskbulge
1094+
Mfac = 1.0;
1095+
1096+
// Hernquist scale a disk basis construction with diskbulge
1097+
HERNA = 0.10;
1098+
10711099
// Width for erf truncation for EOF conditioning density (ignored if zero)
10721100
rwidth = 0.0;
10731101

@@ -1124,16 +1152,19 @@ namespace BasisClasses
11241152
if (conf["dmodel" ]) dmodel = conf["dmodel" ].as<bool>();
11251153

11261154
if (conf["aratio" ]) aratio = conf["aratio" ].as<double>();
1127-
if (conf["hratio" ]) aratio = conf["hratio" ].as<double>();
1128-
if (conf["dweight" ]) aratio = conf["dweight" ].as<double>();
1129-
if (conf["rwidth" ]) aratio = conf["rwidth" ].as<double>();
1130-
if (conf["ashift" ]) aratio = conf["ashift" ].as<double>();
1155+
if (conf["hratio" ]) hratio = conf["hratio" ].as<double>();
1156+
if (conf["dweight" ]) dweight = conf["dweight" ].as<double>();
1157+
if (conf["Mfac" ]) Mfac = conf["Mfac" ].as<double>();
1158+
if (conf["HERNA" ]) HERNA = conf["HERNA" ].as<double>();
1159+
if (conf["rwidth" ]) rwidth = conf["rwidth" ].as<double>();
1160+
if (conf["ashift" ]) ashift = conf["ashift" ].as<double>();
11311161
if (conf["rfactor" ]) rfactor = conf["rfactor" ].as<double>();
11321162
if (conf["rtrunc" ]) rtrunc = conf["rtrunc" ].as<double>();
11331163
if (conf["pow" ]) ppow = conf["ppow" ].as<double>();
11341164
if (conf["mtype" ]) mtype = conf["mtype" ].as<std::string>();
11351165
if (conf["dtype" ]) dtype = conf["dtype" ].as<std::string>();
11361166
if (conf["vflag" ]) vflag = conf["vflag" ].as<int>();
1167+
if (conf["pyname" ]) pyname = conf["pyname" ].as<std::string>();
11371168

11381169
// Deprecation warning
11391170
if (conf["density" ]) {
@@ -1261,7 +1292,7 @@ namespace BasisClasses
12611292
DTYPE = dtlookup.at(dtype); // key is not in the map.
12621293

12631294
if (myid==0) // Report DiskType
1264-
std::cout << "DiskType is <" << dtype << ">" << std::endl;
1295+
std::cout << "---- DiskType is <" << dtype << ">" << std::endl;
12651296
}
12661297
catch (const std::out_of_range& err) {
12671298
if (myid==0) {
@@ -1273,6 +1304,12 @@ namespace BasisClasses
12731304
throw std::runtime_error("Cylindrical::initialize: invalid DiskType");
12741305
}
12751306

1307+
// Check for and initialize the Python density type
1308+
//
1309+
if (DTYPE == DiskType::python) {
1310+
pyDens = std::make_shared<DiskDensityFunc>(pyname);
1311+
}
1312+
12761313
// Use these user models to deproject for the EOF spherical basis
12771314
//
12781315
if (deproject) {

‎exputil/CMakeLists.txt

+3-2
Original file line numberDiff line numberDiff line change
@@ -35,13 +35,14 @@ set(QPDISTF_SRC QPDistF.cc qld.c)
3535
set(SLEDGE_SRC sledge.f)
3636
set(PARTICLE_SRC Particle.cc ParticleReader.cc header.cc)
3737
set(CUDA_SRC cudaParticle.cu cudaSLGridMP2.cu)
38+
set(PYWRAP_SRC DiskDensityFunc.cc)
3839

3940
set(exputil_SOURCES ${ODE_SRC} ${ROOT_SRC} ${QUAD_SRC}
4041
${RANDOM_SRC} ${UTIL_SRC} ${SPECFUNC_SRC}
4142
${PHASE_SRC} ${SYMP_SRC} ${INTERP_SRC} ${MASSMODEL_SRC}
4243
${ORBIT_SRC} ${BIORTH_SRC} ${POLY_SRC} ${GAUSS_SRC}
4344
${QPDISTF_SRC} ${BESSEL_SRC} ${OPTIMIZATION_SRC}
44-
${SLEDGE_SRC} ${PARTICLE_SRC} ${CUDA_SRC})
45+
${SLEDGE_SRC} ${PARTICLE_SRC} ${CUDA_SRC} ${PYWRAP_SRC})
4546

4647
set(common_INCLUDE_DIRS $<INSTALL_INTERFACE:include>
4748
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include/>
@@ -53,7 +54,7 @@ set(common_INCLUDE_DIRS $<INSTALL_INTERFACE:include>
5354

5455
set(common_LINKLIB ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp
5556
${VTK_LIBRARIES} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES}
56-
${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB})
57+
${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB} pybind11::embed)
5758

5859
if(ENABLE_CUDA)
5960
list(APPEND common_LINKLIB CUDA::toolkit CUDA::cudart)

‎exputil/DiskDensityFunc.cc

+31
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
#include <DiskDensityFunc.H>
2+
3+
namespace py = pybind11;
4+
5+
DiskDensityFunc::DiskDensityFunc(const std::string& modulename,
6+
const std::string& funcname)
7+
: funcname(funcname)
8+
{
9+
// Check if a Python interpreter exists
10+
if (Py_IsInitialized() == 0) {
11+
py::initialize_interpreter();
12+
// Mark the interpreter as started by this instance
13+
started = true;
14+
}
15+
16+
// Bind the disk_density function from Python
17+
disk_density =
18+
py::reinterpret_borrow<py::function>
19+
( py::module::import(modulename.c_str()).attr(funcname.c_str()) );
20+
}
21+
22+
DiskDensityFunc::~DiskDensityFunc()
23+
{
24+
// Only end the interpreter if it was started by this instance
25+
if (started) py::finalize_interpreter();
26+
}
27+
28+
double DiskDensityFunc::operator() (double R, double z, double phi)
29+
{
30+
return disk_density(R, z, phi).cast<double>();
31+
}

‎include/DiskDensityFunc.H

+61
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
#ifndef _DiskDensityFunc_H_
2+
#define _DiskDensityFunc_H_
3+
4+
#include <functional>
5+
6+
#include <pybind11/pybind11.h>
7+
#include <pybind11/embed.h>
8+
9+
/** Python disk-density function wrapper
10+
11+
The wrapper class will take a user supplied Python module that
12+
must supply the 'disk-density' function
13+
14+
An example disk-density Python module:
15+
--------------------------cut here--------------------------------
16+
import math
17+
18+
def disk_density(R, z, phi):
19+
"""Computes the disk density at a given radius, height, and
20+
azimuth. This would be the user's target density for basis
21+
creation.
22+
23+
"""
24+
a = 1.0 # Scale radius
25+
h = 0.2 # Scale height
26+
f = math.exp(-math.fabs(z)/h) # Prevent overflows
27+
sech = 2.0*f / (1.0 + f*f) #
28+
return math.exp(-R/a)*sech*sech/(4*math.pi*h*a*a)
29+
30+
--------------------------cut here--------------------------------
31+
32+
*/
33+
class __attribute__((visibility("default")))
34+
DiskDensityFunc
35+
{
36+
private:
37+
38+
//! The disk-density function
39+
pybind11::function disk_density;
40+
41+
//! Interpreter started?
42+
bool started = false;
43+
44+
//! Name of density function
45+
std::string funcname;
46+
47+
public:
48+
49+
//! Constructor
50+
DiskDensityFunc(const std::string& modulename,
51+
const std::string& funcname="disk_density");
52+
53+
//! Destructor
54+
~DiskDensityFunc();
55+
56+
//! Evaluate the density
57+
double operator() (double R, double z, double phi);
58+
59+
};
60+
61+
#endif

‎include/DiskModels.H

+55-2
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
#define _DISK_MODELS_H
33

44
#include <EmpCylSL.H>
5-
5+
#include <DiskDensityFunc.H>
66

77
//! A Ferrers Ellipsoid + Evacuated Exponential Disc (semi-realistic
88
//! bar+disk model)
@@ -223,5 +223,58 @@ public:
223223

224224
};
225225

226-
#endif
227226

227+
//! The usual exponential disk + a Hernquist bulge
228+
class diskbulge : public EmpCylSL::AxiDisk
229+
{
230+
private:
231+
double a, h, as, Mfac, d1, d2;
232+
233+
public:
234+
235+
diskbulge(double a, double h, double as, double Mfac, double M=1) :
236+
a(a), h(h), as(as), Mfac(Mfac), EmpCylSL::AxiDisk(M, "diskbulge")
237+
{
238+
params.push_back(a);
239+
params.push_back(h);
240+
params.push_back(as);
241+
params.push_back(Mfac);
242+
}
243+
244+
double operator()(double R, double z, double phi=0.)
245+
{
246+
double sh = 1.0/cosh(z/h);
247+
d1 = 0.25*M*Mfac/(M_PI*a*a*h) * exp(-R/a) * sh * sh;
248+
249+
double rr = pow(pow(R, 2) + pow(z,2), 0.5);
250+
d2 = M*(1-Mfac)*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0);
251+
252+
return d1 + d2;
253+
254+
}
255+
256+
};
257+
258+
259+
//! A user-defined Python model
260+
class PyModel : public EmpCylSL::AxiDisk
261+
{
262+
private:
263+
264+
std::shared_ptr<DiskDensityFunc> pyDens;
265+
266+
public:
267+
268+
PyModel(std::string& pyname)
269+
{
270+
pyDens = std::make_shared<DiskDensityFunc>(pyname);
271+
}
272+
273+
double operator()(double R, double z, double phi=0.)
274+
{
275+
return (*pyDens)(R, z, phi);
276+
}
277+
278+
};
279+
280+
#endif

0 commit comments

Comments
 (0)