SM_LandauLifschitzSystem.h

#ifndef SM_LandauLifschitzSystem_H
#define SM_LandauLifschitzSystem_H
 
//base class
#include "SM_System.h"
 
//landau lifschtz function
#include "SM_LandauLifschitzFunction.h"
 
//stochastic output data header
#include "SM_StochasticOutput.h"
 
//stochastic functions per thread header
#include "SM_MultiStochasticFunctions.h"
 
//network
#include "SM_Network.h"
 
//material
#include "SM_Material.h"
 
//time
#include "SM_TimeStepper.h"
 
//operator
#include "SM_Operator.h"
 
//noise header
#include "SM_NoiseRateFunction.h"
 
//Temperature <-> noise conversion header
#include "SM_TemperatureNoiseRate.h"
 
//array manipulating
#include <valarray>
 
class SM_LandauLifschitzSystem : public SM_System  {
   
    //attributes
private :
 
    typedef SM_LandauLifschitzSystem Self;
    typedef SM_System SuperClass;
    
    //associations
    
 
    //the landau lishitz function
    CORE_UniquePointer<SM_LandauLifschitzFunction> mLLFunction;
   
    //time stepper
    SM_TimeStepper mTimeStepper;
 
    //total magnetic field of the system of size mParticles number x SM_Constants::DIM at time t
    SM_RealField mHt;
        
    //variation of magnetic moments direction
    SM_RealField mDeltaSt;
    
    //the noise rate functions
    CORE_UniquePointer<SM_NoiseRateFunctionInterface> mNoiseRateFunction;
    
    //noise rate at t
    tReal mEpsilon_t;
 
    //current time
    tReal mTime;
    tIndex mTimeIndex;
 
  
 
    //brownian field
    tBoolean mIsBrownianNull=false;
    tBoolean mIsRealBrownianVariable;
    SM_RealField mDeltaWt;
 
    //relaxation
    tReal mRelaxationTolerance;
    tBoolean mIsRelaxed;
    tReal mRelaxationValue;
    
protected:
    // CONSTRUCTORS 
    SM_LandauLifschitzSystem(void) {
       
      
        mIsRealBrownianVariable=true;
 
        mRelaxationTolerance=1.e-8;
        mRelaxationValue=0;
        mIsRelaxed=false;
        mEpsilon_t=0;
 
        mTime=0;
        mTimeIndex=0;
    }
    
    
    // DESTRUCTORS 
    virtual ~SM_LandauLifschitzSystem(void) {
    }
 
 
public : 
 
     
    //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
 
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+=(mLLFunction.get()==null)?0:mLLFunction->getContentsMemorySize();
        mem+=mTimeStepper.getContentsMemorySize();
        mem+=mHt.getContentsMemorySize();
        mem+=mDeltaSt.getContentsMemorySize();
        mem+=mDeltaWt.getContentsMemorySize();
        mem+=(mNoiseRateFunction.get()==null)?0:mNoiseRateFunction->getMemorySize();
        return mem;
    }
    
    //SET & GET methods
 
    virtual tBoolean updateInitialMagneticMomentDirections(const SM_RealField& S0)  {
        resetTime();
        return SuperClass::updateInitialMagneticMomentDirections(S0);
    }
    
    virtual tBoolean initializeMagneticMomentDirections() override {
        resetTime();
        return SuperClass::initializeMagneticMomentDirections();
       
    }
   
    
    //time stepper data
    //=================
    inline const SM_TimeStepper& getTimeStepper() const {
        return mTimeStepper;
    }
    inline SM_TimeStepper& getTimeStepper() {
        return mTimeStepper;
    }
 
   
    inline const tReal& getTime() const {
        return mTime;
    }
 
    inline const tIndex& getTimeIndex() const {
        return mTimeIndex;
    }
    inline tIndex& getTimeIndex() {
        return mTimeIndex;
    }
 
    inline tReal& getTime()  {
        return mTime;
    }
 
    inline void resetTime() {
        mTime=0;
        mTimeIndex=0;
    }
    
    //noise rate
    //==========
public:
    inline const tReal& getNoiseRate() const {
        return mEpsilon_t;
    }
 
   
    inline tReal& getNoiseRate() {
        return mEpsilon_t;
    }
    inline void setNoiseRateFunction(CORE_UniquePointer<SM_NoiseRateFunctionInterface>& f) {
        mNoiseRateFunction=std::move(f);
    }
 
    inline tBoolean hasNoiseRateFunction() const {
        return (mNoiseRateFunction.get()!=null);
    }
    inline SM_NoiseRateFunctionInterface& getNoiseRateFunction() {
        return *mNoiseRateFunction.get();
    }
    inline const SM_NoiseRateFunctionInterface& getNoiseRateFunction() const {
        return *mNoiseRateFunction.get();
    }
 
    virtual void setNoiseTemperature(const tReal& T)  override {
        if (!hasNoiseRateFunction()) {
            throw CORE_Exception("stochmagnet/landauLifschitz/core",
                                 "setNoiseTemperature(const tReal& T)",
                                 "no noise rate function");
        }
        getNoiseRateFunction().setNoiseRate(SM_TemperatureNoiseRate::ComputeTemperatureNoiseRate(T,*this));
        getNoiseRateFunction().adimensionize(getMaterial().getCharacteristicField());
        SM_System::setNoiseTemperature(T);
    }
   
    //relaxation
    //===========
public:
    inline void setIsRelaxed(const tBoolean& isRelaxed) {
        mIsRelaxed=isRelaxed;
    }
    inline const tBoolean& isRelaxed() const {
        return mIsRelaxed;
    }
    inline void setRelaxationTolerance(const tReal& eps) {
        mRelaxationTolerance=eps;
    }
    inline const tReal&  getRelaxationTolerance() const {
        return mRelaxationTolerance;
    }
    inline const tReal&  getRelaxationValue() const {
        return mRelaxationValue;
    }
    inline tReal&  getRelaxationValue()  {
        return mRelaxationValue;
    }
 
    //relaxation methods
    //====================
 
    inline tReal computeNetTorque() const {
        
        return computeNetTorque(getNetwork().getParticlesNumber(),//number of particles (not halo particles)
                                getMagneticMomentDirections().getValues(),//S values
                                mHt.getElementsNumber(),//number of elements in H
                                mHt.getValues());//H values
    }
    
    virtual tReal computeNetTorque(const tIndex& nParticles,
                                   const tReal* S,
                                   const tIndex& nHParticles,const tReal* H) const=0;
 
protected:
    tReal computeNetSliceTorque(const tIndex& startS,const tIndex& endS,const tReal* S,
                                const tIndex& startH,const tReal* H) const;
    
    
    //brownian function
protected:
    inline void setIsBrownianNull(const tBoolean& v) {
        mIsBrownianNull=v;
    }
    inline const tBoolean& isBrownianNull() const {
        return mIsBrownianNull;
    }
 
    inline void setIsRealBrownianVariable(const tBoolean& v) {
        mIsRealBrownianVariable=v;
    }
 
    inline const tBoolean& isRealBrownianVariable() const {
        return mIsRealBrownianVariable;
    }
 
    inline const SM_RealField& getBrownianField() const {
        return mDeltaWt;
    }
    inline SM_RealField& getBrownianField() {
        return mDeltaWt;
    }
   
 
    //Landau lifschitz function
    //===========================
public:
    inline const SM_LandauLifschitzFunction& getLandauLifschitzFunction() const {
        return *mLLFunction.get();
        
    }
    inline SM_LandauLifschitzFunction& getLandauLifschitzFunction() {
        return *mLLFunction.get();
        
    }
 
protected:
 
    inline void setLandauLifschitzFunction(CORE_UniquePointer<SM_LandauLifschitzFunction> llF) {
        mLLFunction=std::move(llF);
    }
 
   
public:
  
 
    virtual void normalize(const tBoolean& isAdimensionized) override {
        
        SuperClass::normalize(isAdimensionized);
        
        //normalize time step
        getTimeStepper().adimensionize(getMaterial().getCharacteristicTime());
        
        if (isAdimensionized) {
            
            //adimenesionize LL
            getLandauLifschitzFunction().setGamma(1);
            
            
            //normalize the noise rate
            if (hasNoiseRateFunction()) getNoiseRateFunction().adimensionize(getMaterial().getCharacteristicField());
        }
    }
    
 
public:
 
    inline const SM_RealField& getMagneticField() const {
        return mHt;
    }
    inline const SM_RealField& getMagneticMomentDirectionsVariation() const {
        return mDeltaSt;
    }
    
protected:
    inline SM_RealField& getMagneticField() {
        return mHt;
    }
    
    inline SM_RealField& getMagneticMomentDirectionsVariation() {
        return mDeltaSt;
    }
 
    
 
public:
    virtual void discretize() override {
 
        SuperClass::discretize();
     
    }
 
 
    
 
    //simulation methods
    //===================
    virtual tULLInt getGeneratedRandomNumbers() const override {
        
        tULLInt nRNGPerCore=SM_Constants::DIM;
        nRNGPerCore*=getNetwork().getParticlesNumber();
        return nRNGPerCore;
    }
   
    template<class StochOutputI>
    tBoolean  deterministicRun(const tIndex& steppersNumber,
                               SM_StochasticOutput<StochOutputI>& stochasticOutput);
    
    tBoolean stochasticRun(const tIndex& steppersNumber,
                           SM_MultiStochasticFunctionsInterface& sfs);
    
 
      
    template<class StochOutputImplement>
    tBoolean stochasticRun(const tIndex& steppersNumber,
                           SM_MultiStochasticFunctionsInterface& multiStochasticFunctions,
                           SM_StochasticOutput<StochOutputImplement>& stochasticOutput);
 
 
    tBoolean stochasticJRun(const tIndex& steppersNumber,
                            SM_StochasticFunctionsInterface& stochasticFunctions,
                            const std::function<void(const tReal&,const tReal& ,const tReal&,tReal& )>& F);
    
    template<class StochOutputImpl>
    tBoolean stochasticJRun(const tIndex& steppersNumber,
                            SM_StochasticFunctionsInterface& stochasticFunctions,
                            SM_StochasticOutput<StochOutputImpl>& stochasticOutput,
                            const std::function<void(const tReal&,const tReal& ,const tReal&,tReal& )>& F);
    
   
   
 
 
 
    //DH methods
    //===========
 
    virtual void computeVariationMagneticField(const tReal& dt,SM_RealField& H) const=0;
    
    virtual void computeBrownianVariationMagneticField(SM_MultiStochasticFunctionsInterface& stochFs,
                                                       const tReal& dt,const tReal& bdt,const tReal& epsilon,
                                                       const tBoolean& incdWt,tReal* dWt,
                                                       SM_RealField& H)=0;
 
protected:
    inline void computeVariationMagneticFieldSlice(const tIndex& start,const tIndex& end,const tReal& dt,tReal* H) const {
 
        //no H
        if (H==null) return;
            
        //begin iterator on H
        tReal *iH=H;iH+=start*SM_Constants::DIM;
        //end iterator on H
        const tReal *eH=H;eH+=end*SM_Constants::DIM;
        
        while (iH!=eH) {
            (*iH)*=dt;
            iH++;
        }
    }
    inline void computeBrownianVariationMagneticFieldSlice(SM_StochasticFunctionsInterface& randF,
                                                           const tIndex& start,const tIndex& end,
                                                           const tReal& dt,const tReal& bdt,const tReal& epsilon,
                                                           const tBoolean& incdWt,tReal* dWt,tReal* H) {
 
        if (H==null) return;
        
        //iterator on H
        tReal *iH=H;
        iH+=start*SM_Constants::DIM;
        
        //end iterator on H
        const tReal *eH=H;
        eH+=end*SM_Constants::DIM;
        
        //iterator ondWt
        tReal *idWt=dWt;
        idWt+=incdWt*start*SM_Constants::DIM;
        
        while (iH!=eH) {//loop on H values
            
            //compute the brownian
            (*idWt)=randF.normalRandom();
            (*idWt)*=bdt;
            (*idWt)*=epsilon;
            
            //compute time variation of H
            (*iH)*=dt;
            (*iH)+=(*idWt);
            
            //iterator at next value
            iH++;
            idWt+=incdWt;
        }
    }
 
    
    
    template<class StochFImpl>
    inline void computeBrownianVariationMagneticFieldSliceWithSCStochasticFunctions(SM_StochasticFunctions<StochFImpl>& randF,
                                                                                    tIndex start,tIndex end,
                                                                                    const tReal& dt,const tReal& bdt,const tReal& epsilon,
                                                                                    const tBoolean& incdWt,tReal* dWt,tReal* H) {
 
        if (H==null) return;
        
        //particle index -> memory index
        start*=SM_Constants::DIM;
        end*=SM_Constants::DIM;
        
        //iterator on H
        tReal *iH=H;
        iH+=start;
        
        //end iterator on H
        const tReal *eH=H;
        eH+=end;
        
        //iterator ondWt
        tReal *idWt=dWt;
        idWt+=incdWt*start;
        
        while (iH!=eH) {//loop on H values
            
            //compute the brownian
            (*idWt)=randF.scNormalRandom();
            (*idWt)*=bdt;
            (*idWt)*=epsilon;
            
            //compute time variation of H
            (*iH)*=dt;
            (*iH)+=(*idWt);
            
            //iterator at next value
            iH++;
            idWt+=incdWt;
        }
    } 
    
 
    //next S direction
    //================
 
public:
    
    
    virtual  tBoolean computeMagneticMomentDirectionsAtNextTimeStep(const tReal& dt,const tReal& epsilon_t,
                                                                    const SM_RealField& deltaSt,
                                                                    SM_RealField& St)=0;
    
protected:
    static  void  ComputeMagneticMomentDirectionsSliceAtNextTimeStep(tIndex startS,tIndex endS,
                                                                     const tIndex& nDeltaS,const tReal *deltaS,
                                                                     const tIndex& nS,const tReal* S,
                                                                     tReal *nextS);
    
    static void  ComputeMagneticMomentDirectionsSliceAtNextTimeStep(tIndex startS,tIndex endS,
                                                                    const tIndex& nDeltaS,const tReal *deltaSt,
                                                                    const tReal& alpha,
                                                                    const tIndex& nS,const tReal* S,
                                                                    tReal *nextS);
    
    //string representation
    //=====================
public:
    virtual tString toString() const override;
};
 
#include "SM_LandauLifschitzSystem.hpp"
#endif