tfEngine.h

Go to the documentation of this file.

/*******************************************************************************
 * This file is part of mdcore.
 * Coypright (c) 2010 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
 * Coypright (c) 2017 Andy Somogyi (somogyie at indiana dot edu)
 * Copyright (c) 2022-2024 T.J. Sego
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
 
#ifndef _MDCORE_INCLUDE_TFENGINE_H_
#define _MDCORE_INCLUDE_TFENGINE_H_
 
#include "tf_platform.h"
#include <pthread.h>
#include "tfSpace.h"
#include "tf_cycle.h"
#include "tfBoundaryConditions.h"
#include <tfSubEngine.h>
#include <mutex>
#include <thread>
#include <set>
 
 
#define engine_bonds_chunk               100
#define engine_angles_chunk              100
#define engine_rigids_chunk              50
#define engine_dihedrals_chunk           100
#define engine_exclusions_chunk          100
#define engine_readbuff                  16384
#define engine_maxgpu                    10
#define engine_pshake_steps              20
#define engine_maxKcutoff                2
 
#define engine_split_MPI        1
#define engine_split_GPU        2
 
#define engine_bonded_maxnrthreads       16
#define engine_bonded_nrthreads          ((omp_get_num_threads()<engine_bonded_maxnrthreads)?omp_get_num_threads():engine_bonded_maxnrthreads)
 
#ifdef MDCORE_MAXNRTYPES
#define engine_maxnrtypes                MDCORE_MAXNRTYPES
#else
#define engine_maxnrtypes                128
#endif
 
 
namespace TissueForge { 
 
 
    /* some constants */
    enum EngineFlags {
        engine_flag_none                 = 0,
        engine_flag_static               = 1 << 0,
        engine_flag_localparts           = 1 << 1,
        engine_flag_cuda                 = 1 << 2,
        engine_flag_explepot             = 1 << 3,
        engine_flag_verlet               = 1 << 4,
        engine_flag_verlet_pairwise      = 1 << 5,
        engine_flag_affinity             = 1 << 6,
        engine_flag_prefetch             = 1 << 7,
        engine_flag_verlet_pseudo        = 1 << 8,
        engine_flag_shake                = 1 << 9,
        engine_flag_mpi                  = 1 << 10,
        engine_flag_parbonded            = 1 << 11,
        engine_flag_async                = 1 << 12,
        engine_flag_sets                 = 1 << 13,
        engine_flag_nullpart             = 1 << 14,
        engine_flag_initialized          = 1 << 15,
        engine_flag_velocity_clamp       = 1 << 16,
    };
 
    enum EngineIntegrator {
        FORWARD_EULER,
        RUNGE_KUTTA_4
    };
 
    enum {
        engine_timer_step = 0,
        engine_timer_kinetic, 
        engine_timer_prepare,
        engine_timer_verlet,
        engine_timer_exchange1,
        engine_timer_nonbond,
        engine_timer_bonded,
        engine_timer_bonded_sort,
        engine_timer_bonds,
        engine_timer_angles,
        engine_timer_dihedrals,
        engine_timer_exclusions,
        engine_timer_advance,
        engine_timer_rigid,
        engine_timer_exchange2,
        engine_timer_shuffle,
        engine_timer_cuda_load,
        engine_timer_cuda_unload,
        engine_timer_cuda_dopairs,
        engine_timer_render,
        engine_timer_image_data,
        engine_timer_render_total,
        engine_timer_total,
        engine_timer_last
    };
 
 
    enum {
        ENGINE_TIMER_STEP           = 1 << 0,
        ENGINE_TIMER_PREPARE        = 1 << 1,
        ENGINE_TIMER_VERLET         = 1 << 2,
        ENGINE_TIMER_EXCHANGE1      = 1 << 3,
        ENGINE_TIMER_NONBOND        = 1 << 4,
        ENGINE_TIMER_BONDED         = 1 << 5,
        ENGINE_TIMER_BONDED_SORT    = 1 << 6,
        ENGINE_TIMER_BONDS          = 1 << 7,
        ENGINE_TIMER_ANGLES         = 1 << 8,
        ENGINE_TIMER_DIHEDRALS      = 1 << 9,
        ENGINE_TIMER_EXCLUSIONS     = 1 << 10,
        ENGINE_TIMER_ADVANCE        = 1 << 11,
        ENGINE_TIMER_RIGID          = 1 << 12,
        ENGINE_TIMER_EXCHANGE2      = 1 << 13,
        ENGINE_TIMER_SHUFFLE        = 1 << 14,
        ENGINE_TIMER_CUDA_LOAD      = 1 << 15,
        ENGINE_TIMER_CUDA_UNLOAD    = 1 << 16,
        ENGINE_TIMER_CUDA_DOPAIRS   = 1 << 17,
        ENGINE_TIMER_RENDER         = 1 << 18,
        ENGINE_TIMER_LAST           = 1 << 19
    };
 
    enum {
        // forces that set the persistent_force should
        // update values now. Otherwise, the integrator is
        // probably in a multi-step and should use the saved
        // value
        INTEGRATOR_UPDATE_PERSISTENTFORCE    = 1 << 0, 
 
        // intermediate flux values are being calculated between time steps
        INTEGRATOR_FLUX_SUBSTEP              = 1 << 1
    };
 
 
    struct CustomForce;
 
    typedef struct CAPI_EXPORT engine {
 
        unsigned int flags = 0;
 
        unsigned int integrator_flags;
 
    #ifdef HAVE_CUDA
        unsigned int flags_cuda;
    #endif
 
        struct space s;
 
        long time;
        FPTYPE dt;
 
        FPTYPE temperature;
 
        // Boltzmann constant
        FPTYPE K;
 
        static const int max_type;
        static int nr_types;
 
        static struct ParticleType *types;
 
        struct Potential **p, **p_cluster;
 
        struct Force **forces;
 
        struct Fluxes **fluxes;
 
        int nr_fluxsteps = 1;
        int step_flux = 0;
        FPTYPE dt_flux;
 
        pthread_mutex_t barrier_mutex;
        pthread_cond_t barrier_cond;
        pthread_cond_t done_cond;
        int barrier_count;
 
        int nr_runners;
 
        struct runner *runners;
 
        struct queue *queues;
        int nr_queues;
 
        int nodeID;
        int nr_nodes;
 
        struct engine_comm *send, *recv;
 
        std::set<unsigned int> pids_avail;
 
        struct Bond *bonds;
 
        int nr_bonds;
 
        int nr_active_bonds;
 
        int bonds_size;
 
        // mutex for anything that modifies the *number* of bonds.
        std::mutex bonds_mutex;
 
        struct exclusion *exclusions;
 
        int nr_exclusions, exclusions_size;
 
        struct rigid *rigids;
 
        int *part2rigid;
 
        int nr_rigids, rigids_size, nr_constr, rigids_local, rigids_semilocal;
 
        FPTYPE tol_rigid;
 
        struct Angle *angles;
 
        int nr_angles;
 
        int nr_active_angles;
 
        int angles_size;
 
        struct Dihedral *dihedrals;
 
        int nr_dihedrals;
 
        int nr_active_dihedrals;
 
        int dihedrals_size;
 
    #ifdef WITH_MPI
        MPI_Comm comm;
        pthread_mutex_t xchg_mutex;
        pthread_cond_t xchg_cond;
        short int xchg_started, xchg_running;
        pthread_t thread_exchg;
        pthread_mutex_t xchg2_mutex;
        pthread_cond_t xchg2_cond;
        short int xchg2_started, xchg2_running;
        pthread_t thread_exchg2;
    #endif
 
    #ifdef HAVE_CUDA
        void *sortlists_cuda[ engine_maxgpu ];
        int nr_pots_cuda, nr_pots_cluster_cuda, *pind_cuda[ engine_maxgpu ], *pind_cluster_cuda[ engine_maxgpu ], *offsets_cuda[ engine_maxgpu ];
        int nr_devices, devices[ engine_maxgpu ], nr_queues_cuda;
        float *forces_cuda[ engine_maxgpu ];
        void *parts_pos_cuda[engine_maxgpu];
        void *parts_vel_cuda[engine_maxgpu];
        void *parts_datai_cuda[engine_maxgpu];
        void *part_states_cuda[engine_maxgpu];
        int *part_species_flags_cuda[engine_maxgpu];
        void *parts_pos_cuda_local;
        void *parts_vel_cuda_local;
        void *parts_datai_cuda_local;
        void *part_states_cuda_local;
        int *part_species_flags_cuda_local;
        int *cells_cuda_local[ engine_maxgpu];
        int cells_cuda_nr[ engine_maxgpu ];
        int *counts_cuda[ engine_maxgpu ], *counts_cuda_local[ engine_maxgpu ];
        int *ind_cuda[ engine_maxgpu ], *ind_cuda_local[ engine_maxgpu ];
        struct task_cuda *tasks_cuda[ engine_maxgpu ];
        int *taboo_cuda[ engine_maxgpu ];
        int nrtasks_cuda[ engine_maxgpu ];
        void *streams[ engine_maxgpu ];
        int nr_blocks[engine_maxgpu], nr_threads[engine_maxgpu];
        int nr_fluxes_cuda, *fxind_cuda[engine_maxgpu];
        void **fluxes_cuda[engine_maxgpu];
        float *fluxes_next_cuda[engine_maxgpu];
        bool bonds_cuda = false;
        bool angles_cuda = false;
    #endif
 
        ticks timers[engine_timer_last];
 
        struct engine_set *sets;
        int nr_sets;
        
        FPTYPE wall_time;
        
        // bitmask of timers to show in performance counter output.
        uint32_t timers_mask;
        
        long timer_output_period;
 
        std::vector<CustomForce*> custom_forces;
 
        FPTYPE particle_max_dist_fraction;
 
        FPTYPE computed_volume;
 
        EngineIntegrator integrator;
 
        BoundaryConditions boundary_conditions;
 
        std::vector<SubEngine*> subengines;
        
        
        BoundaryConditionsArgsContainer *_init_boundary_conditions;
        int _init_cells[3];
    } engine;
 
 
    typedef struct engine_set {
 
        /* Counts of the different interaction types. */
        int nr_bonds, nr_angles, nr_dihedrals, nr_exclusions, weight;
 
        /* Lists of ID of the relevant bonded types. */
        struct Bond *bonds;
        struct Angle *angles;
        struct Dihedral *dihedrals;
        struct exclusion *exclusions;
 
        /* Nr of sets with which this set conflicts. */
        int nr_confl;
 
        /* IDs of the sets with which this set conflicts. */
        int *confl;
 
    } engine_set;
 
 
    typedef struct engine_comm {
 
        /* Size and count of the cellids. */
        int count, size;
 
        int *cellid;
 
    } engine_comm;
 
 
    /* associated functions */
 
    CAPI_FUNC(HRESULT) engine_addpot(struct engine *e, struct Potential *p, int i, int j);
 
    CAPI_FUNC(HRESULT) engine_addfluxes(struct engine *e, struct Fluxes *f, int i, int j);
 
    Fluxes *engine_getfluxes(struct engine *e, int i, int j);
 
    CAPI_FUNC(HRESULT) engine_add_singlebody_force(struct engine *e, struct Force *p, int typeId);
 
    CAPI_FUNC(int) engine_angle_alloc(struct engine *e, Angle **out);
 
    CAPI_FUNC(HRESULT) engine_angle_eval(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_barrier(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_bond_eval(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_bonded_eval(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_bonded_eval_sets(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_bonded_sets(struct engine *e, int max_sets);
 
    CAPI_FUNC(int) engine_dihedral_alloc(struct engine *e, Dihedral **out);
 
    CAPI_FUNC(HRESULT) engine_dihedral_eval(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exclusion_add(struct engine *e, int i, int j);
 
    CAPI_FUNC(HRESULT) engine_exclusion_eval(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exclusion_shrink(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_finalize(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_flush_ghosts(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_flush(struct engine *e);
 
    CAPI_FUNC(int) engine_gettype(struct engine *e, char *name);
 
    CAPI_FUNC(int) engine_gettype2(struct engine *e, char *name2);
 
    int engine_bond_alloc (struct engine *e, struct Bond **result);
 
    CAPI_FUNC(struct ParticleType*) engine_type(int id);
 
    CAPI_FUNC(HRESULT) engine_addpart(
        struct engine *e, 
        struct Particle *p,
        FPTYPE *x, 
        struct Particle **result
    );
 
    CAPI_FUNC(HRESULT) engine_addparts(struct engine *e, int nr_parts, struct Particle **parts, FPTYPE **x);
 
    CAPI_FUNC(int) engine_addtype(
        struct engine *e, 
        FPTYPE mass, 
        FPTYPE charge,
        const char *name, 
        const char *name2
    );
 
    CAPI_FUNC(HRESULT) engine_init(
        struct engine *e, 
        const FPTYPE *origin, 
        const FPTYPE *dim, 
        int *cells,
        FPTYPE cutoff, 
        BoundaryConditionsArgsContainer *boundaryConditions, 
        int max_type, 
        unsigned int flags, 
        unsigned int nr_fluxsteps
    );
 
    CAPI_FUNC(HRESULT) engine_reset(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_load_ghosts(
        struct engine *e, 
        FPTYPE *x, 
        FPTYPE *v, 
        int *type, 
        int *pid,
        int *vid, FPTYPE *q, 
        unsigned int *flags, 
        int N
    );
 
    CAPI_FUNC(HRESULT) engine_load(
        struct engine *e, 
        FPTYPE *x, 
        FPTYPE *v, 
        int *type, 
        int *pid, 
        int *vid,
        FPTYPE *charge, unsigned int *flags, int N
    );
 
    CAPI_FUNC(HRESULT) engine_nonbond_eval(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_fluxonly_eval(struct engine* e);
 
    CAPI_FUNC(HRESULT) engine_rigid_add(struct engine *e, int pid, int pjd, FPTYPE d);
 
    CAPI_FUNC(HRESULT) engine_rigid_eval(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_rigid_sort(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_rigid_unsort(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_shuffle(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_split_bisect(struct engine *e, int N, int particle_flags);
 
    CAPI_FUNC(HRESULT) engine_split(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_start(struct engine *e, int nr_runners, int nr_queues);
 
    CAPI_FUNC(HRESULT) engine_step(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_timers_reset(struct engine *e);
 
    CAPI_FUNC(int) engine_unload_marked(
        struct engine *e, 
        FPTYPE *x, 
        FPTYPE *v, 
        int *type, 
        int *pid,
        int *vid, 
        FPTYPE *q, 
        unsigned int *flags, 
        FPTYPE *epot, 
        int N
    );
 
    CAPI_FUNC(int) engine_unload_strays(
        struct engine *e, 
        FPTYPE *x, 
        FPTYPE *v, 
        int *type, 
        int *pid,
        int *vid, 
        FPTYPE *q, 
        unsigned int *flags, 
        FPTYPE *epot, 
        int N
    );
 
    CAPI_FUNC(int) engine_unload(
        struct engine *e, 
        FPTYPE *x, 
        FPTYPE *v, 
        int *type, 
        int *pid, 
        int *vid,
        FPTYPE *charge, 
        unsigned int *flags, 
        FPTYPE *epot, 
        int N
    );
 
    CAPI_FUNC(HRESULT) engine_verlet_update(struct engine *e);
 
    CAPI_FUNC(int) engine_next_partid(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_next_partids(struct engine *e, int nr_ids, int *ids);
 
    HRESULT engine_force_prep(struct engine *e);
 
    HRESULT engine_force(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_del_particle(struct engine *e, int pid);
 
    // keep track of how frequently step is called, get average
    // steps per second, averaged over past 10 steps.
    CAPI_FUNC(FPTYPE) engine_steps_per_second();
 
 
    CAPI_FUNC(void) engine_dump();
 
    #define ENGINE_DUMP(msg) {std::cout<<msg<<std::endl; engine_dump();}
 
    CAPI_FUNC(FPTYPE) engine_kinetic_energy(struct engine *e);
 
    CAPI_FUNC(FPTYPE) engine_temperature(struct engine *e);
 
    #ifdef WITH_MPI
 
    CAPI_FUNC(HRESULT) engine_init_mpi(
        struct engine *e, 
        const FPTYPE *origin, 
        const FPTYPE *dim, 
        FPTYPE *L,
        FPTYPE cutoff, 
        unsigned int period, 
        int max_type, 
        unsigned int flags, 
        MPI_Comm comm,
        int rank
    );
 
    CAPI_FUNC(HRESULT) engine_exchange(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exchange_async(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exchange_async_run(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exchange_incomming(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exchange_rigid(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exchange_rigid_async(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exchange_rigid_async_run(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exchange_rigid_wait(struct engine *e);
 
    CAPI_FUNC(HRESULT) engine_exchange_wait(struct engine *e);
    #endif
 
    #if defined(HAVE_CUDA)
 
 
    namespace cuda { 
 
 
        CAPI_FUNC(HRESULT) engine_nonbond_cuda(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_load(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_load_parts(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_unload_parts(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_load_pots(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_unload_pots(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_refresh_particles(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_allocate_particle_states(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_finalize_particle_states(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_refresh_particle_states(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_refresh_pots(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_load_fluxes(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_unload_fluxes(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_refresh_fluxes(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_boundary_conditions_refresh(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_queues_finalize(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_finalize(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_refresh(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_cuda_setthreads(struct engine *e, int id, int nr_threads);
 
        CAPI_FUNC(HRESULT) engine_cuda_setblocks(struct engine *e, int id, int nr_blocks);
 
        CAPI_FUNC(HRESULT) engine_cuda_setdevice(struct engine *e, int id);
 
        CAPI_FUNC(HRESULT) engine_cuda_setdevices(struct engine *e, int nr_devices, int *ids);
 
        CAPI_FUNC(HRESULT) engine_cuda_cleardevices(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_split_gpu(struct engine *e, int N, int flags);
 
        CAPI_FUNC(HRESULT) engine_toCUDA(struct engine *e);
 
        CAPI_FUNC(HRESULT) engine_fromCUDA(struct engine *e);
 
    }
    
    #endif
 
    #ifdef WITH_METIS
    CAPI_FUNC(int) engine_split_METIS(struct engine *e, int N, int flags);
    #endif
 
    FVector3 engine_origin();
    FVector3 engine_dimensions();
    FVector3 engine_center();
 
    // CAPI_DATA(engine) _Engine;
    extern CAPI_EXPORT engine _Engine;
 
    CPPAPI_FUNC(struct engine*) engine_get();
 
};
 
 
// inline definitions...
 
#include "tfParticle.h"
 
 
namespace TissueForge { 
 
 
    inline Particle *Particle::particle(int i) {
        return _Engine.s.partlist[this->parts[i]];
    };
 
    inline FVector3 Particle::global_position() {
        FVector3 position;
        space_getpos(&_Engine.s, this->id, position.data());
        return position;
    }
 
    inline void Particle::set_global_position(const FVector3& pos) {
        FPTYPE x[] = {pos.x(), pos.y(), pos.z()};
        space_setpos(&_Engine.s, this->id, x);
    }
 
    inline Particle *ParticleHandle::part() {
        return _Engine.s.partlist[this->id];
    };
 
    inline ParticleType *ParticleHandle::type() {
        return &_Engine.types[this->getTypeId()];
    }
 
    inline Particle *Particle_FromId(int id) {
        return _Engine.s.partlist[id];
    }
 
    inline Particle *ParticleType::particle(int i) {
        return _Engine.s.partlist[this->parts.parts[i]];
    }
 
};
 
#endif // _MDCORE_INCLUDE_TFENGINE_H_