expansion.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_OPTICAL_MODAL_EXPANSION_H
#define PLASK__SOLVER_OPTICAL_MODAL_EXPANSION_H
 
#include <plask/plask.hpp>
 
#ifdef OPENMP_FOUND
#   include <omp.h>
#endif
 
#include "solverbase.hpp"
#include "matrices.hpp"
#include "meshadapter.hpp"
#include "temp_matrix.hpp"
 
namespace plask { namespace optical { namespace modal {
 
struct PLASK_SOLVER_API Expansion {
 
    enum Component {
        E_UNSPECIFIED = 0,  
        E_TRAN = 1,         
        E_LONG = 2          
    };
 
    WhichField which_field;
    InterpolationMethod field_interpolation;
 
    ModalBase* solver;
 
    dcomplex k0;
 
    double lam0;
 
    LazyData<double> temperature;
 
    bool gain_connected;
 
    bool epsilon_connected;
 
    LazyData<Tensor2<double>> gain;
 
    LazyData<Tensor3<dcomplex>> epsilons;
 
    LazyData<double> carriers;
 
    Expansion(ModalBase* solver): solver(solver), k0(NAN), lam0(NAN) {}
 
    virtual ~Expansion() {}
 
    TempMatrix getTempMatrix() {
        size_t N = matrixSize();
        return temporary.get(N, N);
    }
 
  private:
    dcomplex glambda;
 
  protected:
 
    TempMatrixPool temporary;
 
    virtual void beforeLayersIntegrals(dcomplex PLASK_UNUSED(lam), dcomplex PLASK_UNUSED(glam)) {}
 
    virtual void afterLayersIntegrals() {
        temperature.reset();
        gain.reset();
        carriers.reset();
    }
 
    virtual void layerIntegrals(size_t layer, double lam, double glam) = 0;
 
  public:
 
    virtual void beforeGetEpsilon() {
        computeIntegrals();
    }
 
    virtual void afterGetEpsilon() {}
 
    double getLam0() const { return lam0; }
    void setLam0(double lam) {
        if (lam != lam0 && !(isnan(lam0) && isnan(lam))) {
            lam0 = lam;
            solver->recompute_integrals = true;
            solver->clearFields();
        }
    }
    void clearLam0() {
        if (!isnan(lam0)) {
            lam0 = NAN;
            solver->recompute_integrals = true;
            solver->clearFields();
        }
    }
 
    dcomplex getK0() const { return k0; }
    void setK0(dcomplex k) {
        if (k != k0) {
            k0 = k;
            if (k0 == 0.) k0 = 1e-12;
            if (isnan(lam0)) solver->recompute_integrals = true;
            solver->clearFields();
        }
    }
 
    void computeIntegrals() {
        dcomplex lambda = 2e3*PI/k0;
        if (solver->recompute_integrals) {
            dcomplex lam;
            if (!isnan(lam0)) {
                lam = lam0;
                glambda = (solver->always_recompute_gain)? lambda : lam;
            } else{
                lam = glambda = lambda;
            }
            size_t nlayers = solver->lcount;
            std::exception_ptr error;
            beforeLayersIntegrals(lam, glambda);
            PLASK_OMP_PARALLEL_FOR
            for (plask::openmp_size_t l = 0; l < nlayers; ++l) {
                if (error) continue;
                try {
                    layerIntegrals(l, real(lam), real(glambda));
                } catch(...) {
                    #pragma omp critical
                    error = std::current_exception();
                }
            }
            afterLayersIntegrals();
            if (error) std::rethrow_exception(error);
            solver->recompute_integrals = false;
            solver->recompute_gain_integrals = false;
        } else if (solver->recompute_gain_integrals ||
                   (solver->always_recompute_gain && !is_zero(lambda - glambda))) {
            dcomplex lam = isnan(lam0)? lambda : solver->lam0;
            glambda = solver->always_recompute_gain? lambda : lam;
            std::vector<size_t> recomputed_layers;
            size_t nlayers = solver->lcount;
            recomputed_layers.reserve(nlayers);
            for (size_t l = 0; l != nlayers; ++l) if (solver->lgained[l] || solver->lcomputed[l])
                recomputed_layers.push_back(l);
            std::exception_ptr error;
            beforeLayersIntegrals(lam, glambda);
            PLASK_OMP_PARALLEL_FOR
            for (plask::openmp_size_t l = 0; l < recomputed_layers.size(); ++l) {
                if (error) continue;
                try {
                    layerIntegrals(recomputed_layers[l], real(lam), real(glambda));
                } catch(...) {
                    #pragma omp critical
                    error = std::current_exception();
                }
            }
            afterLayersIntegrals();
            if (error) std::rethrow_exception(error);
            solver->recompute_gain_integrals = false;
        }
    }
 
    virtual bool diagonalQE(size_t PLASK_UNUSED(l)) const { return false; }
 
    virtual size_t matrixSize() const = 0;
 
    virtual void getMatrices(size_t layer, cmatrix& RE, cmatrix& RH) = 0;
 
    virtual LazyData<Tensor3<dcomplex>> getMaterialEps(size_t lay,
                                                       const shared_ptr<const typename LevelsAdapter::Level> &level,
                                                       InterpolationMethod interp) = 0;
    virtual void getDiagonalEigenvectors(cmatrix& Te, cmatrix& Te1, const cmatrix& RE, const cdiagonal& gamma);
 
    void initField(WhichField which, InterpolationMethod method) {
        which_field = which;
        field_interpolation = method;
        prepareField();
    }
 
    virtual double integrateField(WhichField field, size_t layer, const cmatrix& TE, const cmatrix& TH,
                                  const std::function<std::pair<dcomplex,dcomplex>(size_t, size_t)>& vertical) = 0;
 
    virtual double integratePoyntingVert(const cvector& E, const cvector& H) = 0;
 
    inline double getModeFlux(size_t n, const cmatrix& TE, const cmatrix& TH) {
        const cvector E(const_cast<dcomplex*>(TE.data()) + n*TE.rows(), TE.rows()),
                      H(const_cast<dcomplex*>(TH.data()) + n*TH.rows(), TH.rows());
        return integratePoyntingVert(E, H);
    }
 
  protected:
 
    virtual void prepareField() {}
 
  public:
    virtual void cleanupField() {}
 
    virtual LazyData<Vec<3,dcomplex>> getField(size_t l,
                                               const shared_ptr<const typename LevelsAdapter::Level>& level,
                                               const cvector& E,
                                               const cvector& H) = 0;
};
 
 
}}} // namespace plask
 
#endif // PLASK__SOLVER_OPTICAL_MODAL_EXPANSION_H