lab1.cpp

#include <iostream>
#include <vector>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <omp.h>
#include <chrono>

using namespace std;

int main()
{
    auto start = std::chrono::high_resolution_clock::now();
    srand(42);

    const int numParticles = 10000000;
    const int gridRows = 50000;
    const int gridCols = 50000;

    vector<double> velocity(numParticles),
        pressure(numParticles),
        energy(numParticles);

// 1) init velocity & pressure
#pragma omp parallel for
    for (int i = 0; i < numParticles; i++)
    {
        velocity[i] = i * 1.0;
        pressure[i] = (numParticles - i) * 1.0;
    }

// 2) compute energy
#pragma omp parallel for
    for (int i = 0; i < numParticles; i++)
    {
        energy[i] = velocity[i] + pressure[i];
    }

// 3) work loop + update velocity
#pragma omp parallel for
    for (int i = 0; i < numParticles; i++)
    {
        double work = 0.0;
        int loops = (i % 10) * 10 + 10; // 10,20,…,100
        for (int j = 0; j < loops; j++)
        {
            work += sin(j * 0.01);
        }
        velocity[i] = sin(energy[i]) + log1p(fabs(work));
    }

    // 4) fieldSum over 2D grid
    double fieldSum = 0.0;
#pragma omp parallel for collapse(2) reduction(+ : fieldSum)
    for (int r = 0; r < gridRows; r++)
    {
        for (int c = 0; c < gridCols; c++)
        {
            fieldSum += sqrt(r * 2.0) + log1p(c * 2.0);
        }
    }

    // 5) atomicFlux
    double atomicFlux = 0.0;
#pragma omp parallel for reduction(+ : atomicFlux)
    for (int i = 0; i < numParticles; i++)
    {
        atomicFlux += velocity[i] * 1e-6;
    }

    // 6) criticalFlux (remains serial)
    double criticalFlux = 0.0;
    vector<double> A(numParticles), B(numParticles);
#pragma omp parallel for
    for (int i = 0; i < numParticles; i++)
    {
        double tempVal = sqrt(fabs(energy[i])) / 100.0;
        double extraVal = log1p(fabs(velocity[i])) * 0.01;
        A[i] = tempVal + extraVal;
        B[i] = sqrt(tempVal) - extraVal;
    }

    for (int i = 0; i < numParticles; i++)
    {
        if (criticalFlux < 500.0)
        {
            criticalFlux += A[i];
        }
        else
        {
            criticalFlux += B[i];
        }
    }

    // 7) final sums
    double sumVelocity = 0.0, sumPressure = 0.0, sumEnergy = 0.0;
#pragma omp parallel for reduction(+ : sumVelocity, sumPressure, sumEnergy)
    for (int i = 0; i < numParticles; i++)
    {
        sumVelocity += velocity[i];
        sumPressure += pressure[i];
        sumEnergy += energy[i];
    }

    // print results
    cout << "=== result ===\n";
    cout << "fieldValue      = " << fieldSum << "\n";
    cout << "energy[0]       = " << energy[0] << "\n";
    cout << "Sum(velocity)   = " << sumVelocity << "\n";
    cout << "Sum(pressure)   = " << sumPressure << "\n";
    cout << "Sum(energy)     = " << sumEnergy << "\n";
    cout << "atomicFlux      = " << atomicFlux << "\n";
    cout << "criticalFlux    = " << criticalFlux << "\n";

    auto end = std::chrono::high_resolution_clock::now();
    chrono::duration<double, milli> elapsed = end - start;
    cout << "程式執行時間: " << elapsed.count() << " 毫秒\n";

    return 0;
}