solverbase.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_SOLVERBASE_H
#define PLASK__SOLVER_OPTICAL_MODAL_SOLVERBASE_H
 
#include <plask/plask.hpp>
#include "rootdigger.hpp"
#include "transfer.hpp"
 
#undef interface
 
namespace plask { namespace optical { namespace modal {
 
struct PML {
    dcomplex factor;  
    double size;      
    double dist;      
    double order;     
    PML() : factor(1., 0.), size(1.), dist(0.5), order(1) {}
    PML(dcomplex factor, double size, double dist, double order) : factor(factor), size(size), dist(dist), order(order) {}
};
 
struct PLASK_SOLVER_API ModalBase {
    enum Emission {
        EMISSION_UNSPECIFIED = 0,  
        EMISSION_TOP,              
        EMISSION_BOTTOM,           
        EMISSION_FRONT,            
        EMISSION_BACK              
    };
 
    Emission emission;
 
    Transfer::Determinant determinant_type;
 
  protected:
    std::unique_ptr<RootDigger> getRootDigger(const RootDigger::function_type& func, const char* name);
 
    Transfer::Method transfer_method;
 
    void readRootDiggerConfig(XMLReader& reader) {
        root.tolx = reader.getAttribute<double>("tolx", root.tolx);
        root.tolf_min = reader.getAttribute<double>("tolf-min", root.tolf_min);
        root.tolf_max = reader.getAttribute<double>("tolf-max", root.tolf_max);
        root.maxstep = reader.getAttribute<double>("maxstep", root.maxstep);
        root.maxiter = reader.getAttribute<int>("maxiter", root.maxiter);
        root.alpha = reader.getAttribute<double>("alpha", root.alpha);
        root.lambda_min = reader.getAttribute<double>("lambd", root.lambda_min);
        root.initial_dist = reader.getAttribute<dcomplex>("initial-range", root.initial_dist);
        root.method = reader.enumAttribute<RootDigger::Method>("method")
                          .value("brent", RootDigger::ROOT_BRENT)
                          .value("broyden", RootDigger::ROOT_BROYDEN)
                          .value("muller", RootDigger::ROOT_MULLER)
                          .get(root.method);
        reader.requireTagEnd();
    }
 
  public:
    std::unique_ptr<Transfer> transfer;
 
    void initTransfer(Expansion& expansion, bool reflection);
 
    shared_ptr<OrderedAxis> vbounds;
 
    shared_ptr<OrderedAxis> verts;
 
    size_t lcount;
 
    std::vector<bool> lgained;
 
    std::vector<bool> lcomputed;
 
    std::vector<std::size_t> stack;
 
    std::ptrdiff_t interface;
 
    double interface_position;
 
    double lam0;
 
    dcomplex k0;
 
    PML vpml;
 
    RootDigger::Params root;
 
    bool recompute_integrals;
 
    bool recompute_gain_integrals;
 
    bool always_recompute_gain;
 
  protected:
    bool group_layers;
 
    double max_temp_diff;
 
    double temp_dist;
 
    double temp_layer;
 
    void scaleIncidentVector(cvector& incident, size_t layer, double size_factor);
 
  public:
    ModalBase()
        : emission(EMISSION_UNSPECIFIED),
          determinant_type(Transfer::DETERMINANT_EIGENVALUE),
          transfer_method(Transfer::METHOD_AUTO),
          interface(-1),
          interface_position(NAN),
          lam0(NAN),
          k0(NAN),
          vpml(dcomplex(1., -2.), 2.0, 10., 0),
          recompute_integrals(true),
          always_recompute_gain(true),
          group_layers(true),
          max_temp_diff(NAN),
          temp_dist(0.5),
          temp_layer(0.05) {}
 
    virtual ~ModalBase() {}
 
    double getLam0() const { return lam0; }
    void setLam0(double lam) {
        lam0 = lam;
        if (!isnan(lam) && (isnan(real(k0)) || isnan(imag(k0)))) k0 = 2e3 * PI / lam;
    }
    void clearLam0() { lam0 = NAN; }
 
    dcomplex getK0() const { return k0; }
 
    void setK0(dcomplex k) {
        k0 = k;
        if (k0 == 0.) k0 = 1e-12;
    }
 
    dcomplex getLam() const { return 2e3 * PI / k0; }
 
    void setLam(dcomplex lambda) { k0 = 2e3 * PI / lambda; }
 
    void clearFields() {
        if (transfer) transfer->fields_determined = Transfer::DETERMINED_NOTHING;
    }
 
    size_t getLayerFor(double& h) const {
        size_t n = vbounds->findUpIndex(h + 1e-15);
        if (n == 0)
            h -= vbounds->at(0);
        else
            h -= vbounds->at(n - 1);
        return n;
    }
 
    virtual std::string getId() const = 0;
 
    virtual bool initCalculation() = 0;
 
    const std::vector<std::size_t>& getStack() const { return stack; }
 
    virtual void computeIntegrals() = 0;
 
    virtual void clearModes() = 0;
 
    virtual bool setExpansionDefaults(bool with_k0 = true) = 0;
 
    void ensureInterface() {
        if (interface == -1) throw BadInput(this->getId(), "no interface position set");
        if (interface == 0 || interface >= std::ptrdiff_t(stack.size()))
            throw BadInput(this->getId(), "wrong interface position {0} (min: 1, max: {1})", interface, stack.size() - 1);
    }
 
    virtual Expansion& getExpansion() = 0;
 
    dvector getIncidentFluxes(const cvector& incident, Transfer::IncidentDirection side);
 
    dvector getReflectedFluxes(const cvector& incident, Transfer::IncidentDirection side);
 
    dvector getTransmittedFluxes(const cvector& incident, Transfer::IncidentDirection side);
 
    cvector getReflectedCoefficients(const cvector& incident, Transfer::IncidentDirection side);
 
    cvector getTransmittedCoefficients(const cvector& incident, Transfer::IncidentDirection side);
 
    double getReflection(const cvector& incident, Transfer::IncidentDirection side) {
        double R = 0.;
        for (double r : getReflectedFluxes(incident, side)) R += r;
        return R;
    }
 
    double getTransmission(const cvector& incident, Transfer::IncidentDirection side) {
        double T = 0.;
        for (double t : getTransmittedFluxes(incident, side)) T += t;
        return T;
    }
 
#ifndef NDEBUG
    void getMatrices(size_t layer, cmatrix& RE, cmatrix& RH);
#endif
};
 
}}}  // namespace plask::optical::modal
 
#endif  // PLASK__SOLVER_OPTICAL_MODAL_SOLVERBASE_H