#include <cstdlib>
#include <iostream>

#include "palabos2D.h"
#include "palabos2D.hh"

using namespace plb;
using namespace std;

// 将 double 类型定义为一个新的类型别名 T。可以使用 T 代替 double，从而提高代码的可读性和可维护性。
typedef double T; 

#define DESCRIPTOR descriptors::ShanChenD2Q9Descriptor //定义二维9速度ShanChen描述符

// 定义一个模板类，该模板有两个参数，T 是一个类型参数。Descriptor 是一个模板类，需要一个类型参数 U。
template <typename T, template <typename U> class Descriptor> 
class RandomInitializer : public BoxProcessingFunctional2D_L<T, Descriptor> {
public:
    RandomInitializer(T rho0_, T maxRho_, Box2D const &boundingBox_, uint32_t const *seed_) : // 初始化函数
        rho0(rho0_), maxRho(maxRho_), nY(boundingBox_.getNy()), seed(seed_) {};
    virtual void process(Box2D domain, BlockLattice2D<T, Descriptor> &lattice)
    {
        Dot2D location = lattice.getLocation();
        sitmo::prng_engine eng(*seed);
        plint rng_index = 0;
        for (plint iX = domain.x0; iX <= domain.x1; ++iX) {
            plint globalX = nY * (iX + location.x);
            for (plint iY = domain.y0; iY <= domain.y1; ++iY) {
                plint globalY = iY + location.y + globalX;
                PLB_ASSERT(globalY >= rng_index);
                if (globalY > rng_index) {
                    eng.discard(globalY - rng_index);
                    rng_index = globalY;
                }
                T randNumber = (T)eng() / (T)sitmo::prng_engine::max();
                ++rng_index;

                T rho = rho0 + randNumber * maxRho;
                Array<T, 2> zeroVelocity(0., 0.);

                iniCellAtEquilibrium(lattice.get(iX, iY), rho, zeroVelocity); // 随机二维密度

                lattice.get(iX, iY).setExternalField(
                    Descriptor<T>::ExternalField::densityBeginsAt,
                    Descriptor<T>::ExternalField::sizeOfDensity, &rho);
                lattice.get(iX, iY).setExternalField(
                    Descriptor<T>::ExternalField::momentumBeginsAt,
                    Descriptor<T>::ExternalField::sizeOfMomentum, &zeroVelocity[0]);
            }
        }
    };
    virtual RandomInitializer<T, Descriptor> *clone() const
    {
        return new RandomInitializer<T, Descriptor>(*this);
    };
    virtual void getTypeOfModification(std::vector<modif::ModifT> &modified) const
    {
        modified[0] = modif::staticVariables;
    };

private:
    T rho0, maxRho;
    plint nY;
    uint32_t const *seed;
};

int main(int argc, char *argv[]) // 主程序
{
    plbInit(&argc, &argv); // 初始化函数，用于设置和配置 Palabos 环境
    global::directories().setOutputDir("./tmp/"); // 设置输出文件夹

    uint32_t seed = 1;

    // For the choice of the parameters G, rho0, and psi0, we refer to the book
    //   Michael C. Sukop and Daniel T. Thorne (2006),
    //   Lattice Boltzmann Modeling; an Introduction for Geoscientists and Engineers.
    //   Springer-Verlag Berlin/Heidelberg.

    const T omega = 1.0; // 松弛参数？
    const int nx = 400; // x方向格子数目
    const int ny = 400; // y方向格子数目
    const T G = -120; // G为相互作用强度
#ifndef PLB_REGRESSION
    const int maxIter = 100001;
    const int saveIter = 100;
#else
    const int maxIter = 1001;
#endif
    const int statIter = 100;

    const T rho0 = 200.0; // 初始密度
    const T deltaRho = 1.0; // 用在了随机初始化密度上的参数
    const T psi0 = 4.0; // 粒子间势初始值

    MultiBlockLattice2D<T, DESCRIPTOR> lattice( // 拥有外力（色散力？）的BGKdynamics 二维晶格
        nx, ny, new ExternalMomentBGKdynamics<T, DESCRIPTOR>(omega)); // 为了避免计算变量两次，确保在碰撞过程中从外部标量中读取宏观变量

    lattice.periodicity().toggleAll(true); // 所有方向采用周期性边界

    // 使用随机初始条件，激活相分离。
    applyProcessingFunctional(
        new RandomInitializer<T, DESCRIPTOR>(rho0, deltaRho, lattice.getBoundingBox(), &seed),
        lattice.getBoundingBox(), lattice);

    // 添加数据处理器，实现Shan/Chen相互作用势。
    plint processorLevel = 1;
    integrateProcessingFunctional(
        new ShanChenSingleComponentProcessor2D<T, DESCRIPTOR>( // 二维单相ShanChen
            G, new interparticlePotential::PsiShanChen94<T>(psi0, rho0)), // 采用PsiShanChen94(Psi0,rho0)粒子间势
        lattice.getBoundingBox(), lattice, processorLevel);

    lattice.initialize(); // 初始化格子

    pcout << "Starting simulation" << endl;
    for (int iT = 0; iT < maxIter; ++iT) {
        if (iT % statIter == 0) { // 输出到屏幕
            std::unique_ptr<MultiScalarField2D<T> > rho(computeDensity(lattice)); // 计算密度场
            pcout << iT << ": Average rho fluid one = " << computeAverage(*rho) << endl; // 平均密度
            pcout << "Minimum density: " << computeMin(*rho) << endl; // 最小密度
            pcout << "Maximum density: " << computeMax(*rho) << endl; // 最大密度
        }
#ifndef PLB_REGRESSION
        if (iT % saveIter == 0) { // 保存图片
            ImageWriter<T>("leeloo").writeScaledGif(
                createFileName("rho", iT, 6), *computeDensity(lattice));
        }
#endif

        lattice.collideAndStream(); // 每个算例里都会有这一行代码, 碰撞和流动
    }
}