SM_StochasticEnergy.h

#ifndef SM_StochasticEnergy_H
#define SM_StochasticEnergy_H
 
#include "SM_StochasticOutput.h"
 
#include "SM_Beam.h"
#include "SM_LandauLifschitzSystem.h"
#include "CORE_IO.h"
 
class SM_StochasticEnergy  : public SM_StochasticOutput<SM_StochasticEnergy>  {
  
    // ATTRIBUTES
  
public:
   
  
   
private:
 
 
    //type class
    typedef SM_StochasticEnergy SelfClass;
    typedef SM_StochasticOutput<SelfClass> SuperClass;
    
    //number of time steps
    tIndex mTimeStepsNumber;
 
    //number of energies depening on number of operators
    tUCInt mEnergiesNumber;
    
  
    tIndex mEpsilonIndex;//iterator on mEpsilons
    
    std::valarray<tReal> mTimes;//values of times of size mTimeStepsNumber
    std::valarray<tReal> mEpsilons;//Epsilon at t of size mTimeStepsNumber;
    
    std::valarray<tReal> mEnergies;//Energies at t by operator of size mTimeStepsNumber x mEnergiesNumber;
    std::valarray<tReal> mEnergiesAtS;//Energies at simulation s of size mTimeStepsNumber x mEnergiesNumber; 
    std::valarray<tReal> mEs;//energy of the system at t
    
    std::ofstream *mFile;
public:
    // METHODS
  
    // CONSTRUCTORS 
  
    SM_StochasticEnergy(void) {
        mTimeStepsNumber=0;
        mEpsilonIndex=0;
        mFile=NULL;
    }
 
  
  
    // DESTRUCTORS 
public:
    
    virtual ~SM_StochasticEnergy(void) {
    }
 
 
 
public:
    //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
    
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+= mEnergies.size()*sizeof(tReal);
        mem+= mEpsilons.size()*sizeof(tReal);
        mem+= mTimes.size()*sizeof(tReal);
        mem+= mEnergiesAtS.size()*sizeof(tReal);
        mem+= mEs.size()*sizeof(tReal);
        return mem;
    }
    
    inline static CORE_UniquePointer<SM_StochasticEnergy> New() {
        return CORE_UniquePointer<SM_StochasticEnergy>(new SM_StochasticEnergy(),
                                                       SM_StochasticEnergy::Delete());
    }
 
    //SET & GET Methods
    //=================
 
 
    
    virtual void adimensionize(const SM_Material& material) override {
    }
 
    inline void initialize(const SM_Beam& beam) {
 
        
        mEpsilonIndex=0;
        
        //energies number
        mEnergiesNumber=beam.getSystem().getOperatorsNumber()+1;
        
        
        mTimeStepsNumber=beam.getStepsNumber()+1;
        mTimes.resize(mTimeStepsNumber);
        mTimes=0;
        
        
        //epsilon(t)
        mEpsilons.resize(mTimeStepsNumber);
        
        //E(t)
        mEnergies.resize(mTimeStepsNumber*mEnergiesNumber);
        //E(t) at simulation s
        mEnergiesAtS.resize(mTimeStepsNumber*mEnergiesNumber);
        
        //E enregy of the system
        mEs.resize(mEnergiesNumber);
        
    }
 
    //implemented templated methods
    //=============================
    
    tBoolean open(const SM_Beam& beam) {
 
        //verify the system is LL 
        const SM_LandauLifschitzSystem *llSystem=dynamic_cast< const SM_LandauLifschitzSystem *>(&beam.getSystem());
        if (llSystem==null) return false;
 
        if (getIndex()==0) {
            
            //initialize the data from beam
            initialize(beam);
            
            //number of raws values number is null
            tIndex &nRaws=getPackedRawValuesNumber();
            nRaws=0;
        
            tString fn=CORE_IO::GetAbsolutePath(getOutputPath()+"/"+getPrefix()+".me").c_str();
            if (mFile!=NULL) {
                mFile->close();
                delete mFile;
                mFile=NULL;
            }
            if ((mEpsilonIndex==0)&&(llSystem!=null)) {
                
                mFile=new std::ofstream(fn,std::ios::out);
                if ((mFile==null) || !(*mFile)) {
                    std::cout<<"Fatal error "<<fn<<" can not be written \n";
                    return false;
                }
                
                tString str=beam.toString();
                functions_string::replaceAll("\n","\n#",str);
                (*mFile)<<"#"<<str<<"\n";
                (*mFile)<<"#time epsilon ";                 
                (*mFile)<<"\t Etotal";
                std::vector<tString> names;
                llSystem->getOperatorsName(names);
                for(const auto& name:names) (*mFile)<<"\t E"<<name;
                (*mFile)<<"\n";
                  
            }
        }
       
        //reset the energies to 0
        mEnergies=0;
        return true;
    }
 
      
    inline tBoolean open(const tIndex& s,const SM_System& system) {
        mEnergiesAtS=0;//energies at simulation s of size mEnergiesNumber x mTimesSteps
        return true;
    }
 
    
    tBoolean  store(const SM_System& system) {
        
        const SM_LandauLifschitzSystem *llSystem=dynamic_cast< const SM_LandauLifschitzSystem *>(&system);
        //compute the energies of the systems
        llSystem->computeEnergies(llSystem->getTimeIndex(),mEs);
 
        if (getIndex()==0) {
            //time 
            const tReal& time=llSystem->getTime();
            
            //epsilon
            const tReal& epsilon_t=llSystem->getNoiseRate();
            
            if (mTimeStepsNumber==1) {
                //time step
                tReal dt=llSystem->getTimeStepper().getTimeStep()/2.;
                if (fabs(mTimes[0]-time)<dt) {
                    
                    const tReal *Es=&mEs[0];
                    //register energies
                    tReal *EAtS=&mEnergiesAtS[0];
                    for (tUCInt k=0;k<mEnergiesNumber;k++) {
                        (*EAtS)+=(*Es);
                        Es++;
                        EAtS++;
                    }
                    //register epsilon
                    mEpsilons[0]=epsilon_t;
                }
                
            } else {
                //save energies at time
                
                //compute the energies of the systems
                
                const tReal *Es=&mEs[0];
                
                //register the time
                const tIndex& timeIndex=llSystem->getTimeIndex();
                mTimes[timeIndex]=time;
                //register epsilon
                mEpsilons[timeIndex]=epsilon_t;
                
                //register the energies
                tReal *EAtS=&mEnergiesAtS[timeIndex*mEnergiesNumber];
                for (tUCInt k=0;k<mEnergiesNumber;k++) {
                    (*EAtS)+=(*Es);
                    Es++;
                    EAtS++;
                }
            }
        }
        
        return true;
    }
    
  
  
    inline tBoolean close(const tIndex& s,const SM_System& system,const tBoolean& hasSucceeded) {
        
        if (hasSucceeded) {
            if (getIndex()==0) {
                getPackedRawValuesNumber()++;
                
                //register the data for simlation s
                tReal *Es=&mEnergies[0];
                const tReal *EAtS=&mEnergiesAtS[0];
                const tReal *EAtSe=EAtS+mEnergiesAtS.size();
                while (EAtS!=EAtSe) {
                (*Es)+=(*EAtS);
                Es++;
                EAtS++;
                }
            }
        }
        return true;
    }
 
    inline tBoolean close(const SM_Beam& beam) {
        if (getIndex()==0) {
            const tIndex &nRaws=getPackedRawValuesNumber();
            //E(t)=\frac{1}{S} \sum_{s=0}^{s=S-1} E(t,epsilon_t)
            tReal *Es=&mEnergies[0];
            const tReal *Ese=Es+mEnergies.size();
            while (Es!=Ese) {
                (*Es)/=nRaws;
                Es++;
            }
            
            //save the file
            if (mFile!=NULL) {
                const tReal *iEpsilon=&mEpsilons[0];
                const tReal *iEnergies=&mEnergies[0];
                const tReal *iT=&mTimes[0];
                tIndex t;
                tUCInt k;
                (*mFile) <<std::setprecision(12);
                for (t=0;t<mTimeStepsNumber;t++) {
                    (*mFile)<<(*iT);
                    iT++;
                    
                    (*mFile)<<"\t"<<(*iEpsilon);
                    iEpsilon++;
                    
                    for (k=0;k<mEnergiesNumber;k++) {
                        (*mFile)<<"\t"<<(*iEnergies);
                        iEnergies++;
                    }
                    (*mFile)<<"\n";
                    
                    
                }
                
                mFile->close();
                delete mFile;
                mFile=NULL;
            }
            
            //next epsilon index
            mEpsilonIndex++;
        }
        
        return true;
    }
 
    
  
    
    virtual void closePackedSimulations(const SM_Beam& beam,
                                        const tIndex& rawValuesNumber,
                                        const std::valarray<tReal>& rawValues) override {
       
        
    }
 
    virtual tString toString() const override {
        std::stringstream ret;
        ret<<SuperClass::toString();
        return ret.str();
    }
    
public:
 
   
public:
    // OTHERS methods
  
    
 
 
};
 
#endif