/home/julian/svn/gpucpumd/mdgpu2/main.cu

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
#include <assert.h>
#include <cstdlib>
#include <limits>
#include <cuda.h>
#include <ctime>
#include "gpu_forces.cuh"
#include "gpu_integrator.cuh"
#include "gpu_tools.cuh"
#include "define.h"
#include "data.h"
#include <cstring>
#include "tests.cuh"


int main(int argc, char *argv[])
{

    /*###########################################
    TEST AREA
    *///#########################################
    //int returnValue=0;
    //outputSteppingTest();
    //LJ_force_test();
    //summation_test();
    //temperature_scale_test();
    //scaling();
    //drift();
    //test2d();
    //test3d();
    //threadtest(argc, argv);
    //precision_test(argc, argv);
    //returnValue = forceSpeed(argc, argv);
    //return returnValue;
    /*##########################################
    END TEST AREA
    *///########################################


    /*##########################################
    Production code
    *///########################################
    if( argc !=15 )
    {
        std::cout<<"Argument parsing failed!\nPlease edit 'example.sh' for your needs!\n";
        std::string exeName = argv[0];
        std::cout<<exeName<<"\n";
        std::cout<<"Number of given arguments: "<<argc<<"\n";
        bool tmp=generateStartScript(exeName);
        return 1;
    }
    std::string inputfile=argv[1]; //inputdata file name
    std::string basename = argv[2]; //output base name
    int dim = atoi(argv[3]); //dimension of the simulation
    int N = atoi(argv[4]); //number of particles given as command line argument
    int numNN = atoi(argv[5]); //maximum number of coupled nearest neighbors
    float T = atof(argv[6]); //Temperature
    float k = atof(argv[7]); //spring constant
    float boxSize = atof(argv[8]); //edge length of the simulation box
    float dt = atof(argv[9]); //timestep
    int NOffset = atoi(argv[10]);//simulation step count offset
    int NSim = atoi(argv[11]);//final simulation step count
    int NOut = atoi(argv[12]);//output intervall
    int NAdjust =       atoi(argv[13]);//NVT velocity scaling and drift removal intervall
    int gpuDevNo = atoi(argv[14]); //use specific gpu device

    bool success = true;//state flag
    bool linLogEqual = false;
    bool NthreadEqN = true;

    //set device
    success=setCudaDevice(gpuDevNo);
    if(!success)
    {
        return 1;
    }

    //generate script to continue finished simulation
    //success = generateRestartScript(argv[0], basename, dim, N, numNN, T, k, boxSize, dt, NOffset+NSim, NOffset+2*NSim, NOut, NAdjust, gpuDevNo);

    //static simulation parameters
    float rCut = 3. * pow(2., 1./6.);
    float rCutSq = rCut * rCut;
    float rVer = 1.5 * rCut;
    float rVerSq = rVer * rVer;
    int numPointsPerDecade = 20; 
    //cuda parameters
    int threads=128;
    int blocks=N/threads + (N%threads == 0?0:1);
    if(blocks >= 65536)
    {
        std::cout<<"Maximum block number (65536) exeeded: "<<blocks<<"\n";
        return 1;
    }
    if( blocks * threads != N)
    {
        NthreadEqN = false;
    }
    dim3 dimGrid(blocks);
    dim3 dimBlock(threads);

    //reduction tuning parameters
    int redM = 128;
    int redN = 128;
    //derived parameters
    double logInc = (log(10)-log(1))/static_cast<double>(numPointsPerDecade);//logarithmic increment
    int numLogOutputPoints = static_cast<int>( (log(NSim) - log(NOffset) )/logInc);
    int LogOutputCounter = 1;
    int LogOutputOffset = NOffset==0 ? 1 : NOffset;
    int nextLogOutput = static_cast<int>(exp(logInc * LogOutputCounter)*LogOutputOffset);
    std::string backup="_bkp_"+basename;

    //basic array setup
    FLOAT_ARRAY_TYPE *hPos; //particle positions on host
    FLOAT_ARRAY_TYPE *hVel; //particle velocities on host
    FLOAT_ARRAY_TYPE *dPos; //particle positions on device
    FLOAT_ARRAY_TYPE *dVel; //particle velocities on device
    FLOAT_ARRAY_TYPE *velOnDt; //particle velocities on device at dt
    FLOAT_ARRAY_TYPE *dForce; //forces action on particles on device
    float *dEPot; //potentials at each particles on device
    float *reduction; //array used in reduction kernel
    FLOAT_ARRAY_TYPE *compReduction; //array used in component wise reduction
    int *hCnn; //coupled nearest neighbor list on host
    int *dCnn; //coupled nearest neighbor list on device

    //generate verlet list data structures in gpu memory
    int *vCount; //holds number of neigbor list entries for each particle
    int *vLists; //holds in each row the verlet list for one particle
    size_t vListsPitch; //receives pitch value from cudaMallocPitch call
    int NVList = static_cast<int>(0.75 * pow( 2.0 * rVer, 3.));//maximum number of entries in each verlet list
    //std::cout<<"NVList = "<<NVList<<"\n";
    float *dMax; //holds displacement of each particle since last verlet list update
    bool *dListUpdate;//flag for list update on device
    bool hListUpdate=false;//flag for list update on host
    float maxDisplacement=(rVer - rCut) * 0.5; //maximum allowed displacement before next verlet list update
    //std::cout<<"max allowed displacement: "<<maxDisplacement<<"\n";
    //std::cout<<"list update: "<<hListUpdate<<"\n";


    //allocate memory
    //host
    hPos = new FLOAT_ARRAY_TYPE[N];
    hVel = new FLOAT_ARRAY_TYPE[N];
    hCnn = new int[numNN*N];

    //read in data
    success = loadInitData(inputfile, hPos, hVel, hCnn, N, numNN);
    if(!success)
    {
        std::cout<<"Data could not be read!\n";
        delete[] hPos;
        delete[] hVel;
        delete[] hCnn;
        return 1;
    }

    //save inital setup
    success = writeData(basename, hPos, hVel, hCnn, N, numNN, NOffset);

    //allocate device memory
    size_t size=N*sizeof(FLOAT_ARRAY_TYPE);

    cudaMalloc((void **) &dPos, size);
    checkCudaError("Mem Alloc dPos");
    success=memInfo();
    cudaMalloc((void **) &dVel, size);
    checkCudaError("Mem Alloc dVel");
    cudaMalloc((void **) &velOnDt, size);
    checkCudaError("Mem Alloc velOnDt");
    cudaMalloc((void **) &dForce, size);
    checkCudaError("Mem Alloc dForce");
    cudaMalloc((void **) &dEPot, N* sizeof(float) );
    checkCudaError("Mem Alloc dEPot");
    cudaMalloc((void**) &reduction, redM * sizeof(float));
    checkCudaError("Mem Alloc reduction");
    cudaMalloc((void**) &compReduction, redM * sizeof(float3));
    checkCudaError("Mem Alloc compReduction");

    size_t hCnn_pitch=N * sizeof(int); //pitch of host data
    size_t dCnnPitch=0; //receives pitch value from cudaMallocPitch call
    cudaMallocPitch((void**)&dCnn, &dCnnPitch, N * sizeof(int), numNN);
    checkCudaError("Pitch Mem Alloc dCnn");

    //verlet list allocations
    cudaMalloc((void **) &vCount, N * sizeof(int));
    checkCudaError("Mem Alloc vCount");
    cudaMallocPitch((void**)&vLists, &vListsPitch, N * sizeof(int), NVList);
    checkCudaError("Pitch Mem Alloc vLists");
    cudaMalloc((void **) &dListUpdate, sizeof(bool));
    checkCudaError("Mem Alloc dListUpdate");
    cudaMalloc((void **) &dMax, N * sizeof(float));
    checkCudaError("Mem Alloc dMax");
    success=memInfo();

    //copy initial data to device
    cudaMemcpy(dPos, hPos, sizeof(FLOAT_ARRAY_TYPE)*N, cudaMemcpyHostToDevice);
    checkCudaError("Mem Cpy: dPos, hPos");
    cudaMemcpy(dVel, hVel, sizeof(FLOAT_ARRAY_TYPE)*N, cudaMemcpyHostToDevice);
    checkCudaError("Mem Cpy: dVel, hVel");
    cudaMemcpy2D(dCnn, dCnnPitch, hCnn, hCnn_pitch, N * sizeof(int), numNN, cudaMemcpyHostToDevice);
    checkCudaError("Mem Cpy 2D: dCnn, hCnn");
    cudaMemcpy(dListUpdate, &hListUpdate, sizeof(bool), cudaMemcpyHostToDevice);


    //execute one kernel before timing begins
    gpuResetForces<<<dimGrid, dimBlock>>>(dForce, N);
    checkCudaError("Kernel error: gpuResetForces");

    //#############################################################
    //#### start timimg
    //cuda timing
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    float time;
    //host timing
    clock_t cStart, cEnd;
    cudaEventRecord(start, 0);
    cStart=clock();

    //generate first time verlet list
    if(NthreadEqN)
    {
        gpuGenerateVerletListSmem<<<dimGrid, dimBlock, dimBlock.x*sizeof(FLOAT_ARRAY_TYPE)>>>(dPos, vLists, vListsPitch, vCount, dMax, rVerSq, N, boxSize, dListUpdate);
        checkCudaError("Kernel error: gpuGenerateVerletListSmem");
    }
    else
    {
        gpuGenerateVerletListVar1<<<dimGrid, dimBlock>>>(dPos, vLists, vListsPitch, vCount, dMax, rVerSq, N, boxSize, dListUpdate);
        checkCudaError("Kernel error: gpuGenerateVerletList");
    }


    //rescale temperature
    rescaleT(N, dim, redM, redN, T, threads, blocks, dVel, reduction);

    //remove drift
    removeDrift(N, redM, redN, threads, blocks, dVel, compReduction);

    //do initial velocity backstep
    gpuResetForces<<<dimGrid, dimBlock>>>(dForce, N);
    checkCudaError("Kernel error: gpuResetForces");

    //Verlet List kernel
    gpuLJForcesVerletList<<< dimGrid, dimBlock>>>(dPos, dForce, N, vLists, vListsPitch, vCount, boxSize, rCutSq);
    checkCudaError("Kernel error: gpuLJForcesVerletList");

    gpuHarmonicBondForces<<< dimGrid, dimBlock>>>(dPos, dForce, dCnn, dCnnPitch, N, numNN, k, boxSize);
    checkCudaError("Kernel error: gpuHarmonicBondForces");
    gpuInitialBackstep<<<dimGrid, dimBlock>>>(dVel, dForce, N, dt);
    checkCudaError("Kernel error: gpuInitialBackstep");

    std::cout.setf(std::ios_base::scientific);
    std::cout.precision(std::numeric_limits< float >::digits10);

    //main simulation loop
    int counter = NOffset + 1;
    for(int i=NOffset+NOut; i<=NSim; i=i+NOut)
    {
        for(int j=NAdjust; j<= NOut; j=j+NAdjust)
        {
            //adjust velocities and remove drift
            rescaleT(N, dim, redM, redN, T, threads, blocks, dVel, reduction);
            removeDrift(N, redM, redN, threads, blocks, dVel, compReduction);

            for(int mAdj=1; mAdj<=NAdjust; ++mAdj)
            {
                //check if neigbor list update is needed
                cudaMemcpy(&hListUpdate, dListUpdate, sizeof(bool), cudaMemcpyDeviceToHost);
                checkCudaError("Mem Cpy: &hListUpdate, dListUpdate");
                if(hListUpdate == true)
                {
                    //std::cout<<"VL update @: "<<counter<<"\n";
                    if(NthreadEqN)
                    {
                        gpuGenerateVerletListSmem<<<dimGrid, dimBlock, dimBlock.x*sizeof(FLOAT_ARRAY_TYPE)>>>(dPos, vLists, vListsPitch, vCount, dMax, rVerSq, N, boxSize, dListUpdate);
                        checkCudaError("Kernel error: gpuGenerateVerletListSmem");
                    }
                    else
                    {
                        gpuGenerateVerletListVar1<<<dimGrid, dimBlock>>>(dPos, vLists, vListsPitch, vCount, dMax, rVerSq, N, boxSize, dListUpdate);
                        checkCudaError("Kernel error: gpuGenerateVerletList");
                    }
                    hListUpdate = false;
                }

                //calculate forces
                gpuResetForces<<<dimGrid, dimBlock>>>(dForce, N);
                checkCudaError("Kernel error: gpuResetForces");

                //Verlet List kernel
                gpuLJForcesVerletList<<< dimGrid, dimBlock>>>(dPos, dForce, N, vLists, vListsPitch, vCount, boxSize, rCutSq);
                checkCudaError("Kernel error: gpuLJForcesVerletList");
                gpuHarmonicBondForces<<< dimGrid, dimBlock>>>(dPos, dForce, dCnn, dCnnPitch, N, numNN, k, boxSize);
                checkCudaError("Kernel error: gpuHarmonicBondForces");
                //integrate
                gpuIntegratorLeapFrog<<< dimGrid, dimBlock>>>(dPos, dVel, dForce, N, dt, boxSize, dMax, maxDisplacement, dListUpdate);
                checkCudaError("Kernel error: gpuIntegratorLeapFrog");

                if(nextLogOutput == counter)
                {
                    if(i == nextLogOutput)
                    {
                        linLogEqual = true;
                    }
                    gpuResetForces<<<dimGrid, dimBlock>>>(dForce, N);
                    checkCudaError("Kernel error: gpuResetForces");
                    //Verlet List kernel
                    gpuLJForcesVerletList<<< dimGrid, dimBlock>>>(dPos, dForce, N, vLists, vListsPitch, vCount, boxSize, rCutSq);;
                    checkCudaError("Kernel error: gpuLJForcesVerletList");
                    gpuHarmonicBondForces<<< dimGrid, dimBlock>>>(dPos, dForce, dCnn, dCnnPitch, N, numNN, k, boxSize);
                    checkCudaError("Kernel error: gpuHarmonicBondForces");
                    gpuGetVelOnTime<<<dimGrid, dimBlock>>>(dVel, velOnDt, dForce, N, dt);
                    checkCudaError("Kernel error: gpuGetVelOnTime");
                    //transfer x-v-data to host
                    cudaMemcpy(hPos, dPos, sizeof(FLOAT_ARRAY_TYPE)*N, cudaMemcpyDeviceToHost);
                    checkCudaError("Mem Cpy: hPos, dPos");
                    cudaMemcpy(hVel, velOnDt, sizeof(FLOAT_ARRAY_TYPE)*N, cudaMemcpyDeviceToHost);
                    checkCudaError("Mem Cpy: hVel, velOnDt");
                    //write output
                    success = writeData(basename, hPos, hVel, hCnn, N, numNN, counter);
                    int tmp;
                    do
                    {
                        LogOutputCounter += 1;
                        tmp = static_cast<int>(exp(logInc * LogOutputCounter)*LogOutputOffset);
                    }
                    while (tmp == nextLogOutput);
                    nextLogOutput = tmp;
                }
                counter++;
            }//end adjustment loop
        }//end output loop
        //write output
        //get velocities on time slice (due to leap frog shift)
        if(linLogEqual == false)
        {
            gpuResetForces<<<dimGrid, dimBlock>>>(dForce, N);
            checkCudaError("Kernel error: gpuResetForces");

            //Verlet List kernel
            gpuLJForcesVerletList<<< dimGrid, dimBlock>>>(dPos, dForce, N, vLists, vListsPitch, vCount, boxSize, rCutSq);
            checkCudaError("Kernel error: gpuLJForcesVerletList");
            gpuHarmonicBondForces<<< dimGrid, dimBlock>>>(dPos, dForce, dCnn, dCnnPitch, N, numNN, k, boxSize);
            checkCudaError("Kernel error: gpuHarmonicBondForces");
            gpuGetVelOnTime<<<dimGrid, dimBlock>>>(dVel, velOnDt, dForce, N, dt);
            checkCudaError("Kernel error: gpuGetVelOnTime");

            //transfer x-v-data to host
            cudaMemcpy(hPos, dPos, sizeof(FLOAT_ARRAY_TYPE)*N, cudaMemcpyDeviceToHost);
            checkCudaError("Mem Cpy: hPos, dPos");
            cudaMemcpy(hVel, velOnDt, sizeof(FLOAT_ARRAY_TYPE)*N, cudaMemcpyDeviceToHost);
            checkCudaError("Mem Cpy: hVel, velOnDt");
        }
        //write output
        success = writeData(backup, hPos, hVel, hCnn, N, numNN, i);
        linLogEqual = false;
    }//end simulation loop

    //#### stop timing
    cudaThreadSynchronize();//wait until all threads are finished
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);//blocking till stop has been recorded
    cEnd=clock();
    std::cout<<"########### CPU TIME ###########\n";
    std::cout<<"start: "<<cStart<<"\nend: "<<cEnd<<"\n";
    std::cout<<"clocks per second: "<<CLOCKS_PER_SEC<<"\n";
    std::cout<<"cpu time: "<<(cEnd-cStart)/static_cast<double>(CLOCKS_PER_SEC)<<"s\n\n";
    cudaEventElapsedTime( &time, start, stop);
    std::cout<<"########## CUDA TIMER #########\n";
    std::cout<<"elapsed time: "<<time<<"ms\n";
    //#############################################################
    //free memory
    delete[] hPos;
    delete[] hVel;
    delete[] hCnn;
    cudaFree(dPos);
    cudaFree(dVel);
    cudaFree(velOnDt);
    cudaFree(dForce);
    cudaFree(dEPot);
    cudaFree(reduction);
    cudaFree(compReduction);
    cudaFree(dCnn);
    cudaFree(vCount);
    cudaFree(vLists);
    cudaFree(dMax);
    cudaFree(dListUpdate);
    success=memInfo();

    return 0;
}

---