freecarrier.hpp

Go to the documentation of this file.

/*
 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
 * Copyright (c) 2022 Lodz University of Technology
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 */
#ifndef PLASK__SOLVER__GAIN_FREECARRIER_FREECARRIER_H
#define PLASK__SOLVER__GAIN_FREECARRIER_FREECARRIER_H
 
#include <plask/plask.hpp>
 
#include <boost/math/tools/roots.hpp>
using boost::math::tools::toms748_solve;
 
namespace plask { namespace gain { namespace freecarrier {
 
template <typename BaseT> struct GainSpectrum;
 
extern OmpNestedLock gain_omp_lock;
 
template <typename BaseT> struct PLASK_SOLVER_API FreeCarrierGainSolver : public BaseT {
    enum WhichLevel : size_t { EL = EnergyLevels::ELECTRONS, HH = EnergyLevels::HEAVY_HOLES, LH = EnergyLevels::LIGHT_HOLES };
 
    typedef typename BaseT::SpaceType GeometryType;
 
    enum { DIM = GeometryType::DIM };
 
    struct ActiveRegionParams;
 
    struct Level {
        double E;           
        Tensor2<double> M;  
        double thickness;   
 
        Level(double E, const Tensor2<double>& M) : E(E), M(M) {}
 
        Level(double E, const Tensor2<double>& M, double t) : E(E), M(M), thickness(t) {}
 
        // TODO Better effective thickness computation (based on wavefunction)?
        Level(double E, const Tensor2<double>& M, WhichLevel which, const ActiveRegionParams& params);
    };
 
    struct ActiveRegionInfo {
        shared_ptr<StackContainer<DIM>> layers;  
        Vec<DIM> origin;                         
 
        ActiveRegionInfo(Vec<DIM> origin) : layers(plask::make_shared<StackContainer<DIM>>()), origin(origin) {}
 
        size_t size() const { return layers->getChildrenCount(); }
 
        shared_ptr<Material> getLayerMaterial(size_t n) const {
            auto block = static_cast<Block<DIM>*>(static_cast<Translation<DIM>*>(layers->getChildNo(n).get())->getChild().get());
            if (auto m = block->singleMaterial()) return m;
            throw plask::Exception("freeCarrierGainSolver requires solid layers.");
        }
 
        typename Primitive<DIM>::Box getLayerBox(size_t n) const {
            return static_cast<GeometryObjectD<DIM>*>(layers->getChildNo(n).get())->getBoundingBox() + origin;
        }
 
        bool isQW(size_t n) const { return static_cast<Translation<DIM>*>(layers->getChildNo(n).get())->getChild()->hasRole("QW"); }
 
        typename Primitive<DIM>::Box getBoundingBox() const { return layers->getBoundingBox() + origin; }
 
 
        bool contains(const Vec<DIM>& point) const { return getBoundingBox().contains(point); }
 
        bool inQW(const Vec<DIM>& point) const {
            if (!contains(point)) return false;
            assert(layers->getChildForHeight(point.vert() - origin.vert()));
            return layers->getChildForHeight(point.vert() - origin.vert())->getChild()->hasRole("QW");
        }
 
        double averageNr(double lam, double T, double conc = 0.) const {
            double nr = 0.;
            for (size_t i = 0; i != materials.size(); ++i)
                if (isQW(i)) nr += thicknesses[i] * materials[i]->Nr(lam, T, conc).real();
            return nr / totalqw;
        }
 
        std::vector<shared_ptr<Material>> materials;  
        std::vector<double> thicknesses;              
        std::vector<size_t> wells;                    
 
        double total;    
        double totalqw;  
        double bottom;   
        double top;      
 
        enum ConsideredHoles : unsigned {
            NO_HOLES = 0,
            HEAVY_HOLES = 1,
            LIGHT_HOLES = 2,
            BOTH_HOLES = 3
        } holes;  
 
        void summarize(const FreeCarrierGainSolver<BaseT>* solver);
    };
 
    struct PLASK_SOLVER_API ActiveRegionParams {
        const ActiveRegionInfo& region;
        std::vector<double> U[3];           
        std::vector<Tensor2<double>> M[3];  
        double Mt;                          
 
        std::vector<Level> levels[3];  
        double Eg;                     
        size_t nhh,                    
            nlh;                       
 
        ActiveRegionParams(const FreeCarrierGainSolver* solver,
                           const ActiveRegionInfo& region,
                           double T,
                           bool quiet = false,
                           double mt = 0.);
 
        ActiveRegionParams(const FreeCarrierGainSolver* solver, const ActiveRegionInfo& region, bool quiet = false, double mt = 0.)
            : ActiveRegionParams(solver, region, solver->T0, quiet, mt) {}
 
        explicit ActiveRegionParams(const FreeCarrierGainSolver* solver, const ActiveRegionParams& ref, double T, bool quiet = true)
            : ActiveRegionParams(solver, ref.region, T, quiet, ref.Mt) {
            nhh = ref.nhh;
            nlh = ref.nlh;
            for (size_t which = 0; which < 3; ++which) {
                double shift = delta(WhichLevel(which), ref);
                levels[which].reserve(ref.levels[which].size());
                for (Level level : ref.levels[which]) levels[which].emplace_back(level.E + shift, level.M, level.thickness);
            }
        }
 
        double sideU(WhichLevel which) const { return 0.5 * (U[which][0] + U[which][U[which].size() - 1]); }
 
        Tensor2<double> sideM(WhichLevel which) const { return 0.5 * (M[which][0] + M[which][M[which].size() - 1]); }
 
        double delta(WhichLevel which, const ActiveRegionParams& ref) const {
            assert(U[which].size() == ref.U[which].size());
            double delta = 0;
            for (size_t i = 0; i < U[which].size(); ++i) {
                delta += U[which][i] - ref.U[which][i];
            }
            return delta / double(U[which].size());
        }
    };
 
    std::vector<ActiveRegionInfo> regions;
 
    ReceiverFor<Temperature, GeometryType> inTemperature;
 
    ReceiverFor<CarriersConcentration, GeometryType> inCarriersConcentration;
 
    ReceiverFor<BandEdges, GeometryType> inBandEdges;
 
    ReceiverFor<FermiLevels, GeometryType> inFermiLevels;
 
    typename ProviderFor<Gain, GeometryType>::Delegate outGain;
 
    typename ProviderFor<EnergyLevels, GeometryType>::Delegate outEnergyLevels;
 
    FreeCarrierGainSolver(const std::string& name = "");
 
    virtual ~FreeCarrierGainSolver();
 
    void loadConfiguration(plask::XMLReader& reader, plask::Manager& manager) override;
 
  protected:
    shared_ptr<Material> substrateMaterial;
 
    bool explicitSubstrate = false;
 
    virtual void detectActiveRegions() = 0;
 
    double level(WhichLevel which, double E, const ActiveRegionParams& params, size_t start, size_t stop) const;
 
    double level(WhichLevel which, double E, const ActiveRegionParams& params) const {
        return level(which, E, params, 0, params.region.materials.size() - 1);
    }
 
    double level(WhichLevel which, double E, const ActiveRegionParams& params, size_t well) const {
        return level(which, E, params, params.region.wells[well], params.region.wells[well + 1]);
    }
 
    void estimateWellLevels(WhichLevel which, ActiveRegionParams& params, size_t qw) const;
 
    void estimateAboveLevels(WhichLevel which, ActiveRegionParams& params) const;
 
    template <class F>
    std::pair<double, double> fermi_bracket_and_solve(F f, double guess, double step, boost::uintmax_t& max_iter) const {
        double a = guess - 0.5 * step, b = guess + 0.5 * step;
        double fa = f(a), fb = f(b);
        if (fa == 0.) return std::make_pair(a, a);
        if (fb == 0.) return std::make_pair(b, b);
        boost::uintmax_t count = max_iter - 1;
        if ((fa < 0.) == (fa < fb)) {
            while ((fa < 0.) == (fb < 0.)) {
                if (count == 0) return std::make_pair(a, b);
                a = b;
                fa = fb;
                b += step;
                fb = f(b);
                if (fb == 0.) return std::make_pair(b, b);
                --count;
            }
        } else {
            while ((fb < 0.) == (fa < 0.)) {
                if (count == 0) return std::make_pair(a, b);
                b = a;
                fb = fa;
                a -= step;
                fa = f(a);
                if (fa == 0.) return std::make_pair(a, a);
                --count;
            }
        }
        auto res = toms748_solve(
            f, a, b, fa, fb, [this](double l, double r) { return r - l < levelsep; }, count);
        max_iter += count;
        return res;
    }
 
#ifndef NDEBUG
  public:
#endif
    double getN(double F, double T, const ActiveRegionParams& params) const;
 
    double getP(double F, double T, const ActiveRegionParams& params) const;
 
  protected:
    double lifetime;    
    double matrixelem;  
 
    double T0;  
 
  protected:
    double levelsep;  
 
    bool strained;  
 
    void estimateLevels();
 
    void onInitialize() override;
 
    void onInvalidate() override;
 
    void onInputChange(ReceiverBase&, ReceiverBase::ChangeReason) {
        outGain.fireChanged();  // the input changed, so we inform the world that everybody should get the new gain
    }
 
    virtual const LazyData<Tensor2<double>> getGainData(Gain::EnumType what,
                                                        const shared_ptr<const MeshD<DIM>>& dst_mesh,
                                                        double wavelength,
                                                        InterpolationMethod interp = INTERPOLATION_DEFAULT) = 0;
 
    virtual const LazyData<std::vector<double>> getEnergyLevels(EnergyLevels::EnumType which,
                                                                const shared_ptr<const MeshD<DIM>>& dst_mesh,
                                                                InterpolationMethod interp = INTERPOLATION_DEFAULT) = 0;
 
  public:
    std::vector<ActiveRegionParams> params0;  
 
    bool quick_levels;  
 
#ifndef NDEBUG
    double detEl(double E, const ActiveRegionParams& params, size_t well = 0) {
        if (well)
            return level(EL, E, params, well - 1);
        else
            return level(EL, E, params);
    }
    double detHh(double E, const ActiveRegionParams& params, size_t well = 0) {
        if (well)
            return level(HH, E, params, well - 1);
        else
            return level(HH, E, params);
    }
    double detLh(double E, const ActiveRegionParams& params, size_t well = 0) {
        if (well)
            return level(LH, E, params, well - 1);
        else
            return level(LH, E, params);
    }
#endif
 
    shared_ptr<Material> getSubstrate() const { return substrateMaterial; }
 
    void setSubstrate(shared_ptr<Material> material) {
        bool invalid = substrateMaterial != material;
        substrateMaterial = material;
        explicitSubstrate = bool(material);
        if (invalid) this->invalidate();
    }
 
    void findFermiLevels(double& Fc, double& Fv, double n, double T, const ActiveRegionParams& params) const;
 
    Tensor2<double> getGain0(double hw, double Fc, double Fv, double T, double nr, const ActiveRegionParams& params) const;
 
    Tensor2<double> getGain(double hw, double Fc, double Fv, double T, double nr, const ActiveRegionParams& params) const;
 
    double getT0() const { return T0; }
    void setT0(double T) {
        T0 = T;
        this->invalidate();
    }
 
    double getLifeTime() const { return lifetime; }
    void setLifeTime(double iLifeTime) { lifetime = iLifeTime; }
 
    double getMatrixElem() const { return matrixelem; }
    void setMatrixElem(double iMatrixElem) {
        matrixelem = iMatrixElem;
        this->invalidate();
    }
 
    bool getStrained() const { return strained; }
    void setStrained(bool value) {
        strained = value;
        this->invalidate();
    }
 
    friend struct GainSpectrum<BaseT>;
 
    shared_ptr<GainSpectrum<BaseT>> getGainSpectrum(const Vec<DIM>& point) {
        this->initCalculation();
        return make_shared<GainSpectrum<BaseT>>(this, point);
    }
 
    typedef GainSpectrum<BaseT> GainSpectrumType;
};
 
template <typename BaseT> struct GainSpectrum {
    typedef typename FreeCarrierGainSolver<BaseT>::ActiveRegionParams ActiveRegionParams;
 
    FreeCarrierGainSolver<BaseT>* solver;               
    plask::optional<Vec<BaseT::SpaceType::DIM>> point;  
 
    size_t reg;     
    double temp;    
    double conc;    
    double Fc, Fv;  
    std::unique_ptr<ActiveRegionParams> params;
 
    GainSpectrum(FreeCarrierGainSolver<BaseT>* solver, const Vec<BaseT::SpaceType::DIM> point) : solver(solver), point(point) {
        for (size_t i = 0; i != solver->regions.size(); ++i) {
            if (solver->regions[i].contains(point)) {
                reg = i;
                solver->inTemperature.changedConnectMethod(this, &GainSpectrum::onChange);
                solver->inCarriersConcentration.changedConnectMethod(this, &GainSpectrum::onChange);
                temp = solver->inTemperature(plask::make_shared<const OnePointMesh<BaseT::SpaceType::DIM>>(point))[0];
                conc = solver->inCarriersConcentration(CarriersConcentration::PAIRS,
                                                       plask::make_shared<const OnePointMesh<BaseT::SpaceType::DIM>>(point))[0];
                updateParams();
                return;
            };
        }
        throw BadInput(solver->getId(), "point {0} does not belong to any active region", point);
    }
 
    GainSpectrum(FreeCarrierGainSolver<BaseT>* solver, double T, double n, size_t reg)
        : solver(solver), temp(T), conc(n), reg(reg) {
        updateParams();
    }
 
    void updateParams() {
        params.reset(new ActiveRegionParams(solver, solver->params0[reg], temp));
        Fc = Fv = NAN;
        solver->findFermiLevels(Fc, Fv, conc, temp, *params);
    }
 
    void onChange(ReceiverBase&, ReceiverBase::ChangeReason) {
        temp = solver->inTemperature(plask::make_shared<const OnePointMesh<BaseT::SpaceType::DIM>>(*point))[0];
        conc = solver->inCarriersConcentration(CarriersConcentration::PAIRS,
                                               plask::make_shared<const OnePointMesh<BaseT::SpaceType::DIM>>(*point))[0];
        updateParams();
    }
 
    ~GainSpectrum() {
        solver->inTemperature.changedDisconnectMethod(this, &GainSpectrum::onChange);
        solver->inCarriersConcentration.changedDisconnectMethod(this, &GainSpectrum::onChange);
    }
 
    Tensor2<double> getGain(double wavelength) const {
        double nr = solver->regions[reg].averageNr(wavelength, temp, conc);
        return solver->getGain(phys::h_eVc1e9 / wavelength, Fc, Fv, temp, nr, *params);
    }
};
 
}}}  // namespace plask::gain::freecarrier
 
#endif  // PLASK__SOLVER__GAIN_FREECARRIER_FREECARRIER_H