Skip to content

Ability to scale meshes in AEGIS #44

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

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

Ability to scale meshes in AEGIS #44

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
8ec75eb
Testing scaling meshes after read in
Waqar-ukaea Feb 21, 2025
cc1d91d
Added floating point tolerance to bounds checking for particle in mag…
Waqar-ukaea Feb 24, 2025
8e454bf
Refactored mesh_scale to function and various fixes
Waqar-ukaea Feb 24, 2025
4a78605
added mast config file
Waqar-ukaea Feb 24, 2025
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
2 changes: 1 addition & 1 deletion .vscode/launch.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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": [],
Expand Down
7 changes: 6 additions & 1 deletion .vscode/settings.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,12 @@
"*.in": "cpp",
"slist": "cpp",
"*.txx": "cpp",
"*.ipp": "cpp"
"*.ipp": "cpp",
"compare": "cpp",
"concepts": "cpp",
"numbers": "cpp",
"ranges": "cpp",
"span": "cpp"
},

}
2 changes: 1 addition & 1 deletion .vscode/tasks.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
2 changes: 1 addition & 1 deletion 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 RELEASE)
set(CMAKE_BUILD_TYPE DEBUG)


# Import VTK
Expand Down
2 changes: 1 addition & 1 deletion inres1/aegis_settings.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
43 changes: 43 additions & 0 deletions inres1/mast.json
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
}
}

12 changes: 7 additions & 5 deletions src/aegis_lib/EquilData.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,7 @@ EquilData::read_eqdsk(std::string filename)

// Fix for ITER eqdsks

OVERRIDE_ITER = true;
OVERRIDE_ITER = false;

if (OVERRIDE_ITER == true)
{
Expand Down Expand Up @@ -829,10 +829,12 @@ EquilData::b_field(const std::vector<double> & 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);
Expand Down
62 changes: 51 additions & 11 deletions src/aegis_lib/ParticleSimulation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -381,6 +381,21 @@ ParticleSimulation::read_params(const std::shared_ptr<JsonHandler> & configFile)
vectorOfTargetSurfs.push_back(i);
}
}

auto _str = aegisParams.get_optional<std::string>("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<std::string>("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 ";
Expand Down Expand Up @@ -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
Expand All @@ -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());
Expand Down Expand Up @@ -518,27 +557,28 @@ ParticleSimulation::loop_over_particles(int startIndex, int endIndex)
void
ParticleSimulation::setup_sources()
{

std::vector<double> triangleCoords(9);
std::vector<double> 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<moab::EntityHandle> 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;
Expand All @@ -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());
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
}

Expand Down
15 changes: 13 additions & 2 deletions src/aegis_lib/ParticleSimulation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,13 @@

using namespace moab;

enum class MeshUnits
{
UNSET,
M,
CM,
MM
};

enum class meshWriteOptions
{
Expand Down Expand Up @@ -131,7 +138,8 @@ class ParticleSimulation : public AegisBase
void write_particle_launch_positions(std::vector<double> &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
Expand All @@ -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<TriangleSource> listOfTriangles;
int totalNumberOfFacets = 0;
Expand Down Expand Up @@ -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;
Expand Down
3 changes: 1 addition & 2 deletions src/aegis_lib/Source.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,8 +102,7 @@ Sources::set_heatflux_params(const std::shared_ptr<EquilData> & 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++)
Expand Down