Program Listing for File NuMagSANSlib_MagData.h

↰ Return to documentation for file (src/InputData/MagData/NuMagSANSlib_MagData.h)

// File         : NuMagSANSlib_MagData.h
// Author       : Michael Philipp ADAMS, M.Sc.
// Company      : University of Luxembourg
// Department   : Department of Physics and Materials Sciences
// Group        : NanoMagnetism Group
// Group Leader : Prof. Andreas Michels
// Version      : 26 November 2024
// OS           : Linux Ubuntu
// Language     : CUDA C++

#include <iostream>
#include <fstream>
#include <sstream>
#include <sys/stat.h>
#include <sys/types.h>
#include <math.h>
#include <string>
#include <vector>
#include <stdlib.h>
#include <time.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdexcept>
#include <math.h>
#include <chrono>
#include <dirent.h>
#include <unistd.h>

using namespace std;

struct MagnetizationData {

    float* RotMat; // rotation matrix R = [R0, R3, R6; R1, R4, R7; R2, R5, R8]
    float* x;      // position data x
    float* y;      // position data y
    float* z;      // position data z
    float* mx;     // magnetization data mx
    float* my;     // magnetization data my
    float* mz;     // magnetization data mz

    unsigned long int* K;               // number of objects (e.g., orientations)
    unsigned long int* TotalAtomNumber; // total atom number

    unsigned int* NumberOfElements;
    unsigned int* NumberOfNonZeroMoments;
    unsigned long int* N_act; // number of atoms in k-th object
    unsigned long int* N_cum; // cummulated number of atoms from object to object
    unsigned long int* N_avg; // average number of atoms per object
};

void allocate_MagnetizationDataRAM(MagnetizationData* MagData, MagDataProperties* MagDataProp, InputFileData* InputData,
                                   int MagData_File_Index) {

    bool ExcludeZeroMoments_flag = InputData->ExcludeZeroMoments_flag;
    unsigned long int K = MagDataProp->Number_Of_SubFolders;
    unsigned long int TotalAtomNumber = 0;
    if (ExcludeZeroMoments_flag) {
        TotalAtomNumber = MagDataProp->TotalNZMAtomNumber[MagData_File_Index - 1];
    } else {
        TotalAtomNumber = MagDataProp->TotalAtomNumber[MagData_File_Index - 1];
    }

    MagData->x = (float*)malloc(TotalAtomNumber * sizeof(float));
    MagData->y = (float*)malloc(TotalAtomNumber * sizeof(float));
    MagData->z = (float*)malloc(TotalAtomNumber * sizeof(float));
    MagData->mx = (float*)malloc(TotalAtomNumber * sizeof(float));
    MagData->my = (float*)malloc(TotalAtomNumber * sizeof(float));
    MagData->mz = (float*)malloc(TotalAtomNumber * sizeof(float));
    MagData->K = (unsigned long int*)malloc(sizeof(unsigned long int));
    MagData->TotalAtomNumber = (unsigned long int*)malloc(sizeof(unsigned long int));
    MagData->RotMat = (float*)malloc(9 * sizeof(float));
    MagData->NumberOfElements = (unsigned int*)malloc(K * sizeof(unsigned int));
    MagData->NumberOfNonZeroMoments = (unsigned int*)malloc(K * sizeof(unsigned int));
    MagData->N_cum = (unsigned long int*)malloc(K * sizeof(unsigned long int));
    MagData->N_act = (unsigned long int*)malloc(K * sizeof(unsigned long int));
    MagData->N_avg = (unsigned long int*)malloc(sizeof(unsigned long int));

    if (!MagData->x || !MagData->y || !MagData->z || !MagData->mx || !MagData->my || !MagData->mz || !MagData->N_avg ||
        !MagData->K || !MagData->TotalAtomNumber || !MagData->RotMat) {
        perror("Memory allocation failed");
        exit(EXIT_FAILURE);
    }

    *MagData->N_avg = 0;
    *MagData->K = K;
    *MagData->TotalAtomNumber = TotalAtomNumber;

    for (unsigned long int i = 0; i < TotalAtomNumber; i++) {
        MagData->x[i] = 0.0;
        MagData->y[i] = 0.0;
        MagData->z[i] = 0.0;
        MagData->mx[i] = 0.0;
        MagData->my[i] = 0.0;
        MagData->mz[i] = 0.0;
    }

    for (unsigned int i = 0; i < K; i++) {
        MagData->NumberOfElements[i] = MagDataProp->NumberOfElements[i][MagData_File_Index - 1];
        MagData->NumberOfNonZeroMoments[i] = MagDataProp->NumberOfNonZeroMoments[i][MagData_File_Index - 1];
        MagData->N_cum[i] = 0;
        MagData->N_act[i] = 0;
    }
}

void allocate_MagnetizationDataGPU(MagnetizationData* MagData, MagnetizationData* MagData_gpu, int MagData_File_Index) {

    unsigned long int TotalAtomNumber = *MagData->TotalAtomNumber;

    cudaMalloc(&MagData_gpu->x, TotalAtomNumber * sizeof(float));
    cudaMalloc(&MagData_gpu->y, TotalAtomNumber * sizeof(float));
    cudaMalloc(&MagData_gpu->z, TotalAtomNumber * sizeof(float));
    cudaMalloc(&MagData_gpu->mx, TotalAtomNumber * sizeof(float));
    cudaMalloc(&MagData_gpu->my, TotalAtomNumber * sizeof(float));
    cudaMalloc(&MagData_gpu->mz, TotalAtomNumber * sizeof(float));
    cudaMalloc(&MagData_gpu->K, sizeof(unsigned long int));
    cudaMalloc(&MagData_gpu->TotalAtomNumber, sizeof(unsigned long int));
    cudaMalloc(&MagData_gpu->RotMat, 9 * sizeof(float));
    cudaMalloc(&MagData_gpu->NumberOfElements, *MagData->K * sizeof(unsigned int));
    cudaMalloc(&MagData_gpu->NumberOfNonZeroMoments, *MagData->K * sizeof(unsigned int));
    cudaMalloc(&MagData_gpu->N_cum, *MagData->K * sizeof(unsigned long int));
    cudaMalloc(&MagData_gpu->N_act, *MagData->K * sizeof(unsigned long int));
    cudaMalloc(&MagData_gpu->N_avg, sizeof(unsigned long int));

    cudaMemcpy(MagData_gpu->N_avg, MagData->N_avg, sizeof(unsigned long int), cudaMemcpyHostToDevice);
    cudaMemcpy(MagData_gpu->K, MagData->K, sizeof(unsigned long int), cudaMemcpyHostToDevice);
    cudaMemcpy(MagData_gpu->TotalAtomNumber, MagData->TotalAtomNumber, sizeof(unsigned long int),
               cudaMemcpyHostToDevice);
    cudaMemcpy(MagData_gpu->RotMat, MagData->RotMat, 9 * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(MagData_gpu->NumberOfElements, MagData->NumberOfElements, *MagData->K * sizeof(unsigned int),
               cudaMemcpyHostToDevice);
    cudaMemcpy(MagData_gpu->NumberOfElements, MagData->NumberOfElements, *MagData->K * sizeof(unsigned int),
               cudaMemcpyHostToDevice);
    cudaMemcpy(MagData_gpu->N_cum, MagData->N_cum, *MagData->K * sizeof(unsigned long int), cudaMemcpyHostToDevice);
    cudaMemcpy(MagData_gpu->N_act, MagData->N_act, *MagData->K * sizeof(unsigned long int), cudaMemcpyHostToDevice);

    LogSystem::write("");
    // copy data from Host to Device
    LogSystem::write("copy data from RAM to GPU...");
    cudaMemcpy(MagData_gpu->x, MagData->x, TotalAtomNumber * sizeof(float), cudaMemcpyHostToDevice);
    LogSystem::write("  x done...");
    cudaMemcpy(MagData_gpu->y, MagData->y, TotalAtomNumber * sizeof(float), cudaMemcpyHostToDevice);
    LogSystem::write("  y done...");
    cudaMemcpy(MagData_gpu->z, MagData->z, TotalAtomNumber * sizeof(float), cudaMemcpyHostToDevice);
    LogSystem::write("  z done...");
    cudaMemcpy(MagData_gpu->mx, MagData->mx, TotalAtomNumber * sizeof(float), cudaMemcpyHostToDevice);
    LogSystem::write("  mx done...");
    cudaMemcpy(MagData_gpu->my, MagData->my, TotalAtomNumber * sizeof(float), cudaMemcpyHostToDevice);
    LogSystem::write("  my done...");
    cudaMemcpy(MagData_gpu->mz, MagData->mz, TotalAtomNumber * sizeof(float), cudaMemcpyHostToDevice);
    LogSystem::write("  mz done...");
    LogSystem::write("");
    LogSystem::write("data transfer finished...");
    LogSystem::write("");

    cudaMalloc(&MagData_gpu, sizeof(MagnetizationData));
    cudaMemcpy(MagData_gpu, MagData, sizeof(MagnetizationData), cudaMemcpyHostToDevice);
}

void read_MagnetizationData(MagnetizationData* MagData, MagDataProperties* MagDataProp, InputFileData* InputData,
                            int MagData_File_Index) {

    LogSystem::write("");
    LogSystem::write("read magnetization data to RAM...");

    unsigned long int K = *MagData->K;
    bool ExcludeZeroMoments_flag = InputData->ExcludeZeroMoments_flag;

    for (int i = 0; i < 9; i++) {
        MagData->RotMat[i] = InputData->RotMat[i];
    }

    string filename;
    unsigned long int n = 0;
    float x_buf, y_buf, z_buf, mx_buf, my_buf, mz_buf;
    ifstream fin;
    float x_mean = 0.0;
    float y_mean = 0.0;
    float z_mean = 0.0;

    unsigned long int N_cum = 0; // cummulative counter
    unsigned long int N_act = 0; // actual number of atoms per file

    for (unsigned long int k = 1; k <= K; k++) {

        if (ExcludeZeroMoments_flag) {
            N_act = MagData->NumberOfNonZeroMoments[k - 1];
        } else {
            N_act = MagData->NumberOfElements[k - 1];
        }
        MagData->N_act[k - 1] = N_act;

        filename = MagDataProp->GlobalFolderPath + "/" + MagDataProp->SubFolderNames_Nom + "_" + to_string(k) + "/" +
                   MagDataProp->SubFolder_FileNames_Nom + "_" + to_string(MagData_File_Index) + "." +
                   MagDataProp->SubFolder_FileNames_Type;

        LogSystem::write(filename);

        fin.open(filename);
        n = 0;
        x_mean = 0.0;
        y_mean = 0.0;
        z_mean = 0.0;

        // read in the data and calculate the spatial average position
        if (ExcludeZeroMoments_flag) {
            while (fin >> x_buf >> y_buf >> z_buf >> mx_buf >> my_buf >> mz_buf) {
                if (mx_buf != 0.0 || my_buf != 0.0 || mz_buf != 0.0) {
                    MagData->x[n + N_cum] = x_buf * InputData->XYZ_Unit_Factor;
                    MagData->y[n + N_cum] = y_buf * InputData->XYZ_Unit_Factor;
                    MagData->z[n + N_cum] = z_buf * InputData->XYZ_Unit_Factor;
                    MagData->mx[n + N_cum] = mx_buf;
                    MagData->my[n + N_cum] = my_buf;
                    MagData->mz[n + N_cum] = mz_buf;

                    x_mean += MagData->x[n + N_cum] / N_act;
                    y_mean += MagData->y[n + N_cum] / N_act;
                    z_mean += MagData->z[n + N_cum] / N_act;

                    n += 1;
                }
            }
        } else {
            while (fin >> x_buf >> y_buf >> z_buf >> mx_buf >> my_buf >> mz_buf) {

                MagData->x[n + N_cum] = x_buf * InputData->XYZ_Unit_Factor;
                MagData->y[n + N_cum] = y_buf * InputData->XYZ_Unit_Factor;
                MagData->z[n + N_cum] = z_buf * InputData->XYZ_Unit_Factor;
                MagData->mx[n + N_cum] = mx_buf;
                MagData->my[n + N_cum] = my_buf;
                MagData->mz[n + N_cum] = mz_buf;

                x_mean += MagData->x[n + N_cum] / N_act;
                y_mean += MagData->y[n + N_cum] / N_act;
                z_mean += MagData->z[n + N_cum] / N_act;

                n += 1;
            }
        }

        fin.close();

        // centering of the xyz-data set
        for (int l = 0; l < N_act; l++) {
            MagData->x[l + N_cum] = MagData->x[l + N_cum] - x_mean;
            MagData->y[l + N_cum] = MagData->y[l + N_cum] - y_mean;
            MagData->z[l + N_cum] = MagData->z[l + N_cum] - z_mean;
        }

        // update of the cummulative counter
        N_cum += N_act;
        if (k < K) {
            MagData->N_cum[k] = N_cum;
        }
        LogSystem::write("N_act: " + std::to_string(N_act) + ", " + "N_cum: " + std::to_string(N_cum));
    }

    *MagData->N_avg = (unsigned long int)(((float)N_cum) / ((float)K));
    LogSystem::write("N_avg: " + to_string(*MagData->N_avg));
    LogSystem::write("read (x, y, z, mx, my, mz) data to RAM finished...");
}

void replication_read_MagnetizationData(MagnetizationData* MagData, MagDataProperties* MagDataProp,
                                        InputFileData* InputData, int MagData_File_Index) {

    LogSystem::write("");
    LogSystem::write("read replicated magnetization data to RAM...");

    unsigned long int K = *MagData->K;
    bool ExcludeZeroMoments_flag = InputData->ExcludeZeroMoments_flag;

    for (int i = 0; i < 9; i++) {
        MagData->RotMat[i] = InputData->RotMat[i];
    }

    string filename;
    unsigned long int n = 0;
    float x_buf, y_buf, z_buf, mx_buf, my_buf, mz_buf;
    ifstream fin;
    float x_mean = 0.0;
    float y_mean = 0.0;
    float z_mean = 0.0;

    unsigned long int N_cum = 0; // cummulative counter
    unsigned long int N_act = 0; // actual number of atoms per replicated object

    if (ExcludeZeroMoments_flag) {
        N_act = MagData->NumberOfNonZeroMoments[0];
    } else {
        N_act = MagData->NumberOfElements[0];
    }

    filename = MagDataProp->GlobalFolderPath + "/" + MagDataProp->SubFolderNames_Nom + "_1" + "/" +
               MagDataProp->SubFolder_FileNames_Nom + "_" + to_string(MagData_File_Index) + "." +
               MagDataProp->SubFolder_FileNames_Type;

    LogSystem::write(filename);

    std::vector<float> x_template(N_act);
    std::vector<float> y_template(N_act);
    std::vector<float> z_template(N_act);
    std::vector<float> mx_template(N_act);
    std::vector<float> my_template(N_act);
    std::vector<float> mz_template(N_act);

    fin.open(filename);
    n = 0;
    x_mean = 0.0;
    y_mean = 0.0;
    z_mean = 0.0;

    if (ExcludeZeroMoments_flag) {
        while (fin >> x_buf >> y_buf >> z_buf >> mx_buf >> my_buf >> mz_buf) {
            if (mx_buf != 0.0 || my_buf != 0.0 || mz_buf != 0.0) {
                x_template[n] = x_buf * InputData->XYZ_Unit_Factor;
                y_template[n] = y_buf * InputData->XYZ_Unit_Factor;
                z_template[n] = z_buf * InputData->XYZ_Unit_Factor;
                mx_template[n] = mx_buf;
                my_template[n] = my_buf;
                mz_template[n] = mz_buf;

                x_mean += x_template[n] / N_act;
                y_mean += y_template[n] / N_act;
                z_mean += z_template[n] / N_act;

                n += 1;
            }
        }
    } else {
        while (fin >> x_buf >> y_buf >> z_buf >> mx_buf >> my_buf >> mz_buf) {
            x_template[n] = x_buf * InputData->XYZ_Unit_Factor;
            y_template[n] = y_buf * InputData->XYZ_Unit_Factor;
            z_template[n] = z_buf * InputData->XYZ_Unit_Factor;
            mx_template[n] = mx_buf;
            my_template[n] = my_buf;
            mz_template[n] = mz_buf;

            x_mean += x_template[n] / N_act;
            y_mean += y_template[n] / N_act;
            z_mean += z_template[n] / N_act;

            n += 1;
        }
    }

    fin.close();

    for (unsigned long int l = 0; l < N_act; l++) {
        x_template[l] = x_template[l] - x_mean;
        y_template[l] = y_template[l] - y_mean;
        z_template[l] = z_template[l] - z_mean;
    }

    for (unsigned long int k = 1; k <= K; k++) {

        MagData->N_act[k - 1] = N_act;

        for (unsigned long int l = 0; l < N_act; l++) {
            MagData->x[l + N_cum] = x_template[l];
            MagData->y[l + N_cum] = y_template[l];
            MagData->z[l + N_cum] = z_template[l];
            MagData->mx[l + N_cum] = mx_template[l];
            MagData->my[l + N_cum] = my_template[l];
            MagData->mz[l + N_cum] = mz_template[l];
        }

        // update of the cummulative counter
        N_cum += N_act;
        if (k < K) {
            MagData->N_cum[k] = N_cum;
        }
        LogSystem::write("replicated object " + std::to_string(k) + ": N_act: " + std::to_string(N_act) + ", " +
                         "N_cum: " + std::to_string(N_cum));
    }

    *MagData->N_avg = (unsigned long int)(((float)N_cum) / ((float)K));
    LogSystem::write("N_avg: " + to_string(*MagData->N_avg));
    LogSystem::write("read replicated (x, y, z, mx, my, mz) data to RAM finished...");
}

void init_MagnetizationData(MagnetizationData* MagData, MagnetizationData* MagData_gpu, MagDataProperties* MagDataProp,
                            InputFileData* InputData, int MagData_File_Index) {

    allocate_MagnetizationDataRAM(MagData, MagDataProp, InputData, MagData_File_Index);
    if (InputData->MagData_ReplicationImport_flag) {
        replication_read_MagnetizationData(MagData, MagDataProp, InputData, MagData_File_Index);
    } else {
        read_MagnetizationData(MagData, MagDataProp, InputData, MagData_File_Index);
    }
    allocate_MagnetizationDataGPU(MagData, MagData_gpu, MagData_File_Index);
}

void free_MagnetizationData(MagnetizationData* MagData, MagnetizationData* MagData_gpu) {

    LogSystem::write("free MagData...");

    //  cudaError_t err;

    free(MagData->x);
    free(MagData->y);
    free(MagData->z);
    free(MagData->mx);
    free(MagData->my);
    free(MagData->mz);
    free(MagData->K);
    free(MagData->TotalAtomNumber);
    free(MagData->RotMat);
    free(MagData->NumberOfElements);
    free(MagData->NumberOfNonZeroMoments);
    free(MagData->N_cum);
    free(MagData->N_act);
    free(MagData->N_avg);

    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->x);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->y);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->z);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->mx);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->my);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->mz);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->K);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->TotalAtomNumber);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->RotMat);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->NumberOfElements);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->NumberOfNonZeroMoments);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->N_cum);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->N_act);
    cudaDeviceSynchronize();
    cudaFree(MagData_gpu->N_avg);
    cudaDeviceSynchronize();

    LogSystem::write("free MagData finished...");
}

void disp_MagnetizationData(MagnetizationData* MagData) {
    for (int k = 0; k < *MagData->TotalAtomNumber; k++) {
        cout << MagData->x[k] << " " << MagData->y[k] << " " << MagData->z[k] << " " << MagData->mx[k] << " "
             << MagData->my[k] << " " << MagData->mz[k] << "\n";
    }
}