SM_BeamHysteresis.h

#ifndef SM_BeamHysteresis_H
#define SM_BeamHysteresis_H
 
//inherits header
#include "SM_BeamCycle.h"
 
//zeeman operator header
#include "SM_ZeemanOperator.h"
 
//OpenMP header
#include "openMP.h"
 
//ll system header
#include "SM_LandauLifschitzSystem.h"
 
//associated stochastic output
#include "SM_StochasticHysteresis.h"
 
//chrono
#include "CORE_Chrono.h"
 
class SM_BeamHysteresis : public SM_BeamCycle<SM_BeamHysteresis> {
    
    //attributes
private :
    
    //class type
    typedef  SM_BeamHysteresis SelfClass;
    typedef  SM_BeamCycle<SelfClass> SuperClass;
    
    std::array<tReal,7> mHextRange;
 
    //number of steps for a complete cycle of hysteresis
    tInteger mCycleStepsNumber;
    
    std::array<tReal,SM_Constants::DIM> mHextMax;
    std::array<tReal,SM_Constants::DIM> mHextMin;
    std::array<tReal,SM_Constants::DIM> mHext;
    tInteger mHextNumber;
    
    std::array<tReal,SM_Constants::DIM> mDHext;//direction of Hext
    
    tReal mAHext;//amplitude of Hext
    tReal mDH;//variation of hext
    tReal mH;//value of projection of Hext on D
    tReal mHysValue;//hysteresis value
    
    tCInt mHysteresisDirection;//direction of hysteresis
    tRelativeInteger mHysteresisIndex;//index of the hysteresis in direction
 
    //backup
    tInteger mBackupIndex;//index of backup
    tInteger mBackupStep;//cycle steps number from last backup
    tInteger mBackupsNumber;//number of backup files
    tInteger mBackupStepsNumber;//number of cycle steps between 2 backups
 
    //critical point
    tBoolean mAreCriticalPointsSaved;
    tInteger mCriticalPointIndex;
    tBoolean mIsCriticalPoint;
    
    //restoring
    tBoolean mIsRestoring;
 
    //hysteresis file
    tString mHysFileName;
 
  
 
    CORE_Chrono mTimer;
    
    SM_ZeemanOperator* mZeemanOperator;//zeeman operator
 
    
public:
    // CONSTRUCTORS 
    SM_BeamHysteresis(void) {
        mAHext=0;
        mDH=0;
        mH=0;
        mHextNumber=10;
        mHysteresisDirection=0;
        mHysteresisIndex=0;
        mHysFileName="";
        mZeemanOperator=null;
        mBackupIndex=0;
        mBackupsNumber=2;
        mBackupStep=0;
        mBackupStepsNumber=1;
        memset(mHextMin.data(),0,sizeof(tReal)*SM_Constants::DIM);
        memset(mHextMax.data(),0,sizeof(tReal)*SM_Constants::DIM);
        memset(mHext.data(),0,sizeof(tReal)*SM_Constants::DIM);
        memset(mDHext.data(),0,sizeof(tReal)*SM_Constants::DIM);
        mZeemanOperator=null;
        mAreCriticalPointsSaved=true;
        mCriticalPointIndex=0;
        mIsCriticalPoint=false;
        mCycleStepsNumber=0;
    }
  
    // DESTRUCTORS 
    virtual ~SM_BeamHysteresis(void) {
    }
    
 
 
public : 
    // CREATE class
   
 
    //SET & GET methods
 
   
 
    
      //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+=4*SM_Constants::DIM*sizeof(tReal);//mHextMax+mHextMin+mHext+mDHext
        mem+=mHysFileName.length()*sizeof(tChar);
        mem+=mHextRange.size()*sizeof(tReal);
        return mem;
    }
    
    
    static CORE_UniquePointer<SelfClass> New() {
        return CORE_UniquePointer<SelfClass>(new SelfClass(),
                                             SelfClass::Delete()); 
        
    }
public:
  
public:
 
    inline void setAreCriticalPointsSaved(const tBoolean& flag) {
        mAreCriticalPointsSaved=flag;
    }
    
    inline void setBackupsNumber(const tInteger& nBackups) {
        mBackupsNumber=nBackups;
    }
    
    
    inline void setBackupStepsNumber(const tInteger& s) {
        mBackupStepsNumber=s;
    }
    
    
    inline void setHextRange(const std::array<tReal,7>& range) {
        mHextRange=range;
        
    }
    
    
    inline void setHysteresisFileName(const tString& hysteresisFile) {
        mHysFileName=hysteresisFile;
        tIndex iExt=mHysFileName.rfind(".hys");
        if (iExt==tString::npos) {
            mHysFileName+=".hys";
        }
    }
    
    
public:
    
    inline tBoolean open(SM_Beam& beam,SM_StochasticOutputComponent& stochasticOutput,
                         const tInteger& outputId,
                         tIndex& iCycleStep) {
 
        //std::cout<<"begin SM_BeamHysteresis::open("<<outputId<<")\n";
 
        mCycleStepsNumber=0;
        
        //system
        SM_System& system=beam.getSystem();
        
        //compute Hext Min & Hext Max
        const tReal *iRange=mHextRange.data();
        for(auto& Hk:mHextMin) {
            Hk=(*iRange);
            iRange++;
        }
        for(auto& Hk:mHextMax) {
            Hk=(*iRange);
            iRange++;
        }
        //number of Hext
        mHextNumber=mHextRange[6];
        
        //characteristic field
        const tReal& cField=system.getMaterial().getCharacteristicField();
        
        //adimensionize Hext values
        for(auto& Hk:mHextMax) {
            Hk/=cField;
        }
        for(auto& Hk:mHextMin) {
            Hk/=cField;
        }
        
        //D=(Hmax-Hmin)
        //hmax=|Hmax-Hmin|
        //D=Hmax
        memcpy(mDHext.data(),mHextMax.data(),SM_Constants::DIM*sizeof(tReal));//dest,source,size
        mAHext=0;//amplitude of Hext
        //D=Hmax-Hmin
        //hmax=|Hmax-Hmin|
        const tReal *iHext=mHextMin.data();
        for(auto& Dk:mDHext) {
            Dk-=(*iHext);
            mAHext+=Dk*Dk;
            iHext++;
        }
        mAHext=sqrt(mAHext);
        
        //h variation
        mDH=mAHext/mHextNumber;
        
        //normalize D
        for(auto& Dk:mDHext) Dk/=mAHext;
        
        //get zeeman operator
        mZeemanOperator=dynamic_cast<SM_ZeemanOperator*>(system.getOperator("Zeeman")); 
        if (mZeemanOperator==null) {
            std::cout<<"fatal error : use-zeeman must be true in options\n";
            return false;
            
        }
        //set the adimensionized Hext
        memcpy(mHext.data(),mHextMax.data(),sizeof(tReal)*mHextMax.size());//dest,source,size
        mZeemanOperator->setHext(mHext);
        
        //number of cycles steps
        mCycleStepsNumber=2*mHextNumber+1;
        
        
        //set S0 to Hext
        tIndex nParticles=system.getInitialMagneticMomentDirections().getElementsNumber();
        system.setInitialMagneticMomentDirections(nParticles,mDHext);
        
        //set the stochastic data
        SM_StochasticHysteresis* hysOutput=dynamic_cast<SM_StochasticHysteresis*>(&stochasticOutput);
        if (hysOutput==null) {
            throw CORE_Exception("stochMagnet/main/cycle",
                                 "SM_BeamHysteresis::open(...)",
                                 "Fatal Error: bad stochastic output for this cycle. Aborted");
        }
        hysOutput->setHysteresisDirection(mDHext);
 
        mBackupStep=0;
        mBackupIndex=0;
        
        if (isRestoring() && CORE_IO::Exists(mHysFileName))  {
            //restore from last hysteresis step
            if (!restoreFromLastBackup(system,iCycleStep,outputId)) {
                throw CORE_Exception("stochMagnet/main/cycle",
                                     "SM_BeamHysteresis::open(...)",
                                     "Fatal Error: impossible to restore hysterisis cycle. Aborted");
            }
            
        } else {
            iCycleStep=0;
            //not restoring hysteresis: begin at hysteresis first step
            if (outputId==0) {
                std::ofstream file(mHysFileName.c_str(),std::ios::out);
                //std::cout<<"create Hysteresis file:"<<mHysFileName<<"\n";
                if (!file) {
                    std::cout<<"impossible to open Hysteresis file : "<<mHysFileName<<"\n";
                    return false;
                }
                
                //print header in hyteressis file
                printHeader(file,beam);
                
                //close the file
                file.close();
            }
            
            //decreasing hext cycle
            mHysteresisDirection=-1;
            mHysteresisIndex=0;
            mCriticalPointIndex=0;
            mIsCriticalPoint=false;
            mHysValue=1;
        }
        
       
        //start the timer
        mTimer.start();
        
        //std::cout<<"end SM_BeamHysteresis::open("<<outputId<<")\n";
        return true;
    }
    
    inline void  nextCycleStep(const tIndex& iCycleStep,SM_Beam& beam,const tInteger& outputId,const SM_StochasticOutputComponent& output) {
        //std::cout<<"begin SM_BeamHysteresis::nextCycleStep("<<outputId<<")\n";
        //system
        SM_System& system=beam.getSystem();
        SM_LandauLifschitzSystem* llSystem=dynamic_cast<SM_LandauLifschitzSystem*>(&system);
        
        //S
        const SM_RealField& S=system.getMagneticMomentDirections();
        
        //compute the hysteresis value
        tReal oldHysValue=mHysValue;
        mHysValue=output.getOutputValues()[0];
        
        //critical point when hys value changes of sign
        if (!mIsCriticalPoint) mIsCriticalPoint=(mHysValue*oldHysValue<=0);
        
        //mH=<mHext,mDHext>
        mH=0;
        const tReal *iHext=mHext.data();
        for(const auto&  Dk:mDHext) {
            mH+=(*iHext)*Dk;
            iHext++;
        }
 
        tReal netTorque=llSystem->computeNetTorque();
        
        //add the value into file
        if (outputId==0) {
            
            std::ofstream file(mHysFileName.c_str(),std::ios::app);
            file<<((int)mHysteresisDirection)<<"\t"<<mHysteresisIndex<<"\t"<<std::setprecision(12)<<mH<<"\t"<<system.getStepIndex();
            file<<"\t"<<std::setprecision(12)<<mHysValue;
            if (llSystem!=null) {
                file<<"\t"<<std::setprecision(12)<<llSystem->getTime();
                file<<"\t"<<std::setprecision(12)<<netTorque;
            }
            file<<"\t"<<(mTimer.stop()/1000);//chorno is stop and restart
            file<<"\n";
            file.close();
            mBackupStep++;
            if (mIsCriticalPoint && mAreCriticalPointsSaved) {
                tString SFileName=getSFile()+std::to_string(mCriticalPointIndex)+".cpt";
                std::stringstream SHeader;
                SHeader<<"dir="<<((int)mHysteresisDirection)<<"\n";
                SHeader<<"i="<<((int)mHysteresisIndex)<<"\n";
                SHeader<<"h="<<std::setprecision(12)<<mH<<"\n";
                SHeader<<"criticalPointIndex="<<mCriticalPointIndex<<"\n";
                S.saveToFile(SHeader.str(),SFileName,12);
                mCriticalPointIndex++;
                mIsCriticalPoint=false;
            }
             
            //save S for future restore
            if ((mBackupStep==mBackupStepsNumber)) {
                tString SFileName=getSFile()+std::to_string(mBackupIndex)+".back";
                std::stringstream SHeader;
                SHeader<<"dir="<<((int)mHysteresisDirection)<<"\n";
                SHeader<<"i="<<((int)mHysteresisIndex)<<"\n";
                SHeader<<"h="<<std::setprecision(12)<<mH<<"\n";
                SHeader<<"criticalPointIndex="<<mCriticalPointIndex<<"\n";
                if(S.saveToFile(SHeader.str(),SFileName,12)) {
                    mBackupIndex++;
                    mBackupIndex%=mBackupsNumber;
                }
                mBackupStep=0;
            }
           
          
        }
        
        
        //next cycle
        moveToNextCycle(iCycleStep,S,system,outputId);
        
        //std::cout<<"end SM_BeamHysteresis::nextCycleStep("<<outputId<<")\n";
    }
    
private:
    inline void moveToNextCycle(const tIndex& iCycleStep,const SM_RealField& S0,SM_System& system,const tInteger& outputId) {
        
        if (iCycleStep==mHextNumber) {
            //mHext=mHextMin
            //next cycle is in the other dierction
            mHysteresisDirection*=-1;
        }
        mHysteresisIndex-=mHysteresisDirection;
        
        
        
        //next H
        tReal Hold=mH;
        mH+=mDH*mHysteresisDirection;
 
        //critical point when H changes of sign
        if (!mIsCriticalPoint) mIsCriticalPoint=(mH*Hold<=0);
        
        //next Hext=mH.DHext
        const tReal *iD=mDHext.data();
        for(auto& Hk:mHext) {
            Hk=(*iD)*mH;
            iD++;
        }
        
        //set the adimensionized Hext
        mZeemanOperator->setHext(mHext);
        
        //set the next initial direction to last magnetic moment direction computed
        //std::cout<<"S set to "<<S0<<" for id:"<<outputId<<" for cycle step "<<iCycleStep<<"\n";
        system.updateInitialMagneticMomentDirections(S0);
        
    }
    
    inline  tBoolean restoreFromLastBackup(SM_System& system,tIndex& iCycleIndex,const tInteger& outputId) {
        //std::cout<<"begin SM_BeamHysteresis::restoreFromLastBacku(system,"<<iCycleIndex<<")\n";
        std::valarray<tReal> hValues;
        std::valarray<tReal> hysValues;
        tInteger nHysValues;
        if (!LoadFromFile(mHysFileName,mDHext,mAHext,mDH,mCycleStepsNumber,
                          mHysteresisDirection,mHysteresisIndex,mH,
                          nHysValues,hValues,hysValues)) return false;
        mHysValue=hysValues[nHysValues-1];
        
        //set the number of hysteresis values
        mHextNumber=(mCycleStepsNumber-1)/2;
        
        
        //load S at last computed hysteresis values
        SM_RealField S;
        tString SFile;
        tRelativeInteger d=0;//direction of hysteresis
        tRelativeInteger i=-1;//index of the hysteresis in direction
        tString line;
        int nMaxReadData=2;
        int nReadData=0;
        for (mBackupIndex=0;mBackupIndex<mBackupsNumber;mBackupIndex++) {
            SFile=getSFile()+std::to_string(mBackupIndex)+".back";
            std::ifstream file(SFile.c_str(),std::ios::in);
            nReadData=0;
            while ( (!file.eof()) && (std::getline(file,line)) && (nReadData<nMaxReadData)) {
                nReadData+=ReadKeyValue(line,"#dir",d);
                nReadData+=ReadKeyValue(line,"#i",i);
            }
            if ( (((tCInt)d)==mHysteresisDirection) && (i==mHysteresisIndex)) {
                //read critical point
                nMaxReadData++;
                while ( (!file.eof()) && (std::getline(file,line)) && (nReadData<nMaxReadData)) {
                    nReadData+=ReadKeyValue(line,"#criticalPointIndex",mCriticalPointIndex);
                }
                break;
            }
        }
        if (mBackupIndex==mBackupsNumber) {
            std::cout<<"can not find the S backup file from backup files of type "<<SFile<<" and index in [0,"<<mBackupsNumber<<"[ with hysteresis direction \n";
            std::cout<<std::to_string((int)mHysteresisDirection)+" and index ";
            std::cout<<std::to_string((int)mHysteresisIndex)+" ";
            return false;
        }
       
    
        //load S
        S.loadFromFile(SFile);
        
        //next hysteresis step
        iCycleIndex=mHysteresisIndex;
        if (mHysteresisDirection==1) {
            iCycleIndex=mHextNumber;
            iCycleIndex+=mHextNumber;
            iCycleIndex-=mHysteresisIndex;
        }
        moveToNextCycle(iCycleIndex,S,system,outputId);
        iCycleIndex++;
 
        //std::cout<<"end SM_BeamHysteresis::restoreFromLastBacku(system,"<<iCycleIndex<<")\n";
        return true;
    }
    
  
    inline tString getSFile() const {
        tString SFileName=mHysFileName;
        tIndex iExt=SFileName.rfind(".hys");
        SFileName=CORE_IO::GetAbsolutePath(SFileName.substr(0,iExt))+"-S";
        return SFileName;
    }
    
public :
    inline static tBoolean LoadFromFile(const tString& fileName,
                                        std::array<tReal,SM_Constants::DIM>& direction,
                                        tReal& amplitude,
                                        tReal &dh,
                                        tInteger& nCycleSteps,
                                        tCInt& hysDirection,
                                        tRelativeInteger& hysStep,
                                        tReal& h,
                                        tInteger& nHysValues,std::valarray<tReal>& hs,std::valarray<tReal>& hys) {
        //restoring hysteresis
        tString line;
        std::ifstream file(fileName.c_str(),std::ios::in);
        if (!file) {
            std::cout<<"file  "<<fileName<<" does not exist\n";
            return false;
        }
        
        //load the hysteresis data
        int nReadData=0;
      
        std::vector<tString> keys={
            "Direction of Hext",
            "Amplitude of Hext",
            "Variation of Hext",
            "number of steps",
        };
        int nMaxReadData=keys.size();
       
        std::vector<tString>::const_iterator iKey;
        while ( (!file.eof()) && (std::getline(file,line)) && (nReadData<nMaxReadData)) {
            iKey=keys.begin();
            nReadData+=ReadKeyValue(line,*iKey,direction);iKey++;
            nReadData+=ReadKeyValue(line,*iKey,amplitude);iKey++;
            nReadData+=ReadKeyValue(line,*iKey,dh);iKey++;
            nReadData+=ReadKeyValue(line,*iKey,nCycleSteps);iKey++;
        }
        if (nReadData!=nMaxReadData) {
            std::cout<<"error in reading the data from hysteresis file "<<fileName<<" nReadData:"<<nReadData<<" != "<<nMaxReadData<<"\n";
            for(const auto& key:keys) {
                std::cout<<"\t "<<key<<"\n";
            }
            return false;
        }
        hs.resize(nCycleSteps);
        hys.resize(nCycleSteps);
        tReal *ihs=&hs[0];
        tReal* ihys=&hys[0];
        tIndex step;
        //ignore the comments
        while ( (!file.eof()) && (std::getline(file,line)) && (line[0]=='#'));
        //read the hyteresis file
        nHysValues=0;
        do {
            std::stringstream sstream;
            sstream<<line;
            sstream>>step;hysDirection=step;
            sstream>>hysStep;
            sstream>>h;(*ihs)=h;ihs++;
            sstream>>step;
            sstream>>(*ihys);ihys++;
            nHysValues++;
            
        } while ( (!file.eof()) && (std::getline(file,line)) && (line.length()>0));
        
        file.close();
        return true;
    }
    
       
public:
    inline tBoolean  isCycleFinished(const tIndex& iCycle,const SM_Beam& beam) {
        return (iCycle==mCycleStepsNumber);
    }
    
    virtual tString toString() const override {
        std::stringstream ret;
        ret<<SuperClass::toString()<<"\n";
        ret<<"\t Direction of Hext="<<functions_array::toString(mDHext)<<"\n";
        ret<<std::setprecision(12)<<"\t Amplitude of Hext="<<mAHext<<"\n";
        ret<<std::setprecision(12)<<"\t Variation of Hext="<<mDH<<"\n";
        ret<<std::setprecision(12)<<"\t Adimensionized Max of Hext="<<functions_array::toString(mHextMax)<<"\n";
        ret<<"\t number of steps:"<<mCycleStepsNumber<<"\n";
        return ret.str();
    }
 
    inline tBoolean close(const tIndex& iCycle,const SM_Beam& beam) {
        return true;
    }
private:
    inline void printHeader(std::ofstream& file,const SM_Beam& beam) const {
        
        file<<"# Hysteresis cycle Beam generated by soft:"<<CORE_Run::GetSoftName()<<" Version:"<<CORE_Run::GetVersion()<<"\n";
        tString str="Hysteresis Cycle:"+this->toString();
        functions_string::replaceAll("\n","\n#",str);
        file<<"#"<<str<<"\n";
        
        str="beam:"+beam.toString();
        functions_string::replaceAll("\n","\n#",str);
        file<<"#"<<str<<"\n";
 
        file<<"# [...] : for Landau Lifschitz system only \n";
        file<<"#<hysteresis direction> <hysteresis index> <H> <Step Index> <HysValue> [<time> <netTorque>]  <duration in s>\n";
    }
 
 
 
    template<typename T,size_t D>
    static tBoolean ReadKeyValue(tString& line,const tString& key,std::array<T,D>& values) {
        auto *iValues=values.data();
        size_t s=0;
        T v;
        tIndex iKey=line.find(key);
        if (iKey!=tString::npos) {//key
            line=line.substr(iKey+key.length());
            iKey=line.find("=");
            if (iKey==tString::npos)  iKey=line.find(":");
            if (iKey!=tString::npos) {//key = 
                line=line.substr(iKey+1);
                iKey=line.find("{");
                if (iKey==tString::npos)  iKey=line.find("[");
                if (iKey!=tString::npos) {//key = {
                    line=line.substr(iKey+1);
                    iKey=line.rfind("}");
                    if (iKey==tString::npos)  iKey=line.find("]");
                    if (iKey!=tString::npos) {//key = { }
                        
                        //read the vector
                        line=line.substr(0,iKey);
                        iKey=line.rfind(",");
                        while (iKey!=tString::npos) {
                            std::stringstream vss;
                            vss<<std::setprecision(12)<<line.substr(0,iKey);
                            vss>>v;
                            if (s<D) {
                                (*iValues)=v;
                                iValues++;
                                s++;
                            }
                            line=line.substr(iKey+1);
                            iKey=line.rfind(",");
                        }
                        std::stringstream vss;
                        vss<<std::setprecision(12)<<line;
                        if (s<D) {
                            (*iValues)=v;
                            iValues++;
                            s++;
                        }
                        return true;
                    }//end reading key value
                }//end reading {
            }//end reading =
        }//end reading key
        return false;
    }//end method
    
    template<typename T>
        static tBoolean ReadKeyValue(tString& line,const tString& key,T& value) {
        tIndex iKey=line.find(key);
        if (iKey!=tString::npos) {//key
            line=line.substr(iKey+key.length());
            iKey=line.find("=");
            if (iKey==tString::npos)  iKey=line.find(":");
            if (iKey!=tString::npos) {//key = 
                line=line.substr(iKey+1);
                std::stringstream vss;
                vss<<std::setprecision(12)<<line;
                vss>> value;
                return true;
            }//end =
        }//end key
        return false;
    }//end method
    
    
};
    
 
#endif