SM_StochasticTrajectory.h

#ifndef SM_StochasticTrajectory_H
#define SM_StochasticTrajectory_H
 
//super class header
#include "SM_StochasticOutput.h"
 
 
//Beam
#include "SM_Beam.h"
 
//LL System to save
#include "SM_LandauLifschitzSystem.h"
 
//IO header
#include "CORE_IO.h"
 
class SM_StochasticTrajectory  : public SM_StochasticOutput< SM_StochasticTrajectory>  {
  
    // ATTRIBUTES
  
public:
   
  
   
private:
 
    //type class
    typedef SM_StochasticTrajectory SelfClass;
    typedef SM_StochasticOutput<SelfClass> SuperClass;
 
    //particles selection to save the trajectory
    std::vector<tIndex> mParticles;
    
    //simulations selection to save the trajectory of each point stored in mParticles
    std::vector<tIndex> mSimulations;
 
    //iterator on simulation indices
    std::vector<tIndex>::const_iterator mISimulation;
 
    
    //trajectory file 
    std::ofstream *mFile;
 
    //energies of the system
    std::valarray<tReal> mEnergies;
    
    //energy variation
    tReal mT_old;
    tReal mE_old;
    
public:
    // METHODS
  
    // CONSTRUCTORS 
  
    SM_StochasticTrajectory(void) {
        mFile=NULL;
        mParticles.resize(1);
        mParticles[0]=0;
        mT_old=0;
        mE_old=0;
      
    }
 
  
  
    // DESTRUCTORS 
public:
    
    virtual ~SM_StochasticTrajectory(void) {
    }
 
 
 
public:
    //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
    
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+=mSimulations.size()*sizeof(tIndex);
        mem+=mParticles.size()*sizeof(tIndex);
        mem+=mEnergies.size()*sizeof(tReal);
        return mem;
    }
 
    inline static CORE_UniquePointer<SM_StochasticTrajectory> New() {
        return CORE_UniquePointer<SM_StochasticTrajectory>(new SM_StochasticTrajectory(),SM_StochasticTrajectory::Delete());
    }
    inline static CORE_UniquePointer<SM_StochasticTrajectory> New(const SM_StochasticTrajectory& c) {
        CORE_UniquePointer<SM_StochasticTrajectory> p=New();
        p->copy(c);
        return p;
    }
 
    //SET & GET Methods
    //=================
    virtual void copy(const SM_StochasticOutputComponent& c) override {
        SuperClass::copy(c);
        const SM_StochasticTrajectory *cc=dynamic_cast<const SM_StochasticTrajectory* >(&c);
        if (cc!=null) {
            setSimulations(cc->getSimulations());
            setParticles(cc->getParticles());
        }
    }
    
    virtual void adimensionize(const SM_Material& material) override {
    }
    
    inline void setSimulations(std::initializer_list<tIndex>&& simulationIndices) {
        mSimulations=simulationIndices;
        mISimulation=mSimulations.begin();
    }
    inline void setSimulations(const std::vector<tIndex>& simulationIndices) {
        mSimulations=simulationIndices;
        mISimulation=mSimulations.begin();
    }
   
 
    inline const std::vector<tIndex>& getSimulations() const {
        return mSimulations;
    }
    inline void setParticles(std::initializer_list<tIndex>&& particlesIndices) {
        mParticles=particlesIndices;
    }
    inline void setParticles(const std::vector<tIndex>& particlesIndices) {
        mParticles=particlesIndices;
    }
    inline const std::vector<tIndex>& getParticles() const {
        return mParticles;
    }
 
    
    //implemented templated methods
    //===================
    
    inline tBoolean open(const SM_Beam& beam) {
        //sort the index of simulation
        std::sort(mSimulations.begin(),mSimulations.end());
        //sort the index of particles
        std::sort(mParticles.begin(),mParticles.end());
        //iterator on simulationd
        mISimulation=mSimulations.begin();
        //last registered total energy
        mE_old=0;
        //last registered time
        mT_old=0;
        
        return (dynamic_cast< const SM_LandauLifschitzSystem *>(&beam.getSystem())!=null);
    }
    
    inline tBoolean open(const tIndex& s,const SM_System& system) {
        //reset the file
      
        if (mFile!=null) {
            mFile->close();
            delete mFile;
            mFile=NULL;
        }
        
        
        tBoolean ok=true;
        while (mISimulation!=mSimulations.end()) {//loop on simulations index
            if (s<(*mISimulation)) {
                break;//simulations is not saved
            }
            if (s==(*mISimulation)) {
                if ((getIndex()==0) || (getRootIndex()==-1)) {
                    mFile=new std::ofstream(CORE_IO::GetAbsolutePath(getOutputPath()+"/"
                                                                     +getPrefix()+"-s"+std::to_string(s)
                                                                     +"-i"+std::to_string(getIndex())
                                                                     +".txt").c_str(),
                                            std::ios::out);
                    (*mFile)<<"#trajectory file of simulation at index "<<s<<"\n";
                    tString str=system.toString();
                    functions_string::replaceAll("\n","\n#",str);
                    (*mFile) <<"# X: when the particle with index  does not exists \n";
                    (*mFile)<<"#System:"<<str<<"\n";
                    (*mFile)<<"#timeIndex ";
                    (*mFile)<<"\t time ";
                    (*mFile)<<"\t dt ";
                    (*mFile)<<"\t epsilon ";
                    //direction of particles
                    for(const auto& p : mParticles) {
                        (*mFile)<<"\t S_"<<p<<".x";
                        (*mFile)<<"\t S_"<<p<<".y";
                        (*mFile)<<"\t S_"<<p<<".z";
                    }
                    //total energy
                    (*mFile)<<"\t Etotal";
                    //energies
                    std::vector<tString> names;
                    system.getOperatorsName(names);
                    for(const auto& name:names) (*mFile)<<"\t E"<<name;
                    
                    mEnergies.resize(names.size()+1);
                    
                    //variation of total energy
                    (*mFile)<<"\t DEtotal/Dt";
                    
                    //net torque
                    (*mFile)<<"\t netTorque";
                    
                    //end
                    (*mFile)<<"\n";
                }
 
                mT_old=0;
                mE_old=0;
                break;
            }
 
            //next index imulation to save
            mISimulation++;
        }
        return ok;
    }
   
 
 
    inline tBoolean  store(const SM_System& system) {
        //std::cout<<"begin SM_StochTrajectory::store("<<getIndex()<<")\n";
        //get LL system
        const SM_LandauLifschitzSystem * llSystem=dynamic_cast< const SM_LandauLifschitzSystem *>(&system);
            
        //compute the current energies
        llSystem->computeEnergies(llSystem->getTimeIndex(),mEnergies);
        tReal netTorque=llSystem->computeNetTorque();
            
        //write  mu at time step i
        if (mFile!=NULL) {
 
 
            //number of particles of the system
            tIndex nParticles=system.getNetwork().getParticlesNumber();
            
            //get S at time
            const SM_RealField& S=system.getMagneticMomentDirections();
            //get dimension of S
            tDimension k;
            //iterator on S at particle p
            const tReal *iSp=null;
            //time of the mu
            (*mFile)<<llSystem->getTimeIndex();
            (*mFile)<<"\t"<<std::setprecision(12)<<llSystem->getTime();
            (*mFile)<<"\t"<<std::setprecision(12)<<llSystem->getTimeStepper().getTimeStep();
            (*mFile)<<"\t"<<std::setprecision(12)<<llSystem->getNoiseRate();
            //print the particles
            for(const auto& p : mParticles) {
                if (p<nParticles) {
                    iSp=S(p);
                    for (k=0;k<SM_Constants::DIM;k++) {
                        (*mFile)<<"\t"<<std::setprecision(12)<<(*iSp);
                        iSp++;
                    }
                } else {
                    for (k=0;k<SM_Constants::DIM;k++) {
                        (*mFile)<<"\t\t 0.";
                    }
                }
                
            }
 
     
            //total energy
            tReal dE=mEnergies[0];
            //current time
            tReal dT=llSystem->getTime();
            //print the energy in file
            for(const auto& E : mEnergies) (*mFile)<<"\t"<<E;
 
            //compute time variation
            dT-=mT_old;
            //compute energy variation
            dE-=mE_old;
            //save the current energy and time for the next computation
            mE_old=mEnergies[0];
            mT_old=llSystem->getTime();
            //print the variation of energy with respect of time
            if (dT!=0) (*mFile)<<"\t"<<std::setprecision(12)<<(dE/dT);
            else (*mFile)<<"\t"<<std::setprecision(12)<<dE;
 
            //compute the current net torque
            (*mFile)<<"\t"<<std::setprecision(12)<<netTorque;
 
            
            (*mFile)<<"\n";
            
        }
        //std::cout<<"end SM_StochTrajectory::store("<<getIndex()<<")\n";
        return true;
    }
    inline tBoolean close(const tIndex& s,const SM_System& system,const tBoolean& hasSucceeded) {
        if (mFile!=NULL) {
            mFile->close();
            delete mFile;
            mFile=NULL;
        }
        return true;
    }
    inline tBoolean close(const SM_Beam& beam) {
        return true;
    }
 
    virtual void closePackedSimulations(const SM_Beam& beam,
                                        const tIndex& rawValuesNumber,
                                        const std::valarray<tReal>& rawValues) override {
      
        
    }
   
    
public:
 
   
public:
    // OTHERS methods
  
    
 
 
};
 
#endif