openmc-dev / openmc / 15799682238

#include "openmc/random_ray/flat_source_domain.h"

#include "openmc/cell.h"

#include "openmc/constants.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 "openmc/weight_windows.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};

double FlatSourceDomain::diagonal_stabilization_rho_ {1.0};

std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>

  FlatSourceDomain::mesh_domain_map_;

{

  // 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 source regions.

  // Initialize materials

      }

    }

  }

  // 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 and iteration volume tallies to zero

  }

  }

{

  }

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

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

{

// Reset all source regions to zero (important for void regions)

  }

  // Add scattering + fission source

    }

        double sigma_s =

      }

    }

  }

  // 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 sr, double volume, int g)

{

    }

  } else {

  }

{

{

// 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

    }

    // Set the SR to small status if its expected number of hits

    // per iteration is less than 1.5

    } else {

    }

    // 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.

        } else {

        }

      }

      // Halt if NaN implosion is detected

        fatal_error("A source region scalar flux is not finite. "

                    "This indicates a numerical instability in the "

                    "simulation. Consider increasing ray density or adjusting "

                    "the source region mesh.");

      }

    }

  }

  // 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;

    }

      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

            // 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.

      // For non-void regions, we store the external source pre-divided by

      // sigma_t. We need to multiply non-void regions back up by sigma_t

      // to get the total source strength in the expected units.

      }

    }

  }

  // 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.

      double flux =

      // Determine numerical score value

        case SCORE_TOTAL:

          if (material != MATERIAL_VOID) {

            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.

            }

          }

        }

      }

    }

  }

  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) reduction(min : min_weight)

    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;

          p.r_last() = sample;

          p.E() = 1.0;

          p.E_last() = 1.0;

          p.u() = {1.0, 0.0, 0.0};

          bool found = exhaustive_find_cell(p);

          if (!found) {

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

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

            weight_windows[z * Ny * Nx + y * Nx + x] = 0.0;

            continue;

          }

          int i_cell = p.lowest_coord().cell;

          int64_t sr = source_region_offsets_[i_cell] + p.cell_instance();

          if (RandomRay::mesh_subdivision_enabled_) {

            int mesh_idx = base_source_regions_.mesh(sr);

            int mesh_bin;

            if (mesh_idx == C_NONE) {

              mesh_bin = 0;

            } else {

              mesh_bin = model::meshes[mesh_idx]->get_bin(p.r());

            }

            SourceRegionKey sr_key {sr, mesh_bin};

            auto it = source_region_map_.find(sr_key);

            if (it != source_region_map_.end()) {

              sr = it->second;

            } else {

              sr = -1;

            }

          }

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

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

          if (variance_reduction::weight_windows.size() == 1) {

            WeightWindow ww =

              variance_reduction::weight_windows[0]->get_weight_window(p);

            float weight = ww.lower_weight;

            weight_windows[z * Ny * Nx + y * Nx + x] = weight;

            if (weight < min_weight)

              min_weight = weight;

          }

        }

      }

    }

    double source_normalization_factor =

    // Open file for writing

    // Write vtk metadata

    // Plot multigroup flux data

        }

          }

        }

      }

    }

    // Slightly negative fluxes can be normal when sampling corners of linear

    // source regions. However, very common and high magnitude negative fluxes

    // may indicate numerical instability.

                          "(minumum found = {:.2e} maximum_found = {:.2e})",

    }

    // Plot FSRs

    }

    // Plot Materials

    }

    // Plot fission source

            }

          }

        }

      }

    } else {

            // External sources are already divided by sigma_t, so we need to

            // multiply it back to get the true external source.

            }

          }

        }