fadvance_core.cpp

/*
# =============================================================================
# Copyright (c) 2016 - 2021 Blue Brain Project/EPFL
#
# See top-level LICENSE file for details.
# =============================================================================.
*/

#include <cstdlib>
#include <functional>

#include "coreneuron/coreneuron.hpp"
#include "coreneuron/nrnconf.h"
#include "coreneuron/apps/corenrn_parameters.hpp"
#include "coreneuron/sim/multicore.hpp"
#include "coreneuron/mpi/nrnmpi.h"
#include "coreneuron/sim/fast_imem.hpp"
#include "coreneuron/gpu/nrn_acc_manager.hpp"
#include "coreneuron/io/reports/nrnreport.hpp"
#include "coreneuron/network/netcvode.hpp"
#include "coreneuron/network/netpar.hpp"
#include "coreneuron/network/partrans.hpp"
#include "coreneuron/utils/nrnoc_aux.hpp"
#include "coreneuron/utils/progressbar/progressbar.h"
#include "coreneuron/utils/profile/profiler_interface.h"
#include "coreneuron/io/nrn2core_direct.h"
#include "coreneuron/io/nrn_checkpoint.hpp"

// Do an auto checkpoint only if execution lasted longer than this var (secs)
#define CHECKPOINT_MIN_RUNTIME (4 * 3600)  // 4h

namespace coreneuron {

extern corenrn_parameters corenrn_param;
static void* nrn_fixed_step_thread(NrnThread*);
static void* nrn_fixed_step_group_thread(NrnThread*, int, int, int&);
static bool nrn_auto_checkpoint();
static time_t sim_start_time;


void dt2thread(double adt) { /* copied from nrnoc/fadvance.c */
    if (adt != nrn_threads[0]._dt) {
        for (int i = 0; i < nrn_nthread; ++i) {
            NrnThread* nt = nrn_threads + i;
            nt->_t = t;
            nt->_dt = dt;
            if (secondorder) {
                nt->cj = 2.0 / dt;
            } else {
                nt->cj = 1.0 / dt;
            }
        }
    }
}

void nrn_fixed_step_minimal() { /* not so minimal anymore with gap junctions */
    if (t != nrn_threads->_t) {
        dt2thread(-1.);
    } else {
        dt2thread(dt);
    }
    nrn_thread_table_check();
    nrn_multithread_job(nrn_fixed_step_thread);
    if (nrn_have_gaps) {
        nrnmpi_v_transfer();
        nrn_multithread_job(nrn_fixed_step_lastpart);
    }
#if NRNMPI
    if (nrn_threads[0]._stop_stepping) {
        nrn_spike_exchange(nrn_threads);
    }
#endif

#if defined(ENABLE_BIN_REPORTS) || defined(ENABLE_SONATA_REPORTS)
    nrn_flush_reports(nrn_threads[0]._t);
#endif
    t = nrn_threads[0]._t;
}

/* better cache efficiency since a thread can do an entire minimum delay
integration interval before joining
*/
/// --> Coreneuron
static progressbar* progress;

void initialize_progress_bar(int nstep) {
    if (nrnmpi_myid == 0 && !corenrn_param.is_quiet()) {
        printf("\n");
        progress = progressbar_new(" psolve", nstep);
    }
}

void update_progress_bar(int step, double time) {
    if (nrnmpi_myid == 0 && !corenrn_param.is_quiet()) {
        progressbar_update(progress, step, time);
    }
}

void finalize_progress_bar() {
    if (nrnmpi_myid == 0 && !corenrn_param.is_quiet()) {
        progressbar_finish(progress);
    }
}

void nrn_fixed_single_steps_minimal(int total_sim_steps, double tstop) {
    const int progressbar_update_interval = 5;
    static int current_steps = 0;
    initialize_progress_bar(total_sim_steps);
#if NRNMPI
    double updated_tstop = tstop - dt;
    nrn_assert(nrn_threads->_t <= tstop);
    // It may very well be the case that we do not advance at all
    while (nrn_threads->_t <= updated_tstop) {
#else
    double updated_tstop = tstop - .5 * dt;
    while (nrn_threads->_t < updated_tstop) {
#endif
        nrn_fixed_step_minimal();
        if (stoprun) {
            nrn_auto_checkpoint();
            break;
        }
        current_steps++;
        if (!(current_steps % progressbar_update_interval)) {
            update_progress_bar(current_steps, nrn_threads[0]._t);
        }
    }
    finalize_progress_bar();
}

void nrn_fixed_step_group_minimal(int total_sim_steps) {
    static int current_steps = 0;
    dt2thread(dt);
    nrn_thread_table_check();
    int step_group_n = total_sim_steps;
    int step_group_begin = 0;
    int step_group_end = 0;
    initialize_progress_bar(step_group_n);

    while (step_group_end < step_group_n) {
        nrn_multithread_job(nrn_fixed_step_group_thread,
                            step_group_n,
                            step_group_begin,
                            step_group_end);
#if NRNMPI
        nrn_spike_exchange(nrn_threads);
#endif

#if defined(ENABLE_BIN_REPORTS) || defined(ENABLE_SONATA_REPORTS)
        nrn_flush_reports(nrn_threads[0]._t);
#endif
        if (stoprun) {
            nrn_auto_checkpoint();
            break;
        }
        current_steps++;
        step_group_begin = step_group_end;
        update_progress_bar(step_group_end, nrn_threads[0]._t);
    }
    t = nrn_threads[0]._t;
    finalize_progress_bar();
}

static void* nrn_fixed_step_group_thread(NrnThread* nth,
                                         int step_group_max,
                                         int step_group_begin,
                                         int& step_group_end) {
    nth->_stop_stepping = 0;
    for (int i = step_group_begin; i < step_group_max; ++i) {
        nrn_fixed_step_thread(nth);
        if (nth->_stop_stepping) {
            if (nth->id == 0) {
                step_group_end = i + 1;
            }
            nth->_stop_stepping = 0;
            return nullptr;
        }
    }
    if (nth->id == 0) {
        step_group_end = step_group_max;
    }
    return nullptr;
}

void update(NrnThread* _nt) {
    double* vec_v = &(VEC_V(0));
    double* vec_rhs = &(VEC_RHS(0));
    int i2 = _nt->end;
#if defined(_OPENACC)
    int stream_id = _nt->stream_id;
#endif

    /* do not need to worry about linmod or extracellular*/
    if (secondorder) {
        // clang-format off

        #pragma acc parallel loop present(          \
            vec_v[0:i2], vec_rhs[0:i2])             \
            if (_nt->compute_gpu) async(stream_id)
        // clang-format on
        for (int i = 0; i < i2; ++i) {
            vec_v[i] += 2. * vec_rhs[i];
        }
    } else {
        // clang-format off

        #pragma acc parallel loop present(              \
                vec_v[0:i2], vec_rhs[0:i2])             \
                if (_nt->compute_gpu) async(stream_id)
        // clang-format on
        for (int i = 0; i < i2; ++i) {
            vec_v[i] += vec_rhs[i];
        }
    }

    // update_matrix_to_gpu(_nt);

    if (_nt->tml) {
        assert(_nt->tml->index == CAP);
        nrn_cur_capacitance(_nt, _nt->tml->ml, _nt->tml->index);
    }
    if (nrn_use_fast_imem) {
        nrn_calc_fast_imem(_nt);
    }
}

void nonvint(NrnThread* _nt) {
    if (nrn_have_gaps) {
        Instrumentor::phase p("gap-v-transfer");
        nrnthread_v_transfer(_nt);
    }
    errno = 0;

    Instrumentor::phase_begin("state-update");
    for (auto tml = _nt->tml; tml; tml = tml->next)
        if (corenrn.get_memb_func(tml->index).state) {
            mod_f_t s = corenrn.get_memb_func(tml->index).state;
            std::string ss("state-");
            ss += nrn_get_mechname(tml->index);
            {
                Instrumentor::phase p(ss.c_str());
                (*s)(_nt, tml->ml, tml->index);
            }
#ifdef DEBUG
            if (errno) {
                hoc_warning("errno set during calculation of states", nullptr);
            }
#endif
        }
    Instrumentor::phase_end("state-update");
}

void nrn_ba(NrnThread* nt, int bat) {
    for (auto tbl = nt->tbl[bat]; tbl; tbl = tbl->next) {
        mod_f_t f = tbl->bam->f;
        int type = tbl->bam->type;
        Memb_list* ml = tbl->ml;
        (*f)(nt, ml, type);
    }
}

void nrncore2nrn_send_init() {
    if (nrn2core_trajectory_values_ == nullptr) {
        // standalone execution : no callbacks
        return;
    }
    // if per time step transfer, need to call nrn_record_init() in NEURON.
    // if storing full trajectories in CoreNEURON, need to initialize
    // vsize for all the trajectory requests.
    (*nrn2core_trajectory_values_)(-1, 0, nullptr, 0.0);
    for (int tid = 0; tid < nrn_nthread; ++tid) {
        NrnThread& nt = nrn_threads[tid];
        if (nt.trajec_requests) {
            nt.trajec_requests->vsize = 0;
        }
    }
}

void nrncore2nrn_send_values(NrnThread* nth) {
    if (nrn2core_trajectory_values_ == nullptr) {
        // standalone execution : no callbacks
        return;
    }

    TrajectoryRequests* tr = nth->trajec_requests;
    if (tr) {
        // \todo Check if user has requested voltages for this NrnThread object.
        //       Currently we are updating voltages if there is any trajectory
        //       requested by NEURON.
        update_voltage_from_gpu(nth);

        if (tr->varrays) {  // full trajectories into Vector data
            double** va = tr->varrays;
            int vs = tr->vsize++;
            assert(vs < tr->bsize);
            for (int i = 0; i < tr->n_trajec; ++i) {
                va[i][vs] = *(tr->gather[i]);
            }
        } else if (tr->scatter) {  // scatter to NEURON and notify each step.
            nrn_assert(nrn2core_trajectory_values_);
            for (int i = 0; i < tr->n_trajec; ++i) {
                *(tr->scatter[i]) = *(tr->gather[i]);
            }
            (*nrn2core_trajectory_values_)(nth->id, tr->n_pr, tr->vpr, nth->_t);
        }
    }
}

static void* nrn_fixed_step_thread(NrnThread* nth) {
    /* check thresholds and deliver all (including binqueue)
       events up to t+dt/2 */
    Instrumentor::phase_begin("timestep");

    {
        Instrumentor::phase p("deliver_events");
        deliver_net_events(nth);
    }

    nth->_t += .5 * nth->_dt;

    if (nth->ncell) {
#if defined(_OPENACC)
        int stream_id = nth->stream_id;
        /*@todo: do we need to update nth->_t on GPU: Yes (Michael, but can launch kernel) */
        // clang-format off

        #pragma acc update device(nth->_t) if (nth->compute_gpu) async(stream_id)
        #pragma acc wait(stream_id)
// clang-format on
#endif
        fixed_play_continuous(nth);

        {
            Instrumentor::phase p("setup_tree_matrix");
            setup_tree_matrix_minimal(nth);
        }

        {
            Instrumentor::phase p("matrix-solver");
            nrn_solve_minimal(nth);
        }

        {
            Instrumentor::phase p("second_order_cur");
            second_order_cur(nth, secondorder);
        }

        {
            Instrumentor::phase p("update");
            update(nth);
        }
    }
    if (!nrn_have_gaps) {
        nrn_fixed_step_lastpart(nth);
    }
    Instrumentor::phase_end("timestep");
    return nullptr;
}

void* nrn_fixed_step_lastpart(NrnThread* nth) {
    nth->_t += .5 * nth->_dt;

    if (nth->ncell) {
#if defined(_OPENACC)
        int stream_id = nth->stream_id;
        /*@todo: do we need to update nth->_t on GPU */
        // clang-format off

        #pragma acc update device(nth->_t) if (nth->compute_gpu) async(stream_id)
        #pragma acc wait(stream_id)
// clang-format on
#endif

        fixed_play_continuous(nth);
        nonvint(nth);
        nrncore2nrn_send_values(nth);
        nrn_ba(nth, AFTER_SOLVE);
        nrn_ba(nth, BEFORE_STEP);
    } else {
        nrncore2nrn_send_values(nth);
    }

    {
        Instrumentor::phase p("deliver_events");
        nrn_deliver_events(nth); /* up to but not past texit */
    }

    return nullptr;
}

/**
 * \brief Does a checkpoint of the simulation in enough time has passed
 * \return True if a checkpoint was performed. False otherwise (not enough elapsed time)
 */
static bool nrn_auto_checkpoint() {
    time_t cur_time = time(NULL);
    int elapsed_secs = difftime(sim_start_time, cur_time);
    if (elapsed_secs < CHECKPOINT_MIN_RUNTIME) {
        return false;
    }
    // Write to tmp location first because allocated time may not be enough to complete
    const auto ckpt_tmp = corenrn_param.outpath + "/_corenrn_ckpt_dirty",
               ckpt_dir = corenrn_param.outpath + "/_corenrn_ckpt";
    Instrumentor::phase p("Checkpointing");
    write_checkpoint(nrn_threads, nrn_nthread, ckpt_tmp.c_str());
    system(("/bin/rm -rf '" + ckpt_dir + "'; " + "/bin/mv '" + ckpt_tmp + "' '" + ckpt_dir + "'")
               .c_str());
    return true;
}


}  // namespace coreneuron