diag.c

Go to the documentation of this file.

 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "ascot5.h"
#include "simulate.h"
#include "B_field.h"
#include "diag.h"
#include "diag/diag_orb.h"
#include "diag/dist_5D.h"
#include "diag/dist_6D.h"
#include "diag/dist_rho5D.h"
#include "diag/dist_rho6D.h"
#include "diag/dist_com.h"
#include "diag/diag_transcoef.h"
#include "particle.h"
 
void diag_arraysum(size_t start, size_t stop, real* array1, real* array2);
 
int diag_init_offload(diag_offload_data* data, real** offload_array, int Nmrk){
    /* Determine how long array we need and allocate it */
    size_t n = 0;
 
    if(data->dist5D_collect) {
        data->offload_dist5D_index = n;
        n += (size_t)(data->dist5D.n_r)     * (size_t)(data->dist5D.n_phi)
           * (size_t)(data->dist5D.n_z)     * (size_t)(data->dist5D.n_ppara)
           * (size_t)(data->dist5D.n_pperp) * (size_t)(data->dist5D.n_time)
           * (size_t)(data->dist5D.n_q);
    }
 
    if(data->dist6D_collect) {
        data->offload_dist6D_index = n;
        n += (size_t)(data->dist6D.n_r)    * (size_t)(data->dist6D.n_phi)
           * (size_t)(data->dist6D.n_z)    * (size_t)(data->dist6D.n_pr)
           * (size_t)(data->dist6D.n_pphi) * (size_t)(data->dist6D.n_pz)
           * (size_t)(data->dist6D.n_time) * (size_t)(data->dist6D.n_q);
    }
 
    if(data->distrho5D_collect) {
        data->offload_distrho5D_index = n;
        n += (size_t)(data->distrho5D.n_rho) * (size_t)(data->distrho5D.n_theta)
           * (size_t)(data->distrho5D.n_phi) * (size_t)(data->distrho5D.n_ppara)
           * (size_t)(data->distrho5D.n_pperp)
           * (size_t)(data->distrho5D.n_time) * (size_t)(data->distrho5D.n_q);
    }
 
    if(data->distrho6D_collect) {
        data->offload_distrho6D_index = n;
        n += (size_t)(data->distrho6D.n_rho) * (size_t)(data->distrho6D.n_theta)
           * (size_t)(data->distrho6D.n_phi) * (size_t)(data->distrho6D.n_pr)
           * (size_t)(data->distrho6D.n_pphi) * (size_t)(data->distrho6D.n_pz)
           * (size_t)(data->distrho6D.n_time) * (size_t)(data->distrho6D.n_q);
    }
 
    if(data->distCOM_collect) {
        data->offload_distCOM_index = n;
        n += (size_t)(data->distCOM.n_mu) * (size_t)(data->distCOM.n_Ekin)
            * (size_t)(data->distCOM.n_Ptor);
    }
 
    data->offload_dist_length = n;
 
    if(data->diagorb_collect) {
        data->offload_diagorb_index = n;
        data->diagorb.Nmrk = Nmrk;
 
        switch(data->diagorb.record_mode) {
 
            case simulate_mode_fo:
                data->diagorb.Nfld = DIAG_ORB_FOFIELDS;
                break;
 
            case simulate_mode_gc:
                data->diagorb.Nfld = DIAG_ORB_GCFIELDS;
                break;
 
            case simulate_mode_ml:
                data->diagorb.Nfld = DIAG_ORB_MLFIELDS;
                break;
 
            case simulate_mode_hybrid:
                data->diagorb.Nfld = DIAG_ORB_HYBRIDFIELDS;
                break;
 
        }
 
        if(data->diagorb.mode == DIAG_ORB_POINCARE) {
            n += (size_t)(data->diagorb.Nfld+2)
                * (size_t)(data->diagorb.Nmrk) * (size_t)(data->diagorb.Npnt);
        }
        else if(data->diagorb.mode == DIAG_ORB_INTERVAL) {
            n += (size_t)(data->diagorb.Nfld)
                * (size_t)(data->diagorb.Nmrk) * (size_t)(data->diagorb.Npnt);
        }
    }
    if(data->diagtrcof_collect) {
        data->offload_diagtrcof_index = n;
        data->diagtrcof.Nmrk = Nmrk;
        n += (size_t)(3*data->diagtrcof.Nmrk);
    }
 
    data->offload_array_length = n;
    *offload_array = malloc(n * sizeof(real));
    if(*offload_array == NULL) {
        return 1;
    }
 
    memset(*offload_array, 0, n * sizeof(real));
 
    return 0;
}
 
void diag_free_offload(diag_offload_data* data, real** offload_array) {
    free(*offload_array);
    *offload_array = NULL;
}
 
void diag_init(diag_data* data, diag_offload_data* offload_data,
               real* offload_array) {
    data->diagorb_collect   = offload_data->diagorb_collect;
    data->dist5D_collect    = offload_data->dist5D_collect;
    data->dist6D_collect    = offload_data->dist6D_collect;
    data->distrho5D_collect = offload_data->distrho5D_collect;
    data->distrho6D_collect = offload_data->distrho6D_collect;
    data->distCOM_collect   = offload_data->distCOM_collect;
    data->diagtrcof_collect = offload_data->diagtrcof_collect;
 
    if(data->dist5D_collect) {
        dist_5D_init(&data->dist5D, &offload_data->dist5D,
                     &offload_array[offload_data->offload_dist5D_index]);
    }
 
    if(data->dist6D_collect) {
        dist_6D_init(&data->dist6D, &offload_data->dist6D,
                     &offload_array[offload_data->offload_dist6D_index]);
    }
 
    if(data->distrho5D_collect) {
        dist_rho5D_init(&data->distrho5D, &offload_data->distrho5D,
                        &offload_array[offload_data->offload_distrho5D_index]);
    }
 
    if(data->distrho6D_collect) {
        dist_rho6D_init(&data->distrho6D, &offload_data->distrho6D,
                        &offload_array[offload_data->offload_distrho6D_index]);
    }
 
    if(data->distCOM_collect) {
        dist_COM_init(&data->distCOM, &offload_data->distCOM,
                        &offload_array[offload_data->offload_distCOM_index]);
    }
 
    if(data->diagorb_collect) {
        diag_orb_init(&data->diagorb, &offload_data->diagorb,
                      &offload_array[offload_data->offload_diagorb_index]);
    }
 
    if(data->diagtrcof_collect) {
        diag_transcoef_init(
            &data->diagtrcof, &offload_data->diagtrcof,
            &offload_array[offload_data->offload_diagtrcof_index]);
    }
}
 
void diag_free(diag_data* data) {
    if(data->diagorb_collect) {
        diag_orb_free(&data->diagorb);
    }
    if(data->diagtrcof_collect) {
        diag_transcoef_free(&data->diagtrcof);
    }
}
 
void diag_update_fo(diag_data* data, B_field_data* Bdata, particle_simd_fo* p_f,
                    particle_simd_fo* p_i) {
    if(data->diagorb_collect) {
        diag_orb_update_fo(&data->diagorb, p_f, p_i);
    }
    if(data->dist5D_collect) {
        dist_5D_update_fo(&data->dist5D, p_f, p_i);
    }
    if(data->dist6D_collect) {
        dist_6D_update_fo(&data->dist6D, p_f, p_i);
    }
    if(data->distrho5D_collect) {
        dist_rho5D_update_fo(&data->distrho5D, p_f, p_i);
    }
    if(data->distrho6D_collect) {
        dist_rho6D_update_fo(&data->distrho6D, p_f, p_i);
    }
    if(data->distCOM_collect){
        dist_COM_update_fo(&data->distCOM, Bdata, p_f, p_i);
    }
    if(data->diagtrcof_collect){
        diag_transcoef_update_fo(&data->diagtrcof, p_f, p_i);
    }
}
 
void diag_update_gc(diag_data* data, B_field_data* Bdata, particle_simd_gc* p_f,
                    particle_simd_gc* p_i) {
    if(data->diagorb_collect) {
        diag_orb_update_gc(&data->diagorb, p_f, p_i);
    }
 
    if(data->dist5D_collect){
        dist_5D_update_gc(&data->dist5D, p_f, p_i);
    }
 
    if(data->dist6D_collect){
        dist_6D_update_gc(&data->dist6D, p_f, p_i);
    }
 
    if(data->distrho5D_collect){
        dist_rho5D_update_gc(&data->distrho5D, p_f, p_i);
    }
 
    if(data->distrho6D_collect){
        dist_rho6D_update_gc(&data->distrho6D, p_f, p_i);
    }
 
    if(data->distCOM_collect){
        dist_COM_update_gc(&data->distCOM, Bdata, p_f, p_i);
    }
 
    if(data->diagtrcof_collect){
        diag_transcoef_update_gc(&data->diagtrcof, p_f, p_i);
    }
}
 
void diag_update_ml(diag_data* data, particle_simd_ml* p_f,
                    particle_simd_ml* p_i) {
    if(data->diagorb_collect) {
        diag_orb_update_ml(&data->diagorb, p_f, p_i);
    }
 
    if(data->diagtrcof_collect){
        diag_transcoef_update_ml(&data->diagtrcof, p_f, p_i);
    }
}
 
void diag_sum(diag_offload_data* data, real* array1, real* array2) {
    if(data->diagorb_collect) {
        size_t arr_start = data->offload_diagorb_index;
        size_t arr_length = (size_t)(data->diagorb.Nfld)
            * (size_t)(data->diagorb.Nmrk) * (size_t)(data->diagorb.Npnt);
 
        memcpy(&(array1[arr_start+arr_length]),
               &(array2[arr_start]),
               arr_length*sizeof(real));
    }
 
    if(data->diagtrcof_collect) {
        size_t arr_start = data->offload_diagtrcof_index;
        size_t arr_length = (size_t)(3 * data->diagtrcof.Nmrk);
 
        memcpy(&(array1[arr_start+arr_length]),
               &(array2[arr_start]),
               arr_length*sizeof(real));
    }
 
    if(data->dist5D_collect){
        size_t start = data->offload_dist5D_index;
        size_t stop = start + (size_t)(data->dist5D.n_r)
            * (size_t)(data->dist5D.n_z) * (size_t)(data->dist5D.n_ppara)
            * (size_t)(data->dist5D.n_pperp) * (size_t)(data->dist5D.n_time)
            * (size_t)(data->dist5D.n_q);
        diag_arraysum(start, stop, array1, array2);
    }
 
    if(data->dist6D_collect){
        size_t start = data->offload_dist6D_index;
        size_t stop = start + (size_t)(data->dist6D.n_r)
            * (size_t)(data->dist6D.n_phi)
            * (size_t)(data->dist6D.n_z) * (size_t)(data->dist6D.n_pr)
            * (size_t)(data->dist6D.n_pphi) * (size_t)(data->dist6D.n_pz)
            * (size_t)(data->dist6D.n_time) * (size_t)(data->dist6D.n_q);
        diag_arraysum(start, stop, array1, array2);
    }
 
    if(data->distrho5D_collect){
        size_t start = data->offload_distrho5D_index;
        size_t stop = start + (size_t)(data->distrho5D.n_rho)
            * (size_t)(data->distrho5D.n_theta)
            * (size_t)(data->distrho5D.n_phi)
            * (size_t)(data->distrho5D.n_ppara)
            * (size_t)(data->distrho5D.n_pperp)
            * (size_t)(data->distrho5D.n_time)
            * (size_t)(data->distrho5D.n_q);
        diag_arraysum(start, stop, array1, array2);
    }
 
    if(data->distrho6D_collect){
        size_t start = data->offload_distrho6D_index;
        size_t stop = start + (size_t)(data->distrho6D.n_rho)
            * (size_t)(data->distrho6D.n_theta)
            * (size_t)(data->distrho6D.n_phi)  * (size_t)(data->distrho6D.n_pr)
            * (size_t)(data->distrho6D.n_pphi) * (size_t)(data->distrho6D.n_pz)
            * (size_t)(data->distrho6D.n_time) * (size_t)(data->distrho6D.n_q);
        diag_arraysum(start, stop, array1, array2);
    }
 
    if(data->distCOM_collect){
        size_t start = data->offload_distCOM_index;
        size_t stop = start + (size_t)(data->distCOM.n_mu)
            * (size_t)(data->distCOM.n_Ekin) * (size_t)(data->distCOM.n_Ptor);
        diag_arraysum(start, stop, array1, array2);
    }
}
 
void diag_arraysum(size_t start, size_t stop, real* array1, real* array2) {
    for(size_t i = start; i < stop; i++) {
        array1[i] += array2[i];
    }
}