Mesh tally amalgamation #2

magnoxemo · 2025-01-09T20:32:00Z

Element Amalgamation - Post-Processing

This PR is for doing the element amalgamation in the post-processing. Once OpenMC is done with its simulation, Cardinal will read the mesh data. We can get the element_id and use it to get the current MOOSE element:
libMesh::Elem* elem_ptr = _mesh_filter->get_elem(elem_id);

Then we can check if that element is marked for AMALGAMATION:
if (elem_ptr->refinement_flag() == Elem::AMALGAMATE).

If the element is marked for AMALGAMATION, we have to find the cluster (for now, elements could be clustered only with the neighboring elements).

Tally score of the main element:
Real total_tally_of_the_cluster = tally_vals[local_score](ext_bin * _mesh_filter->n_bins() + e);
Volume of the main element:
Real total_volume_of_the_cluster = elem_ptr->volume();

We iterate through the neighbor elements (todo: need to check if the mapping is correct or not), and if they are also marked for AMALGAMATION, we increment total_tally_of_the_cluster and total_volume_of_the_cluster. Finally, the volumetric mean is calculated as:
Real power_fraction = elem_ptr->volume() * total_tally_of_the_cluster / total_volume_of_the_cluster;.

If the if (elem_ptr->refinement_flag() != Elem::AMALGAMATE), then presume the normal operation.

confusion:
total_tally_of_the_cluster += tally_vals[local_score](ext_bin * _mesh_filter->n_bins() + side_id*e); //need some clarification on that (not sure but ig side_id and e they aren't of the same data type)
@gonuke

magnoxemo · 2025-01-10T17:07:09Z

@lewisgross1296

…::MarkerValue. '

nuclearkevin

Took a look - great work! I've left a few comments which can be divided into C++ / programming style comments and some functional comments.

The only major functional issue I can see so far is in AmalgamationMarker which is that you're accessing _u_neighbor, which is not initialized by the base classes (QuadraturePointMarker and it's ancestors) as they don't participate in element neighbor loops. This is an issue with the marker system (there are no markers that loop over neighbors) which isn't easily solvable. I'll look into adding a new marker base class that addresses this since it's something that's come up a few times now as we've been working with AMR in MOOSE.

nuclearkevin · 2025-01-14T14:46:51Z

include/tallies/MeshTally.h

@@ -23,6 +23,7 @@

 #include "openmc/tallies/filter_mesh.h"

+#include "libmesh/mesh.h"


It's generally bad form to include a header file inside another header unless absolutely necessary. In this case, #include/libmesh/mesh.h should be in MeshTally.C as it's only used there.

Also, are you sure you need to include libmesh/mesh.h? We include libmesh/replicated_mesh in MeshTally.C, which should provide the same functionality.

nuclearkevin · 2025-01-14T14:50:19Z

src/markers/AmalgamationMarker.C

+AmalgamationMarker::AmalgamationMarker(const InputParameters& parameters)
+  : QuadraturePointMarker(parameters)
+{
+  _tolerance = getParam<Real>("tolerance");


I'd recommend setting tolerance in the initializer list; it's the preferred approach in C++.

nuclearkevin · 2025-01-14T14:55:36Z

src/markers/AmalgamationMarker.C

+  Real sum=0;
+  for (unsigned int qp=0;qp<_qrule->n_points();qp++)
+    {
+        sum=sum+_u_neighbor[qp];


You can use sum += _u_neighbor[qp]; instead.

I also don't think _u_neighbor is well posed in QuadraturePointMarker as this base class doesn't participate in element neighbour loops (one of my many gripes with the Indicator / Marker system). This isn't something that's easy to fix, I'll look into it and get back to you on this thought.

nuclearkevin · 2025-01-14T14:59:18Z

src/markers/AmalgamationMarker.C

+  // qurle->n_points() should never be zero. But do I need to use if statement to check? 
+  local_avg=sum/_qrule->n_points();
+  //local_avg it could be zero. So ig it needs to be checked 
+  if (local_avg==0)


Comparing a double (local_avg) to 0 will almost never return true due to the imprecise nature of floating point arithmetic. I'd use a fuzzy equals instead.

nuclearkevin · 2025-01-14T14:59:42Z

src/tallies/MeshTally.C

+              const unsigned int n_sides = elem_ptr->n_sides();
+              for (unsigned int side_id = 0; side_id < n_sides; ++side_id)
+              {
+                const libMesh::Elem* neighbor_ptr = elem_ptr->neighbor(side_id); //maybe wrong


This is correct.

nuclearkevin · 2025-01-14T15:02:20Z

src/tallies/MeshTally.C

+                {
+                    if (neighbor_ptr->refinement_flag()==Elem::AMALGAMATE)
+                    {
+                        total_tally_of_the_cluster += tally_vals[local_score](ext_bin * _mesh_filter->n_bins() + side_id + e);  //need some clarification on that


This will not work: you'd need an inverse map which goes from element DoF IDs to tally bin IDs (a _total_to_active_mapping). You'd then do:

auto neighbor_bin = _total_to_active[neighbor_ptr->id()]; total_tally_of_the_cluster += tally_vals[local_score](ext_bin * _mesh_filter->n_bins() + neighbor_bin);

…d only work given that neighbor_ids is subset of the element_ids

lewisgross1296 · 2025-01-16T15:51:28Z

include/markers/AmalgamationMarker.h

+  const Real _tolerance;
+  const Real TOLERENCE=1E-20; // to check if local_avg isn't zero.


I think this naming choice might be a little unintuitive. What is the difference between _tolerance and TOLERANCE? Hard to know from just looking, though you added a comment.

I'm also not sure comparing something to 1E-20 is the best way to check if something is non-zero. IIRC floating point precision means anything 1E-16 or smaller is basically zero. There might be some kind of C++ library function to check if a float is non-zero, but I don't know one off the top of my head.

libMesh has a built-in tolerance which I would recommend using instead of defining another - you can access it with libMesh::TOLERANCE.

lewisgross1296 · 2025-01-16T15:52:43Z

src/markers/AmalgamationMarker.C

+  Loop over the local elements and store the solution in a Real var*/
+  Real local_avg=0;
+  Real sum=0;
+  for (unsigned int qp=0;qp<_qrule->n_points();qp++)


Is this clang-formatted?

lewisgross1296 · 2025-01-16T15:57:19Z

src/markers/AmalgamationMarker.C

+    {
+        sum+=_u[qp];
+    }
+  // qurle->n_points() should never be zero. But do I need to use if statement to check?


You probably want to remove this comment before merging into cardinal. It looks like many places in MOOSE use _qrule -> n_points() without checking if it's zero.

lewisgross1296 · 2025-01-16T16:14:02Z

src/tallies/MeshTally.C

+      //error: 'class openmc::MeshFilter' has no member named 'get_elem'
+      //openmc::mesh_filter doesn't have anything like get_elem


Do you mean to have this in here? Let's learn about git add -p

lewisgross1296 · 2025-01-16T16:15:16Z

src/tallies/MeshTally.C

+                    //TODO
+					//I have to do an inverse map which goes from element DoF IDs to tally bin IDs
+                    //(a _total_to_active_mapping


Is this still a TODO or is that what _total_to_active does?

lewisgross1296 · 2025-01-16T16:18:36Z

include/tallies/MeshTally.h

@@ -113,4 +114,5 @@ class MeshTally : public TallyBase
  std::unique_ptr<libMesh::ReplicatedMesh> _libmesh_mesh_copy;
  /// A mapping between the elements in '_libmesh_mesh_copy' and the elements in the MooseMesh.
  std::vector<unsigned int> _active_to_total_mapping;
+  std::unordered_map<unsigned int, unsigned int>_total_to_active_mapping; //hash map for inverse mapping


Might be worth having a conversation about what this is named. @nuclearkevin what are your thoughts? I understand it's supposed to be the inverse mapping of element IDs to bins, but I'm not sure _total_to_active_mapping is a very intuitive name for what it does.

The original intent of the name for _active_to_total_mapping (and now _total_to_active_mapping) was to indicate that we're mapping from active elements to all elements. This uses libMesh's AMR terminology, where active elements are the elements participating in the solve (in the case of MeshTally the elements getting tallied on) and all elements includes the parent/child elements in the refinement hierarchy.

It's not the most intuitive naming scheme and I'm open to changing it to _bin_to_element_mapping and _element_to_bin_mapping (especially since we've added block restriction to MeshTally which uses _active_to_total_mapping) in an attempt to make the purpose of these maps clearer.

@magnoxemo we've merged the proposed name change into Cardinal (_active_to_total_mapping is now _bin_to_total_mapping). You should be good to change the name of _total_to_active_mapping to _element_to_bin_mapping.

magnoxemo and others added 8 commits December 30, 2024 14:39

setting environment to check if the test passes

bf85eb0

adding AmalgamationMarker

a45d17e

fixing a typo

61f3462

modifying the MeshTally.C file

0c94fab

fixing a typo

a003c24

Update .gitmodules

8ce26e6

fixing a typo

419e9d6

typo

6b16aa7

magnoxemo added 2 commits January 12, 2025 02:51

changing the mapping

9b10f5f

debuging the marker

7a0f55b

magnoxemo marked this pull request as draft January 12, 2025 20:21

magnoxemo added 3 commits January 13, 2025 10:52

marker fix

69cc3e6

marker test is failing error: 'AMALGAMATE' is not a member of 'Marker…

5c54c98

…::MarkerValue. '

update the elem_ptr

f810cc1

nuclearkevin reviewed Jan 14, 2025

View reviewed changes

magnoxemo added 2 commits January 14, 2025 10:49

implementing the suggestions from Kevin

41ce934

implementing hash map for the _total_to_active_maping. But this shoul…

513be2a

…d only work given that neighbor_ids is subset of the element_ids

magnoxemo requested review from gonuke and lewisgross1296 January 14, 2025 22:35

magnoxemo added 3 commits January 15, 2025 09:55

fixing typo

827cc1c

fixing typo

d30dc2c

fixed some bugs. Tests are passing but with a warning

aeb019b

lewisgross1296 reviewed Jan 16, 2025

View reviewed changes

applying clang format

a21f13d

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Mesh tally amalgamation #2

Mesh tally amalgamation #2

magnoxemo commented Jan 9, 2025 •

edited

Loading

magnoxemo commented Jan 10, 2025

nuclearkevin left a comment

nuclearkevin Jan 14, 2025

nuclearkevin Jan 14, 2025

nuclearkevin Jan 14, 2025

nuclearkevin Jan 14, 2025

nuclearkevin Jan 14, 2025

nuclearkevin Jan 14, 2025

lewisgross1296 Jan 16, 2025

nuclearkevin Jan 16, 2025

lewisgross1296 Jan 16, 2025

lewisgross1296 Jan 16, 2025

lewisgross1296 Jan 16, 2025

lewisgross1296 Jan 16, 2025

lewisgross1296 Jan 16, 2025

nuclearkevin Jan 16, 2025 •

edited

Loading

nuclearkevin Jan 17, 2025

magnoxemo Jan 17, 2025

		@@ -23,6 +23,7 @@

		#include "openmc/tallies/filter_mesh.h"

		#include "libmesh/mesh.h"

		const Real _tolerance;
		const Real TOLERENCE=1E-20; // to check if local_avg isn't zero.

		//error: 'class openmc::MeshFilter' has no member named 'get_elem'
		//openmc::mesh_filter doesn't have anything like get_elem

Mesh tally amalgamation #2

Are you sure you want to change the base?

Mesh tally amalgamation #2

Conversation

magnoxemo commented Jan 9, 2025 • edited Loading

magnoxemo commented Jan 10, 2025

nuclearkevin left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nuclearkevin Jan 16, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

magnoxemo commented Jan 9, 2025 •

edited

Loading

nuclearkevin Jan 16, 2025 •

edited

Loading