EvolutionSolver.hxx

Go to the documentation of this file.

/*
  This file is part of MURE,
  Copyright (C) 2007-2021 MURE developers.
 
  MURE is free software: you can redistribute it and/or modify
  it under the terms of the GNU Lesser General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
 
  MURE is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU Lesser General Public License for more details.
 
  You should have received a copy of the GNU Lesser General Public License
  along with MURE.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef _EvolutionSolver_HXX_
#define _EvolutionSolver_HXX_
 
#include <nlohmann/json.hpp>
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>
#include <string>
#include <algorithm>
#include <map>
 
#include "MUREGlobal.hxx"
#include "Cell.hxx"
#include "EvolutionControl.hxx"
#include "EvolutiveSystem.hxx"
#include "ControlMaterial.hxx"
#include "Nucleus.hxx"
#include "Reaction.hxx"
#include "libValErr/ValErr.hxx"
 
class ControlMaterial;
class EvolutionControl;
class EvolutiveSystem;
class Material;
class ZAI;
 
using json = nlohmann::ordered_json; // preserve order of insertion
 
using namespace std;
 
 
 
class EvolutionSolver
{
  public:
    EvolutionSolver(); 
    EvolutionSolver(const EvolutionSolver &s); 
    ~EvolutionSolver(); 
 
 
    void Evolution(
        vector < double > T,       
        int Start = 0,             
        bool BrutalStop = false    
    );
 
 
    void UseMultiGroupTallies(bool StdTallyFor238U = false);
 
    bool IsU8StdTally()
    {
        return fU8StdTally;   
    }
    void SetMultiGroupDecadeMultiplicator(int N = 1)
    {
        fMultiGroupDecadeMultiplicator = N;   
    }
    int     GetMultiGroupDecadeMultiplicator()
    {
        return fMultiGroupDecadeMultiplicator;   
    }
    void SetMultiGroupLowerEnergy(double Emin = 1e-4)
    {
        fMultiGroupEmin = Emin; 
    }
    void SetMultiGroupUpperEnergy(double Emax = 20e6)
    {
        fMultiGroupEmax = Emax; 
    }
    double GetMultiGroupLowerEnergy()
    {
        return fMultiGroupEmin;   
    }
    double GetMultiGroupUpperEnergy()
    {
        return fMultiGroupEmax;   
    }
    bool IsMultiGroupTalliesUsed()
    {
        return fMultigroupTallies;   
    }
 
    void SetXe135Equilibrium(double t = 3 * 9.14 * 3600);
    bool IsXe135Equilibrium()
    {
        return fXe135Equilibrium;   
    }
    double GetXe135EquilibriumTime()
    {
        return fXe135EquilibriumTime;   
    }
 
 
    void SetCooling(vector < bool > C);
    void SetCooling(int i, bool flag){fIsCooling[i] = flag ;}
    
    bool IsCooling(int MCRunNumber);    
    void FinishOnMCRun(bool flag = true)
    {
        fFinishOnMCRun = flag;   
    }
 
    void AddCellToGroup(Cell *TheCell, int GroupNumber = 0);
    vector < Cell *>  &GetCellGroup(int GroupNumber)
    {
        return fGroupCellVector[GroupNumber];   
    }
    int NumberOfCellGroup()
    {
        return fGroupCellVector.size();   
    }
 
 
    void SetPower(double Power);    
    double GetPower()
    {
        return fPower;   
    }
    vector < double > GetPowerVector()
    {
        return fPowerVector;
    }
    void UpdateBurnUp();
    double GetBurnUp();     
 
    void VerifyPower();
 
 
    void SetBrutalStopKValue(double brutalstopvalue)
    {
        fBrutalStopKValue = brutalstopvalue;
        fBrutalStop = true;
    }
 
 
    void SetEvolutionControl(EvolutionControl *EC);
 
    EvolutionControl *GetEvolutionControl()
    {
        return fEvolutionControl;   
    }
 
    void SetMaterialControl()
    {
        fMaterialControl = true;    
    }
    bool IsMaterialControl()
    {
        return fMaterialControl;    
    }
    ControlMaterial *GetControlMaterial();    
 
    void SetControlRate(Reaction r, ValErr_t rate)
    {
        fControlRate[r] = rate;
    }
 
 
    void SetControlRate(int ReactionCode, ValErr_t rate)
    {
        fControlRate[Reaction(ReactionCode)] = rate;
    }
    ValErr_t GetControlRate(Reaction r)
    {
        return fControlRate[r];   
    }
    ValErr_t GetControlRate(int ReactionCode)
    {
        return fControlRate[Reaction(ReactionCode)];   
    }
 
    void AddControlRodCell(Cell *TheCell)
    {
        fControlRodCell.push_back(TheCell);   
    }
    vector < Cell *> &GetControlRodCell()
    {
        return fControlRodCell;   
    }
 
    double GetInitialHNMass()
    {
        if ( ! fInitHNMassInTon)
        {
            CalculateInitialHNMass();   
        }
        return fInitHNMassInTon;
    }
    void CalculateInitialHNMass();    
 
    ValErr_t GetAbs()
    {
        return fGlobalAbsorptions;   
    }
    ValErr_t GetN2N()
    {
        return fGlobalN2N;   
    }
    ValErr_t GetN3N()
    {
        return fGlobalN3N;   
    }
    ValErr_t GetFiss()
    {
        return fGlobalFissions;   
    }
    ValErr_t GetNuFiss()
    {
        return fGlobalNuFissions;   
    }
    ValErr_t GetNu()
    {
        return fGlobalNu;   
    }
    ValErr_t GetNeutronLosses()
    {
        return fGlobalNeutronLosses;   
    }
    double GetAvgQ()
    {
        return fAvgQ;   
    }
    void SetNeutronLosses(ValErr_t losses)
    {
        fGlobalNeutronLosses = losses;   
    }
    void SetGlobalLossesTallyNumber(int num)
    {
        fGlobalLossesTallyNumber = num;   
    }
    int GetGlobalLossesTallyNumber()
    {
        return fGlobalLossesTallyNumber;   
    }
    void SetGlobalRates(); 
    void IntegrateMassInEvolvingCells();    
    double GetTotalAtoms();
 
 
 
    void SetDepletionSolver(string s = "CRAM")
    {
        fSolverMethod = s;   
    }
    string GetDepletionSolver()
    {
        return fSolverMethod;
    }
 
 
    void SetCRAMOrder(int order)
    {
        fCRAMOrder = order;
    }
    int GetCRAMOrder()
    {
        return fCRAMOrder;   
    }
 
 
    void SetCRAMMethod(string method)
    {
        fCRAMMethod = method;
    }
    string GetCRAMMethod()
    {
        return fCRAMMethod;   
    }
 
 
 
    void SetNSubSteps(int NSubStep = 10)
    {
        fNSubStep = NSubStep;
    }
 
 
    void SetSubDT(double dt = - 1)
    {
        fSubDT = dt;
    }
 
 
    void SetMaxSubDT(double dtmax = - 1)
    {
        fMaxSubDT = dtmax;
    }
 
    int GetNSubSteps()
    {
        return fNSubStep;   
    }
 
    int FindLastMCRunNumber();  
    void SetCurrentMCTime(double time)
    {
        fCurrentMCTime = time;  
    }
    double GetCurrentMCTime()
    {
        return fCurrentMCTime;  
    }
    void SetNextMCTime(double time)
    {
        fNextMCTime = time;     
    }
    double GetNextMCTime()
    {
        return fNextMCTime;     
    }
    void SetCurrentTime(double time)
    {
        fCurrentTime = time;    
    }
    double GetCurrentTime()
    {
        return fCurrentTime;    
    }
    double GetCurrentSubStep()
    {
        return fCurrentSubStep; 
    }
    int GetCurrentSubStepNumber()
    {
        return fCurrentSubStepNumber;   
    }
    int GetMCRunNumber()
    {
        return fMCRunNumber;    
    }
 
    double GetTimeUnitConverter()
    {
        return fTimeUnitConverter;  
    }
    string GetTimeUnit()
    {
        return fTimeUnit;   
    }
    void SetTimeUnit(double t); 
 
 
 
 
    void SetPredictorCorrector(bool middle_step = false);
 
    bool IsPredictorCorrector()
    {
        return fPredictorCorrector; 
    }
    bool IsPredictorStep()
    {
        return fPredictorStep;      
    }
    void SetPredictorStep(bool flag)
    {
        fPredictorStep = flag;      
    }
    double GetPCEvaluationTime()
    {
        return fPCEvaluationTime;   
    }
 
 
    void FindPredictorCorrectorStep();
 
 
 
    void SetPartialMCRun(bool flag = true)
    {
        fPartialMCRun = flag;
    }
    bool IsPartialMCRun()
    {
        return fPartialMCRun;   
    }
    bool IsPartialParticleMCRun()
    {
        return fPartialParticleMCRun;   
    }
    void SetPartialParticleMCRun(bool flag = true)
    {
        fPartialParticleMCRun = flag;
    }
    void SetPartialMCRunDirectory(string name = "PartialRun")
    {
        fPartialMCRunDir = name;    
    }
    string GetPartialMCRunDirectory()
    {
        return fPartialMCRunDir;    
    }
    void SetSourcePartialFactor(int N = 2)
    {
        fSourcePartialFactor = N;   
    }
    int GetSourcePartialFactor()
    {
        return fSourcePartialFactor;    
    }
 
 
    void SetWriteBinaryData(bool flag = true)
    {
        fWriteBinaryData = flag;    
    }
    void SetWriteASCIIData(bool flag = true)
    {
        fWriteASCIIData = flag;     
    }
    void SetWriteJsonFormat(bool flag = true)
    {
        fWrite2Json = flag;   
    }
    bool IsWriteJsonData()
    {
        return fWrite2Json;   
    }
    bool IsWriteBinaryData()
    {
        return fWriteBinaryData;    
    }
    bool IsWriteASCIIData()
    {
        return fWriteASCIIData;     
    }
    void WriteKeffData();       
    void WriteNFissionsData();  
    void ReadEvolutionData(string filename);    
    void ReadEvolutionPowerData(string filename, int version = 1);  
    void WriteEvolutionData(string filename);   
    void ReadBinaryEvolutionData(string filename);  
    void ReadBinaryEvolutionPowerData(string filename, int version = 1);    
    void WriteBinaryEvolutionData(string filename); 
 
    void BuildPerturbativeMaterials();  
 
    void BuildTallies();
 
    void BuildTalliesMultiGroup();
    void BuildTHTalliesMultiGroup();
 
    void UpdateSigmaPhi();
 
    void UpdateSigmaPhiMultiGroup();
    void UpdateTHSigmaPhiMultiGroup();
 
    /*
     Only for "pertubative material" automatically built, the flux is normalized using the TallyNormalizationFactor.
    */
    void UpdateSigma();
 
    void RebuildPerturbativeMaterials();
 
    void RebuildTallies();
    void TransMogrify();            
    void PredictorTransMogrify();   
 
 
    void FitSigmaPhi(bool flag = true)
    {
        fFitSigmaPhi = flag;
    }
 
 
    void FitTallyNormalizationFactor();
    bool IsFitTallyNormalizationFactor()
    {
        return fFittedTNF;  
    }
    void SetFitTallyNormalizationFactor()
    {
        fFittedTNF = true;  
    }
    bool IsFitSigmaPhi()
    {
        return fFitSigmaPhi;
    }
    int GetFitRangeNumber()
    {
        return fFitRange;   
    }
    void SetFitRangeNumber(int n = 4)
    {
        fFitRange = n;  
    }
 
 
    void SmoothSigmaPhi(bool flag = true)
    {
        fSmoothSigmaPhi = flag;
    }
 
    bool IsSmoothSigmaPhi()
    {
        return fSmoothSigmaPhi;
    }
    bool IsFitSigmaPhiPC()
    {
        return fFitSigmaPhiPC;  
    }
    void SetFitSigmaPhiPC(bool flag = true)
    {
        fFitSigmaPhiPC = flag;  
    }
 
 
 
    vector < EvolutiveSystem *> &GetEvolutiveSystemVector()
    {
        return fEvolutiveSystemVector;   
    }
    void AddEvolutiveSystem(EvolutiveSystem *theEvolutiveSystem);    
    vector < EvolutiveSystem *> &GetOutCoreEvolutiveSystemVector()
    {
        return fOutCoreEvolutiveSystemVector;   
    }
    void AddOutCoreEvolutiveSystem(EvolutiveSystem *theEvolutiveSystem);
 
    vector < Material *> &GetPerturbativeMaterial()
    {
        return fPerturbativeMaterial;   
    }
    Material *FindPerturbativeMaterial(Nucleus_ptr TheNucleus); 
    Material *FindPerturbativeMaterial(Material *TheCellMat);   
 
    void AddZAIThatFissions(ZAI *zai);  
    vector < ZAI *> &GetZAIThatFission()
    {
        return fZAIThatFission;   
    }
 
    void UpdateTallyNormalizationFactor();  
    double GetTallyNormalizationFactor();   
 
    void SetNuSigmaFisPhiTotal(double value)
    {
        fNuSigmaFisPhiTotal = value;    
    }
    double GetNuSigmaFisPhiTotal()
    {
        return fNuSigmaFisPhiTotal;     
    }
    void AddToNuSigmaFisPhiTotal(double value)
    {
        fNuSigmaFisPhiTotal += value;   
    }
 
 
    int GetMCFissionsPerSecondSize()
    {
        return fMCFissionsPerSecond.size();
    }
 
    double GetMCFissionsPerSecond(int i)
    {
        return fMCFissionsPerSecond[i];
    }
 
    void SetMCFissionsPerSourceNeutron(int i, double val)
    {
        fMCFissionsPerSecond[i] = val;
    }
 
    void AddToMCFissionsPerSourceNeutron(int i, double val)
    {
        fMCFissionsPerSecond[i] += val;
    }
 
 
    void SetTallyNormalizationFactor(double val);
    void SetPowerNormalizationFactorNuOverKeff(bool flag = true)
    {
        fPowerNormalizationFactorNuOverKeff = flag;
    }
    bool GetPowerNormalizationFactorNuOverKeff()
    {
        return fPowerNormalizationFactorNuOverKeff;
    }
 
    void SetCorrectNormalizationFactor(bool flag = true)
    {
        fCorrectNormalizationFactor = flag;
    }
    bool GetCorrectNormalizationFactor()
    {
        return fCorrectNormalizationFactor;   
    }
    void ReplaceVector();
    
    
    int GetDataWritingVersion(){return fMureDataWritingVersion;}
    
 
 
  protected:
    void InitEvolution(vector < double > &T, int &Start);
 
    double NormalizeMCFissionsPerSecond(double fractional_power = 1.);
 
    void CheckCompositions();   
    void UpdateMaterial();      
 
    void InitTNFFit();  
    void FindExtrapolateSigmaPhi(); 
 
    void validate_json_against_schema(json &data, const json &schema);
 
    vector < EvolutiveSystem *>  fEvolutiveSystemVector;        
    vector < EvolutiveSystem *>  fOutCoreEvolutiveSystemVector;
    vector < Material *>    fPerturbativeMaterial;  
    vector < ZAI *> fZAIThatFission;    
    vector < Cell *> fControlRodCell;   
    vector < vector < Cell *> > fGroupCellVector;   
    EvolutionControl    *fEvolutionControl;         
 
    double  fPower;                 
    bool    fIsCoolingUsed;         
    vector  < bool > fIsCooling;    
    vector  < double > fPowerVector;
    bool    fBrutalStop;            
    double  fBrutalStopKValue;      
 
    bool    fXe135Equilibrium;      
    double  fXe135EquilibriumTime;  
 
    bool    fFromRebuildTallies;    
    bool    fTallyRebuild;          
    bool    fMultigroupTallies;     
    int     fMultiGroupDecadeMultiplicator; 
    double  fMultiGroupEmin;        
    double  fMultiGroupEmax;        
    bool    fU8StdTally;            
    bool    fFinishOnMCRun;         
    double  fTallyNormalizationFactor;  
    bool    fCorrectNormalizationFactor;
    bool    fPowerNormalizationFactorNuOverKeff;
 
    int     fMCRunNumber;       
    int     fCurrentSubStepNumber;  
    double  fCurrentSubStep;        
    double  fCurrentMCTime; 
    double  fNextMCTime;        
    double  fCurrentTime;       
    int     fNSubStep;              
    double  fSubDT;             
    double  fMaxSubDT;          
 
    double  fDeltaK;            
    double  fConstantK;         
    bool    fMaterialControl;   
    double  fInitHNMassInTon;   
    bool    fHNMassInTonInitialized;    
    double fBurnup;
    double fAvgQ;
 
    map < Reaction, ValErr_t, EqualReaction >   fControlRate ;
    ValErr_t fGlobalAbsorptions;        
    ValErr_t fGlobalFissions;           
    ValErr_t fGlobalNuFissions;         
    ValErr_t fGlobalNu;                 
    ValErr_t fGlobalN2N;                
    ValErr_t fGlobalN3N;                
    ValErr_t fGlobalNeutronLosses;      
    double fNuSigmaFisPhiTotal;         
    vector < double > fMCFissionsPerSecond; 
 
    int fGlobalLossesTallyNumber;       
 
    double  fTimeUnitConverter; 
    string  fTimeUnit;          
 
    bool    fFitSigmaPhi;       
    bool    fSmoothSigmaPhi;    
    int     fFitRange;          
    vector < double >   fFitTNF;        
    vector < double >   fMCTimeForFit;
    double  fFitTNFSlope;           
    double  fInterseptTNF;          
    bool    fFittedTNF;             
 
    bool fPredictorCorrector;   
    double fPCEvaluationTime;   
    bool fPredictorStep;        
    bool fFitSigmaPhiPC;        
 
    bool fPartialMCRun;         
    bool fPartialParticleMCRun; 
    string fPartialMCRunDir;    
    int fSourcePartialFactor;   
 
    bool    fWriteBinaryData;       
    bool    fWriteASCIIData;        
    int fMureDataWritingVersion;    
 
    string fSolverMethod;
    int fCRAMOrder;
    string fCRAMMethod;
 
    bool fWrite2Json; 
    static const json EvolDataSchema;
};
 
#endif