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