sampling.hpp

Go to the documentation of this file.

// Copyright 2015-2018 The ALMA Project Developers
//
// 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.

#pragma once


#include <vector>
#include <utility>
#include <tuple>
#include <map>
#include <Eigen/Dense>
#include <randutils.hpp>
#include <structures.hpp>
#include <qpoint_grid.hpp>
#include <dos.hpp>
#include <processes.hpp>
#include <deviational_particle.hpp>

namespace alma {
inline double random_dt(double w, randutils::mt19937_rng& rng) {
    if (w == 0.)
        throw value_error("invalid scattering rate");
    double u = rng.uniform(0., 1.);
    return -std::log(u) / w;
}


class Grid_distribution {
private:
    std::vector<double> cumulative;
    randutils::mt19937_rng& rng;

public:
    const std::size_t nqpoints;
    const std::size_t nmodes;
    Grid_distribution(const Gamma_grid& grid, randutils::mt19937_rng& _rng)
        : rng(_rng), nqpoints(grid.nqpoints),
          nmodes(grid.get_spectrum_at_q(0).omega.size()) {
    }


    virtual ~Grid_distribution() {
    }


    void fill_cumulative(const std::vector<double>& p);

    std::array<std::size_t, 2> sample() {
        double u = this->rng.uniform(0., 1.);
        auto pos = std::lower_bound(
                       this->cumulative.begin(), this->cumulative.end(), u) -
                   this->cumulative.begin();

        std::array<std::size_t, 2> nruter;
        nruter[0] = pos % this->nmodes;
        nruter[1] = pos / this->nmodes;
        return nruter;
    }
};

class BE_derivative_distribution : public Grid_distribution {
public:
    BE_derivative_distribution(const Gamma_grid& grid,
                               const Eigen::Ref<const Eigen::ArrayXXd>& w,
                               double T,
                               randutils::mt19937_rng& _rng);
};

class Nabla_T_distribution : public Grid_distribution {
public:
    Nabla_T_distribution(const Gamma_grid& grid,
                         const Eigen::Ref<const Eigen::Vector3d>& nablaT,
                         double T,
                         randutils::mt19937_rng& _rng);
};

class Isothermal_wall_distribution : public Grid_distribution {
public:
    Isothermal_wall_distribution(
        const Gamma_grid& grid,
        double Twall,
        double Teq,
        const Eigen::Ref<const Eigen::Vector3d>& normal,
        randutils::mt19937_rng& _rng);
};

class planar_source_distribution : public Grid_distribution {
public:
    planar_source_distribution(const Gamma_grid& grid,
                               double Tref,
                               const Eigen::Ref<const Eigen::Vector3d>& normal,
                               randutils::mt19937_rng& _rng);
};


class Diffuse_mismatch_distribution {
private:
    std::map<char, std::size_t> Nbranches;
    std::map<char, std::size_t> Ntot;
    std::map<char, double> volume;

    using dos_tuple = std::
        tuple<char, std::size_t, std::size_t, double, double, Gaussian_for_DOS>;
    typedef std::tuple<char, std::size_t, double> lookup_entry;
    std::vector<std::vector<lookup_entry>> lookup_incidentA;
    std::vector<std::vector<lookup_entry>> lookup_incidentB;
    randutils::mt19937_rng& rng;
    double Tref;
    std::vector<dos_tuple> get_modes(
        const Gamma_grid& gridA,
        const Gamma_grid& gridB,
        const Eigen::Ref<const Eigen::Vector3d>& normal,
        double scalebroad);

public:
    Diffuse_mismatch_distribution(
        const Gamma_grid& gridA,
        const Crystal_structure& poscarA,
        const Gamma_grid& gridB,
        const Crystal_structure& poscarB,
        const Eigen::Ref<const Eigen::Vector3d>& normal,
        double scalebroad,
        randutils::mt19937_rng& _rng,
        double Tref = 300.0);
    char reemit(char incidence, D_particle& particle);
};

class Elastic_interface_distribution {
private:
    std::size_t nqA;
    std::size_t nqB;
    std::size_t nmodesA;
    std::size_t nmodesB;
    double VA;
    double VB;
    using dos_tuple =
        std::tuple<char, std::size_t, std::size_t, Gaussian_for_DOS>;
    std::vector<std::vector<double>> cumulative;
    std::vector<std::vector<std::size_t>> allowed;
    Eigen::MatrixXd directions;
    randutils::mt19937_rng& rng;
    std::vector<dos_tuple> get_all_modes(const Gamma_grid& gridA,
                                         const Gamma_grid& gridB,
                                         double scalebroad) const;

public:
    Elastic_interface_distribution(const Gamma_grid& gridA,
                                   const Gamma_grid& gridB,
                                   double scalebroad,
                                   randutils::mt19937_rng& _rng);
    char reemit(char incidence,
                const Eigen::Ref<const Eigen::Vector3d>& normal,
                D_particle& particle);
};

// Simple and general class that allows us to simulate completely
// random elastic scattering.
class Elastic_distribution {
private:
    std::size_t nq;
    std::size_t nmodes;
    using dos_tuple = std::tuple<std::size_t, std::size_t, Gaussian_for_DOS>;
    std::vector<std::vector<double>> cumulative;
    std::vector<std::vector<std::size_t>> allowed;
    Eigen::MatrixXd directions;
    randutils::mt19937_rng& rng;
    std::vector<dos_tuple> get_all_modes(const Gamma_grid& grid,
                                         double scalebroad) const;

public:
    Elastic_distribution(const Gamma_grid& grid,
                         double scalebroad,
                         randutils::mt19937_rng& _rng);
    void scatter(D_particle& particle);

    void scatter(const Eigen::Ref<const Eigen::Vector3d>& normal,
                 const int refsign,
                 D_particle& particle);
};
} // namespace alma