Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Optimise/coordtransforms #29

Merged
merged 4 commits into from
Jun 13, 2024
Merged

Optimise/coordtransforms #29

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
bee50da
broken heatflux calculations, working trajectories
Waqar-ukaea Jun 11, 2024
2486fb3
now working. Seems to be much faster
Waqar-ukaea Jun 11, 2024
e202192
updated config file
Waqar-ukaea Jun 12, 2024
dafed40
various optimisations mostly around passing const references
Waqar-ukaea Jun 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ include_directories(externals/json/include)
#include(${CMAKE_SOURCE_DIR}/cmake/SetBuildType.cmake)

# Default to cpmpile with debugging symbols
set(CMAKE_BUILD_TYPE DEBUG)
set(CMAKE_BUILD_TYPE RELEASE)
# set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
# set (CMAKE_LINKER_FLAGS_DEBUG "${CMAKE_LINKER_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
# set(BUILD_TYPE "PROFILING")
Expand Down Expand Up @@ -92,7 +92,7 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/lib)
############################################################
# Set the build type environment variable
############################################################
if(BUILD_TYPE STREQUAL "TEST")
if(CMAKE_BUILD_TYPE STREQUAL "TEST")
set(BUILD_TYPE_COMPILE_FLAGS "-g;-O0;--coverage;")
set(TEST_LIBRARIES "gcov")
message(STATUS "BUILD_TYPE SET TO TEST")
Expand All @@ -102,18 +102,18 @@ elseif(CMAKE_BUILD_TYPE STREQUAL "DEBUG")
set(TEST_LIBRARIES "")
message(STATUS "BUILD_TYPE SET TO DEBUG")

elseif(BUILD_TYPE STREQUAL "RELEASE")
set(BUILD_TYPE_COMPILE_FLAGS "-O2")
elseif(CMAKE_BUILD_TYPE STREQUAL "RELEASE")
set(BUILD_TYPE_COMPILE_FLAGS "-g;-O2")
set(TEST_LIBRARIES "")
message(STATUS "BUILD_TYPE SET TO RELEASE")

elseif(BUILD_TYPE STREQUAL "PROFILING")
elseif(CMAKE_BUILD_TYPE STREQUAL "PROFILING")
set(BUILD_TYPE_COMPILE_FLAGS "-g;-O0;-pg")
set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg")
set (CMAKE_EXE_LINKER_FLAGS_DEBUG "${CMAKE_LINKER_FLAGS_DEBUG} -pg")
message(STATUS "BUILD_TYPE SET TO PROFILING")

elseif(BUILD_TYPE STREQUAL "ADDRESS-SANITIZE")
elseif(CMAKE_BUILD_TYPE STREQUAL "ADDRESS-SANITIZE")
set(BUILD_TYPE_COMPILE_FLAGS "-g;-O0;")
set(TEST_LIBRARIES "")
set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
Expand Down
79 changes: 40 additions & 39 deletions inres1/aegis_settings.json
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,40 +1,41 @@
{
"description": "inres1 test case to compare against SMARDDA",
"aegis_params":{
"DAGMC": "inres1_shad.h5m",
"step_size": 0.01,
"max_steps": 10000,
"launch_pos": "mc",

"monte_carlo_params":{
"number_of_particles_per_facet":5,
"write_particle_launch_positions":false
},

"target_surfs": [1,2,3,4,5,6],
"force_no_deposition": false,
"coordinate_system": "flux",
"execution_type": "dynamic",

"dynamic_batch_params":{
"batch_size":16,
"worker_profiling": false,
"debug":false
}
},
"equil_params":{
"eqdsk": "EQ3.eqdsk",
"power_sol": 7.5e+06,
"lambda_q": 0.012,
"r_outrbdry": 8.07385,
"cenopt": 2,
"psiref": -2.1628794,
"rmove": -0.006,
"draw_equil_rz": false,
"draw_equil_xyz": false,
"print_debug_info": false
},
"vtk_params":{
"draw_particle_tracks": false
}
}
"description": "inres1 test case to compare against SMARDDA",
"aegis_params":{
"DAGMC": "fine_inres1.h5m",
"step_size": 0.01,
"max_steps": 10000,
"launch_pos": "fixed",

"monte_carlo_params":{
"number_of_particles_per_facet":1,
"write_particle_launch_positions":false
},

"target_surfs": [1,2,3,4,5,6],
"force_no_deposition": false,
"coordinate_system": "flux",
"execution_type": "dynamic",

"dynamic_batch_params":{
"batch_size":16,
"worker_profiling": false,
"debug":false
}
},
"equil_params":{
"eqdsk": "EQ3.eqdsk",
"power_sol": 7.5e+06,
"lambda_q": 0.012,
"r_outrbdry": 8.07385,
"cenopt": 2,
"psiref": -2.1628794,
"rmove": -0.006,
"draw_equil_rz": false,
"draw_equil_xyz": false,
"print_debug_info": false
},
"vtk_params":{
"draw_particle_tracks": false
}
}

180 changes: 44 additions & 136 deletions src/aegis_lib/CoordTransform.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,177 +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;
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;
}
34 changes: 28 additions & 6 deletions src/aegis_lib/CoordTransform.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,22 +6,44 @@
#include <cmath>
#include <string>
#include "EquilData.h"
#include "alglib/interpolation.h"


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(double e0, double e1, double e2);
std::vector<double> cart_to_polar(const std::vector<double> &xyz);

// 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(std::vector<double> inputVector);
std::vector<double> polar_to_cart(double e0, double e1, 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);

};

#endif
Loading
Loading