From bee50daef6d07771a97256401c55af651a7f20b6 Mon Sep 17 00:00:00 2001 From: Waqar Butt Date: Tue, 11 Jun 2024 17:15:52 +0100 Subject: [PATCH 1/4] broken heatflux calculations, working trajectories --- CMakeLists.txt | 8 ++++---- inres1/aegis_settings.json | 4 ++-- src/aegis_lib/CoordTransform.cpp | 24 +++++++++++++++++++++++- src/aegis_lib/CoordTransform.h | 6 ++++++ src/aegis_lib/EquilData.cpp | 21 ++++++++++++++------- src/aegis_lib/EquilData.h | 3 +++ src/aegis_lib/Particle.cpp | 10 +++------- src/aegis_lib/ParticleSimulation.cpp | 19 ++++++++----------- src/aegis_lib/Source.cpp | 3 +-- 9 files changed, 64 insertions(+), 34 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d5af8f9..4560fed 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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") diff --git a/inres1/aegis_settings.json b/inres1/aegis_settings.json index 52cd7a5..4ac066b 100644 --- a/inres1/aegis_settings.json +++ b/inres1/aegis_settings.json @@ -7,7 +7,7 @@ "launch_pos": "mc", "monte_carlo_params":{ - "number_of_particles_per_facet":5, + "number_of_particles_per_facet":1, "write_particle_launch_positions":false }, @@ -35,6 +35,6 @@ "print_debug_info": false }, "vtk_params":{ - "draw_particle_tracks": false + "draw_particle_tracks": true } } diff --git a/src/aegis_lib/CoordTransform.cpp b/src/aegis_lib/CoordTransform.cpp index effa031..fec2cc1 100644 --- a/src/aegis_lib/CoordTransform.cpp +++ b/src/aegis_lib/CoordTransform.cpp @@ -31,6 +31,28 @@ CoordTransform::cart_to_polar(std::vector inputVector) return outputVector; } +// cart_to_polar without the calculation of phi +std::vector +CoordTransform::cart_to_polar_xy(std::vector inputVector) +{ + std::vector 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 CoordTransform::polar_to_cart(std::vector inputVector) { @@ -182,4 +204,4 @@ CoordTransform::polar_to_flux(double e0, double e1, double e2, outputVector[2] = phi; return outputVector; -} +} \ No newline at end of file diff --git a/src/aegis_lib/CoordTransform.h b/src/aegis_lib/CoordTransform.h index c43f8b1..d358350 100644 --- a/src/aegis_lib/CoordTransform.h +++ b/src/aegis_lib/CoordTransform.h @@ -6,6 +6,8 @@ #include #include #include "EquilData.h" +#include "alglib/interpolation.h" + namespace CoordTransform { @@ -13,6 +15,7 @@ namespace CoordTransform // Cartesian to polar toroidal (direction = "backwards" for polar->cartesian) std::vector cart_to_polar(std::vector inputVector); +std::vector cart_to_polar_xy(std::vector inputVector); std::vector cart_to_polar(double e0, double e1, double e2); std::vector polar_to_cart(std::vector inputVector); @@ -21,7 +24,10 @@ std::vector polar_to_cart(double e0, double e1, double e2); // Polar toroidal to flux coordinates defined by psi spline std::vector polar_to_flux(std::vector inputVector, const std::shared_ptr& equilibrium); + + std::vector polar_to_flux(double e0, double e1, double e2, const std::shared_ptr& equilibrium); + }; #endif \ No newline at end of file diff --git a/src/aegis_lib/EquilData.cpp b/src/aegis_lib/EquilData.cpp index 7d4e18e..4af4c40 100644 --- a/src/aegis_lib/EquilData.cpp +++ b/src/aegis_lib/EquilData.cpp @@ -746,7 +746,7 @@ EquilData::r_extrema() work1[j] = 2; break; } - zpsi = alglib::spline2dcalc(psiSpline, re, ze); + zpsi = get_psi(re, ze); } } @@ -825,7 +825,7 @@ EquilData::b_field(std::vector position, std::string startingFrom) else { std::vector polarPosition; - polarPosition = CoordTransform::cart_to_polar(position); + polarPosition = CoordTransform::cart_to_polar_xy(position); zr = polarPosition[0]; zz = polarPosition[1]; } @@ -1074,7 +1074,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) { @@ -1305,13 +1305,12 @@ EquilData::psi_limiter(std::vector> vertices) for (auto & i : vertices) { std::vector cartPos = i; - std::vector polarPos = CoordTransform::cart_to_polar(i); + std::vector 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) @@ -1409,7 +1408,7 @@ EquilData::get_midplane_params() bool EquilData::check_if_in_bfield(std::vector xyzPos) { - std::vector polarPos = CoordTransform::cart_to_polar(xyzPos); + std::vector polarPos = CoordTransform::cart_to_polar_xy(xyzPos); double r = polarPos[0]; double z = polarPos[1]; if (r < rmin || r > rmax || z < zmin || z > zmax) @@ -1418,3 +1417,11 @@ EquilData::check_if_in_bfield(std::vector xyzPos) } return true; } + +double +EquilData::get_psi(double r, double z) +{ + auto psi = alglib::spline2dcalc(psiSpline, r, z); + std::cout << psi << std::endl; + return psi; +} \ No newline at end of file diff --git a/src/aegis_lib/EquilData.h b/src/aegis_lib/EquilData.h index 7e6565e..acbde9d 100644 --- a/src/aegis_lib/EquilData.h +++ b/src/aegis_lib/EquilData.h @@ -91,6 +91,9 @@ class EquilData : public AegisBase // Return eqdsk struct eqdskData get_eqdsk_struct(); + // return psi value at given (R,Z) coords + double get_psi(double r, double z); + // Read eqdsk file void read_eqdsk(std::string filename); diff --git a/src/aegis_lib/Particle.cpp b/src/aegis_lib/Particle.cpp index d02de88..9d12f49 100644 --- a/src/aegis_lib/Particle.cpp +++ b/src/aegis_lib/Particle.cpp @@ -73,21 +73,17 @@ ParticleBase::get_xyz_pos() double ParticleBase::get_psi(const std::shared_ptr & equilibrium) { - std::vector fluxPos; std::vector polarPos; double psi; switch (coordSystem) { case coordinateSystem::CARTESIAN: - polarPos = CoordTransform::cart_to_polar(pos); - fluxPos = CoordTransform::polar_to_flux(polarPos, equilibrium); - psi = fluxPos[0]; + polarPos = CoordTransform::cart_to_polar_xy(pos); + psi = equilibrium->get_psi(polarPos[0], polarPos[1]); break; case coordinateSystem::POLAR: - fluxPos = CoordTransform::polar_to_flux(polarPos, equilibrium); - psi = fluxPos[0]; + psi = equilibrium->get_psi(pos[0], pos[1]); } - return psi; } diff --git a/src/aegis_lib/ParticleSimulation.cpp b/src/aegis_lib/ParticleSimulation.cpp index 4481bf7..eb24300 100644 --- a/src/aegis_lib/ParticleSimulation.cpp +++ b/src/aegis_lib/ParticleSimulation.cpp @@ -64,7 +64,6 @@ ParticleSimulation::Execute_dynamic_mpi() MPI_Status mpiStatus; - const int noMoreWork = -1; std::vector particleHeatfluxes(num_particles_launched()); @@ -442,8 +441,7 @@ ParticleSimulation::init_geometry() } dagmcMeshReadTime = MPI_Wtime() - dagmcMeshReadStart; - - // setup B Field data + // setup B Field data double prepSurfacesStart = MPI_Wtime(); std::vector vertexCoordinates; @@ -483,7 +481,7 @@ ParticleSimulation::init_geometry() targetFacets = select_target_surface(); integrator = std::make_unique(targetFacets); prepSurfacesTime = MPI_Wtime() - prepSurfacesStart; - + double setupArrayOfParticlesTimeStart = MPI_Wtime(); setup_sources(); setupArrayOfParticlesTime = MPI_Wtime() - setupArrayOfParticlesTimeStart; @@ -819,7 +817,6 @@ ParticleSimulation::Execute_serial() vtkInterface->init(); moab::ErrorCode rval; - int start = 0; int end = num_particles_launched(); @@ -870,7 +867,6 @@ ParticleSimulation::Execute_serial() vtkInterface->write_multiBlockData("particle_tracks.vtm"); print_particle_stats(particleStats); - } void @@ -882,7 +878,6 @@ ParticleSimulation::Execute_mpi() vtkInterface->init(); MPI_Status mpiStatus; - int totalFacets = target_num_facets(); int nFacetsPerProc = totalFacets / nprocs; int remainder = totalFacets % nprocs; @@ -955,7 +950,6 @@ ParticleSimulation::Execute_mpi() attach_mesh_attribute("Heatflux", targetFacets, rootQvalues); write_out_mesh(meshWriteOptions::BOTH, targetFacets); mainParticleLoopTime = MPI_Wtime() - mainParticleLoopStart; - } mpi_particle_stats(); @@ -1171,9 +1165,12 @@ ParticleSimulation::get_profiling_times() { std::map profilingTimes; profilingTimes.insert(std::make_pair("DAGMC Mesh read runtime = ", dagmcMeshReadTime)); - profilingTimes.insert(std::make_pair("Preparing surfaces for particles runtime = ", prepSurfacesTime)); - profilingTimes.insert(std::make_pair("Pool of particles generation runtime = ", setupArrayOfParticlesTime)); - profilingTimes.insert(std::make_pair("Main particle tracking loop runtime = ", mainParticleLoopTime)); + profilingTimes.insert( + std::make_pair("Preparing surfaces for particles runtime = ", prepSurfacesTime)); + profilingTimes.insert( + std::make_pair("Pool of particles generation runtime = ", setupArrayOfParticlesTime)); + profilingTimes.insert( + std::make_pair("Main particle tracking loop runtime = ", mainParticleLoopTime)); profilingTimes.insert(std::make_pair("Mesh Write out runtime = ", aegisMeshWriteTime)); return profilingTimes; diff --git a/src/aegis_lib/Source.cpp b/src/aegis_lib/Source.cpp index 944ca8e..e0a10b7 100644 --- a/src/aegis_lib/Source.cpp +++ b/src/aegis_lib/Source.cpp @@ -101,7 +101,6 @@ Sources::set_heatflux_params(const std::shared_ptr & equilibrium, std::vector bFieldXYZ; polarPos = CoordTransform::cart_to_polar(launchPos); - fluxPos = CoordTransform::polar_to_flux(polarPos, equilibrium); bField = equilibrium->b_field(polarPos, "forwards"); @@ -114,7 +113,7 @@ Sources::set_heatflux_params(const std::shared_ptr & equilibrium, Bn = product; // store B.n - psi = fluxPos[0]; // store psi at particle start + psi = equilibrium->get_psi(polarPos[0], polarPos[1]); // store psi at particle start psid = psi + equilibrium->psibdry; _totalHeatflux = equilibrium->omp_power_dep(psid, Bn, "exp"); // store Q at particle start } From 2486fb3d13950ddb2789eedec48d45dd7f746e51 Mon Sep 17 00:00:00 2001 From: Waqar Butt Date: Tue, 11 Jun 2024 17:38:38 +0100 Subject: [PATCH 2/4] now working. Seems to be much faster --- CMakeLists.txt | 4 ++-- inres1/aegis_settings.json | 6 +++--- src/aegis_lib/EquilData.cpp | 4 +--- 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4560fed..6a38aab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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") @@ -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") diff --git a/inres1/aegis_settings.json b/inres1/aegis_settings.json index 4ac066b..0dc319f 100644 --- a/inres1/aegis_settings.json +++ b/inres1/aegis_settings.json @@ -4,10 +4,10 @@ "DAGMC": "inres1_shad.h5m", "step_size": 0.01, "max_steps": 10000, - "launch_pos": "mc", + "launch_pos": "fixed", "monte_carlo_params":{ - "number_of_particles_per_facet":1, + "number_of_particles_per_facet":5, "write_particle_launch_positions":false }, @@ -35,6 +35,6 @@ "print_debug_info": false }, "vtk_params":{ - "draw_particle_tracks": true + "draw_particle_tracks": false } } diff --git a/src/aegis_lib/EquilData.cpp b/src/aegis_lib/EquilData.cpp index 4af4c40..2b0875b 100644 --- a/src/aegis_lib/EquilData.cpp +++ b/src/aegis_lib/EquilData.cpp @@ -1421,7 +1421,5 @@ EquilData::check_if_in_bfield(std::vector xyzPos) double EquilData::get_psi(double r, double z) { - auto psi = alglib::spline2dcalc(psiSpline, r, z); - std::cout << psi << std::endl; - return psi; + return -(alglib::spline2dcalc(psiSpline, r, z)); // spline returns sign flipped psi } \ No newline at end of file From e202192ef31dd48c62482271db518db3f7212e7f Mon Sep 17 00:00:00 2001 From: Waqar Butt Date: Wed, 12 Jun 2024 10:53:21 +0100 Subject: [PATCH 3/4] updated config file --- inres1/aegis_settings.json | 79 +++++++++++++++++++------------------- 1 file changed, 40 insertions(+), 39 deletions(-) diff --git a/inres1/aegis_settings.json b/inres1/aegis_settings.json index 0dc319f..bc7a0e4 100644 --- a/inres1/aegis_settings.json +++ b/inres1/aegis_settings.json @@ -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": "fixed", - - "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": "", + + "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 + } + } + From dafed400c6342fb0e94bb2b0de4ca01123536e7f Mon Sep 17 00:00:00 2001 From: Waqar Butt Date: Thu, 13 Jun 2024 16:12:10 +0100 Subject: [PATCH 4/4] various optimisations mostly around passing const references --- inres1/aegis_settings.json | 2 +- src/aegis_lib/CoordTransform.cpp | 204 +++++++------------------------ src/aegis_lib/CoordTransform.h | 30 +++-- src/aegis_lib/EquilData.cpp | 48 ++------ src/aegis_lib/EquilData.h | 10 +- src/aegis_lib/Particle.cpp | 12 +- src/aegis_lib/Particle.h | 4 +- src/aegis_lib/Source.cpp | 2 +- 8 files changed, 96 insertions(+), 216 deletions(-) diff --git a/inres1/aegis_settings.json b/inres1/aegis_settings.json index bc7a0e4..34b08d6 100644 --- a/inres1/aegis_settings.json +++ b/inres1/aegis_settings.json @@ -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, diff --git a/src/aegis_lib/CoordTransform.cpp b/src/aegis_lib/CoordTransform.cpp index fec2cc1..c5690a2 100644 --- a/src/aegis_lib/CoordTransform.cpp +++ b/src/aegis_lib/CoordTransform.cpp @@ -9,199 +9,85 @@ #include "SimpleLogger.h" std::vector -CoordTransform::cart_to_polar(std::vector inputVector) +CoordTransform::cart_to_polar(const std::vector & inputVector) { - std::vector 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 -CoordTransform::cart_to_polar_xy(std::vector inputVector) -{ - std::vector 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 outputVector = {r, z, phi}; return outputVector; } std::vector -CoordTransform::polar_to_cart(std::vector inputVector) +CoordTransform::polar_to_cart(const std::vector & inputVector) { - std::vector 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 outputVector = {x, y, z}; return outputVector; } std::vector -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 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 outputVector = {r, z, phi}; return outputVector; } std::vector -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 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 outputVector = {x, y, z}; return outputVector; } std::vector -CoordTransform::polar_to_flux(std::vector inputVector, +CoordTransform::polar_to_flux(const std::vector & inputVector, const std::shared_ptr & equilibrium) { - std::vector 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 outputVector = {-psi, theta, phi}; return outputVector; } -// TODO -// std::vector -// CoordTransform::flux_to_polar(std::vector inputVector, -// const std::shared_ptr & equilibrium) -// { -// std::vector 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 -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 & equilibrium) { - std::vector 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 outputVector = {-psi, theta, phi}; return outputVector; -} \ No newline at end of file +} diff --git a/src/aegis_lib/CoordTransform.h b/src/aegis_lib/CoordTransform.h index d358350..18746a7 100644 --- a/src/aegis_lib/CoordTransform.h +++ b/src/aegis_lib/CoordTransform.h @@ -14,19 +14,35 @@ namespace CoordTransform // Cartesian to polar toroidal (direction = "backwards" for polar->cartesian) -std::vector cart_to_polar(std::vector inputVector); -std::vector cart_to_polar_xy(std::vector inputVector); -std::vector cart_to_polar(double e0, double e1, double e2); +std::vector cart_to_polar(const std::vector &xyz); -std::vector polar_to_cart(std::vector inputVector); -std::vector polar_to_cart(double e0, double e1, double e2); +// polar (R,Z) coords from (X,Y,Z) +inline std::vector cart_to_polar_xy(const std::vector &xyz) +{ + double r = sqrt(pow(xyz[0], 2) + pow(xyz[1], 2)); + std::vector rz = {r, xyz[2]}; + return rz; +}; + +// get polar angle (PHI) from (X,Y) +inline double cart_to_polar_phi(const std::vector &xyz) +{ + double phi = atan2(-xyz[1], xyz[0]); + return phi; +}; + + +std::vector cart_to_polar(const double &e0, const double &e1, const double &e2); + +std::vector polar_to_cart(const std::vector &inputVector); +std::vector polar_to_cart(const double &e0, const double &e1, const double &e2); // Polar toroidal to flux coordinates defined by psi spline -std::vector polar_to_flux(std::vector inputVector, +std::vector polar_to_flux(const std::vector &inputVector, const std::shared_ptr& equilibrium); -std::vector polar_to_flux(double e0, double e1, double e2, const std::shared_ptr& equilibrium); +std::vector polar_to_flux(const double &e0, const double &e1, const double &e2, const std::shared_ptr& equilibrium); }; diff --git a/src/aegis_lib/EquilData.cpp b/src/aegis_lib/EquilData.cpp index 2b0875b..5291b6b 100644 --- a/src/aegis_lib/EquilData.cpp +++ b/src/aegis_lib/EquilData.cpp @@ -803,7 +803,7 @@ EquilData::rz_splines() // Caculate B field vector (in toroidal polars) at given position std::vector -EquilData::b_field(std::vector position, std::string startingFrom) +EquilData::b_field(const std::vector & position, std::string startingFrom) { std::vector bVector(3); // B in toroidal polars double zr; // local R from position vector supplied @@ -824,8 +824,7 @@ EquilData::b_field(std::vector position, std::string startingFrom) // otherwise transform from cartesian -> polar before calculating B else { - std::vector polarPosition; - polarPosition = CoordTransform::cart_to_polar_xy(position); + auto polarPosition = CoordTransform::cart_to_polar_xy(position); zr = polarPosition[0]; zz = polarPosition[1]; } @@ -851,35 +850,14 @@ EquilData::b_field(std::vector position, std::string startingFrom) } std::vector -EquilData::b_field_cart(std::vector polarBVector, double phi, int normalise) +EquilData::b_field_cart(const std::vector & polarBVector, const double & phi) { - std::vector 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 cartBVector = {zbx, zby, zbz}; return cartBVector; } @@ -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 @@ -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; @@ -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 xyzPos) +EquilData::check_if_in_bfield(const std::vector & xyzPos) { std::vector polarPos = CoordTransform::cart_to_polar_xy(xyzPos); double r = polarPos[0]; @@ -1419,7 +1397,7 @@ EquilData::check_if_in_bfield(std::vector 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 } \ No newline at end of file diff --git a/src/aegis_lib/EquilData.h b/src/aegis_lib/EquilData.h index acbde9d..bf92129 100644 --- a/src/aegis_lib/EquilData.h +++ b/src/aegis_lib/EquilData.h @@ -92,7 +92,7 @@ class EquilData : public AegisBase eqdskData get_eqdsk_struct(); // return psi value at given (R,Z) coords - double get_psi(double r, double z); + double get_psi(const double &r, const double &z); // Read eqdsk file void read_eqdsk(std::string filename); @@ -117,10 +117,10 @@ class EquilData : public AegisBase // Caculate B field vector (in toroidal polars) at given position // set string to "polar" if position vector already in polars - std::vector b_field(std::vector position, std::string startingFrom); + std::vector b_field(const std::vector &position, std::string startingFrom); // Convert B Field vectors to cartesian given polar form and value of angle - std::vector b_field_cart(std::vector polarBVector, double phi, int normalise); + std::vector b_field_cart(const std::vector &polarBVector,const double &phi); // Write out positions and associated BField vectors in cartesian and/or polar toroidal void write_bfield(int phiSamples = 12); @@ -130,7 +130,7 @@ class EquilData : public AegisBase // Determine Rm and Bpm (R and Bpol at omp) void boundary_rb(); - double omp_power_dep(double psi, double bn, std::string formula); + double omp_power_dep(const double &psi, const double &bn, std::string formula); void psi_limiter(std::vector> vertices); @@ -138,7 +138,7 @@ class EquilData : public AegisBase void psiref_override(); - bool check_if_in_bfield(std::vector xyzPos); + bool check_if_in_bfield(const std::vector &xyzPos); std::array get_midplane_params(); // return rInnerMidplane, rOuterMidplane, zMidplane diff --git a/src/aegis_lib/Particle.cpp b/src/aegis_lib/Particle.cpp index 9d12f49..fa52239 100644 --- a/src/aegis_lib/Particle.cpp +++ b/src/aegis_lib/Particle.cpp @@ -91,8 +91,6 @@ ParticleBase::get_psi(const std::shared_ptr & equilibrium) void ParticleBase::set_dir(const std::shared_ptr & equilibrium) { - double norm = 0; // normalisation constant for magnetic field - std::vector normB(3); // normalised magnetic field std::vector polarPos; // polar position std::vector magn(3); @@ -101,7 +99,7 @@ ParticleBase::set_dir(const std::shared_ptr & equilibrium) case coordinateSystem::CARTESIAN: polarPos = CoordTransform::cart_to_polar(pos); magn = equilibrium->b_field(pos, "cart"); // magnetic field - magn = equilibrium->b_field_cart(magn, polarPos[2], 0); + magn = equilibrium->b_field_cart(magn, polarPos[2]); break; case coordinateSystem::POLAR: @@ -112,8 +110,10 @@ ParticleBase::set_dir(const std::shared_ptr & equilibrium) break; } - norm = sqrt(pow(magn[0], 2) + pow(magn[1], 2) + pow(magn[2], 2)); + // normalise particle direction vector + double norm = sqrt(pow(magn[0], 2) + pow(magn[1], 2) + pow(magn[2], 2)); Bfield = magn; + std::vector normB(3); normB[0] = magn[0] / norm; normB[1] = magn[1] / norm; normB[2] = magn[2] / norm; @@ -146,7 +146,7 @@ ParticleBase::align_dir_to_surf(double Bn) // update position vector of particle with set distance travelled void -ParticleBase::update_vectors(double distanceTravelled) +ParticleBase::update_vectors(const double & distanceTravelled) { // This will update the position std::vector newPt(3); @@ -168,7 +168,7 @@ ParticleBase::update_vectors(double distanceTravelled) // update position vector of particle with set distance travelled and update // the particle direction at new position void -ParticleBase::update_vectors(double distanceTravelled, +ParticleBase::update_vectors(const double & distanceTravelled, const std::shared_ptr & equilibrium) { // This will update the position and direction vector of the particle diff --git a/src/aegis_lib/Particle.h b/src/aegis_lib/Particle.h index 24b5345..49c405a 100644 --- a/src/aegis_lib/Particle.h +++ b/src/aegis_lib/Particle.h @@ -73,8 +73,8 @@ class ParticleBase : public AegisBase // return STL vector of unit direction std::vector get_dir(); void align_dir_to_surf(double Bn); // align particle dir to surface normal - void update_vectors(double distanceTravelled); // update position - void update_vectors(double distanceTravelled, const std::shared_ptr &EquData); // overload to update dir as well + void update_vectors(const double &distanceTravelled); // update position + void update_vectors(const double &distanceTravelled, const std::shared_ptr &EquData); // overload to update dir as well void check_if_midplane_crossed(const std::array &midplaneParameters); // check if particle has reached inner or outer midplane and set value of atMidplane void set_intersection_threshold(double distanceThreshold); bool check_if_threshold_crossed(); diff --git a/src/aegis_lib/Source.cpp b/src/aegis_lib/Source.cpp index e0a10b7..24d9151 100644 --- a/src/aegis_lib/Source.cpp +++ b/src/aegis_lib/Source.cpp @@ -104,7 +104,7 @@ Sources::set_heatflux_params(const std::shared_ptr & equilibrium, bField = equilibrium->b_field(polarPos, "forwards"); - bField = equilibrium->b_field_cart(bField, polarPos[2], 0); + bField = equilibrium->b_field_cart(bField, polarPos[2]); double product = 0; for (int i = 0; i < 3; i++) {