solver.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_SOLVER_H
#define PLASK__SOLVER_OPTICAL_MODAL_SOLVER_H
 
#include <plask/plask.hpp>
#include "expansion.hpp"
#include "rootdigger.hpp"
#include "solverbase.hpp"
#include "transfer.hpp"
 
#undef interface
 
namespace plask { namespace optical { namespace modal {
 
template <typename BaseT> class PLASK_SOLVER_API ModalSolver : public BaseT, public ModalBase {
    void onInputChanged(ReceiverBase&, ReceiverBase::ChangeReason) {
        this->clearModes();
        this->recompute_integrals = true;
    }
 
    void onGainChanged(ReceiverBase&, ReceiverBase::ChangeReason reason) {
        if (reason == ReceiverBase::ChangeReason::REASON_VALUE) {
            this->clearModes();
            this->recompute_gain_integrals = true;
        } else {
            this->invalidate();
        }
    }
 
  protected:
    void onInitialize() override { setupLayers(); }
 
    void onGeometryChange(const Geometry::Event& evt) override {
        BaseT::onGeometryChange(evt);
        if (this->geometry) {
            if (evt.flags() == 0) {
                auto objects = this->geometry->getObjectsWithRole("interface");
                if (objects.size() > 1) {
                    Solver::writelog(LOG_WARNING, "More than one object with 'interface' role: interface not set");
                } else if (objects.size() == 1) {
                    setInterfaceOn(objects[0]);
                }
            }
        } else {
            vbounds->clear();
        }
    }
 
    void parseCommonModalConfiguration(XMLReader& reader, Manager& manager);
 
    void setupLayers();
 
    double smooth;
 
  public:
    ReceiverFor<Temperature, typename BaseT::SpaceType> inTemperature;
 
    ReceiverFor<Gain, typename BaseT::SpaceType> inGain;
 
    ReceiverFor<Epsilon, typename BaseT::SpaceType> inEpsilon;
 
    ReceiverFor<CarriersConcentration, typename BaseT::SpaceType> inCarriersConcentration;
 
    typename ProviderFor<Epsilon, typename BaseT::SpaceType>::Delegate outEpsilon;
 
    typename ProviderFor<RefractiveIndex, typename BaseT::SpaceType>::Delegate outRefractiveIndex;
 
    typename ProviderFor<ModeWavelength>::Delegate outWavelength;
 
    typename ProviderFor<ModeLightMagnitude, typename BaseT::SpaceType>::Delegate outLightMagnitude;
 
    typename ProviderFor<ModeLightE, typename BaseT::SpaceType>::Delegate outLightE;
 
    typename ProviderFor<ModeLightH, typename BaseT::SpaceType>::Delegate outLightH;
 
    typename ProviderFor<ModeLightE, typename BaseT::SpaceType>::Delegate outUpwardsLightE;
 
    typename ProviderFor<ModeLightH, typename BaseT::SpaceType>::Delegate outUpwardsLightH;
 
    typename ProviderFor<ModeLightE, typename BaseT::SpaceType>::Delegate outDownwardsLightE;
 
    typename ProviderFor<ModeLightH, typename BaseT::SpaceType>::Delegate outDownwardsLightH;
 
    ModalSolver(const std::string& name = "");
 
    ~ModalSolver();
 
    double getSmooth() const { return smooth; }
 
    void setSmooth(double value) {
        bool changed = smooth != value;
        smooth = value;
        if (changed) this->invalidate();
    }
 
    Transfer::Method getTransferMethod() const { return transfer_method; }
 
    void setTransferMethod(Transfer::Method method) {
        if (method != transfer_method) this->invalidate();
        transfer_method = method;
    }
 
    bool getGroupLayers() const { return group_layers; }
 
    void setGroupLayers(bool value) {
        bool changed = group_layers != value;
        group_layers = value;
        if (changed) this->invalidate();
    }
 
    double getMaxTempDiff() const { return max_temp_diff; }
 
    void setMaxTempDiff(double value) {
        bool changed = max_temp_diff != value;
        max_temp_diff = value;
        if (changed) this->invalidate();
    }
 
    double getTempDist() const { return temp_dist; }
 
    void setTempDist(double value) {
        bool changed = temp_dist != value;
        temp_dist = value;
        if (changed) this->invalidate();
    }
 
    double getTempLayer() const { return temp_layer; }
 
    void setTempLayer(double value) {
        bool changed = temp_layer != value;
        temp_layer = value;
        if (changed) this->invalidate();
    }
 
    std::string getId() const override { return Solver::getId(); }
 
    bool initCalculation() override { return Solver::initCalculation(); }
 
    inline size_t getInterface() {
        Solver::initCalculation();
        return interface;
    }
 
    inline void setInterfaceAt(double pos) {
        if (pos != interface_position) {
            this->invalidate();
            interface_position = pos;
            Solver::writelog(LOG_DEBUG, "Setting interface at position {:g}um", interface_position);
        }
    }
 
    void setInterfaceOn(const shared_ptr<const GeometryObject>& object, const PathHints* path = nullptr) {
        auto boxes = this->geometry->getObjectBoundingBoxes(object, path);
        if (boxes.size() != 1) throw NotUniqueObjectException();
        if (interface_position != boxes[0].lower.vert()) {
            this->invalidate();
            interface_position = boxes[0].lower.vert();
            Solver::writelog(LOG_DEBUG, "Setting interface on an object at position {:g}um", interface_position);
        }
    }
 
    void setInterfaceOn(const shared_ptr<const GeometryObject>& object, const PathHints& path) { setInterfaceOn(object, &path); }
 
    dcomplex getDeterminant() {
        initCalculation();
        ensureInterface();
        if (!transfer) initTransfer(getExpansion(), false);
        return transfer->determinant();
    }
 
    void prepareExpansionIntegrals(Expansion* expansion,
                                   const shared_ptr<MeshD<BaseT::SpaceType::DIM>>& mesh,
                                   dcomplex lam,
                                   dcomplex glam) {
        expansion->temperature = inTemperature(mesh);
        expansion->gain_connected = inGain.hasProvider();
        expansion->epsilon_connected = inEpsilon.hasProvider();
        if (expansion->gain_connected) {
            if (isnan(glam)) glam = lam;
            expansion->gain = inGain(mesh, real(glam));
        }
        if (expansion->epsilon_connected) {
            if (isnan(glam)) glam = lam;
            expansion->epsilons = inEpsilon(mesh, glam);
        }
        expansion->carriers = inCarriersConcentration.hasProvider() ? inCarriersConcentration(CarriersConcentration::MAJORITY, mesh)
                                                                    : LazyData<double>(mesh->size(), 0.);
    }
 
    cvector incidentVector(Transfer::IncidentDirection side, size_t idx, dcomplex lam = NAN);
 
    cvector incidentVector(Transfer::IncidentDirection side, const cvector& incident, dcomplex lam = NAN);
 
  protected:
    // Initialize the incident field vector
    size_t initIncidence(Transfer::IncidentDirection side, dcomplex lam = NAN);
 
    void scaleIncidentVector(cvector& incident, size_t layer);
 
    virtual size_t nummodes() const = 0;
 
    virtual double applyMode(size_t n) = 0;
 
    DataVector<const Tensor3<dcomplex>> getEpsilonProfile(const shared_ptr<const MeshD<BaseT::SpaceType::DIM>>& dst_mesh,
                                                          dcomplex lam,
                                                          InterpolationMethod interp = INTERPOLATION_DEFAULT);
 
    LazyData<dcomplex> getRefractiveIndex(RefractiveIndex::EnumType component,
                                          const shared_ptr<const MeshD<BaseT::SpaceType::DIM>>& dst_mesh,
                                          dcomplex lam,
                                          InterpolationMethod interp = INTERPOLATION_DEFAULT);
 
    template <PropagationDirection part = PROPAGATION_TOTAL>
    LazyData<Vec<3, dcomplex>> getLightE(size_t num,
                                         shared_ptr<const MeshD<BaseT::SpaceType::DIM>> dst_mesh,
                                         InterpolationMethod method);
 
    template <PropagationDirection part = PROPAGATION_TOTAL>
    LazyData<Vec<3, dcomplex>> getLightH(size_t num,
                                         shared_ptr<const MeshD<BaseT::SpaceType::DIM>> dst_mesh,
                                         InterpolationMethod method);
 
    LazyData<double> getLightMagnitude(size_t num,
                                       shared_ptr<const MeshD<BaseT::SpaceType::DIM>> dst_mesh,
                                       InterpolationMethod method);
 
    virtual double getWavelength(size_t n) = 0;
};
 
}}}  // namespace plask::optical::modal
 
#endif  // PLASK__SOLVER_OPTICAL_MODAL_SOLVER_H