Skip to content

Commit 5af5071

Browse files
author
Martin D. Weinberg
committed
Merge branch 'getAccelArray' into biorthFactory
2 parents e4bc79a + 4b0a991 commit 5af5071

7 files changed

Lines changed: 203 additions & 89 deletions

File tree

‎expui/BiorthBasis.H‎

Lines changed: 82 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -84,6 +84,18 @@ namespace BasisClasses
8484
//! Subspace index
8585
virtual const std::string harmonic() = 0;
8686

87+
//@{
88+
//! Compute acceleration helper
89+
virtual void computeAccel(double x, double y, double z,
90+
std::vector<double>& vstore) = 0;
91+
92+
//! C++ local storage for getAccel
93+
std::vector<double> vstore {0.0, 0.0, 0.0};
94+
95+
//! C++ local storage for getAccelArray
96+
std::vector<std::vector<double>> vstore_array;
97+
//@}
98+
8799
public:
88100

89101
//! Constructor from YAML node
@@ -172,9 +184,46 @@ namespace BasisClasses
172184

173185
//! Evaluate acceleration in Cartesian coordinates in centered
174186
//! coordinate system
175-
virtual std::vector<double> getAccel(double x, double y, double z) = 0;
187+
virtual std::vector<double>& getAccel(double x, double y, double z)
188+
{
189+
computeAccel(x, y, z, vstore);
190+
return vstore;
191+
}
192+
193+
//! Evaluate acceleration in Cartesian coordinates in centered
194+
//! coordinate system for a collections of points
195+
virtual RowMatrixXd getAccelArray(Eigen::VectorXd& x,
196+
Eigen::VectorXd& y,
197+
Eigen::VectorXd& z)
198+
{
199+
// Sanity check
200+
if (x.size() != y.size() || x.size() != z.size())
201+
throw std::runtime_error("BiorthBasis::getAccelArray: "
202+
"x, y, z vectors must be of the same size");
203+
204+
int max_threads = 1;
205+
#if defined(_OPENMP)
206+
max_threads = omp_get_max_threads();
207+
#endif
208+
vstore_array.resize(max_threads);
209+
for (auto & v : vstore_array) v.resize(3);
210+
211+
RowMatrixXd ret(x.size(), 3);
212+
std::vector<double> v(3);
213+
int tid = 0; // Default thread ID if no OpenMP
214+
#pragma omp parallel for
215+
for (int i=0; i<x.size(); ++i) {
216+
#if defined(_OPENMP)
217+
tid = omp_get_thread_num();
218+
#endif
219+
computeAccel(x(i), y(i), z(i), vstore_array[tid]);
220+
for (int k=0; k<3; k++) ret(i, k) = vstore_array[tid][k];
221+
}
222+
223+
return ret;
224+
}
176225
};
177-
226+
178227
/**
179228
An abstract spherical basis to evaluate expansion coeffients and
180229
provide potential and density basis fields
@@ -263,6 +312,11 @@ namespace BasisClasses
263312

264313
//@}
265314

315+
//! Evaluate acceleration in Cartesian coordinates in centered
316+
//! coordinate system
317+
virtual void computeAccel(double x, double y, double z,
318+
std::vector<double>& vstore);
319+
266320
public:
267321

268322
//! Constructor from YAML node
@@ -317,9 +371,6 @@ namespace BasisClasses
317371
return ret;
318372
}
319373

320-
//! Evaluate acceleration in Cartesian coordinates in centered
321-
//! coordinate system
322-
virtual std::vector<double> getAccel(double x, double y, double z);
323374
};
324375

325376
/**
@@ -504,6 +555,11 @@ namespace BasisClasses
504555
//! Subspace index
505556
virtual const std::string harmonic() { return "m";}
506557

558+
//! Evaluate acceleration in Cartesian coordinates in centered
559+
//! coordinate system
560+
void computeAccel(double x, double y, double z,
561+
std::vector<double>& vstore);
562+
507563
public:
508564

509565
//! Constructor from YAML node
@@ -555,10 +611,6 @@ namespace BasisClasses
555611
return ret;
556612
}
557613

558-
//! Evaluate acceleration in Cartesian coordinates in centered
559-
//! coordinate system
560-
std::vector<double> getAccel(double x, double y, double z);
561-
562614
};
563615

564616
/**
@@ -659,6 +711,11 @@ namespace BasisClasses
659711
//! Subspace index
660712
virtual const std::string harmonic() { return "m";}
661713

714+
//! Evaluate acceleration in Cartesian coordinates in centered
715+
//! coordinate system
716+
void computeAccel(double x, double y, double z,
717+
std::vector<double>& vstore);
718+
662719
public:
663720

664721
//! Constructor from YAML node
@@ -710,10 +767,6 @@ namespace BasisClasses
710767
return ret;
711768
}
712769

713-
//! Evaluate acceleration in Cartesian coordinates in centered
714-
//! coordinate system
715-
std::vector<double> getAccel(double x, double y, double z);
716-
717770
};
718771

719772
/**
@@ -809,6 +862,11 @@ namespace BasisClasses
809862
//! Subspace index
810863
virtual const std::string harmonic() { return "m";}
811864

865+
//! Evaluate acceleration in Cartesian coordinates in centered
866+
//! coordinate system
867+
void computeAccel(double x, double y, double z,
868+
std::vector<double>& vstore);
869+
812870
public:
813871

814872
//! Constructor from YAML node
@@ -870,10 +928,6 @@ namespace BasisClasses
870928
return ret;
871929
}
872930

873-
//! Evaluate acceleration in Cartesian coordinates in centered
874-
//! coordinate system
875-
std::vector<double> getAccel(double x, double y, double z);
876-
877931
};
878932

879933
/**
@@ -967,6 +1021,11 @@ namespace BasisClasses
9671021
std::tuple<double, double, double, double, double>
9681022
eval(double x, double y, double z);
9691023

1024+
//! Evaluate acceleration in Cartesian coordinates in centered
1025+
//! coordinate system
1026+
void computeAccel(double x, double y, double z,
1027+
std::vector<double>& vstore);
1028+
9701029
public:
9711030

9721031
//! Constructor from YAML node
@@ -1020,10 +1079,6 @@ namespace BasisClasses
10201079
return true;
10211080
}
10221081

1023-
//! Evaluate acceleration in Cartesian coordinates in centered
1024-
//! coordinate system
1025-
std::vector<double> getAccel(double x, double y, double z);
1026-
10271082
};
10281083

10291084
/**
@@ -1101,6 +1156,12 @@ namespace BasisClasses
11011156
//! Readable index name
11021157
virtual const std::string harmonic() { return "n";}
11031158

1159+
//! Evaluate acceleration in Cartesian coordinates in centered
1160+
//! coordinate system
1161+
void computeAccel(double x, double y, double z,
1162+
std::vector<double>& vstore);
1163+
1164+
11041165
public:
11051166

11061167
//! Constructor from YAML node
@@ -1146,10 +1207,6 @@ namespace BasisClasses
11461207
return true;
11471208
}
11481209

1149-
//! Evaluate acceleration in Cartesian coordinates in centered
1150-
//! coordinate system
1151-
std::vector<double> getAccel(double x, double y, double z);
1152-
11531210
};
11541211

11551212

‎expui/BiorthBasis.cc‎

Lines changed: 19 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -692,8 +692,8 @@ namespace BasisClasses
692692
// Return force not potential gradient
693693
}
694694

695-
std::vector<double>
696-
Spherical::getAccel(double x, double y, double z)
695+
void Spherical::computeAccel(double x, double y, double z,
696+
std::vector<double>& acc)
697697
{
698698
// Get polar coordinates
699699
double R = sqrt(x*x + y*y);
@@ -795,7 +795,7 @@ namespace BasisClasses
795795

796796
// Return force not potential gradient
797797
//
798-
return {tpotx, tpoty, potz};
798+
acc = {tpotx, tpoty, potz};
799799
}
800800

801801

@@ -1599,7 +1599,8 @@ namespace BasisClasses
15991599
}
16001600

16011601
// Evaluate in cartesian coordinates
1602-
std::vector<double> Cylindrical::getAccel(double x, double y, double z)
1602+
void Cylindrical::computeAccel(double x, double y, double z,
1603+
std::vector<double>& acc)
16031604
{
16041605
double R = sqrt(x*x + y*y);
16051606
double phi = atan2(y, x);
@@ -1613,7 +1614,7 @@ namespace BasisClasses
16131614
double tpotx = tpotR*x/R - tpotp*y/R ;
16141615
double tpoty = tpotR*y/R + tpotp*x/R ;
16151616

1616-
return {tpotx, tpoty, tpotz};
1617+
acc = {tpotx, tpoty, tpotz};
16171618
}
16181619

16191620
// Evaluate in cylindrical coordinates
@@ -2180,7 +2181,8 @@ namespace BasisClasses
21802181
return {den0, den1, den0+den1, pot0, pot1, pot0+pot1, rpot, zpot, ppot};
21812182
}
21822183

2183-
std::vector<double> FlatDisk::getAccel(double x, double y, double z)
2184+
void FlatDisk::computeAccel(double x, double y, double z,
2185+
std::vector<double>& acc)
21842186
{
21852187
// Get thread id
21862188
int tid = omp_get_thread_num();
@@ -2204,7 +2206,7 @@ namespace BasisClasses
22042206
rpot = -totalMass*R/(r*r2 + 10.0*std::numeric_limits<double>::min());
22052207
zpot = -totalMass*z/(r*r2 + 10.0*std::numeric_limits<double>::min());
22062208

2207-
return {rpot, zpot, ppot};
2209+
acc = {rpot, zpot, ppot};
22082210
}
22092211

22102212
// Get the basis fields
@@ -2268,7 +2270,7 @@ namespace BasisClasses
22682270
double potx = rpot*x/R - ppot*y/R;
22692271
double poty = rpot*y/R + ppot*x/R;
22702272

2271-
return {potx, poty, zpot};
2273+
acc = {potx, poty, zpot};
22722274
}
22732275

22742276

@@ -2955,7 +2957,8 @@ namespace BasisClasses
29552957
}
29562958

29572959

2958-
std::vector<double> CBDisk::getAccel(double x, double y, double z)
2960+
void CBDisk::computeAccel(double x, double y, double z,
2961+
std::vector<double>& acc)
29592962
{
29602963
// Get thread id
29612964
int tid = omp_get_thread_num();
@@ -3020,7 +3023,7 @@ namespace BasisClasses
30203023
double potx = rpot*x/R - ppot*y/R;
30213024
double poty = rpot*y/R + ppot*x/R;
30223025

3023-
return {potx, poty, zpot};
3026+
acc = {potx, poty, zpot};
30243027
}
30253028

30263029

@@ -3427,8 +3430,8 @@ namespace BasisClasses
34273430
}
34283431

34293432

3430-
std::vector<double>
3431-
Slab::getAccel(double x, double y, double z)
3433+
void Slab::computeAccel(double x, double y, double z,
3434+
std::vector<double>& acc)
34323435
{
34333436
// Loop indices
34343437
//
@@ -3501,7 +3504,7 @@ namespace BasisClasses
35013504
}
35023505
}
35033506

3504-
return {accx.real(), accy.real(), accz.real()};
3507+
acc = {accx.real(), accy.real(), accz.real()};
35053508
}
35063509

35073510

@@ -3866,7 +3869,8 @@ namespace BasisClasses
38663869
return {0, den1, den1, 0, pot1, pot1, frcx, frcy, frcz};
38673870
}
38683871

3869-
std::vector<double> Cube::getAccel(double x, double y, double z)
3872+
void Cube::computeAccel(double x, double y, double z,
3873+
std::vector<double>& acc)
38703874
{
38713875
// Get thread id
38723876
int tid = omp_get_thread_num();
@@ -3877,7 +3881,7 @@ namespace BasisClasses
38773881
// Get the basis fields
38783882
auto frc = ortho->get_force(expcoef, pos);
38793883

3880-
return {-frc(0).real(), -frc(1).real(), -frc(2).real()};
3884+
acc = {-frc(0).real(), -frc(1).real(), -frc(2).real()};
38813885
}
38823886

38833887
std::vector<double> Cube::cyl_eval(double R, double z, double phi)

‎expui/expMSSA.H‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -85,7 +85,10 @@ namespace MSSA
8585
int ncomp;
8686

8787
//! Normalization options
88-
bool totVar, totPow;
88+
bool varFlag, powFlag;
89+
90+
//! Normalization values
91+
double totVar, totPow;
8992

9093
//! Toggle for detrending
9194
bool useMean;

0 commit comments

Comments
 (0)