simulate.c

Go to the documentation of this file.

 
#include <string.h>
#include <unistd.h>
#include "endcond.h"
#include "offload.h"
#include "particle.h"
#include "plasma.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 "copytogpu.h"
 
void sim_monitor(char* filename, volatile int* n, volatile int* finished);
 
void simulate(
    int id, int n_particles, particle_state* p, sim_offload_data* sim_offload,
    offload_package* offload_data, real* offload_array, int* int_offload_array,
    real* diag_offload_array) {
 
    // 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.                                          */
    /*                                                                        */
    /**************************************************************************/
    sim_data sim;
    sim_init(&sim, sim_offload);
 
#ifdef GPU
    if(sim_offload->sim_mode != 1) {
        print_err("Only GO mode ported to GPU. Please set SIM_MODE=1.");
        exit(1);
    }
    if(sim_offload->record_mode) {
        print_err("RECORD_MODE=1 not ported to GPU. Please disable it.");
        exit(1);
    }
    if(sim_offload->enable_atomic) {
        print_err("Atomic not yet ported to GPU. Please set ENABLE_ATOMIC=0.");
        exit(1);
    }
    if(sim_offload->enable_mhd) {
        print_err("MHD not yet ported to GPU. Please set ENABLE_MHD=0.");
        exit(1);
    }
    if(sim_offload->diag_offload_data.diagorb_collect) {
        print_err(
            "ENABLE_ORBITWRITE=1 not ported to GPU. Please disable it.");
        exit(1);
    }
    if(sim_offload->diag_offload_data.diagtrcof_collect) {
        print_err(
            "ENABLE_TRANSCOEF=1 not ported to GPU. Please disable it.");
        exit(1);
    }
#endif
    real* ptr; int* ptrint;
    offload_unpack(offload_data, offload_array,
                   sim_offload->B_offload_data.offload_array_length,
                   NULL, 0, &ptr, &ptrint);
    B_field_init(&sim.B_data, &sim_offload->B_offload_data, ptr);
 
    offload_unpack(offload_data, offload_array,
                   sim_offload->E_offload_data.offload_array_length,
                   NULL, 0, &ptr, &ptrint);
    E_field_init(&sim.E_data, &sim_offload->E_offload_data, ptr);
 
    offload_unpack(offload_data, offload_array,
                   sim_offload->plasma_offload_data.offload_array_length,
                   NULL, 0, &ptr, &ptrint);
    plasma_init(&sim.plasma_data, &sim_offload->plasma_offload_data, ptr);
 
    offload_unpack(offload_data, offload_array,
                   sim_offload->neutral_offload_data.offload_array_length,
                   NULL, 0, &ptr, &ptrint);
    neutral_init(&sim.neutral_data, &sim_offload->neutral_offload_data, ptr);
 
    offload_unpack(offload_data, offload_array,
                   sim_offload->wall_offload_data.offload_array_length,
                   int_offload_array,
                   sim_offload->wall_offload_data.int_offload_array_length,
                   &ptr, &ptrint);
    wall_init(&sim.wall_data, &sim_offload->wall_offload_data, ptr, ptrint);
 
    offload_unpack(offload_data, offload_array,
                   sim_offload->boozer_offload_data.offload_array_length,
                   NULL, 0, &ptr, &ptrint);
    boozer_init(&sim.boozer_data, &sim_offload->boozer_offload_data, ptr);
 
    offload_unpack(offload_data, offload_array,
                   sim_offload->mhd_offload_data.offload_array_length,
                   NULL, 0, &ptr, &ptrint);
    mhd_init(&sim.mhd_data, &sim_offload->mhd_offload_data, ptr);
 
    offload_unpack(offload_data, offload_array,
                   sim_offload->asigma_offload_data.offload_array_length,
                   NULL, 0, &ptr, &ptrint);
    asigma_init(&sim.asigma_data, &sim_offload->asigma_offload_data, ptr);
 
    /* Offload complete. Reset struct so it can be reused. */
    offload_data->unpack_pos     = 0;
    offload_data->int_unpack_pos = 0;
 
    diag_init(&sim.diag_data, &sim_offload->diag_offload_data,
              diag_offload_array);
    simulate_copy_to_gpu(&sim);
 
    /**************************************************************************/
    /* 2. Meta data (e.g. random number generator) is initialized.            */
    /*                                                                        */
    /**************************************************************************/
    random_init(&sim.random_data, 0);
 
    /**************************************************************************/
    /* 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());
 
    /**************************************************************************/
    /* 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);
                }
                else {
                    OMP_PARALLEL_CPU_ONLY
                    simulate_gc_fixed(&pq, &sim);
                }
            }
            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(id == 0) {
                char filename[519], outfn[256];
                strcpy(outfn, sim_offload->hdf5_out);
                outfn[strlen(outfn)-3] = '\0';
                sprintf(filename, "%s_%s.stdout", outfn, sim_offload->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(id == 0) {
                    char filename[519], outfn[256];
                    strcpy(outfn, sim_offload->hdf5_out);
                    outfn[strlen(outfn)-3] = '\0';
                    sprintf(filename, "%s_%s.stdout", outfn, sim_offload->qid);
                    sim_monitor(filename, &pq.n, &pq.finished);
                }
            }
        }
#endif
    }
 
    /**************************************************************************/
    /* 7. Simulation data is deallocated.                                     */
    /**************************************************************************/
    free(pq.p);
    diag_free(&sim.diag_data);
 
    print_out(VERBOSE_NORMAL, "Simulation complete.\n");
}
 
void simulate_init_offload(sim_offload_data* sim) {
    if(sim->disable_gctransform) {
        gctransform_setorder(0);
    }
    asigma_extrapolate(sim->enable_atomic==2);
}
 
void sim_init(sim_data* sim, sim_offload_data* offload_data) {
    sim->sim_mode             = offload_data->sim_mode;
    sim->enable_ada           = offload_data->enable_ada;
    sim->record_mode          = offload_data->record_mode;
 
    sim->fix_usrdef_use       = offload_data->fix_usrdef_use;
    sim->fix_usrdef_val       = offload_data->fix_usrdef_val;
    sim->fix_gyrodef_nstep    = offload_data->fix_gyrodef_nstep;
 
    sim->ada_tol_orbfol       = offload_data->ada_tol_orbfol;
    sim->ada_tol_clmbcol      = offload_data->ada_tol_clmbcol;
    sim->ada_max_drho         = offload_data->ada_max_drho;
    sim->ada_max_dphi         = offload_data->ada_max_dphi;
 
    sim->enable_orbfol        = offload_data->enable_orbfol;
    sim->enable_clmbcol       = offload_data->enable_clmbcol;
    sim->enable_mhd           = offload_data->enable_mhd;
    sim->enable_atomic        = offload_data->enable_atomic;
    sim->disable_gctransform  = offload_data->disable_gctransform;
    sim->disable_energyccoll  = offload_data->disable_energyccoll;
    sim->disable_pitchccoll   = offload_data->disable_pitchccoll;
    sim->disable_gcdiffccoll  = offload_data->disable_gcdiffccoll;
    sim->reverse_time         = offload_data->reverse_time;
 
    sim->endcond_active       = offload_data->endcond_active;
    sim->endcond_lim_simtime  = offload_data->endcond_lim_simtime;
    sim->endcond_max_mileage  = offload_data->endcond_max_mileage;
    sim->endcond_max_cputime  = offload_data->endcond_max_cputime;
    sim->endcond_min_rho      = offload_data->endcond_min_rho;
    sim->endcond_max_rho      = offload_data->endcond_max_rho;
    sim->endcond_min_ekin     = offload_data->endcond_min_ekin;
    sim->endcond_min_thermal  = offload_data->endcond_min_thermal;
    sim->endcond_max_tororb   = offload_data->endcond_max_tororb;
    sim->endcond_max_polorb   = offload_data->endcond_max_polorb;
    sim->endcond_torandpol    = offload_data->endcond_torandpol;
 
    mccc_init(&sim->mccc_data, !sim->disable_energyccoll,
              !sim->disable_pitchccoll, !sim->disable_gcdiffccoll);
 
}
 
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);
}