Skip to content

Commit

Permalink
various optimisations mostly around passing const references
Browse files Browse the repository at this point in the history
  • Loading branch information
@Waqar-ukaea
Waqar-ukaea committed Jun 13, 2024
1 parent e202192 commit dafed40
Show file tree
Hide file tree
Showing 8 changed files with 96 additions and 216 deletions.
2 changes: 1 addition & 1 deletion inres1/aegis_settings.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
"DAGMC": "fine_inres1.h5m",
"step_size": 0.01,
"max_steps": 10000,
"launch_pos": "",
"launch_pos": "fixed",

"monte_carlo_params":{
"number_of_particles_per_facet":1,
Expand Down
204 changes: 45 additions & 159 deletions src/aegis_lib/CoordTransform.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,199 +9,85 @@
#include "SimpleLogger.h"

std::vector<double>
CoordTransform::cart_to_polar(std::vector<double> inputVector)
CoordTransform::cart_to_polar(const std::vector<double> & inputVector)
{
std::vector<double> outputVector(3);
double r; // polar r
double phi; // polar phi
double x; // cart x
double y; // cart y
double z; // cart z

x = inputVector[0];
y = inputVector[1];
z = inputVector[2];

r = sqrt(pow(x, 2) + pow(y, 2)); // calculate
phi = atan2(-y, x); // calculate phi

outputVector[0] = r;
outputVector[1] = z;
outputVector[2] = phi;
return outputVector;
}

// cart_to_polar without the calculation of phi
std::vector<double>
CoordTransform::cart_to_polar_xy(std::vector<double> inputVector)
{
std::vector<double> outputVector(3);
double r; // polar r
double phi; // polar phi
double x; // cart x
double y; // cart y
double z; // cart z

x = inputVector[0];
y = inputVector[1];
z = inputVector[2];

r = sqrt(pow(x, 2) + pow(y, 2)); // calculate

outputVector[0] = r;
outputVector[1] = z;
double x = inputVector[0];
double y = inputVector[1];
double z = inputVector[2];
double r = sqrt(pow(x, 2) + pow(y, 2)); // calculate
double phi = atan2(-y, x); // calculate phi
std::vector<double> outputVector = {r, z, phi};
return outputVector;
}

std::vector<double>
CoordTransform::polar_to_cart(std::vector<double> inputVector)
CoordTransform::polar_to_cart(const std::vector<double> & inputVector)
{
std::vector<double> outputVector(3);
double r; // polar r
double phi; // polar phi
double x; // cart x
double y; // cart y
double z; // cart z

r = inputVector[0];
z = inputVector[1];
phi = inputVector[2];

x = r * cos(phi); // calculate x
y = -r * sin(phi); // calculate y

outputVector[0] = x;
outputVector[1] = y;
outputVector[2] = z;

double r = inputVector[0];
double z = inputVector[1];
double phi = inputVector[2];
double x = r * cos(phi); // calculate x
double y = -r * sin(phi); // calculate y
std::vector<double> outputVector = {x, y, z};
return outputVector;
}

std::vector<double>
CoordTransform::cart_to_polar(double e0, double e1, double e2)
CoordTransform::cart_to_polar(const double & e0, const double & e1, const double & e2)
{
std::vector<double> outputVector(3);
double r; // polar r
double phi; // polar phi
double x; // cart x
double y; // cart y
double z; // cart z
x = e0;
y = e1;
z = e2;

r = sqrt(pow(x, 2) + pow(y, 2)); // calculate
phi = atan2(-y, x); // calculate phi

outputVector[0] = r;
outputVector[1] = z;
outputVector[2] = phi;
double x = e0;
double y = e1;
double z = e2;
double r = sqrt(pow(x, 2) + pow(y, 2)); // calculate
double phi = atan2(-y, x); // calculate phi
std::vector<double> outputVector = {r, z, phi};
return outputVector;
}

std::vector<double>
CoordTransform::polar_to_cart(double e0, double e1, double e2)
CoordTransform::polar_to_cart(const double & e0, const double & e1, const double & e2)
{
std::vector<double> outputVector(3);
double r; // polar r
double phi; // polar phi
double x; // cart x
double y; // cart y
double z; // cart z

r = e0;
z = e1;
phi = e2;

x = r * cos(phi); // calculate x

y = -r * sin(phi); // calculate y

outputVector[0] = x;
outputVector[1] = y;
outputVector[2] = z;

double r = e0;
double z = e1;
double phi = e2;
double x = r * cos(phi); // calculate x
double y = -r * sin(phi); // calculate y
std::vector<double> outputVector = {x, y, z};
return outputVector;
}

std::vector<double>
CoordTransform::polar_to_flux(std::vector<double> inputVector,
CoordTransform::polar_to_flux(const std::vector<double> & inputVector,
const std::shared_ptr<EquilData> & equilibrium)
{
std::vector<double> outputVector(3);
double r; // local polar r
double z; // local polar z
double phi; // local polar phi
double psi; // local psi
double theta; // local theta

r = inputVector[0];
z = inputVector[1];
phi = inputVector[2];

psi = alglib::spline2dcalc(equilibrium->psiSpline, r, z); // spline interpolation of psi(R,Z)

theta = atan2(z - equilibrium->zcen, r - equilibrium->rcen);
double r = inputVector[0];
double z = inputVector[1];
double phi = inputVector[2];
double psi =
alglib::spline2dcalc(equilibrium->psiSpline, r, z); // spline interpolation of psi(R,Z)
double theta = atan2(z - equilibrium->zcen, r - equilibrium->rcen);
if (theta < -M_PI_2)
{
theta = 2 * M_PI + theta;
}

outputVector[0] = -psi;
outputVector[1] = theta;
outputVector[2] = phi;

std::vector<double> outputVector = {-psi, theta, phi};
return outputVector;
}

// TODO
// std::vector<double>
// CoordTransform::flux_to_polar(std::vector<double> inputVector,
// const std::shared_ptr<EquilData> & equilibrium)
// {
// std::vector<double> outputVector(3);
// double r; // local polar r
// double z; // local polar z
// double phi; // local polar phi
// double psi; // local psi
// double theta; // local theta

// if (direction == "backwards") // backwards transform (flux -> polar coords)
// {
// psi = inputVector[0];
// theta = inputVector[1];
// phi = inputVector[2];
// }

// return outputVector;
// }

std::vector<double>
CoordTransform::polar_to_flux(double e0, double e1, double e2,
CoordTransform::polar_to_flux(const double & e0, const double & e1, const double & e2,
const std::shared_ptr<EquilData> & equilibrium)
{
std::vector<double> outputVector(3);
double r; // local polar r
double z; // local polar z
double phi; // local polar phi
double psi; // local psi
double theta; // local theta

r = e0;
z = e1;
phi = e2;

psi = alglib::spline2dcalc(equilibrium->psiSpline, r, z); // spline interpolation of psi(R,Z)

theta = atan2(z - equilibrium->zcen, r - equilibrium->rcen);
double r = e0;
double z = e1;
double phi = e2;
double psi =
alglib::spline2dcalc(equilibrium->psiSpline, r, z); // spline interpolation of psi(R,Z)
double theta = atan2(z - equilibrium->zcen, r - equilibrium->rcen);
if (theta < -M_PI_2)
{
theta = 2 * M_PI + theta;
}

outputVector[0] = -psi;
outputVector[1] = theta;
outputVector[2] = phi;

std::vector<double> outputVector = {-psi, theta, phi};
return outputVector;
}
}
30 changes: 23 additions & 7 deletions src/aegis_lib/CoordTransform.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,19 +14,35 @@ namespace CoordTransform


// Cartesian to polar toroidal (direction = "backwards" for polar->cartesian)
std::vector<double> cart_to_polar(std::vector<double> inputVector);
std::vector<double> cart_to_polar_xy(std::vector<double> inputVector);
std::vector<double> cart_to_polar(double e0, double e1, double e2);
std::vector<double> cart_to_polar(const std::vector<double> &xyz);

std::vector<double> polar_to_cart(std::vector<double> inputVector);
std::vector<double> polar_to_cart(double e0, double e1, double e2);
// polar (R,Z) coords from (X,Y,Z)
inline std::vector<double> cart_to_polar_xy(const std::vector<double> &xyz)
{
double r = sqrt(pow(xyz[0], 2) + pow(xyz[1], 2));
std::vector<double> rz = {r, xyz[2]};
return rz;
};

// get polar angle (PHI) from (X,Y)
inline double cart_to_polar_phi(const std::vector<double> &xyz)
{
double phi = atan2(-xyz[1], xyz[0]);
return phi;
};


std::vector<double> cart_to_polar(const double &e0, const double &e1, const double &e2);

std::vector<double> polar_to_cart(const std::vector<double> &inputVector);
std::vector<double> polar_to_cart(const double &e0, const double &e1, const double &e2);

// Polar toroidal to flux coordinates defined by psi spline
std::vector<double> polar_to_flux(std::vector<double> inputVector,
std::vector<double> polar_to_flux(const std::vector<double> &inputVector,
const std::shared_ptr<EquilData>& equilibrium);


std::vector<double> polar_to_flux(double e0, double e1, double e2, const std::shared_ptr<EquilData>& equilibrium);
std::vector<double> polar_to_flux(const double &e0, const double &e1, const double &e2, const std::shared_ptr<EquilData>& equilibrium);

};

Expand Down
48 changes: 13 additions & 35 deletions src/aegis_lib/EquilData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -803,7 +803,7 @@ EquilData::rz_splines()

// Caculate B field vector (in toroidal polars) at given position
std::vector<double>
EquilData::b_field(std::vector<double> position, std::string startingFrom)
EquilData::b_field(const std::vector<double> & position, std::string startingFrom)
{
std::vector<double> bVector(3); // B in toroidal polars
double zr; // local R from position vector supplied
Expand All @@ -824,8 +824,7 @@ EquilData::b_field(std::vector<double> position, std::string startingFrom)
// otherwise transform from cartesian -> polar before calculating B
else
{
std::vector<double> polarPosition;
polarPosition = CoordTransform::cart_to_polar_xy(position);
auto polarPosition = CoordTransform::cart_to_polar_xy(position);
zr = polarPosition[0];
zz = polarPosition[1];
}
Expand All @@ -851,35 +850,14 @@ EquilData::b_field(std::vector<double> position, std::string startingFrom)
}

std::vector<double>
EquilData::b_field_cart(std::vector<double> polarBVector, double phi, int normalise)
EquilData::b_field_cart(const std::vector<double> & polarBVector, const double & phi)
{
std::vector<double> cartBVector(3);
double zbx; // cartesian Bx component
double zby; // cartesian By component
double zbz; // Bz component
double zbr; // polar Br component
double zbt; // polar toroidal Bt component

zbr = polarBVector[0];
zbz = polarBVector[1];
zbt = polarBVector[2];

zbx = zbr * cos(phi) - zbt * sin(phi);
zby = -(zbr * sin(phi) + zbt * cos(phi));

cartBVector[0] = zbx;
cartBVector[1] = zby;
cartBVector[2] = zbz;

if (normalise == 1)
{
double norm;
norm = sqrt(pow(cartBVector[0], 2) + pow(cartBVector[1], 2) + pow(cartBVector[2], 2));
cartBVector[0] = cartBVector[0] / norm;
cartBVector[1] = cartBVector[1] / norm;
cartBVector[2] = cartBVector[2] / norm;
}

double zbr = polarBVector[0];
double zbz = polarBVector[1];
double zbt = polarBVector[2];
double zbx = zbr * cos(phi) - zbt * sin(phi);
double zby = -(zbr * sin(phi) + zbt * cos(phi));
std::vector<double> cartBVector = {zbx, zby, zbz};
return cartBVector;
}

Expand Down Expand Up @@ -958,7 +936,7 @@ EquilData::write_bfield(int phiSamples)
"Fixup eqdsk required");
}
polarB = b_field(polarPos, "polar"); // calculate B(R,Z,phi)
cartB = b_field_cart(polarB, polarPos[2], 0); // transform to B(x,y,z)
cartB = b_field_cart(polarB, polarPos[2]); // transform to B(x,y,z)
cartPos = CoordTransform::polar_to_cart(polarPos); // transform position to cartesian

// write out magnetic field data for plotting
Expand Down Expand Up @@ -1249,7 +1227,7 @@ EquilData::boundary_rb()
}

double
EquilData::omp_power_dep(double psi, double bn, std::string formula)
EquilData::omp_power_dep(const double & psi, const double & bn, std::string formula)
{
double heatFlux;
double exponential, fpfac;
Expand Down Expand Up @@ -1406,7 +1384,7 @@ EquilData::get_midplane_params()

// check if in b_field from polar coords (R,Z)
bool
EquilData::check_if_in_bfield(std::vector<double> xyzPos)
EquilData::check_if_in_bfield(const std::vector<double> & xyzPos)
{
std::vector<double> polarPos = CoordTransform::cart_to_polar_xy(xyzPos);
double r = polarPos[0];
Expand All @@ -1419,7 +1397,7 @@ EquilData::check_if_in_bfield(std::vector<double> xyzPos)
}

double
EquilData::get_psi(double r, double z)
EquilData::get_psi(const double & r, const double & z)
{
return -(alglib::spline2dcalc(psiSpline, r, z)); // spline returns sign flipped psi
}
Loading

0 comments on commit dafed40

Please sign in to comment.