OpenMCCellAverageProblem.h

Go to the documentation of this file.

/********************************************************************/
/*                  SOFTWARE COPYRIGHT NOTIFICATION                 */
/*                             Cardinal                             */
/*                                                                  */
/*                  (c) 2021 UChicago Argonne, LLC                  */
/*                        ALL RIGHTS RESERVED                       */
/*                                                                  */
/*                 Prepared by UChicago Argonne, LLC                */
/*               Under Contract No. DE-AC02-06CH11357               */
/*                With the U. S. Department of Energy               */
/*                                                                  */
/*             Prepared by Battelle Energy Alliance, LLC            */
/*               Under Contract No. DE-AC07-05ID14517               */
/*                With the U. S. Department of Energy               */
/*                                                                  */
/*                 See LICENSE for full restrictions                */
/********************************************************************/

#pragma once

#include "OpenMCProblemBase.h"
#include "SymmetryPointGenerator.h"

#include "TallyBase.h"
#include "FilterBase.h"

#ifdef ENABLE_DAGMC
#include "MoabSkinner.h"
#include "DagMC.hpp"
#endif

class OpenMCVolumeCalculation;

class OpenMCCellAverageProblem : public OpenMCProblemBase
{
public:
  OpenMCCellAverageProblem(const InputParameters & params);
  static InputParameters validParams();

  virtual void initialSetup() override;
  virtual void externalSolve() override;
  virtual void syncSolutions(ExternalProblem::Direction direction) override;
  virtual bool converged(unsigned int) override { return true; }

  void read2DBlockParameters(const std::string name,
                             std::vector<std::vector<SubdomainName>> & names,
                             std::vector<SubdomainID> & flattened_ids);

  void checkBlocksInMesh(const std::string name,
                         const std::vector<SubdomainID> & ids,
                         const std::vector<SubdomainName> & names) const;

  void setupProblem();

  virtual void addExternalVariables() override;

  virtual Real cellVolume(const cellInfo & cell_info) const;

  virtual const OpenMCVolumeCalculation * volumeCalculation() const { return _volume_calc; }

  virtual const std::map<cellInfo, std::vector<unsigned int>> & cellToElem() const
  {
    return _cell_to_elem;
  }

  virtual std::unordered_set<SubdomainID> getCellToElementSub(const cellInfo & info)
  {
    return _cell_to_elem_subdomain.at(info);
  }

  virtual bool hasPointTransformations() const { return _symmetry != nullptr; }

  virtual const std::vector<std::string> & getTallyScores() const { return _all_tally_scores; }

  unsigned int getNumScoringTallies(const std::string & score) const
  {
    return _score_count.count(score) > 0 ? _score_count.at(score) : 0;
  }

  bool hasScore(const std::string & score)
  {
    return std::find(_all_tally_scores.begin(), _all_tally_scores.end(), score) !=
           _all_tally_scores.end();
  }

  const TallyBase * getTally(const std::string & name);

  std::vector<const MooseVariableFE<Real> *> getTallyScoreVariables(const std::string & score,
                                                                    const std::string & tally_name,
                                                                    THREAD_ID tid,
                                                                    const std::string & output = "",
                                                                    bool skip_func_exp = false);

  std::vector<const VariableValue *> getTallyScoreVariableValues(const std::string & score,
                                                                 const std::string & tally_name,
                                                                 THREAD_ID tid,
                                                                 const std::string & output = "",
                                                                 bool skip_func_exp = false);

  std::vector<const VariableValue *>
  getTallyScoreNeighborVariableValues(const std::string & score,
                                      const std::string & tally_name,
                                      THREAD_ID tid,
                                      const std::string & output = "",
                                      bool skip_func_exp = false);

  bool hasOutput(const std::string & score, const std::string & output) const;

  virtual Point transformPoint(const Point & pt) const
  {
    return this->hasPointTransformations() ? _symmetry->transformPoint(pt) : pt;
  }

  virtual void createQRules(QuadratureType type,
                            Order order,
                            Order volume_order,
                            Order face_order,
                            SubdomainID block,
                            bool allow_negative_weights = true) override;

  typedef std::unordered_map<int32_t, std::vector<int32_t>> containedCells;

  int32_t elemToCellIndex(const int & elem_id) const { return elemToCellInfo(elem_id).first; }

  int32_t elemToCellID(const int & elem_id) const { return cellID(elemToCellIndex(elem_id)); }

  int32_t elemToCellInstance(const int & elem_id) const { return elemToCellInfo(elem_id).second; }

  cellInfo elemToCellInfo(const int & elem_id) const { return _elem_to_cell[elem_id]; }

  int32_t cellToMaterialIndex(const cellInfo & cell_info) const
  {
    return _cell_to_material.at(cell_info);
  }

  coupling::CouplingFields cellFeedback(const cellInfo & cell_info) const;

  bool hasDensityFeedback(const cellInfo & cell_info) const
  {
    std::vector<coupling::CouplingFields> phase = {coupling::density,
                                                   coupling::density_and_temperature};
    return std::find(phase.begin(), phase.end(), cellFeedback(cell_info)) != phase.end();
  }

  bool hasTemperatureFeedback(const cellInfo & cell_info) const
  {
    std::vector<coupling::CouplingFields> phase = {coupling::temperature,
                                                   coupling::density_and_temperature};
    return std::find(phase.begin(), phase.end(), cellFeedback(cell_info)) != phase.end();
  }

  bool hasFilter(const std::string & filter_name) const { return _filters.count(filter_name) > 0; }

  std::shared_ptr<FilterBase> & getFilter(const std::string & filter_name)
  {
    return _filters.at(filter_name);
  }

  const std::vector<std::shared_ptr<TallyBase>> & getLocalTallies() const { return _local_tallies; }

  double cellTemperature(const cellInfo & cell_info) const;

  double cellMappedVolume(const cellInfo & cell_info) const;

  void reloadDAGMC();

  void addFilter(const std::string & type,
                 const std::string & name,
                 InputParameters & moose_object_pars);

  std::shared_ptr<TallyBase>
  addTally(const std::string & type, const std::string & name, InputParameters & moose_object_pars);

  Real tallyMultiplier(const std::string & score_name, const Real & local_mean_tally) const;

  const Real getReferenceDensity(const Elem * elem) const;

  template <typename T>
  void checkEmptyVector(const std::vector<T> & vector, const std::string & name) const
  {
    if (vector.empty())
      mooseError(name + " cannot be empty!");
  }

  int fixedPointIteration() const { return _fixed_point_iteration; }

  bool hasAdaptivity() const { return _has_adaptivity; }

  bool hasSkinner() const { return _using_skinner; }

  static constexpr int32_t UNMAPPED{-1};

  static constexpr int DIMENSION{3};

  virtual MooseMesh & getMooseMesh();

  virtual const MooseMesh & getMooseMesh() const;

  const bool & useDisplaced() const { return _use_displaced; }

protected:
  void reinitCouplingAndApplyFeedback();

  virtual void critSearchStep() override;

  unsigned int getCellLevel(const Point & c) const;

  void
  readBlockVariables(const std::string & param,
                     const std::string & default_name,
                     std::map<std::string, std::vector<SubdomainName>> & vars_to_specified_blocks,
                     std::vector<SubdomainID> & specified_blocks);

  bool cellHasIdenticalFill(const cellInfo & cell_info) const;

  containedCells shiftCellInstances(const cellInfo & cell_info) const;

  bool cellMapsToSubdomain(const cellInfo & cell_info,
                           const std::unordered_set<SubdomainID> & id) const;

  cellInfo firstContainedMaterialCell(const cellInfo & cell_info) const;

  containedCells containedMaterialCells(const cellInfo & cell_info) const;

  void updateOpenMCGeometry();

  std::vector<std::string> getMaterialInEachSubdomain() const;

  Point transformPointToOpenMC(const Point & pt) const;

  void
  printTrisoHelp(const std::chrono::time_point<std::chrono::high_resolution_clock> & start) const;

  void printAuxVariableIO();

  std::vector<int32_t> materialsInCells(const containedCells & contained_cells) const;

  void subdomainsToMaterials();

  std::set<SubdomainID> coupledSubdomains() const;

  template <typename T>
  void gatherCellSum(std::vector<T> & local, std::map<cellInfo, T> & global) const;

  template <typename T>
  void gatherCellVector(std::vector<T> & local, std::vector<unsigned int> & n_local, std::map<cellInfo, std::vector<T>> & global);

  coupling::CouplingFields elemFeedback(const Elem * elem) const;

  void getTallyTriggerParameters(const InputParameters & params);

  void readBlockParameters(const std::string name, std::unordered_set<SubdomainID> & blocks);

  void cacheContainedCells();

  void setContainedCells(const cellInfo & cell_info,
                         const Point & hint,
                         std::map<cellInfo, containedCells> & map);

  void checkContainedCellsStructure(const cellInfo & cell_info,
                                    containedCells & reference,
                                    containedCells & compare) const;

  void setMinimumVolumeQRules(Order & volume_order, const std::string & type);

  std::string printNewline() const
  {
    if (_verbose)
      return "\n";
    else
      return "";
  }

  void storeElementPhase();

  void getCellMappedPhase();

  void checkCellMappedPhase();

  void getCellMappedSubdomains();

  void computeCellMappedVolumes();

  void initializeElementToCellMapping();

  void mapElemsToCells();

  void validateLocalTallies();

  void initializeTallies();

  void resetTallies();

  void getMaterialFills();

  void getPointInCell();

  std::map<cellInfo, Real> computeVolumeWeightedCellInput(
      const std::map<SubdomainID, std::pair<unsigned int, std::string>> & var_num,
      const std::vector<coupling::CouplingFields> * phase = nullptr,
      const std::map<SubdomainID, Real> * scaling = nullptr) const;

  void sendTemperatureToOpenMC() const;

  void sendDensityToOpenMC() const;

  void latticeOuterCheck(const Point & c, int level) const;

  void latticeOuterError(const Point & c, int level) const;

  bool findCell(const Point & point);

  void compareContainedCells(std::map<cellInfo, containedCells> & reference,
                             std::map<cellInfo, containedCells> & compare) const;

  virtual std::vector<int32_t> getMappedTallyIDs() const override;

  NumericVector<Number> & _serialized_solution;

  const bool & _output_cell_mapping;

  const coupling::OpenMCInitialCondition _initial_condition;

  const relaxation::RelaxationEnum _relaxation;

  const trigger::TallyTriggerTypeEnum _k_trigger;

  unsigned int _cell_level;

  const bool & _export_properties;

  const bool _using_skinner;

  bool _need_to_reinit_coupling;

  const bool _has_identical_cell_fills;

  const bool & _check_identical_cell_fills;

  const bool & _assume_separate_tallies;

  bool _map_density_by_cell;

  const bool _specified_density_feedback;

  const bool _specified_temperature_feedback;

  bool _has_cell_tallies = false;

  bool _needs_to_map_cells;

  std::map<std::string, std::shared_ptr<FilterBase>> _filters;

  std::vector<std::shared_ptr<TallyBase>> _local_tallies;

  std::vector<std::string> _all_tally_scores;

  std::map<std::string, unsigned int> _score_count;

  std::vector<std::vector<unsigned int>> _tally_var_ids;

  std::vector<std::vector<std::vector<unsigned int>>> _tally_ext_var_ids;

  std::vector<SubdomainID> _density_blocks;

  std::vector<SubdomainID> _temp_blocks;

  std::unordered_set<SubdomainID> _identical_cell_fill_blocks;

  std::vector<cellInfo> _elem_to_cell{};

  std::map<cellInfo, coupling::CouplingFields> _cell_phase;

  int _n_moose_density_elems;

  int _n_moose_temp_elems;

  int _n_moose_temp_density_elems;

  int _n_moose_none_elems;

  int _n_mapped_temp_elems;

  int _n_mapped_density_elems;

  int _n_mapped_temp_density_elems;

  int _n_mapped_none_elems;

  Real _uncoupled_volume;

  bool _material_cells_only{true};

  std::map<cellInfo, std::vector<unsigned int>> _cell_to_elem;

  std::map<cellInfo, std::vector<unsigned int>> _local_cell_to_elem;

  std::map<cellInfo, std::unordered_set<SubdomainID>> _cell_to_elem_subdomain;

  std::map<SubdomainID, std::set<int32_t>> _subdomain_to_material;

  std::map<cellInfo, Point> _cell_to_point;

  std::map<cellInfo, Real> _cell_to_elem_volume;

  std::map<cellInfo, Real> _cell_volume;

  std::map<cellInfo, int32_t> _cell_to_material;

  std::map<cellInfo, containedCells> _cell_to_contained_material_cells;

  std::map<cellInfo, int32_t> _cell_to_n_contained;

  static bool _first_transfer;

  static bool _printed_initial;

  static bool _printed_triso_warning;

  openmc::Particle _particle;

  unsigned int _n_particles_1;

  std::map<std::string, std::vector<SubdomainName>> _temp_vars_to_blocks;

  std::map<std::string, std::vector<SubdomainName>> _density_vars_to_blocks;

  OpenMCVolumeCalculation * _volume_calc;

  const SymmetryPointGenerator * _symmetry;

  std::map<cellInfo, int> _n_temp;

  std::map<cellInfo, int> _n_rho;

  std::map<cellInfo, int> _n_temp_rho;

  std::map<cellInfo, int> _n_none;

  std::shared_ptr<TallyBase> _source_rate_norm_tally;

  int _source_rate_score = -1;

#ifdef ENABLE_DAGMC
  MoabSkinner * _skinner = nullptr;

  std::shared_ptr<moab::DagMC> _dagmc = nullptr;
#endif

  long unsigned int _n_openmc_cells;

  int32_t _dagmc_universe_id;

  bool _dagmc_root_universe = true;

  int32_t _cell_using_dagmc_universe_id;

  const int32_t _initial_num_openmc_surfaces;

  static constexpr Real EV_TO_JOULE = 1.6022e-19;

  static constexpr Real ZERO_TALLY_THRESHOLD = 1e-12;

private:
  void dufekGudowskiParticleUpdate();

  std::vector<int32_t> _flattened_ids;

  std::vector<int32_t> _flattened_instances;

  containedCells _instance_offsets;

  std::map<cellInfo, int32_t> _n_offset;

  cellInfo _first_identical_cell;

  std::vector<int32_t> _first_identical_cell_materials;

  bool _use_displaced;

  std::map<SubdomainID, std::pair<unsigned int, std::string>> _subdomain_to_temp_vars;

  std::map<SubdomainID, std::pair<unsigned int, std::string>> _subdomain_to_density_vars;

  std::map<SubdomainID, Real> _subdomain_to_ref_density;
};