openmc-dev / openmc / 12605079372

#include "openmc/random_ray/flat_source_domain.h"

#include "openmc/cell.h"

#include "openmc/eigenvalue.h"

#include "openmc/geometry.h"

#include "openmc/material.h"

#include "openmc/message_passing.h"

#include "openmc/mgxs_interface.h"

#include "openmc/output.h"

#include "openmc/plot.h"

#include "openmc/random_ray/random_ray.h"

#include "openmc/simulation.h"

#include "openmc/tallies/filter.h"

#include "openmc/tallies/tally.h"

#include "openmc/tallies/tally_scoring.h"

#include "openmc/timer.h"

#include <cstdio>

namespace openmc {

//==============================================================================

// FlatSourceDomain implementation

//==============================================================================

// Static Variable Declarations

RandomRayVolumeEstimator FlatSourceDomain::volume_estimator_ {

  RandomRayVolumeEstimator::HYBRID};

bool FlatSourceDomain::volume_normalized_flux_tallies_ {false};

bool FlatSourceDomain::adjoint_ {false};

{

  // Count the number of source regions, compute the cell offset

  // indices, and store the material type The reason for the offsets is that

  // some cell types may not have material fills, and therefore do not

  // produce FSRs. Thus, we cannot index into the global arrays directly

    } else {

    }

  }

  // Initialize cell-wise arrays

  // Initialize element-wise arrays

    // If in eigenvalue mode, set starting flux to guess of unity

  } else {

    // If in fixed source mode, set starting flux to guess of zero

    // and initialize external source arrays

  }

  // Initialize material array

      }

    }

  }

  // Sanity check

  }

  // Initialize tally volumes

      //  Get the shape of the 3D result tensor

      // Create a new 2D tensor with the same size as the first

      // two dimensions of the 3D tensor

    }

  }

  // Compute simulation domain volume based on ray source

{

  // Reset scalar fluxes, iteration volume tallies, and region hit flags to

  // zero

{

  }

// Compute new estimate of scattering + fission sources in each source region

// based on the flux estimate from the previous iteration.

{

  // Add scattering source

        double sigma_s =

      }

    }

  }

  // Add fission source

      }

    }

  }

  // Add external source if in fixed source mode

    }

  }

// Normalizes flux and updates simulation-averaged volume estimate

  double total_active_distance_per_iteration)

{

// Normalize scalar flux to total distance travelled by all rays this iteration

  }

// Accumulate cell-wise ray length tallies collected this iteration, then

// update the simulation-averaged cell-wise volume estimates

  }

  int64_t idx, double volume, int material, int g)

{

{

}

{

// Combine transport flux contributions and flat source contributions from the

// previous iteration to generate this iteration's estimate of scalar flux.

{

    // Increment the number of hits if cell was hit this iteration

    }

    // Check if an external source is present in this source region

    bool external_source_present =

    // The volume treatment depends on the volume estimator type

    // and whether or not an external source is present in the cell.

    double volume;

      } else {

      }

    default:

      fatal_error("Invalid volume estimator type");

    }

      // There are three scenarios we need to consider:

        // 1. If the FSR was hit this iteration, then the new flux is equal to

        // the flat source from the previous iteration plus the contributions

        // from rays passing through the source region (computed during the

        // transport sweep)

        // 2. If the FSR was not hit this iteration, but has been hit some

        // previous iteration, then we need to make a choice about what

        // to do. Naively we will usually want to set the flux to be equal

        // to the reduced source. However, in fixed source problems where

        // there is a strong external source present in the cell, and where

        // the cell has a very low cross section, this approximation will

        // cause a huge upward bias in the flux estimate of the cell (in these

        // conditions, the flux estimate can be orders of magnitude too large).

        // Thus, to avoid this bias, if any external source is present

        // in the cell we will use the previous iteration's flux estimate. This

        // injects a small degree of correlation into the simulation, but this

        // is going to be trivial when the miss rate is a few percent or less.

          set_flux_to_old_flux(idx);

        } else {

        }

      }

      // If the FSR was not hit this iteration, and it has never been hit in

      // any iteration (i.e., volume is zero), then we want to set this to 0

      // to avoid dividing anything by a zero volume. This happens implicitly

      // given that the new scalar flux arrays are set to zero each iteration.

    }

  }

  // Return the number of source regions that were hit this iteration

}

// Generates new estimate of k_eff based on the differences between this

// iteration's estimate of the scalar flux and the last iteration's estimate.

{

  // Vector for gathering fission source terms for Shannon entropy calculation

    // If simulation averaged volume is zero, don't include this cell

      continue;

    }

    }

    // Compute total fission rates in FSR

    // Accumulate totals

    // Store total fission rate in the FSR for Shannon calculation

  }

  // defining an inverse sum for better performance

    // Only if FSR has non-negative and non-zero fission source

      // Normalize to total weight of bank sites. p_i for better performance

      // Sum values to obtain Shannon entropy.

    }

  }

  // Adds entropy value to shared entropy vector in openmc namespace.

// This function is responsible for generating a mapping between random

// ray flat source regions (cell instances) and tally bins. The mapping

// takes the form of a "TallyTask" object, which accounts for one single

// score being applied to a single tally. Thus, a single source region

// may have anywhere from zero to many tally tasks associated with it ---

// meaning that the global "tally_task" data structure is in 2D. The outer

// dimension corresponds to the source element (i.e., each entry corresponds

// to a specific energy group within a specific source region), and the

// inner dimension corresponds to the tallying task itself. Mechanically,

// the mapping between FSRs and spatial filters is done by considering

// the location of a single known ray midpoint that passed through the

// FSR. I.e., during transport, the first ray to pass through a given FSR

// will write down its midpoint for use with this function. This is a cheap

// and easy way of mapping FSRs to spatial tally filters, but comes with

// the downside of adding the restriction that spatial tally filters must

// share boundaries with the physical geometry of the simulation (so as

// not to subdivide any FSR). It is acceptable for a spatial tally region

// to contain multiple FSRs, but not the other way around.

// TODO: In future work, it would be preferable to offer a more general

// (but perhaps slightly more expensive) option for handling arbitrary

// spatial tallies that would be allowed to subdivide FSRs.

// Besides generating the mapping structure, this function also keeps track

// of whether or not all flat source regions have been hit yet. This is

// required, as there is no guarantee that all flat source regions will

// be hit every iteration, such that in the first few iterations some FSRs

// may not have a known position within them yet to facilitate mapping to

// spatial tally filters. However, after several iterations, if all FSRs

// have been hit and have had a tally map generated, then this status will

// be passed back to the caller to alert them that this function doesn't

// need to be called for the remainder of the simulation.

{

  // Tracks if we've generated a mapping yet for all source regions.

// Attempt to generate mapping for all source regions

    // If this source region has not been hit by a ray yet, then

    // we aren't going to be able to map it, so skip it.

      all_source_regions_mapped = false;

      continue;

    }

    // A particle located at the recorded midpoint of a ray

    // crossing through this source region is used to estabilish

    // the spatial location of the source region

    // Loop over energy groups (so as to support energy filters)

      // Set particle to the current energy

      // If this task has already been populated, we don't need to do

      // it again.

        continue;

      }

      // Loop over all active tallies. This logic is essentially identical

      // to what happens when scanning for applicable tallies during

      // MC transport.

        // Initialize an iterator over valid filter bin combinations.

        // If there are no valid combinations, use a continue statement

        // to ensure we skip the assume_separate break below.

        // Loop over filter bins.

          // Loop over scores

               score_index++) {

            // If a valid tally, filter, and score combination has been found,

            // then add it to the list of tally tasks for this source element.

            // Also add this task to the list of volume tasks for this source

            // region.

          }

        }

      }

      // Reset all the filter matches for the next tally event.

    }

// Set the volume accumulators to zero for all tallies

{

    }

  }

// In fixed source mode, due to the way that volumetric fixed sources are

// converted and applied as volumetric sources in one or more source regions,

// we need to perform an additional normalization step to ensure that the

// reported scalar fluxes are in units per source neutron. This allows for

// direct comparison of reported tallies to Monte Carlo flux results.

// This factor needs to be computed at each iteration, as it is based on the

// volume estimate of each FSR, which improves over the course of the

// simulation

{

  // If we are not in fixed source mode, then there are no external sources

  // so no normalization is needed.

  }

  // Step 1 is to sum over all source regions and energy groups to get the

  // total external source strength in the simulation.

    }

  }

  // Step 2 is to determine the total user-specified external source strength

  }

  // The correction factor is the ratio of the user-specified external source

  // strength to the simulation external source strength.

    user_external_source_strength / simulation_external_source_strength;

}

// Tallying in random ray is not done directly during transport, rather,

// it is done only once after each power iteration. This is made possible

// by way of a mapping data structure that relates spatial source regions

// (FSRs) to tally/filter/score combinations. The mechanism by which the

// mapping is done (and the limitations incurred) is documented in the

// "convert_source_regions_to_tallies()" function comments above. The present

// tally function simply traverses the mapping data structure and executes

// the scoring operations to OpenMC's native tally result arrays.

{

  // Reset our tally volumes to zero

  double source_normalization_factor =

// We loop over all source regions and energy groups. For each

// element, we check if there are any scores needed and apply

// them.

    // The fsr.volume_ is the unitless fractional simulation averaged volume

    // (i.e., it is the FSR's fraction of the overall simulation volume). The

    // simulation_volume_ is the total 3D physical volume in cm^3 of the

    // entire global simulation domain (as defined by the ray source box).

    // Thus, the FSR's true 3D spatial volume in cm^3 is found by multiplying

    // its fraction of the total volume by the total volume. Not important in

    // eigenvalue solves, but useful in fixed source solves for returning the

    // flux shape with a magnitude that makes sense relative to the fixed

    // source strength.

      // Determine numerical score value

        double score;

        case SCORE_TOTAL:

          score = flux * volume * sigma_t_[material * negroups_ + g];

          break;

        case SCORE_EVENTS:

          score = 1.0;

          break;

        default:

          fatal_error("Invalid score specified in tallies.xml. Only flux, "

                      "total, fission, nu-fission, and events are supported in "

                      "random ray mode.");

          break;

        }

        // Apply score to the appropriate tally bin

#pragma omp atomic

          score;

      }

    }

    // For flux tallies, the total volume of the spatial region is needed

    // for normalizing the flux. We store this volume in a separate tensor.

    // We only contribute to each volume tally bin once per FSR.

#pragma omp atomic

            volume;

        }

      }

    }

  } // end FSR loop

  // Normalize any flux scores by the total volume of the FSRs scoring to that

  // bin. To do this, we loop over all tallies, and then all filter bins,

  // and then scores. For each score, we check the tally data structure to

  // see what index that score corresponds to. If that score is a flux score,

  // then we divide it by volume.

            }

          }

        }

      }

    }

  }

{

#ifdef OPENMC_MPI

  // If we only have 1 MPI rank, no need

  // to reduce anything.

    return;

  // The "position_recorded" variable needs to be allreduced (and maxed),

  // as whether or not a cell was hit will affect some decisions in how the

  // source is calculated in the next iteration so as to avoid dividing

  // by zero. We take the max rather than the sum as the hit values are

  // expected to be zero or 1.

    MPI_INT, MPI_MAX, mpi::intracomm);

  // The position variable is more complicated to reduce than the others,

  // as we do not want the sum of all positions in each cell, rather, we

  // want to just pick any single valid position. Thus, we perform a gather

  // and then pick the first valid position we find for all source regions

  // that have had a position recorded. This operation does not need to

  // be broadcast back to other ranks, as this value is only used for the

  // tally conversion operation, which is only performed on the master rank.

  // While this is expensive, it only needs to be done for active batches,

  // and only if we have not mapped all the tallies yet. Once tallies are

  // fully mapped, then the position vector is fully populated, so this

  // operation can be skipped.

  // First, we broadcast the fully mapped tally status variable so that

  // all ranks are on the same page

  // Then, we perform the gather of position data, if needed

    // Master rank will gather results and pick valid positions

      // Initialize temporary vector for receiving positions

      }

      // Copy master rank data into gathered vector for convenience

      // Receive all data into gather vector

          0, mpi::intracomm, MPI_STATUS_IGNORE);

      }

      // Scan through gathered data and pick first valid cell posiiton

            }

          }

        }

      }

      // Other ranks just send in their data

        mpi::intracomm);

    }

  }

  // For the rest of the source region data, we simply perform an all reduce,

  // as these values will be needed on all ranks for transport during the

  // next iteration.

    MPI_SUM, mpi::intracomm);

    MPI_DOUBLE, MPI_SUM, mpi::intracomm);

#endif

  Position r, int64_t sr, int g) const

{

}

// Outputs all basic material, FSR ID, multigroup flux, and

// fission source data to .vtk file that can be directly

// loaded and displayed by Paraview. Note that .vtk binary

// files require big endian byte ordering, so endianness

// is checked and flipped if necessary.

{

  // Rename .h5 plot filename(s) to .vtk filenames

  }

  // Print header information

  // Outer loop over plots

    // Get handle to OpenMC plot object and extract params

    // Random ray plots only support voxel plots

                          "is allowed in random ray mode.",

        p));

                          "is allowed in random ray mode.",

        p));

    }

    // Perform sanity checks on file size

              "slow.");

    }

    // Relate voxel spatial locations to random ray source regions

#pragma omp parallel for collapse(3)

    for (int z = 0; z < Nz; z++) {

      for (int y = 0; y < Ny; y++) {

        for (int x = 0; x < Nx; x++) {

          Position sample;

          sample.z = ll.z + z_delta / 2.0 + z * z_delta;

          sample.y = ll.y + y_delta / 2.0 + y * y_delta;

          sample.x = ll.x + x_delta / 2.0 + x * x_delta;

          Particle p;

          p.r() = sample;

          bool found = exhaustive_find_cell(p);

          int i_cell = p.lowest_coord().cell;

          int64_t source_region_idx =

            source_region_offsets_[i_cell] + p.cell_instance();

          voxel_indices[z * Ny * Nx + y * Nx + x] = source_region_idx;

          voxel_positions[z * Ny * Nx + y * Nx + x] = sample;

        }

      }

    }

    double source_normalization_factor =

    // Open file for writing

    // Write vtk metadata

    // Plot multigroup flux data

      }

    }

    // Plot FSRs

    }

    // Plot Materials

    }

    // Plot fission source

      }

    }

  }

}

  Discrete* discrete, double strength_factor, int64_t source_region)

{

  }