diff --git a/.vscode/launch.json b/.vscode/launch.json index 3f04732..5d83244 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -125,7 +125,7 @@ "type": "cppdbg", "request": "launch", "program": "${workspaceRoot}/bin/aegis", - "args": ["aegis_settings.json"], + "args": ["mast.json"], "stopAtEntry": false, "cwd": "${workspaceRoot}/inres1", "environment": [], diff --git a/.vscode/settings.json b/.vscode/settings.json index f13b1c5..296542c 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -83,7 +83,12 @@ "*.in": "cpp", "slist": "cpp", "*.txx": "cpp", - "*.ipp": "cpp" + "*.ipp": "cpp", + "compare": "cpp", + "concepts": "cpp", + "numbers": "cpp", + "ranges": "cpp", + "span": "cpp" }, } \ No newline at end of file diff --git a/.vscode/tasks.json b/.vscode/tasks.json index d882b49..b97d208 100644 --- a/.vscode/tasks.json +++ b/.vscode/tasks.json @@ -7,7 +7,7 @@ "options": { "cwd": "${workspaceRoot}/bld" }, - "command": "cmake ../ -DDAGMC_DIR=/home/waqar/dagmc_bld/DAGMC" + "command": "cmake ../ -DDAGMC_DIR=/home/waqar/aegis-deps/DAGMC" }, { "label": "clang-format-files", diff --git a/CMakeLists.txt b/CMakeLists.txt index 5a7ad7f..731ca67 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 RELEASE) +set(CMAKE_BUILD_TYPE DEBUG) # Import VTK diff --git a/inres1/aegis_settings.json b/inres1/aegis_settings.json index 4a65175..c837a95 100644 --- a/inres1/aegis_settings.json +++ b/inres1/aegis_settings.json @@ -4,7 +4,7 @@ "DAGMC": "inres1_shad.h5m", "step_size": 0.001, "max_steps": 100000, - "launch_pos": "random", + "launch_pos": "fixed", "monte_carlo_params":{ "number_of_particles_per_facet":10, diff --git a/inres1/mast.json b/inres1/mast.json new file mode 100644 index 0000000..4120f1f --- /dev/null +++ b/inres1/mast.json @@ -0,0 +1,43 @@ +{ + "description": "inres1 test case to compare against SMARDDA", + "aegis_params":{ + "DAGMC": "/home/waqar/moose-workshop/mast-dagmc.h5m", + "step_size": 0.001, + "max_steps": 1000000, + "launch_pos": "fixed", + + "monte_carlo_params":{ + "number_of_particles_per_facet":1, + "write_particle_launch_positions":true + }, + + "target_surfs": [79436,79437,79438,79439,79440,79441,79442,79443], + "mesh_units": "cm", + "write_mesh": "target", + "force_no_deposition": false, + "coordinate_system": "xyz", + "execution_type": "dynamic", + + "dynamic_batch_params":{ + "batch_size":16, + "worker_profiling": false, + "debug":false + } + }, + "equil_params":{ + "eqdsk": "/home/waqar/moose-workshop/equil_MAST.eqdsk", + "power_sol": 3.2e+06, + "r_outrbdry": 8.07385, + "lambda_q": 6.187, + "cenopt": 2, + "psiref": -2.1628794, + "rmove": 0.0, + "draw_equil_rz": false, + "draw_equil_xyz": false, + "print_debug_info": false + }, + "vtk_params":{ + "draw_particle_tracks": true + } + } + diff --git a/src/aegis_lib/EquilData.cpp b/src/aegis_lib/EquilData.cpp index 5291b6b..1e05d67 100644 --- a/src/aegis_lib/EquilData.cpp +++ b/src/aegis_lib/EquilData.cpp @@ -223,7 +223,7 @@ EquilData::read_eqdsk(std::string filename) // Fix for ITER eqdsks - OVERRIDE_ITER = true; + OVERRIDE_ITER = false; if (OVERRIDE_ITER == true) { @@ -829,10 +829,12 @@ EquilData::b_field(const std::vector & position, std::string startingFro zz = polarPosition[1]; } - if (zr < rmin || zr > rmax || zz < zmin || zz > zmax) - { - return {}; - } + double tolerance = 1e-9; + // check if out of bounds + if (zr < rmin - tolerance) return {}; + if (zr > rmax + tolerance) return {}; + if (zz < zmin - tolerance) return {}; + if (zz > zmax + tolerance) return {}; // evaluate psi and psi derivs at given (R,Z) coords alglib::spline2ddiff(psiSpline, zr, zz, zpsi, zdpdr, zdpdz, null); diff --git a/src/aegis_lib/ParticleSimulation.cpp b/src/aegis_lib/ParticleSimulation.cpp index 454a3cd..ad773a1 100644 --- a/src/aegis_lib/ParticleSimulation.cpp +++ b/src/aegis_lib/ParticleSimulation.cpp @@ -381,6 +381,21 @@ ParticleSimulation::read_params(const std::shared_ptr & configFile) vectorOfTargetSurfs.push_back(i); } } + + auto _str = aegisParams.get_optional("mesh_units"); + if (_str) { // Check if the optional has a value and then set to enum + if (*_str == "m") meshDimension = MeshUnits::M; + else if (*_str == "cm") meshDimension = MeshUnits::CM; + else if (*_str == "mm") meshDimension = MeshUnits::MM; + } + + _str = aegisParams.get_optional("write_mesh"); + if (_str) { // Check if the optional has a value and then set to enum + if (*_str == "target") writeMesh = meshWriteOptions::TARGET; + else if (*_str == "full") writeMesh = meshWriteOptions::FULL; + else if (*_str == "partial") writeMesh = meshWriteOptions::PARTIAL; + } + } std::string launchType = "Launching particles from triangles at "; @@ -415,6 +430,28 @@ ParticleSimulation::select_coordinate_system() } } + +void ParticleSimulation::scale_mesh(MeshUnits from) +{ + float scalingFactor = 1.0f; + switch (from) + { + case MeshUnits::CM: + scalingFactor = 0.01f; + break; + + case MeshUnits::MM: + scalingFactor = 0.001f; + break; + } + + for (auto &i:nodeCoords) + { + i = i*scalingFactor; + } + DAG->moab_instance()->set_coords(&nodesList[0], nodesList.size(), nodeCoords.data()); +} + // initialise CAD geometry for AEGIS and magnetic field equilibrium for CAD // Geometry void @@ -432,6 +469,8 @@ ParticleSimulation::init_geometry() nodeCoords.resize(nodesList.size() * 3); DAG->moab_instance()->get_coords(&nodesList[0], nodesList.size(), nodeCoords.data()); + if (meshDimension != MeshUnits::UNSET) scale_mesh (meshDimension); // scale mesh to Metres + totalNumberOfFacets = facetsList.size(); implComplementVol = DAG->entity_by_index(3, volsList.size()); @@ -518,27 +557,28 @@ ParticleSimulation::loop_over_particles(int startIndex, int endIndex) void ParticleSimulation::setup_sources() { - std::vector triangleCoords(9); std::vector trianglePtsA(3), trianglePtsB(3), trianglePtsC(3); + + // Preallocate exactly num_particles_launched() elements listOfParticles.reserve(num_particles_launched()); - for (const auto & facetEH : targetFacets) - { + + // Index-based modification + size_t particleIndex = 0; + for (const auto &facetEH : targetFacets) { std::vector triangleNodes; DAG->moab_instance()->get_adjacencies(&facetEH, 1, 0, false, triangleNodes); DAG->moab_instance()->get_coords(&triangleNodes[0], triangleNodes.size(), triangleCoords.data()); - for (int j = 0; j < 3; j++) - { + for (int j = 0; j < 3; j++) { trianglePtsA[j] = triangleCoords[j]; trianglePtsB[j] = triangleCoords[j + 3]; trianglePtsC[j] = triangleCoords[j + 6]; } TriangleSource triangle(trianglePtsA, trianglePtsB, trianglePtsC, facetEH, particleLaunchPos); - if (!equilibrium->check_if_in_bfield(triangle.launch_pos())) - { + if (!equilibrium->check_if_in_bfield(triangle.launch_pos())) { log_string(LogLevel::INFO, "Triangle start outside of magnetic field. Skipping to next triangle..."); continue; @@ -548,8 +588,8 @@ ParticleSimulation::setup_sources() double heatflux = equilibrium->omp_power_dep(psid, triangle.BdotN(), "exp"); triangle.update_heatflux(heatflux); double heatfluxPerParticle = heatflux / nParticlesPerFacet; - for (int i = 0; i < nParticlesPerFacet; ++i) - { + + for (int i = 0; i < nParticlesPerFacet; ++i) { auto launchPos = (particleLaunchPos == "fixed") ? triangle.centroid() : triangle.random_pt(); ParticleBase particle(coordinateSystem::CARTESIAN, launchPos, heatfluxPerParticle, facetEH); particle.align_dir_to_surf(triangle.BdotN()); @@ -831,7 +871,7 @@ ParticleSimulation::Execute_serial() // attach_mesh_attribute("Psi Start", targetFacets, psiStart); // attach_mesh_attribute("BdotN", targetFacets, bnList); // attach_mesh_attribute("Normal", targetFacets, normalList); - write_out_mesh(meshWriteOptions::BOTH); + write_out_mesh(writeMesh); mainParticleLoopTime = MPI_Wtime() - mainParticleLoopStart; // write out data and print final @@ -922,7 +962,7 @@ ParticleSimulation::Execute_mpi() } attach_mesh_attribute("Heatflux", targetFacets, rootQvalues); - write_out_mesh(meshWriteOptions::BOTH, targetFacets); + write_out_mesh(writeMesh, targetFacets); mainParticleLoopTime = MPI_Wtime() - mainParticleLoopStart; } diff --git a/src/aegis_lib/ParticleSimulation.h b/src/aegis_lib/ParticleSimulation.h index 89a4c9f..2d2c98a 100644 --- a/src/aegis_lib/ParticleSimulation.h +++ b/src/aegis_lib/ParticleSimulation.h @@ -59,6 +59,13 @@ using namespace moab; +enum class MeshUnits +{ + UNSET, + M, + CM, + MM +}; enum class meshWriteOptions { @@ -131,7 +138,8 @@ class ParticleSimulation : public AegisBase void write_particle_launch_positions(std::vector &particleHeatfluxes); void mesh_coord_transform(coordinateSystem coordSys); void select_coordinate_system(); - + void scale_mesh(MeshUnits from); + // Simulation config parameters std::string exeType = "dynamic"; std::string dagmcInputFile; // no defaults because @@ -150,6 +158,9 @@ class ParticleSimulation : public AegisBase bool noMidplaneTermination = false; int nParticlesPerFacet = 1; // Number of particles launched per facet bool writeParticleLaunchPos = false; + MeshUnits meshDimension = MeshUnits::UNSET; + meshWriteOptions writeMesh = meshWriteOptions::TARGET; + std::vector listOfTriangles; int totalNumberOfFacets = 0; @@ -178,7 +189,7 @@ class ParticleSimulation : public AegisBase moab::EntityHandle implComplementVol; double nextSurfDist = 0.0; // distance to next surface moab::EntityHandle intersectedFacet; - int rayOrientation = 1; // rays are fired along surface normals + int rayOrientation = -1; // rays are fired along surface normals double trackLength = 0.0; double startTime; diff --git a/src/aegis_lib/Source.cpp b/src/aegis_lib/Source.cpp index 24d9151..cf84110 100644 --- a/src/aegis_lib/Source.cpp +++ b/src/aegis_lib/Source.cpp @@ -102,8 +102,7 @@ Sources::set_heatflux_params(const std::shared_ptr & equilibrium, polarPos = CoordTransform::cart_to_polar(launchPos); - bField = equilibrium->b_field(polarPos, "forwards"); - + bField = equilibrium->b_field(polarPos, "polar"); bField = equilibrium->b_field_cart(bField, polarPos[2]); double product = 0; for (int i = 0; i < 3; i++)