velociraptor_interface.c 10.3 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2018 James Willis (james.s.willis@durham.ac.uk)
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 *
 * 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/>.
 *
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

23
/* Some standard headers. */
24
25
26
#include <errno.h>
#include <unistd.h>

27
28
29
/* This object's header. */
#include "velociraptor_interface.h"

30
31
32
/* Local includes. */
#include "common_io.h"

33
34
35
36
37
38
/**
 * @brief Initialise VELOCIraptor with input and output file names along with cosmological info needed to run.
 *
 * @param e The #engine.
 *
 */
39
void velociraptor_init(struct engine *e) {
40

41
#ifdef HAVE_VELOCIRAPTOR
James Willis's avatar
James Willis committed
42
    struct space *s = e->s;
43
44
45
46
    struct cosmoinfo cosmo_info;
    struct unitinfo unit_info;
    struct siminfo sim_info;
    
47
    /* Set cosmological constants. */
48
    cosmo_info.atime = e->cosmology->a;
49
50
51
52
    cosmo_info.littleh = e->cosmology->h;
    cosmo_info.Omega_m = e->cosmology->Omega_m;
    cosmo_info.Omega_b = e->cosmology->Omega_b;
    cosmo_info.Omega_Lambda = e->cosmology->Omega_lambda;
53
    cosmo_info.Omega_cdm = e->cosmology->Omega_m - e->cosmology->Omega_b;
54
    cosmo_info.w_de = e->cosmology->w;
55
56
57
58
59
60
61
62

    message("Scale factor: %e", cosmo_info.atime);
    message("Little h: %e", cosmo_info.littleh);
    message("Omega_m: %e", cosmo_info.Omega_m);
    message("Omega_b: %e", cosmo_info.Omega_b);
    message("Omega_Lambda: %e", cosmo_info.Omega_Lambda);
    message("Omega_cdm: %e", cosmo_info.Omega_cdm);
    message("w_de: %e", cosmo_info.w_de);
63
64

    /* Set unit conversions. */
65
66
67
68
69
70
    unit_info.lengthtokpc = 1.0;
    unit_info.velocitytokms = 1.0;
    unit_info.masstosolarmass = 1.0;    
    unit_info.energyperunitmass = 1.0;
    unit_info.gravity = e->physical_constants->const_newton_G;
    unit_info.hubbleunit = e->cosmology->H; /* TODO: double check this. */
71

72
73
74
    message("Length conversion factor: %e", unit_info.lengthtokpc);
    message("Velocity conversion factor: %e", unit_info.velocitytokms);
    message("Mass conversion factor: %e", unit_info.masstosolarmass);
James Willis's avatar
James Willis committed
75
    message("Potential conversion factor: %e", unit_info.energyperunitmass);
76
77
78
    message("G: %e", unit_info.gravity);
    message("H: %e", unit_info.hubbleunit);

79
80
    /* TODO: Find the total number of DM particles when running with star particles and BHs. */
    const int total_nr_dmparts = e->total_nr_gparts - e->total_nr_parts;
81
    
82
    /* Set simulation information. */
James Willis's avatar
James Willis committed
83
84
85
86
    if(e->s->periodic) {
        sim_info.period = unit_info.lengthtokpc * s->dim[0]; /* Physical size of box in VELOCIraptor units (kpc). */
    }
    else sim_info.period = 0.0;
87
    sim_info.zoomhigresolutionmass = -1.0; /* Placeholder. */
88
    sim_info.interparticlespacing = sim_info.period / cbrt(total_nr_dmparts);
89
90
    if(e->policy & engine_policy_cosmology) sim_info.icosmologicalsim = 1;
    else sim_info.icosmologicalsim = 0;
91
92
93
94
    sim_info.spacedimension[0] = unit_info.lengthtokpc * s->dim[0];
    sim_info.spacedimension[1] = unit_info.lengthtokpc * s->dim[1];
    sim_info.spacedimension[2] = unit_info.lengthtokpc * s->dim[2];
    sim_info.numcells = s->nr_cells;
95
96
97
98
99
100
101
102
103
 
    sim_info.cellwidth[0] = unit_info.lengthtokpc * s->cells_top[0].width[0];
    sim_info.cellwidth[1] = unit_info.lengthtokpc * s->cells_top[0].width[1];
    sim_info.cellwidth[2] = unit_info.lengthtokpc * s->cells_top[0].width[2];

    sim_info.icellwidth[0] = s->iwidth[0] / unit_info.lengthtokpc;
    sim_info.icellwidth[1] = s->iwidth[1] / unit_info.lengthtokpc;
    sim_info.icellwidth[2] = s->iwidth[2] / unit_info.lengthtokpc;
   
104
    /* Allocate and populate top-level cell locations. */
105
    if (posix_memalign((void **)&(e->cell_loc), 32,
106
107
108
109
                       s->nr_cells * sizeof(struct cell_loc)) != 0)
        error("Failed to allocate top-level cell locations for VELOCIraptor.");

    for(int i=0; i<s->nr_cells; i++) {
110
111
112
        e->cell_loc[i].loc[0] = unit_info.lengthtokpc * s->cells_top[i].loc[0];
        e->cell_loc[i].loc[1] = unit_info.lengthtokpc * s->cells_top[i].loc[1];
        e->cell_loc[i].loc[2] = unit_info.lengthtokpc * s->cells_top[i].loc[2];
113
114
    }

115
    sim_info.cell_loc = e->cell_loc;
116

James Willis's avatar
James Willis committed
117
    char configfilename[PARSER_MAX_LINE_SIZE], outputFileName[PARSER_MAX_LINE_SIZE + 128];
118
    parser_get_param_string(e->parameter_file, "StructureFinding:config_file_name", configfilename);
James Willis's avatar
James Willis committed
119
    snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s.VELOCIraptor", e->stfBaseName);
120
121
    
    message("Config file name: %s", configfilename);
122
123
124
125
    message("Period: %e", sim_info.period);
    message("Zoom high res mass: %e", sim_info.zoomhigresolutionmass);
    message("Inter-particle spacing: %e", sim_info.interparticlespacing);
    message("Cosmological: %d", sim_info.icosmologicalsim);
126
127
    message("Space dimensions: (%e,%e,%e)", sim_info.spacedimension[0], sim_info.spacedimension[1], sim_info.spacedimension[2]);
    message("No. of top-level cells: %d", sim_info.numcells);
128
    message("Top-level cell locations range: (%e,%e,%e) -> (%e,%e,%e)", sim_info.cell_loc[0].loc[0], sim_info.cell_loc[0].loc[1], sim_info.cell_loc[0].loc[2], sim_info.cell_loc[sim_info.numcells - 1].loc[0], sim_info.cell_loc[sim_info.numcells - 1].loc[1], sim_info.cell_loc[sim_info.numcells - 1].loc[2]);
129

James Willis's avatar
James Willis committed
130
    /* Initialise VELOCIraptor. */
131
    if(!InitVelociraptor(configfilename, outputFileName, cosmo_info, unit_info, sim_info)) error("Exiting. VELOCIraptor initialisation failed.");
132
133
#else
  error("SWIFT not configure to run with VELOCIraptor.");
134
#endif /* HAVE_VELOCIRAPTOR */
135

136
137
138
139
140
141
142
143
}

/**
 * @brief Run VELOCIraptor with current particle data.
 *
 * @param e The #engine.
 *
 */
144
void velociraptor_invoke(struct engine *e) {
145

146
#ifdef HAVE_VELOCIRAPTOR
James Willis's avatar
James Willis committed
147
    struct space *s = e->s;
148
    struct gpart *gparts = s->gparts;
149
    struct part *parts = s->parts;
150
151
    const size_t nr_gparts = s->nr_gparts;
    const size_t nr_hydro_parts = s->nr_parts;
James Willis's avatar
James Willis committed
152
    const int nr_cells = s->nr_cells;
James Willis's avatar
James Willis committed
153
    int *cell_node_ids = NULL;
154
155
156

    /* Allow thread to run on any core for the duration of the call to VELOCIraptor so that 
     * when OpenMP threads are spawned they can run on any core on the processor. */
157
    const int nr_cores = sysconf(_SC_NPROCESSORS_ONLN);
158
159
160
161
162
163
164
165
166
167
168
169
    cpu_set_t cpuset;
    pthread_t thread = pthread_self();

    /* Set affinity mask to include all cores on the CPU for VELOCIraptor. */
    CPU_ZERO(&cpuset);
    for (int j = 0; j < nr_cores; j++)
        CPU_SET(j, &cpuset);

    pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset);

    ticks tic = getticks();

170
171
172
173
174
175
    /* Allocate and populate array of cell node IDs. */
    if (posix_memalign((void **)&cell_node_ids, 32,
                       nr_cells * sizeof(int)) != 0)
        error("Failed to allocate list of cells node IDs for VELOCIraptor.");
    
    for(int i=0; i<nr_cells; i++) cell_node_ids[i] = s->cells_top[i].nodeID;    
James Willis's avatar
James Willis committed
176

James Willis's avatar
James Willis committed
177
    message("MPI rank %d sending %zu gparts to VELOCIraptor.", engine_rank, nr_gparts);
James Willis's avatar
James Willis committed
178
    
179
    /* Append base name with either the step number or time depending on what format is specified in the parameter file. */
James Willis's avatar
James Willis committed
180
    char outputFileName[PARSER_MAX_LINE_SIZE + 128];
181
    if(e->stf_output_freq_format == IO_STF_OUTPUT_FREQ_FORMAT_STEPS) {
James Willis's avatar
James Willis committed
182
        snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04i.VELOCIraptor", e->stfBaseName,
183
             e->step);
184
185
    }
    else if(e->stf_output_freq_format == IO_STF_OUTPUT_FREQ_FORMAT_TIME) {
James Willis's avatar
James Willis committed
186
        snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04e.VELOCIraptor", e->stfBaseName,
187
188
             e->time);
    }
189

190
    /* Allocate and populate an array of swift_vel_parts to be passed to VELOCIraptor. */
James Willis's avatar
James Willis committed
191
    struct swift_vel_part *swift_parts = NULL;
192
193
194
195
196
197
198
199

    if (posix_memalign((void **)&swift_parts, part_align,
                       nr_gparts * sizeof(struct swift_vel_part)) != 0)
        error("Failed to allocate array of particles for VELOCIraptor.");
    
    bzero(swift_parts, nr_gparts * sizeof(struct swift_vel_part));

    const float energy_scale = 1.0;
200
201
    const float a2 = e->cosmology->a * e->cosmology->a;
    
202
    message("Energy scaling factor: %f", energy_scale);
203
    message("a^2: %f", a2);
204
205
   
    /* Convert particle properties into VELOCIraptor units */ 
James Willis's avatar
James Willis committed
206
    for(size_t i=0; i<nr_gparts; i++) {
207
208
209
210
211
212
213
214
      swift_parts[i].x[0] = gparts[i].x[0];
      swift_parts[i].x[1] = gparts[i].x[1];
      swift_parts[i].x[2] = gparts[i].x[2];
      swift_parts[i].v[0] = gparts[i].v_full[0] / a2;
      swift_parts[i].v[1] = gparts[i].v_full[1] / a2;
      swift_parts[i].v[2] = gparts[i].v_full[2] / a2;
      swift_parts[i].mass = gparts[i].mass;
      swift_parts[i].potential = gparts[i].potential;
215
216
      swift_parts[i].type = gparts[i].type;
      
217
      /* Set gas particle IDs from their hydro counterparts and set internal energies. */
218
219
220
221
222
223
224
225
226
227
228
      if(gparts[i].type == swift_type_gas) {  
        swift_parts[i].id = parts[-gparts[i].id_or_neg_offset].id;
        swift_parts[i].u = hydro_get_physical_internal_energy(&parts[-gparts[i].id_or_neg_offset], e->cosmology) * energy_scale;
      }
      else {
        swift_parts[i].id = gparts[i].id_or_neg_offset;
        swift_parts[i].u = 0.f;
      }

    }

229
    /* Call VELOCIraptor. */
230
    if(!InvokeVelociraptor(nr_gparts, nr_hydro_parts, swift_parts, cell_node_ids, outputFileName)) error("Exiting. Call to VELOCIraptor failed on rank: %d.", e->nodeID);
231

232
233
234
    /* Reset the pthread affinity mask after VELOCIraptor returns. */
    pthread_setaffinity_np(thread, sizeof(cpu_set_t), engine_entry_affinity());
    
235
    /* Free cell node ids after VELOCIraptor has copied them. */
236
    free(cell_node_ids);
237
    free(swift_parts);
238
    
239
    message("VELOCIraptor took %.3f %s on rank %d.",
240
            clocks_from_ticks(getticks() - tic), clocks_getunit(), engine_rank);
241
242
#else
    error("SWIFT not configure to run with VELOCIraptor.");
243
#endif /* HAVE_VELOCIRAPTOR */
244
}