TimeEvolver.h

/* Copyright 2016 Kristofer Björnson
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
 
 
#ifndef COM_DAFER45_TBTK_TIME_EVOLVER
#define COM_DAFER45_TBTK_TIME_EVOLVER
 
#include "TBTK/Solver/Diagonalizer.h"
#include "TBTK/Model.h"
#include "TBTK/UnitHandler.h"
 
#include <complex>
#include <vector>
 
namespace TBTK{
namespace Solver{
 
class TimeEvolver : public Solver{
public:
    TimeEvolver();
 
    virtual ~TimeEvolver();
 
    void run();
 
    void setCallback(bool (*callback)(TimeEvolver *timeEvolver));
 
    void setMaxSCIterations(int maxIterations);
 
    void setNumTimeSteps(int numTimeSteps);
 
    void setTimeStep(double dt);
 
    void setNumberOfParticles(int numberOfParticles);
 
    int getNumberOfParticles();
 
    void fixParticleNumber(bool particleNumberIsFixed);
 
    double getOccupancy(int state);
 
    Diagonalizer *getDiagonalizer();
 
    double getEigenValue(int state);
 
    const std::complex<double> getAmplitude(int state, const Index &index);
 
    enum class DecayMode{None, Instantly, Interpolate, Custom};
 
    void setDecayMode(DecayMode decayMode);
 
    class DecayHandler{
    public:
        virtual void decay(
            TimeEvolver *timeEvolver,
            double *occupancy,
            double *eigenValues,
            std::complex<double> **eigenVectorsMap
        ) = 0;
    private:
    };
 
    void setDecayHandler(DecayHandler *decayHandler);
 
    int getCurrentTimeStep();
 
    void setOrthogonalityCheckInterval(int orthogonalityCheckInterval);
 
    double getOrthogonalityError();
private:
    Diagonalizer dSolver;
 
    double *eigenValues;
 
    std::complex<double> *eigenVectors;
 
    std::complex<double> **eigenVectorsMap;
 
    double* occupancy;
 
    int numberOfParticles;
 
    bool particleNumberIsFixed;
 
    DecayMode decayMode;
 
    DecayHandler *decayHandler;
 
    bool (*callback)(TimeEvolver *timeEvolver);
 
    int numTimeSteps;
 
    double dt;
 
    int currentTimeStep;
 
    static std::vector<TimeEvolver*> timeEvolvers;
 
    static std::vector<Diagonalizer*> dSolvers;
 
    class SelfConsistencyCallback :
        public Diagonalizer::SelfConsistencyCallback
    {
    public:
        virtual bool selfConsistencyCallback(Diagonalizer &diagonalizer);
    };
//  static bool selfConsistencyCallback(Diagonalizer &dSolver);
    static SelfConsistencyCallback selfConsistencyCallback;
 
    void onDiagonalizationFinished();
 
    void sort();
 
    void updateOccupancy();
 
    void decayInstantly();
 
    void decayInterpolate();
 
    double orthogonalityError;
 
    int orthogonalityCheckInterval;
 
    void calculateOrthogonalityError();
};
 
inline void TimeEvolver::setCallback(
    bool (*callback)(TimeEvolver *timeEvolver)
){
    this->callback = callback;
}
 
inline void TimeEvolver::setMaxSCIterations(int maxIterations){
    dSolver.setMaxIterations(maxIterations);
}
 
inline void TimeEvolver::setNumTimeSteps(int numTimeSteps){
    this->numTimeSteps = numTimeSteps;
}
 
inline void TimeEvolver::setTimeStep(double dt){
    this->dt = dt;
}
 
inline void TimeEvolver::setNumberOfParticles(int numberOfParticles){
    this->numberOfParticles = numberOfParticles;
}
 
inline int TimeEvolver::getNumberOfParticles(){
    return numberOfParticles;
}
 
inline void TimeEvolver::fixParticleNumber(bool particleNumberIsFixed){
    this->particleNumberIsFixed = particleNumberIsFixed;
}
 
inline Diagonalizer* TimeEvolver::getDiagonalizer(){
    return &dSolver;
}
 
inline double TimeEvolver::getEigenValue(int state){
    return eigenValues[state];
}
 
inline double TimeEvolver::getOccupancy(int state){
    return occupancy[state];
}
 
inline const std::complex<double> TimeEvolver::getAmplitude(
    int state,
    const Index &index
){
    return eigenVectorsMap[state][getModel().getBasisIndex(index)];
}
 
inline void TimeEvolver::setDecayMode(DecayMode decayMode){
    this->decayMode = decayMode;
}
 
inline void TimeEvolver::setDecayHandler(DecayHandler *decayHandler){
    this->decayHandler = decayHandler;
}
 
inline int TimeEvolver::getCurrentTimeStep(){
    return currentTimeStep;
}
 
inline void TimeEvolver::setOrthogonalityCheckInterval(
    int orthogonalityCheckInterval
){
    this->orthogonalityCheckInterval = orthogonalityCheckInterval;
}
 
inline double TimeEvolver::getOrthogonalityError(){
    return orthogonalityError;
}
 
};  //End of namespace Solver
};  //End of namespace TBTK
 
#endif