Skip to content

Conversation

@Fuad-HH
Copy link
Contributor

@Fuad-HH Fuad-HH commented Feb 27, 2025

As a continuation of standardizing the search method (next steps of #143), this wraps the search related particle state arrays inside this class. It includes the following tasks:

  • Create ParticleTracer class
  • Write test case
  • Add consequent search in the test case
  • remove finishedUnmoved
  • Move particle_done out as a parameter
  • Make maxLoop a variable in ParticleTracer class
  • Add a test case with GMSH mesh to test node orientation check or correction
  • Rename trace_particle_through_mesh function parameters
  • Clean up: hardcoded parameters, comments, etc.

- passing all tests except smoke_test_particle
It is giving a warning for MPI oversubscribing
and "Structure not initalized at Particle"
- passing all tests
instead of checking if tol was passed using negative values,
it is made optional. In the future, tol will be calculated
outside according to the standardization plan and for this
reason, we don't need it to be passed as reference or pointer
it is passed as a parameter rather than computing it every time
the search is called.
it's an initial step to make the function standard with
proper docuemntation. [important] there are some breaking
changes and the function is under heavy development
the functor decides what happens to the particles at the boundary
next_element parameter holds the next element when a particle
intersects a face. This is added to both the trace_particle_through_mesh
and the functor. The functor may or may not use the value, but there's
a high chance to be used.

[Important] The test is failing. This is because although the particles
doesn't intersect face 0 (it is along the way of the ray going forward)
the closeness algorithm finds it as the intersected face.
- works for particles going to a new element
- or remains in the same

- [to be fixed] particles leaking out
next_element becomes -1 when particle leaks
fix search_mesh in accordance with changes in
trace_particle_through_mesh

[future fixes]
- interPoints should give the face intersection location when leaking
  out instead of the destination
- same test case repurposed from search_mesh test case
@Fuad-HH
Copy link
Contributor Author

Fuad-HH commented Feb 27, 2025

@jacobmerson, Could you please take a look? Does the class look like you suggested?

- it creates confusion without solving anything significant
- tested
- it should save the allocation time for every call of search
Comment on lines 50 to 56
Omega_h::Write<Omega_h::LO> GetElemIds() { return elem_ids_; }

Omega_h::Write<Omega_h::LO> GetInterFaces() { return inter_faces_; }

Omega_h::Write<Omega_h::Real> GetInterPoints() { return inter_points_; }

Omega_h::Write<Omega_h::LO> GetLastExits() { return last_exits_; }
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Couple comments on this:

  • Please make any words fully spelled out as it is easier to understand. i.e. GetElemIds -> GetElementIds, GetInterFaces -> GetIntersectionFaces, etc.
  • The second question is more philosophical. Why should these be in the public interface? Do we expect the user of this class to need to modify this?
  • Should be marked [[nodiscard]] and const and should not return Write instead, should return a Read, or better a span<const T>

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have changed them to Read. It makes more sense here since we don't want the user to change the internal state of the particles. However, I think there can be cases where the user may want to know the states after the search, for example, through which face particles leaked. It can be argued that these can be done inside the functor, but it will provide an easier option to read them whereas the functor can be difficult to understand for a new user of this class.

Comment on lines +32 to +34
auto particle_orig = ptcls_->template get<0>();
auto particle_dest = ptcls_->template get<1>();
auto particle_ids = ptcls_->template get<2>();
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should provide indirection here so the user can decide which IDs correspond to the locations. This will also make it easier if we try to move the other particle arrays into the particle datastructure.

The way to do this is as follows:

  1. Create a template struct called ParticleDataStructureIndexMap we may change name later. This should take a single template argument which will be the particle datastructure.
  2. Overload the map for your particle datastructure with constexpr int members called OriginalPosition, DestinationPosition, ParticleID. For your default particle datastructure they will be set to 0, 1, 2.
  3. Have your class take an additional template argument which will be this map. The default type should be ParticleDataStructureIndexMap<ParticleType>

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you mean something like this?

template<typename ParticleType>
struct ParticleDataStructureIndexMap {
    static constexpr int OriginalPosition       = 0;
    static constexpr int DestinationPosition    = 1;
    static constexpr int ParticleID             = 2;
};

template<typename ParticleType, typename Func, typename ParticleIDMap = ParticleDataStructureIndexMap<ParticleType>>
class ParticleTracer {
...

If this is what you mean,

  1. Why does the struct need the ParticleType template?
  2. Can't it be just an enum class?
  3. If a user wants to specify the map themselves, they have to specify all three of them which can be difficult.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this is what I mean. This extra indirection let's us decouple the position of the information in the particle struct from the the items. I.e., if the user wants to put the location in different positions they can easily do that.

I'm not 100% sure that we want to include the template on the particle structure since you could have the same particle structure type with different semantic meaning. But the goal there is that the compiler will "intelligently" pick the right index based on the particle structure the user passes in.

Omega_h::Write<Omega_h::LO> ptcl_done_;


bool validate_internal_data_sizes() {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what is this function doing? Needs a comment that describes what this check does.

return true;
}

void reset_particle_done() {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The name done is still a bit confusing to me. What about particle_requires_tracing or particle_at_final_location?

Comment on lines +329 to +330
//o::Real quality = -1;
//o::LO bestFace = -1;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove commented code.

Comment on lines 347 to 351
//if (dproj > -tol && (quality < 0 || closeness < quality) && lastExit[ptcl] == -1) {
// quality = closeness;
// bestFace = face_id;
// xPoints[3*ptcl] = xpts[0]; xPoints[3*ptcl + 1] = xpts[1]; xPoints[3*ptcl + 2] = xpts[2];
//}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove commented code

Comment on lines 355 to 358
//if (lastExit[ptcl] == -1) {
// printf("!!!!!!! Best face %d\n", bestFace);
// lastExit[ptcl] = bestFace;
//}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove commented code.

xPoints[3*ptcl+2] = x_ps_tgt(ptcl,2);
}
ptcl_done[ptcl] = (lastExit[ptcl] == -1);
//ptcl_done[ptcl] = (lastExit[ptcl] == -1);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove commented code

}
else {
elem_ids[pid] = exposed ? -1 : elem_ids[pid]; //leaves domain if exposed
auto checkIntersectionwGeomBoundary = PS_LAMBDA(const int e, const int pid, const int mask) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rename this to intersectsWithGeometryInterface when you write functions that will return a bool (or array of bool) put a question mark after the name in your head. If the answer to that question matches the result of your boolean then it will be clear what the function is doing.

With the current name, if you say "check the intersection with the geometry boundary?` the name does not give any clues about what true means and what false means.

- added update particle position function for the class
- constructor gets picparts instead of the mesh
- rebuild: update particle positions and migrate (based on the passed
  flag)
- major bug
- Related test cases are fixed too
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants