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/spline interpolation #28

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Optimise/spline interpolation #28

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
1e858d9
splines for psi derivatives ends up slower than spline differentiation
Waqar-ukaea Jun 12, 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
80 changes: 41 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,42 @@
{
"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": "inres1_shad.h5m",
"step_size": 0.01,
"max_steps": 10000,
"launch_pos": "mc",

"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": "cart",
"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,
"spline_type":"cubic",
"draw_equil_rz": false,
"draw_equil_xyz": false,
"print_debug_info": false
},
"vtk_params":{
"draw_particle_tracks": false
}
}

24 changes: 23 additions & 1 deletion src/aegis_lib/CoordTransform.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,28 @@ CoordTransform::cart_to_polar(std::vector<double> inputVector)
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;
return outputVector;
}

std::vector<double>
CoordTransform::polar_to_cart(std::vector<double> inputVector)
{
Expand Down Expand Up @@ -182,4 +204,4 @@ CoordTransform::polar_to_flux(double e0, double e1, double e2,
outputVector[2] = phi;

return outputVector;
}
}
6 changes: 6 additions & 0 deletions src/aegis_lib/CoordTransform.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,16 @@
#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_xy(std::vector<double> inputVector);
std::vector<double> cart_to_polar(double e0, double e1, double e2);

std::vector<double> polar_to_cart(std::vector<double> inputVector);
Expand All @@ -21,7 +24,10 @@ std::vector<double> polar_to_cart(double e0, double e1, double e2);
// Polar toroidal to flux coordinates defined by psi spline
std::vector<double> polar_to_flux(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);

};

#endif
97 changes: 75 additions & 22 deletions src/aegis_lib/EquilData.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,9 @@ EquilData::read_params(const std::shared_ptr<JsonHandler> & configFile)
drawEquRZ = equilParams.get_optional<bool>("draw_equil_rz").value_or(drawEquRZ);
drawEquXYZ = equilParams.get_optional<bool>("draw_equil_xyz").value_or(drawEquXYZ);
debug = equilParams.get_optional<bool>("debug").value_or(debug);

auto splineTypeStr = equilParams.get_optional<std::string>("spline_type").value();
_splineType = (splineTypeStr == "linear") ? SplineType::LINEAR : SplineType::CUBIC;
}
else
{
Expand Down Expand Up @@ -438,7 +441,7 @@ EquilData::init_interp_splines()

// set_rsig();

// loop over Z creating spline knots
// loop over R creating spline knots
for (int i = 1; i < nw; i++)
{
rPts[i] = rPts[i - 1] + dr;
Expand Down Expand Up @@ -481,24 +484,62 @@ EquilData::init_interp_splines()
psi_1dgrid.setcontent(nw, psi1dpts.data());
f_grid.setcontent(nw, eqdsk.fpol.data());

// for (int i=0 ; i<nh; i++)
// {
// psi1d_out << psi_1dgrid[i] << std::endl;
// }

// Construct the spline interpolant for flux function psi(R,Z)
alglib::spline2dbuildbicubicv(r_grid, nw, z_grid, nh, psi_grid, 1, psiSpline);
switch (_splineType)
{
case SplineType::CUBIC:
alglib::spline2dbuildbicubicv(r_grid, nw, z_grid, nh, psi_grid, 1, psiSpline);
alglib::spline1dbuildcubic(psi_1dgrid, f_grid, fSpline);

log_string(LogLevel::WARNING,
"Splines initialised as cubic splines (cubic splines are default, use "
"spline_type:'linear' to switch to linear splines)");
break;

case SplineType::LINEAR:
alglib::spline2dbuildbilinearv(r_grid, nw, z_grid, nh, psi_grid, 1, psiSpline); // psi spline
alglib::spline1dbuildlinear(psi_1dgrid, f_grid, fSpline); // I(psi) aka f(psi)

// Construct the spline interpolant for toroidal component flux function
// I(psi) aka f(psi)
alglib::spline1dbuildcubic(psi_1dgrid, f_grid, fSpline);
log_string(LogLevel::WARNING, "Splines intialised as linear splines");
break;
}

int count2 = 0;
std::vector<double> dpsidrPts(nw * nh);
std::vector<double> dpsidzPts(nw * nh);
double psi, dpsidr, dpsidz, dpsi2drdz;
for (int j = 0; j < nh; j++) // flatten Psi(R,Z) array into Psi(index)
{
for (int i = 0; i < nw; i++)
{
alglib::spline2ddiff(psiSpline, rPts[i], zPts[j], psi, dpsidr, dpsidz, dpsi2drdz);
dpsidrPts[count2] = dpsidr;
dpsidzPts[count2] = dpsidz;
count2 += 1;
}
}

alglib::real_1d_array dpsidr_grid;
alglib::real_1d_array dpsidz_grid;
dpsidr_grid.setcontent(nw * nh, dpsidrPts.data());
dpsidz_grid.setcontent(nw * nh, dpsidzPts.data());

switch (_splineType)
{
case SplineType::CUBIC:
alglib::spline2dbuildbicubicv(r_grid, nw, z_grid, nh, dpsidr_grid, 1, dpsidrSpline);
alglib::spline2dbuildbicubicv(r_grid, nw, z_grid, nh, dpsidz_grid, 1, dpsidzSpline);
break;

case SplineType::LINEAR:
alglib::spline2dbuildbilinearv(r_grid, nw, z_grid, nh, dpsidr_grid, 1, dpsidrSpline);
alglib::spline2dbuildbilinearv(r_grid, nw, z_grid, nh, dpsidz_grid, 1, dpsidzSpline);
break;
}

// Alglib::spline2ddiff function can return a value of spline at (R,Z) as well
// as derivatives. I.e no need to have a separate spline for each derivative
// dPsidR and dPsidZ

// create 1d spline for f(psi) aka eqdsk.fpol
// alglib::spline1dbuildlinear(f_grid,)
}

// Write out psi(R,Z) data for gnuplotting
Expand Down Expand Up @@ -746,7 +787,7 @@ EquilData::r_extrema()
work1[j] = 2;
break;
}
zpsi = alglib::spline2dcalc(psiSpline, re, ze);
zpsi = get_psi(re, ze);
}
}

Expand Down Expand Up @@ -825,7 +866,7 @@ EquilData::b_field(std::vector<double> position, std::string startingFrom)
else
{
std::vector<double> polarPosition;
polarPosition = CoordTransform::cart_to_polar(position);
polarPosition = CoordTransform::cart_to_polar_xy(position);
zr = polarPosition[0];
zz = polarPosition[1];
}
Expand All @@ -836,7 +877,10 @@ EquilData::b_field(std::vector<double> position, std::string startingFrom)
}

// evaluate psi and psi derivs at given (R,Z) coords
alglib::spline2ddiff(psiSpline, zr, zz, zpsi, zdpdr, zdpdz, null);
zpsi = alglib::spline2dcalc(psiSpline, zr, zz);
zdpdr = alglib::spline2dcalc(dpsidrSpline, zr, zz);
zdpdz = alglib::spline2dcalc(dpsidzSpline, zr, zz);
// alglib::spline2ddiff(psiSpline, zr, zz, zpsi, zdpdr, zdpdz, null);

// evaluate I aka f at psi (I aka f is the flux function)
zf = alglib::spline1dcalc(fSpline, zpsi);
Expand Down Expand Up @@ -1074,7 +1118,7 @@ EquilData::boundary_rb()
// extremal Z reached
break;
}
zpsi = alglib::spline2dcalc(psiSpline, re, ze);
zpsi = get_psi(re, ze);

if (rsig > 0) // psiaxis < psiqbdry (psi increasing outwards)
{
Expand Down Expand Up @@ -1115,6 +1159,7 @@ EquilData::boundary_rb()

// check monotone
alglib::spline2ddiff(psiSpline, re, ze, zpsi, zdpdr, zdpdz, zddpdz);

zdpdsr = zdpdr * zcostheta + zdpdz * zsintheta;
if (idplset == 1)
{
Expand Down Expand Up @@ -1305,13 +1350,12 @@ EquilData::psi_limiter(std::vector<std::vector<double>> vertices)
for (auto & i : vertices)
{
std::vector<double> cartPos = i;
std::vector<double> polarPos = CoordTransform::cart_to_polar(i);
std::vector<double> polarPos = CoordTransform::cart_to_polar_xy(i);

zr = polarPos[0];
zz = polarPos[1];

zpsi = alglib::spline2dcalc(this->psiSpline, zr,
zz); // spline interpolation of psi(R,Z)
zpsi = get_psi(zr, zz);

ztheta = atan2(zz - zcen, zr - rcen);
if (ztheta < -M_PI_2)
Expand All @@ -1338,7 +1382,10 @@ EquilData::psi_limiter(std::vector<std::vector<double>> vertices)

double zdpdr;
double zdpdz;
alglib::spline2ddiff(psiSpline, zr, zz, zpsi, zdpdr, zdpdz, null);
// alglib::spline2ddiff(psiSpline, zr, zz, zpsi, zdpdr, zdpdz, null);
zpsi = alglib::spline2dcalc(psiSpline, zr, zz);
zdpdr = alglib::spline2dcalc(dpsidrSpline, zr, zz);
zdpdz = alglib::spline2dcalc(dpsidzSpline, zr, zz);

double zbpbdry, zbtotbdry;
zbpbdry = (1 / zr) * sqrt(std::max(0.0, (pow(zdpdr, 2) + pow(zdpdz, 2))));
Expand Down Expand Up @@ -1409,7 +1456,7 @@ EquilData::get_midplane_params()
bool
EquilData::check_if_in_bfield(std::vector<double> xyzPos)
{
std::vector<double> polarPos = CoordTransform::cart_to_polar(xyzPos);
std::vector<double> polarPos = CoordTransform::cart_to_polar_xy(xyzPos);
double r = polarPos[0];
double z = polarPos[1];
if (r < rmin || r > rmax || z < zmin || z > zmax)
Expand All @@ -1418,3 +1465,9 @@ EquilData::check_if_in_bfield(std::vector<double> xyzPos)
}
return true;
}

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