Solver.h

Go to the documentation of this file.

#ifndef SOLVER_H_
#define SOLVER_H_

#ifdef __cplusplus
#define _USE_MATH_DEFINES
#ifdef SWIG
#include "Python.h"
#endif
#include "constants.h"
#include "Timer.h"
#include "Quadrature.h"
#include "TrackGenerator3D.h"
#include "Cmfd.h"
#include "Progress.h"
#include "ExpEvaluator.h"
#include "segmentation_type.h"
#include <math.h>
#endif

#define _scalar_flux(r,e) (_scalar_flux[(r)*_num_groups + (e)])

#define _old_scalar_flux(r,e) (_old_scalar_flux[(r)*_num_groups + (e)])

#define _reference_flux(r,e) (_reference_flux[(r)*_num_groups + (e)])

#define reference_flux(r,e) (reference_flux[(r)*_num_groups + (e)])

#define _stabilizing_flux(r,e) (_stabilizing_flux[(r)*_num_groups + (e)])

#define _reduced_sources(r,e) (_reduced_sources[(r)*_num_groups + (e)])

#define _boundary_flux(i,j,pe) (_boundary_flux[(i)*2*_fluxes_per_track \
                                                + (j)*_fluxes_per_track \
                                               + (pe)])
#define _start_flux(i,j,pe) (_start_flux[(i)*2*_fluxes_per_track \
                                                + (j)*_fluxes_per_track \
                                               + (pe)])

#define _fixed_sources(r,e) (_fixed_sources[(r)*_num_groups + (e)])

#define fission_sources(r,e) (fission_sources[(r)*_num_groups + (e)])

#define scatter_sources(r,e) (scatter_sources[(r)*_num_groups + (e)])


enum solverMode {

  FORWARD,

  ADJOINT,
};


enum residualType {

  SCALAR_FLUX,

  FISSION_SOURCE,

  TOTAL_SOURCE,
};


enum stabilizationType {

  DIAGONAL,

  YAMAMOTO,

  GLOBAL,
};


class Solver {

protected:

  int _num_azim;

  int _num_groups;

  long _num_FSRs;

  long _num_fissionable_FSRs;

  FP_PRECISION* _FSR_volumes;

  Material** _FSR_materials;

  Material* _chi_spectrum_material;

  TrackGenerator* _track_generator;

  Geometry* _geometry;

  int _num_materials;

  Quadrature* _quad;

  int _num_polar;

  int _fluxes_per_track;

  Track** _tracks; //FIXME

  std::vector<int> _limit_xs_materials;

  int _reset_iteration;

  bool _limit_xs;

  std::map<int, Material*> _original_materials;

  std::map<int, Material*> _limit_materials;

  int*** _tracks_per_stack;

#ifndef THREED

  bool _SOLVE_3D;
#endif

  solverMode _solver_mode;

  bool _is_restart;

  bool _user_fluxes;

  bool _fixed_sources_on;

  bool _correct_xs;

  bool _stabilize_transport;

  bool _verbose;

  bool _calculate_initial_spectrum;

  double _initial_spectrum_thresh;

  bool _load_initial_FSR_fluxes;

  bool _calculate_residuals_by_reference;

  bool _negative_fluxes_allowed;

  std::string _initial_FSR_fluxes_file;

  std::string _reference_file;

  logLevel _xs_log_level;

  segmentationType _segment_formation;

  long _tot_num_tracks;

  float* _boundary_flux;
  float* _start_flux;

  float* _boundary_leakage;

  FP_PRECISION* _scalar_flux;

  FP_PRECISION* _old_scalar_flux;

  FP_PRECISION* _reference_flux;

  FP_PRECISION* _stabilizing_flux;

  FP_PRECISION* _fixed_sources;

  std::vector<FP_PRECISION*> _groupwise_scratch;
  double* _regionwise_scratch;

  std::map< std::pair<int, int>, FP_PRECISION > _fix_src_FSR_map;

  std::map< std::pair<Cell*, int>, FP_PRECISION > _fix_src_cell_map;

  std::map< std::pair<Material*, int>, FP_PRECISION > _fix_src_material_map;

  FP_PRECISION* _reduced_sources;

  double _k_eff;

  bool _keff_from_fission_rates;

  int _num_iterations;

  double _converge_thresh;

  double _stabilization_factor;

  stabilizationType _stabilization_type;

  ExpEvaluator*** _exp_evaluators;

  int _num_exp_evaluators_azim;

  int _num_exp_evaluators_polar;

  Timer* _timer;

  Cmfd* _cmfd;

  std::string _source_type;

  virtual void initializeFluxArrays() =0;

  virtual void initializeSourceArrays() =0;

  /* Initialize interp. tables, constants for computing source exponentials */
  virtual void initializeExpEvaluators();

  /* Build fission matrices, transpose production for adjoint calculations */
  void initializeMaterials(solverMode mode);

  virtual void initializeFSRs();
  void countFissionableFSRs();
  void checkXS();
  virtual void initializeCmfd();
  void calculateInitialSpectrum(double threshold);

  virtual void zeroTrackFluxes() =0;

  virtual void flattenFSRFluxes(FP_PRECISION value) =0;

  virtual void flattenFSRFluxesChiSpectrum() =0;

  virtual void storeFSRFluxes() =0;

  virtual double normalizeFluxes() =0;

  virtual void computeStabilizingFlux() =0;

  virtual void stabilizeFlux() =0;

  virtual void computeFSRSources(int iteration) =0;

  virtual void computeFSRFissionSources() =0;

  virtual void computeFSRScatterSources() =0;

  virtual double computeResidual(residualType res_type) =0;

  virtual void computeKeff() =0;

  virtual void addSourceToScalarFlux() =0;

  virtual void transportSweep() =0;

  void clearTimerSplits();

  /* Whether to ray-trace and propagate fluxes at the same time */
  bool _OTF_transport;

public:
  Solver(TrackGenerator* track_generator=NULL);
  virtual ~Solver();

  void setGeometry(Geometry* geometry);

  Geometry* getGeometry();
  TrackGenerator* getTrackGenerator();
  FP_PRECISION getFSRVolume(long fsr_id);
  int getNumPolarAngles();
  int getNumIterations();
  double getTotalTime();
  double getKeff();
  double getConvergenceThreshold();
  FP_PRECISION getMaxOpticalLength();
  bool isUsingDoublePrecision();
  bool isUsingExponentialInterpolation();
  bool is3D();

  void initializeSolver(solverMode solver_mode);
  virtual void initializeFixedSources();

  void printFissionRates(std::string fname, int nx, int ny, int nz);
  virtual void printInputParamsSummary();

  void setKeffFromNeutronBalance();
  void setResidualByReference(std::string fname);
  void dumpFSRFluxes(std::string fname);
  void loadInitialFSRFluxes(std::string fname);
  void loadFSRFluxes(std::string fname, bool assign_k_eff=false, double tolerance=0.01);

  double getFlux(long fsr_id, int group);
  virtual void getFluxes(FP_PRECISION* out_fluxes, int num_fluxes) = 0;
  double getFSRSource(long fsr_id, int group);

  void setTrackGenerator(TrackGenerator* track_generator);
  void setConvergenceThreshold(double threshold);
  virtual void setFluxes(FP_PRECISION* in_fluxes, int num_fluxes) = 0;
  virtual void setFixedSourceByFSR(long fsr_id, int group, double source);
  void setFixedSourceByCell(Cell* cell, int group, double source);
  void setFixedSourceByMaterial(Material* material, int group,
                                double source);
  void setMaxOpticalLength(FP_PRECISION max_optical_length);
  void setExpPrecision(double precision);
  void useExponentialInterpolation();
  void useExponentialIntrinsic();
  void setSolverMode(solverMode solver_mode);
  void setRestartStatus(bool is_restart);
  void allowNegativeFluxes(bool negative_fluxes_on);
  void correctXS();
  void stabilizeTransport(double stabilization_factor,
                          stabilizationType stabilization_type=DIAGONAL);
  void setInitialSpectrumCalculation(double threshold);
  void setCheckXSLogLevel(logLevel log_level);
  void setChiSpectrumMaterial(Material* material);
  void resetMaterials(solverMode mode);

  void fissionTransportSweep();
  void scatterTransportSweep();
  void computeFlux(int max_iters=1000, bool only_fixed_source=true);
  void computeSource(int max_iters=1000, double k_eff=1.0,
                     residualType res_type=TOTAL_SOURCE);
  void computeEigenvalue(int max_iters=1000,
                         residualType res_type=FISSION_SOURCE);

  void printBGQMemory();

  virtual void computeFSRFissionRates(double* fission_rates, long num_FSRs,
                                      bool nu = false) =0;

  inline float* getBoundaryFlux(long track_id, bool fwd) {
    return &_boundary_flux(track_id, !fwd, 0);
  }

  void setVerboseIterationReport();
  void printTimerReport();
  FP_PRECISION* getFluxesArray();

  /* Functions to limit cross sections, to attempt to stabilize MOC */
  void limitXS();
  void setLimitingXSMaterials(std::vector<int> material_ids,
                              int reset_iteration);
  void checkLimitXS(int iteration);


  inline void setOTFTransport() {
    _OTF_transport = true;
    log_printf(NORMAL, "Using On-The-Fly transport");
  }
};


#endif /* SOLVER_H_ */