simulate.c

Go to the documentation of this file.

#include <string.h>
#include <unistd.h>
#include "endcond.h"
#include "particle.h"
#include "plasma.h"
#include "wall.h"
#include "boozer.h"
#include "mhd.h"
#include "neutral.h"
#include "B_field.h"
#include "E_field.h"
#include "random.h"
#include "simulate.h"
#include "print.h"
#include "simulate/simulate_ml_adaptive.h"
#include "simulate/simulate_gc_adaptive.h"
#include "simulate/simulate_gc_fixed.h"
#include "simulate/simulate_fo_fixed.h"
#include "simulate/mccc/mccc.h"
#include "gctransform.h"
#include "asigma.h"
#include "rfof.h"
 
void sim_monitor(char* filename, volatile int* n, volatile int* finished);

void simulate(int n_particles, particle_state* p, sim_data* sim) {
 
    // Size = NSIMD on CPU and Size = Total number of particles on GPU
    int n_queue_size;
#ifdef GPU
    n_queue_size = n_particles;
#else
    n_queue_size = NSIMD;
#endif
    /**************************************************************************/
    /* 1. Input offload data is unpacked and initialized by calling           */
    /*    respective init functions.                                          */
    /*                                                                        */
    /**************************************************************************/
 
    simulate_init(sim);
 
#ifdef GPU
    if((sim->sim_mode != 1) && (sim->sim_mode != 2)) {
        print_err("Only GO and GC mode ported to GPU. Please set SIM_MODE=1 or 2.");
        exit(1);
    }
    if(sim->enable_icrh) {
        print_err("ENABLE_ICRH=1 not ported to GPU. Please disable it.");
        exit(1);
    }
    if(sim->record_mode) {
        print_err("RECORD_MODE=1 not ported to GPU. Please disable it.");
        exit(1);
    }
    if(sim->enable_atomic) {
        print_err("Atomic not yet ported to GPU. Please set ENABLE_ATOMIC=0.");
        exit(1);
    }
    if(sim->enable_mhd) {
        print_err("MHD not yet ported to GPU. Please set ENABLE_MHD=0.");
        exit(1);
    }
    if(sim->diag_data.diagorb_collect) {
        print_err(
            "ENABLE_ORBITWRITE=1 not ported to GPU. Please disable it.");
        exit(1);
    }
    if(sim->diag_data.diagtrcof_collect) {
        print_err(
            "ENABLE_TRANSCOEF=1 not ported to GPU. Please disable it.");
        exit(1);
    }
    
#endif
 
    diag_init(&sim->diag_data, n_particles);
    GPU_MAP_TO_DEVICE(sim[0:1])
    B_field_offload(&sim->B_data);
    E_field_offload(&sim->E_data);
    plasma_offload(&sim->plasma_data);
    neutral_offload(&sim->neutral_data);
    wall_offload(&sim->wall_data);
    boozer_offload(&sim->boozer_data);
    mhd_offload(&sim->mhd_data);
    asigma_offload(&sim->asigma_data);
    diag_offload(&sim->diag_data);
 
    /**************************************************************************/
    /* 2. Meta data (e.g. random number generator) is initialized.            */
    /*                                                                        */
    /**************************************************************************/
#ifdef GPU
    /* On GPU we need one RNG state per parallel slot (thread/lane).
     * data_size is a safety factor (number of RNGs used per particle/loop).
     * Total number of RNG states = data_size * n_queue_size.
     */
    int data_size = 5;
    /* Allocate an array of RNG states (one per parallel element) */
    sim->random_data = malloc(data_size * n_queue_size * sizeof(random_data));
    /* Initialize each RNG state with a different seed.
     * Using i ensures independent sequences across threads.
     */
    for(int i = 0; i < data_size * n_queue_size; ++i) {
        random_init(&sim->random_data[i], i);
    }
    /* Transfer RNG states to the GPU device memory */
    GPU_MAP_TO_DEVICE(sim->random_data[0:data_size * n_queue_size])
#else
    sim->random_data = malloc(sizeof(random_data));
    random_init(sim->random_data, 0);
#endif
    /**************************************************************************/
    /* 3. Markers are put into simulation queue.                              */
    /*                                                                        */
    /**************************************************************************/
    particle_queue pq;
 
    pq.n = 0;
    for(int i = 0; i < n_particles; i++) {
        pq.n++;
    }
 
    pq.p = (particle_state**) malloc(pq.n * sizeof(particle_state*));
    pq.finished = 0;
 
    pq.next = 0;
    for(int i = 0; i < n_particles; i++) {
        pq.p[pq.next++] = &p[i];
 
    }
    pq.next = 0;
 
    print_out(VERBOSE_NORMAL, "Simulation begins; %d threads.\n",
              omp_get_max_threads());
    fflush(stdout);
 
    /**************************************************************************/
    /* 4. Threads are spawned. One thread is dedicated for monitoring         */
    /*    progress, if monitoring is active.                                  */
    /*                                                                        */
    /**************************************************************************/
#ifndef GPU
    omp_set_max_active_levels(2);
#endif
#if !defined(GPU) && VERBOSE > 1
    #pragma omp parallel sections num_threads(2)
    {
        #pragma omp section
#endif
        {
            /******************************************************************/
            /* 5. Other threads execute marker simulation using the mode the  */
            /*    user has chosen.                                            */
            /*                                                                */
            /******************************************************************/
            if(pq.n > 0 && (sim->sim_mode == simulate_mode_gc
                        || sim->sim_mode == simulate_mode_hybrid)) {
                if(sim->enable_ada) {
            OMP_PARALLEL_CPU_ONLY
            simulate_gc_adaptive(&pq, sim, n_queue_size);
                }
                else {
                    OMP_PARALLEL_CPU_ONLY
                    simulate_gc_fixed(&pq, sim, n_queue_size);
                }
            }
            else if(pq.n > 0 && sim->sim_mode == simulate_mode_fo) {
                OMP_PARALLEL_CPU_ONLY
                simulate_fo_fixed(&pq, sim, n_queue_size);
            }
            else if(pq.n > 0 && sim->sim_mode == simulate_mode_ml) {
                OMP_PARALLEL_CPU_ONLY
                simulate_ml_adaptive(&pq, sim);
            }
        }
#if !defined(GPU) && VERBOSE > 1
        #pragma omp section
        {
            /* Update progress until simulation is complete.             */
            /* Trim .h5 from filename and replace it with _<QID>.stdout  */
            if(sim->mpi_rank == sim->mpi_root) {
                char filename[519], outfn[256];
                strcpy(outfn, sim->hdf5_out);
                outfn[strlen(outfn)-3] = '\0';
                sprintf(filename, "%s_%s.stdout", outfn, sim->qid);
                sim_monitor(filename, &pq.n, &pq.finished);
            }
        }
    }
#endif
 
    /**************************************************************************/
    /* 6. (If hybrid mode is active) Markers with hybrid end condition active */
    /*    are placed on a new queue, and they have their end condition        */
    /*    deactivated and they are simulated with simulate_fo_fixed.c until   */
    /*    they have met some other end condition. Threads are spawned and     */
    /*    progress is monitored as previously.                                */
    /*                                                                        */
    /**************************************************************************/
    int n_new = 0;
    if(sim->sim_mode == simulate_mode_hybrid) {
 
        /* Determine the number markers that should be run
         * in fo after previous gc simulation */
        for(int i = 0; i < pq.n; i++) {
            if(pq.p[i]->endcond == endcond_hybrid) {
                /* Check that there was no wall between when moving from
                   gc to fo */
                real w_coll;
                int tile = wall_hit_wall(pq.p[i]->r, pq.p[i]->phi, pq.p[i]->z,
                        pq.p[i]->rprt, pq.p[i]->phiprt, pq.p[i]->zprt,
                                         &sim->wall_data, &w_coll);
                if(tile > 0) {
                    pq.p[i]->walltile = tile;
                    pq.p[i]->endcond |= endcond_wall;
                }
                else {
                    n_new++;
                }
            }
        }
    }
    if(n_new > 0) {
 
        /* Reset hybrid marker end condition */
        for(int i = 0; i < pq.n; i++) {
            if(pq.p[i]->endcond & endcond_hybrid) {
                pq.p[i]->endcond ^= endcond_hybrid;
            }
        }
        pq.next = 0;
        pq.finished = 0;
 
#if !defined(GPU) && VERBOSE > 1
        #pragma omp parallel sections num_threads(2)
        {
            #pragma omp section
#endif
            {
                OMP_PARALLEL_CPU_ONLY
                simulate_fo_fixed(&pq, sim, n_queue_size);
            }
#if !defined(GPU) && VERBOSE > 1
            #pragma omp section
            {
                /* Trim .h5 from filename and replace it with _<qid>.stdout */
                if(sim->mpi_rank == sim->mpi_root) {
                    char filename[519], outfn[256];
                    strcpy(outfn, sim->hdf5_out);
                    outfn[strlen(outfn)-3] = '\0';
                    sprintf(filename, "%s_%s.stdout", outfn, sim->qid);
                    sim_monitor(filename, &pq.n, &pq.finished);
                }
            }
        }
#endif
    }
 
    /**************************************************************************/
    /* 7. Simulation data is deallocated.                                     */
    /**************************************************************************/
    free(pq.p);
 
 
    /**************************************************************************/
    /* 8. Execution returns to host where this function was called.           */
    /*                                                                        */
    /**************************************************************************/
 
    diag_onload(&sim->diag_data);
    print_out(VERBOSE_NORMAL, "Simulation complete.\n");
}

void simulate_init(sim_data* sim) {
 
    mccc_init(&sim->mccc_data, !sim->disable_energyccoll,
              !sim->disable_pitchccoll, !sim->disable_gcdiffccoll);
 
    if(sim->disable_gctransform) {
        gctransform_setorder(0);
    }
    asigma_extrapolate(sim->enable_atomic==2);
}

void sim_monitor(char* filename, volatile int* n, volatile int* finished) {
    /* Open a file for writing simulation progress */
    FILE *f = fopen(filename, "w");
    if (f == NULL) {
        print_out(VERBOSE_DEBUG,
                  "Warning. %s could not be opened for progress updates.\n",
                  filename);
        return;
    }
 
    real time_sim_started = A5_WTIME;
    int stopflag = 1; /* Ensures progress is written one last time at 100% */
    int n_temp, finished_temp; /* Use these to store volatile variables so that
                                  their value does not change during one loop */
    while(stopflag) {
        n_temp = *n;
        finished_temp = *finished;
        real fracprog = ((real) finished_temp)/n_temp;
        real timespent = (A5_WTIME)-time_sim_started;
 
        if(n_temp == finished_temp) {
            stopflag = 0;
        }
 
        if(fracprog == 0) {
            fprintf(f, "No marker has finished simulation yet. "
                    "Time spent: %.2f h\n", timespent/3600);
        }
        else {
            fprintf(f, "Progress: %d/%d, %.2f %%. Time spent: %.2f h, "
                    "estimated time to finish: %.2f h\n", finished_temp, n_temp,
                    100*fracprog, timespent/3600,
                    (1/fracprog-1)*timespent/3600);
        }
        fflush(f);
        sleep(A5_PRINTPROGRESSINTERVAL);
    }
 
    fprintf(f, "Simulation finished.\n");
    fclose(f);
}