EMM_LandauLifschitzSystem.h

Go to the documentation of this file.

#ifndef EMM_LandauLifschitzSystem_H
#define EMM_LandauLifschitzSystem_H

#include "CORE_File.h"

#include "EMM_LandauLifschitzFunction.h"
#include "EMM_Stepper.h"
#include "EMM_Grid3D.h"
#include "EMM_Matter.h"
#include "EMM_MatterField.h"
#include "EMM_RealField.h"
#include "EMM_RealArray.h"
#include "EMM_MagneticExcitationOperator.h"

DEFINE_SPTR(EMM_LandauLifschitzSystem);
class EMM_LandauLifschitzSystem : public  EMM_Object {
  
    SP_OBJECT(EMM_LandauLifschitzSystem);
    // ATTRIBUTES
  
public:
   
   
    //retcode for relaxation
    static const tSInt MIN_CODE;
    static const tSInt MAGNETIZATION_FIELD_UPDATE_ERROR;
    static const tSInt OPERATORS_FIELDS_UPDATE_ERROR;
    static const tSInt E_INCREASING;
    static const tSInt E_IS_NAN;
    static const tSInt RELAXATION_IS_REACHED;
    static const tSInt DE_DT_MIN_IS_REACHED;
    static const tSInt MAXIMUM_ITERATION_IS_REACHED;
    static const tSInt M_PARALLEL_H;
    static const tSInt DT_MIN_IS_REACHED;

    static const tFlag OsM=1;//Operataor field and  M are evoluating
    static const tFlag MOs=2;//M and operator field are evoluating
    static const tFlag M=3;//only M is evoluating
    static const tFlag Os=4;//only Operator fields is evoluating
    
private:
    static const tString RET_STRING[];

    
private:

    //max number of operators
    static const tFlag MAX_OPERATORS=5;
private:

    //Backup  data
    //=============
    //the output path for saving file
    tString mOutputPath;
   
    //the prefix for loading/saving the files
    tString mPrefix;


    //backup step
    tUInteger mBackupStep;
    
    //backups number
    tUInteger mBackupsNumber;

   
    //indicates if the backup is enabled
    tBoolean mIsBackupEnabled;

    //log
    //====
      
    //indicates if the log is printed in log file
    tBoolean mIsLogPrinted;
    
    //M evolution data
    //============
    //for initialization of M in a file
    tString mMInitialDataFile;
    
    //for initialization of uniform M 
    SP::EMM_RealField mMInitialData;
    
    //indicate the evolution scheme between operators fields & Magnetization field
    tFlag mEvolutionScheme;
    

    //integration data
    //=================
    //min time step reachabled
    tReal mMinDt;
    //the relaxation criterium
    tReal mMinEnergyVariation;

    //Energy data
    //============
    
    //energy
    tReal mEs[MAX_OPERATORS+1];
    
    //variation of Energy among t
    tReal mdE_dt;
    

    //ASSOCIATIONS
    //============
    
    //Landau Lifschitz function
    SP::EMM_LandauLifschitzFunction mLLFunction;
  
 
    //cells of the problem
    SP::EMM_Grid3D mMesh;
    
    //matter indices  for each cell
    SP::EMM_MatterField mMatters;
    
    //operators of the problem
    SV::EMM_MagneticExcitationOperator mOperators;

    //time step of the system
    SP::EMM_Stepper mStepper;
      
    //the weight of the M at eah cell (default 1. for uniform matters
    SP::EMM_RealArray mSigma;
    
    //M variable: the normalized magnetization field at t=0 into current time step
    SP::EMM_RealField mM;

    //the landau lifschitz function F = dM_dt
    SP::EMM_RealField mF;
    
    //the complete excitation field at t=0 into current time step
    SP::EMM_RealField mH;
    
    //magnetic excitation field for each operator
    mutable SP::EMM_RealField mHop;

    
   

    
    //states wich defines the relaxation process:
    //===========================================
   
    
protected: 
    // METHODS
  
    // CONSTRUCTORS 
  
    EMM_LandauLifschitzSystem(void);
  
  
  
    // DESTRUCTORS 
  
    
    virtual ~EMM_LandauLifschitzSystem(void);   
  
  
public:

    //NEW constructor
protected:
 

    
    //mesh data
    //=============
public:
    inline tBoolean setMesh(SP::EMM_Grid3D mesh) {
        if (mesh.get()!=null) {
            mMesh=mesh;
            return true;
        }
        return false;
    }
    inline const EMM_Grid3D& getMesh() const {
        return *mMesh.get();
    }
    inline EMM_Grid3D& getMesh()  {
        return *mMesh.get();
    }  

    inline const tUIndex& getMagnetizedElementsNumber() const {
        return mMesh->getMagnetizedElementsNumber();
    }
    
    //matter data
    //============
    
    inline const EMM_MatterField& getMatterField() const {
        return *mMatters.get();
    }
    inline EMM_MatterField& getMatterField()  {
        return *mMatters.get();
    }

    //output data

    inline const tString& getPrefix() const {
        return mPrefix;
    }
    inline const tString& getOutputPath() const {
        return mOutputPath;
    }


    // initializer methods
    //====================
    
    virtual void adimensionize();
    

    virtual tBoolean discretize(const tBoolean& isMnormalized);
    
    inline tBoolean discretize() {
        return discretize(false);
    }
    

    //OPERATOR  methods
    //==================
      
public:
    virtual tBoolean addOperator(SP::EMM_MagneticExcitationOperator op) {
        if (op.get()!=null) {
            mOperators.add(op);
            return true;
        }
        return false;
    }
    inline void clearOperators() {
        mOperators.clear();
    }
   
    inline const SV::EMM_MagneticExcitationOperator& getOperators() const {
        return mOperators;
    } 
   
    inline SV::EMM_MagneticExcitationOperator& getOperators()  {
        return mOperators;
    }
    
public:
    inline tUIndex getOperatorsNumber() const {
        return mOperators.getSize();
    }
    inline const EMM_Operator&  getOperator(const tUIndex& index) const {
        if (index>=mOperators.getSize())
            throw EMM_Exception("EMM_LandauLifschitzSystem::getOperator("+CORE_Integer::toString(index)+")",
                                "index out of bounds [0,"+CORE_Integer::toString(mOperators.getSize())+"[");
        const EMM_Operator *op=mOperators[index].get();
        if (op==null) {
            throw EMM_Exception("EMM_LandauLifschitzSystem::getOperator("+CORE_Integer::toString(index)+")",
                                "null operator at index");
        }
        return *op;
    }
    
    inline  EMM_Operator&  getOperator(const tUIndex& index)  {
        if (index>=mOperators.getSize())
            throw EMM_Exception("EMM_LandauLifschitzSystem::getOperator("+CORE_Integer::toString(index)+")",
                                "index out of bounds [0,"+CORE_Integer::toString(mOperators.getSize())+"[");
        EMM_Operator *op=mOperators[index].get();
        if (op==null) {
            throw EMM_Exception("EMM_LandauLifschitzSystem::getOperator("+CORE_Integer::toString(index)+")",
                                "null operator at index");
        }
        return *op;
        
    }
 
    tBoolean isAffine() const;
    
    tBoolean resetToInitialState();

    virtual tBoolean resetOperatorsToInitialState();

   
    
protected:
    tBoolean computeOperatorsFieldsAtTime(const tReal& t,
                                          const tFlag& order,
                                          const EMM_RealArray& sigma,
                                          const EMM_RealField& dM_dt,
                                          const EMM_RealField& M);
    
    tBoolean updateOperatorsAtNextTimeStep(const tReal& dt,
                                           const EMM_RealArray& sigma,
                                           const EMM_RealField& Mt);

    tBoolean initializeOperators(const EMM_RealArray& sigma,
                                 const EMM_RealField& M);
    

    
    
    
    //MAGNETIZATION field methods
    //===========================

public:
    
    inline void computeFunction(const EMM_RealField& M,
                                const EMM_RealField& H,
                                EMM_RealField& F) const {
        mLLFunction->computeFunction(getSigma(),M,H,F);
    }
    void computeGradFunction(const EMM_RealField& M,
                             const EMM_RealField& H,
                             const EMM_RealField& D,
                             const EMM_RealField& gradH,
                             EMM_RealField& gradF) const {
        mLLFunction->computeGradFunction(getSigma(),M,H,D,gradH,gradF);
    }
    void computePartialGradMFunction(const EMM_RealField& M,
                                     const EMM_RealField& H,
                                     const EMM_RealField& D,
                                     EMM_RealField& gradF) const {
        mLLFunction->computePartialGradMFunction(getSigma(),M,H,D,gradF);
    }
  
   

    inline void setInitialMagnetization(const tString& dataFile) {
        mMInitialDataFile=dataFile;
    }
    inline void setInitialMagnetization(const tReal data[3]) {
        mMInitialData->setSize(1);
        EMM_RealField &H=*mMInitialData.get();
        H.setValue(0,data);
        mMInitialDataFile="";
    }
    inline void setInitialMagnetization(const tReal &M0,const tReal& M1,const tReal& M2) {
        mMInitialData->setSize(1);
        EMM_RealField &H=*mMInitialData.get();
        H.setValue(0,(tDimension)0,M0);
        H.setValue(0,(tDimension)1,M1);
        H.setValue(0,(tDimension)2,M2);
        mMInitialDataFile="";
    }
    inline void setInitialMagnetization(const EMM_RealField& data) {
        mMInitialData->copy(data);
        mMInitialDataFile="";
    }

    inline void setEvolutionScheme(const tFlag& f) {
        mEvolutionScheme=f;
    }
    inline const tFlag& getEvolutionScheme() const {
        return mEvolutionScheme;
    }
    inline tBoolean isMagnetizationStatic() const {
        return (mEvolutionScheme==Os);
    }

    inline const EMM_RealField& getMagnetizationField() const {
        return *mM.get();
    }
    
    inline EMM_RealField& getMagnetizationField()  {
        return *mM.get();
    }
     
    inline const EMM_RealField& getMagnetizationFieldTimeDerivative()  const {
        return *mF.get();
    }
    inline EMM_RealField& getMagnetizationFieldTimeDerivative()   {
        return *mF.get();
    }
    
    inline const EMM_RealArray& getSigma() const {
        return *mSigma.get();
    }
    

    inline  void computeMagnetizationFieldTimeDerivative(const EMM_RealField& Mt,
                                                         const EMM_RealField& Ht,
                                                         EMM_RealField& Ft) const {
        computeFunction(Mt,Ht,Ft);
    }


  
    
    //MAGNETIC EXCITATION methods
    //===========================
    
public:
   
    inline EMM_RealField& getWorkingMagneticExcitationField() const {
        return *mHop.get();
    }
    inline const EMM_RealField& getMagneticExcitationField() const {
        return *mH.get();
    }
    inline EMM_RealField& getMagneticExcitationField()  {
        return *mH.get();
    }

    
    void computeMagneticExcitationField(const EMM_RealField& M,
                                        EMM_RealField& H) const;

    void computeMagnetizationExcitationField(const EMM_RealField& M,
                                             EMM_RealField& H) const;
    
    inline const EMM_RealField&  computeMagneticExcitationField()  {
        computeMagneticExcitationField(getMagnetizationField(),
                                       getMagneticExcitationField());
        return getMagneticExcitationField();
    }


    
    tReal computeMagneticExcitationFieldAndEnergies(const EMM_RealField& M,
                                                    EMM_RealField& H,
                                                    tReal Es[]) const;
   
    
   
  

public:
    virtual void computeLandauLifschitzFields(const EMM_RealField& M,
                                              const EMM_RealField& H,
                                              EMM_RealField& F) {
        computeMagnetizationFieldTimeDerivative(M,H,F);
    }
    
    void computeMeanMagnetizationField(tReal* meanM) const;

    void computeMeanField(const EMM_RealField& F,tReal *meanF) const;
    
    
    //ENERGY methods
    //==============
    inline static tUSInt getEnergiesNumber() {
        return MAX_OPERATORS+1;
    }
    inline const tReal& getEnergy() const {
        return *mEs;
    }
    inline const tReal& getEnergyTimeDerivative(tReal& dE_dt) const {
        dE_dt=mdE_dt;
        return *mEs;
    }
   
   
    inline void setMinimumEnergyVariation(const tReal& v) {
        mMinEnergyVariation=v;
    }
    inline const tReal&  getMinimumEnergyVariation() const {
        return mMinEnergyVariation;
    }
    tReal computeEnergy() const;
   
    inline tReal computeEnergyTimeDerivative(const EMM_RealField& F,
                                             const EMM_RealField& H) const {
        return (-2*getMesh().getAdimensionizedVolume()*F.dot(getSigma(),H));
    }
   
  
protected:
    inline tReal* getEnergies() {
        return mEs;
    }
    inline tReal& getEnergyTimeDerivative() {
        return mdE_dt;
    }
    
    //stepper methods
    //=============
public:

    
    inline void setMinimumTimeStep(const tReal& v) {
        mMinDt=v;
    }
    inline const tReal&  getMinimumTimeStep() const {
        return mMinDt;
    }
    virtual const tUInteger& getTimeStepsNumber() const=0;
    
    virtual const tReal& getTime() const=0;
    
    virtual const tReal& getTimeStep() const=0;
 
    
public:
    inline tBoolean setStepper(SP::EMM_Stepper time) {
        //verify time is not null
        if (time.get()==null) {
            resetStepper();
            return false;
        }
        //the time is already set
        if (mStepper.get()==time.get()) return true;
        
        //reset the time relation
        resetStepper();
        //set the time relation
        if (time->setSystem(getThis())) {
            mStepper=time;
            //set the inverse relation
            mStepper->setSystem(getThis());
            return true;
        }
        return false;
    }

    inline void resetStepper() {
        //do nothing if the time is null
        if (mStepper.get()==null) return;
        //store the time of this
        SP::EMM_Stepper oldTime=mStepper;
        //set to null the time of thies
        mStepper.reset();
        //reset the old time inverse time system relation
        oldTime->resetSystem();
        
    }

    
    inline const EMM_Stepper& getStepper() const {
        return *mStepper.get();
    }
    inline EMM_Stepper& getStepper()  {
        return *mStepper.get();
    }

    //RELAXATION methods
    //===================


public:
   
   

    virtual void integrate(const tBoolean& isRestoring,tSInt& retCode)=0;
    
    virtual tString getInformation(const tUSInt& retCode) {
        return RET_STRING[retCode-MIN_CODE];
    }
    
    //BACKUP & LOG methods
    //====================
public:

    inline void setIsLogPrinted(const tBoolean& v) {
        mIsLogPrinted=v;
    }
    inline const tBoolean& isLogPrinted() const {
        return mIsLogPrinted;
    }
    
    void printLog(const tUInteger& timeStep,
                  const tReal& t,
                  const map<tString,tString>& stepperData,
                  const tUIndex& nMagnetizedElements,
                  const tReal Es[],
                  const tReal& dE_dt,
                  const tReal& DeltaE_DeltaT,
                  const tReal *meanM,
                  const tReal *meanH)  const;
    
    inline void setOutputPath(const tString& currentPath,const tString& path,const tString& prefix) {

        //set the output path
        mOutputPath=path;
        if (mOutputPath.rfind(CORE_File::PATH_SEPARATOR)+1!=mOutputPath.length()) {
            mOutputPath+=CORE_File::PATH_SEPARATOR;
        }
         //verify the path terminates by a path separator
        if (!CORE_File::isAbsolutePath(mOutputPath))  {
            //the output path is not absolute
            if (currentPath.rfind(CORE_File::PATH_SEPARATOR)+1!=currentPath.length()) {
                mOutputPath=currentPath+CORE_File::PATH_SEPARATOR+mOutputPath;
            } else {
                mOutputPath=currentPath+mOutputPath;
            }
        }
        
        mPrefix=prefix;
    }
    
    
    inline void setBackupStep(const tUInteger& n) {
        mBackupStep=n;
    }
    inline void setBackupsNumber(const tUInteger& n) {
        if (n==0) mBackupsNumber=1;
        else mBackupsNumber=n;
    }

   
    inline void setIsBackupEnabled(const tBoolean& v) {
        mIsBackupEnabled=v;
    }
     
    inline const tBoolean& isBackupEnabled() const {
        return mIsBackupEnabled;
    }
    
    inline const tUInteger& getBackupsNumber() const {
        return mBackupsNumber;
    }
    inline const tUInteger& getBackupStep() const {
        return mBackupStep;
    }
    
   

    

   
    
public:
    virtual tString toString() const;


};

#endif