src/mdlib/domdec_setup.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  * This file is part of Gromacs        Copyright (c) 1991-2008
   5  * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation; either version 2
  10  * of the License, or (at your option) any later version.
  11  *
  12  * To help us fund GROMACS development, we humbly ask that you cite
  13  * the research papers on the package. Check out http://www.gromacs.org
  14  *
  15  * And Hey:
  16  * Gnomes, ROck Monsters And Chili Sauce
  17  */
  18
  19 #ifdef HAVE_CONFIG_H
  20 #include <config.h>
  21 #endif
  22
  23 #include <stdio.h>
  24 #include "domdec.h"
  25 #include "network.h"
  26 #include "perf_est.h"
  27 #include "physics.h"
  28 #include "smalloc.h"
  29 #include "typedefs.h"
  30 #include "vec.h"
  31 #include "names.h"
  32
  33 /* Margin for setting up the DD grid */
  34 #define DD_GRID_MARGIN_PRES_SCALE 1.05
  35
  36 static int factorize(int n,int **fac,int **mfac)
  37 {
  38         int d,ndiv;
  39
  40         /* Decompose n in factors */
  41         snew(*fac,n/2);
  42         snew(*mfac,n/2);
  43         d = 2;
  44         ndiv = 0;
  45         while (n > 1)
  46                 {
  47                 while (n % d == 0)
  48                         {
  49                         if (ndiv == 0 || (*fac)[ndiv-1] != d)
  50                                 {
  51                                 ndiv++;
  52                                 (*fac)[ndiv-1] = d;
  53                                 }
  54                         (*mfac)[ndiv-1]++;
  55                         n /= d;
  56                         }
  57                 d++;
  58                 }
  59
  60         return ndiv;
  61 }
  62
  63 static gmx_bool fits_pme_ratio(int nnodes,int npme,float ratio)
  64 {
  65     return ((double)npme/(double)nnodes > 0.95*ratio);
  66 }
  67
  68 static gmx_bool fits_pp_pme_perf(FILE *fplog,
  69                              t_inputrec *ir,matrix box,gmx_mtop_t *mtop,
  70                              int nnodes,int npme,float ratio)
  71 {
  72     int ndiv,*div,*mdiv,ldiv;
  73
  74     ndiv = factorize(nnodes-npme,&div,&mdiv);
  75     ldiv = div[ndiv-1];
  76     sfree(div);
  77     sfree(mdiv);
  78     /* The check below gives a reasonable division:
  79      * factor 5 allowed at 5 or more PP nodes,
  80      * factor 7 allowed at 49 or more PP nodes.
  81      */
  82     if (ldiv > 3 + (int)(pow(nnodes-npme,1.0/3.0) + 0.5))
  83     {
  84         return FALSE;
  85     }
  86
  87     /* Does this division gives a reasonable PME load? */
  88     return fits_pme_ratio(nnodes,npme,ratio);
  89 }
  90
  91 static int guess_npme(FILE *fplog,gmx_mtop_t *mtop,t_inputrec *ir,matrix box,
  92                                           int nnodes)
  93 {
  94         float ratio;
  95         int  npme,nkx,nky;
  96         t_inputrec ir_try;
  97
  98         ratio = pme_load_estimate(mtop,ir,box);
  99
 100         if (fplog)
 101     {
 102                 fprintf(fplog,"Guess for relative PME load: %.2f\n",ratio);
 103     }
 104
 105         /* We assume the optimal node ratio is close to the load ratio.
 106          * The communication load is neglected,
 107          * but (hopefully) this will balance out between PP and PME.
 108          */
 109
 110     if (!fits_pme_ratio(nnodes,nnodes/2,ratio))
 111     {
 112         /* We would need more than nnodes/2 PME only nodes,
 113          * which is not possible. Since the PME load is very high,
 114          * we will not loose much performance when all nodes do PME.
 115          */
 116
 117         return 0;
 118     }
 119
 120     /* First try to find npme as a factor of nnodes up to nnodes/3.
 121      * We start with a minimum PME node fraction of 1/16
 122      * and avoid ratios which lead to large prime factors in nnodes-npme.
 123      */
 124     npme = (nnodes + 15)/16;
 125     while (npme <= nnodes/3) {
 126         if (nnodes % npme == 0)
 127         {
 128             /* Note that fits_perf might change the PME grid,
 129              * in the current implementation it does not.
 130              */
 131             if (fits_pp_pme_perf(fplog,ir,box,mtop,nnodes,npme,ratio))
 132                         {
 133                                 break;
 134                         }
 135         }
 136         npme++;
 137     }
 138     if (npme > nnodes/3)
 139     {
 140         /* Try any possible number for npme */
 141         npme = 1;
 142         while (npme <= nnodes/2)
 143         {
 144             /* Note that fits_perf may change the PME grid */
 145             if (fits_pp_pme_perf(fplog,ir,box,mtop,nnodes,npme,ratio))
 146             {
 147                 break;
 148             }
 149             npme++;
 150         }
 151     }
 152     if (npme > nnodes/2)
 153     {
 154         gmx_fatal(FARGS,"Could not find an appropriate number of separate PME nodes. i.e. >= %5f*#nodes (%d) and <= #nodes/2 (%d) and reasonable performance wise (grid_x=%d, grid_y=%d).\n"
 155                   "Use the -npme option of mdrun or change the number of processors or the PME grid dimensions, see the manual for details.",
 156                   ratio,(int)(0.95*ratio*nnodes+0.5),nnodes/2,ir->nkx,ir->nky);
 157         /* Keep the compiler happy */
 158         npme = 0;
 159     }
 160     else
 161     {
 162         if (fplog)
 163         {
 164             fprintf(fplog,
 165                     "Will use %d particle-particle and %d PME only nodes\n"
 166                     "This is a guess, check the performance at the end of the log file\n",
 167                     nnodes-npme,npme);
 168         }
 169         fprintf(stderr,"\n"
 170                 "Will use %d particle-particle and %d PME only nodes\n"
 171                 "This is a guess, check the performance at the end of the log file\n",
 172                 nnodes-npme,npme);
 173     }
 174
 175     return npme;
 176 }
 177
 178 static int lcd(int n1,int n2)
 179 {
 180     int d,i;
 181
 182     d = 1;
 183     for(i=2; (i<=n1 && i<=n2); i++)
 184     {
 185         if (n1 % i == 0 && n2 % i == 0)
 186         {
 187             d = i;
 188         }
 189     }
 190
 191   return d;
 192 }
 193
 194 static int div_up(int n,int f)
 195 {
 196     return (n + f - 1)/f;
 197 }
 198
 199 real comm_box_frac(ivec dd_nc,real cutoff,gmx_ddbox_t *ddbox)
 200 {
 201     int  i,j,k,npp;
 202     rvec bt,nw;
 203     real comm_vol;
 204
 205     for(i=0; i<DIM; i++)
 206     {
 207         bt[i] = ddbox->box_size[i]*ddbox->skew_fac[i];
 208         nw[i] = dd_nc[i]*cutoff/bt[i];
 209     }
 210
 211     npp = 1;
 212     comm_vol = 0;
 213     for(i=0; i<DIM; i++)
 214     {
 215         if (dd_nc[i] > 1)
 216         {
 217             npp *= dd_nc[i];
 218             comm_vol += nw[i];
 219             for(j=i+1; j<DIM; j++)
 220             {
 221                 if (dd_nc[j] > 1)
 222                 {
 223                     comm_vol += nw[i]*nw[j]*M_PI/4;
 224                     for(k=j+1; k<DIM; k++)
 225                     {
 226                         if (dd_nc[k] > 1)
 227                         {
 228                             comm_vol += nw[i]*nw[j]*nw[k]*M_PI/6;
 229                         }
 230                     }
 231                 }
 232             }
 233         }
 234     }
 235
 236     return comm_vol;
 237 }
 238
 239 static gmx_bool inhomogeneous_z(const t_inputrec *ir)
 240 {
 241     return ((EEL_PME(ir->coulombtype) || ir->coulombtype==eelEWALD) &&
 242             ir->ePBC==epbcXYZ && ir->ewald_geometry==eewg3DC);
 243 }
 244
 245 static float comm_cost_est(gmx_domdec_t *dd,real limit,real cutoff,
 246                            matrix box,gmx_ddbox_t *ddbox,
 247                            int natoms,t_inputrec *ir,
 248                            float pbcdxr,
 249                            int npme_tot,ivec nc)
 250 {
 251     ivec npme={1,1,1};
 252     int  i,j,k,nk,overlap;
 253     rvec bt;
 254     float comm_vol,comm_vol_xf,comm_pme,cost_pbcdx;
 255     /* This is the cost of a pbc_dx call relative to the cost
 256      * of communicating the coordinate and force of an atom.
 257      * This will be machine dependent.
 258      * These factors are for x86 with SMP or Infiniband.
 259      */
 260     float pbcdx_rect_fac = 0.1;
 261     float pbcdx_tric_fac = 0.2;
 262
 263     /* Check the DD algorithm restrictions */
 264     if ((ir->ePBC == epbcXY && ir->nwall < 2 && nc[ZZ] > 1) ||
 265         (ir->ePBC == epbcSCREW && (nc[XX] == 1 || nc[YY] > 1 || nc[ZZ] > 1)))
 266     {
 267         return -1;
 268     }
 269
 270     if (inhomogeneous_z(ir) && nc[ZZ] > 1)
 271     {
 272         return -1;
 273     }
 274
 275     /* Check if the triclinic requirements are met */
 276     for(i=0; i<DIM; i++)
 277     {
 278         for(j=i+1; j<ddbox->npbcdim; j++)
 279         {
 280             if (box[j][i] != 0 || ir->deform[j][i] != 0 ||
 281                 (ir->epc != epcNO && ir->compress[j][i] != 0))
 282             {
 283                 if (nc[j] > 1 && nc[i] == 1)
 284                 {
 285                     return -1;
 286                 }
 287             }
 288         }
 289     }
 290
 291     for(i=0; i<DIM; i++)
 292     {
 293         bt[i] = ddbox->box_size[i]*ddbox->skew_fac[i];
 294
 295         /* Without PBC there are no cell size limits with 2 cells */
 296         if (!(i >= ddbox->npbcdim && nc[i] <= 2) && bt[i] < nc[i]*limit)
 297         {
 298             return -1;
 299         }
 300     }
 301
 302     if (npme_tot > 1)
 303     {
 304         /* The following choices should match those
 305          * in init_domain_decomposition in domdec.c.
 306          */
 307         if (nc[XX] == 1 && nc[YY] > 1)
 308         {
 309             npme[XX] = 1;
 310             npme[YY] = npme_tot;
 311         }
 312         else if (nc[YY] == 1)
 313         {
 314             npme[XX] = npme_tot;
 315             npme[YY] = 1;
 316         }
 317         else
 318         {
 319             /* Will we use 1D or 2D PME decomposition? */
 320             npme[XX] = (npme_tot % nc[XX] == 0) ? nc[XX] : npme_tot;
 321             npme[YY] = npme_tot/npme[XX];
 322         }
 323     }
 324
 325     /* When two dimensions are (nearly) equal, use more cells
 326      * for the smallest index, so the decomposition does not
 327      * depend sensitively on the rounding of the box elements.
 328      */
 329     for(i=0; i<DIM; i++)
 330     {
 331         for(j=i+1; j<DIM; j++)
 332         {
 333             /* Check if the box size is nearly identical,
 334              * in that case we prefer nx > ny  and ny > nz.
 335              */
 336             if (fabs(bt[j] - bt[i]) < 0.01*bt[i] && nc[j] > nc[i])
 337             {
 338                 /* The XX/YY check is a bit compact. If nc[YY]==npme[YY]
 339              * this means the swapped nc has nc[XX]==npme[XX],
 340              * and we can also swap X and Y for PME.
 341              */
 342                 /* Check if dimension i and j are equivalent for PME.
 343                  * For x/y: if nc[YY]!=npme[YY], we can not swap x/y
 344                  * For y/z: we can not have PME decomposition in z
 345                  */
 346                 if (npme_tot <= 1 ||
 347                     !((i == XX && j == YY && nc[YY] != npme[YY]) ||
 348                       (i == YY && j == ZZ && npme[YY] > 1)))
 349                 {
 350                     return -1;
 351                 }
 352             }
 353         }
 354     }
 355
 356     /* This function determines only half of the communication cost.
 357      * All PP, PME and PP-PME communication is symmetric
 358      * and the "back"-communication cost is identical to the forward cost.
 359      */
 360
 361     comm_vol = comm_box_frac(nc,cutoff,ddbox);
 362
 363     comm_pme = 0;
 364     for(i=0; i<2; i++)
 365     {
 366         /* Determine the largest volume for PME x/f redistribution */
 367         if (nc[i] % npme[i] != 0)
 368         {
 369             if (nc[i] > npme[i])
 370             {
 371                 comm_vol_xf = (npme[i]==2 ? 1.0/3.0 : 0.5);
 372             }
 373             else
 374             {
 375                 comm_vol_xf = 1.0 - lcd(nc[i],npme[i])/(double)npme[i];
 376             }
 377             comm_pme += 3*natoms*comm_vol_xf;
 378         }
 379
 380         /* Grid overlap communication */
 381         if (npme[i] > 1)
 382         {
 383             nk = (i==0 ? ir->nkx : ir->nky);
 384             overlap = (nk % npme[i] == 0 ? ir->pme_order-1 : ir->pme_order);
 385             comm_pme += npme[i]*overlap*ir->nkx*ir->nky*ir->nkz/nk;
 386         }
 387     }
 388
 389     /* PME FFT communication volume.
 390      * This only takes the communication into account and not imbalance
 391      * in the calculation. But the imbalance in communication and calculation
 392      * are similar and therefore these formulas also prefer load balance
 393      * in the FFT and pme_solve calculation.
 394      */
 395     comm_pme += (npme[YY] - 1)*npme[YY]*div_up(ir->nky,npme[YY])*div_up(ir->nkz,npme[YY])*ir->nkx;
 396     comm_pme += (npme[XX] - 1)*npme[XX]*div_up(ir->nkx,npme[XX])*div_up(ir->nky,npme[XX])*ir->nkz;
 397
 398     /* Add cost of pbc_dx for bondeds */
 399     cost_pbcdx = 0;
 400     if ((nc[XX] == 1 || nc[YY] == 1) || (nc[ZZ] == 1 && ir->ePBC != epbcXY))
 401     {
 402         if ((ddbox->tric_dir[XX] && nc[XX] == 1) ||
 403             (ddbox->tric_dir[YY] && nc[YY] == 1))
 404         {
 405             cost_pbcdx = pbcdxr*pbcdx_tric_fac;
 406         }
 407         else
 408         {
 409             cost_pbcdx = pbcdxr*pbcdx_rect_fac;
 410         }
 411     }
 412
 413     if (debug)
 414     {
 415         fprintf(debug,
 416                 "nc %2d %2d %2d %2d %2d vol pp %6.4f pbcdx %6.4f pme %9.3e tot %9.3e\n",
 417                 nc[XX],nc[YY],nc[ZZ],npme[XX],npme[YY],
 418                 comm_vol,cost_pbcdx,comm_pme,
 419                 3*natoms*(comm_vol + cost_pbcdx) + comm_pme);
 420     }
 421
 422     return 3*natoms*(comm_vol + cost_pbcdx) + comm_pme;
 423 }
 424
 425 static void assign_factors(gmx_domdec_t *dd,
 426                            real limit,real cutoff,
 427                            matrix box,gmx_ddbox_t *ddbox,
 428                            int natoms,t_inputrec *ir,
 429                            float pbcdxr,int npme,
 430                            int ndiv,int *div,int *mdiv,ivec ir_try,ivec opt)
 431 {
 432     int x,y,z,i;
 433     float ce;
 434
 435     if (ndiv == 0)
 436     {
 437         ce = comm_cost_est(dd,limit,cutoff,box,ddbox,
 438                            natoms,ir,pbcdxr,npme,ir_try);
 439         if (ce >= 0 && (opt[XX] == 0 ||
 440                         ce < comm_cost_est(dd,limit,cutoff,box,ddbox,
 441                                            natoms,ir,pbcdxr,
 442                                            npme,opt)))
 443         {
 444             copy_ivec(ir_try,opt);
 445         }
 446
 447         return;
 448     }
 449
 450     for(x=mdiv[0]; x>=0; x--)
 451     {
 452         for(i=0; i<x; i++)
 453         {
 454             ir_try[XX] *= div[0];
 455         }
 456         for(y=mdiv[0]-x; y>=0; y--)
 457         {
 458             for(i=0; i<y; i++)
 459             {
 460                 ir_try[YY] *= div[0];
 461             }
 462             for(i=0; i<mdiv[0]-x-y; i++)
 463             {
 464                 ir_try[ZZ] *= div[0];
 465             }
 466
 467             /* recurse */
 468             assign_factors(dd,limit,cutoff,box,ddbox,natoms,ir,pbcdxr,npme,
 469                            ndiv-1,div+1,mdiv+1,ir_try,opt);
 470
 471             for(i=0; i<mdiv[0]-x-y; i++)
 472             {
 473                 ir_try[ZZ] /= div[0];
 474             }
 475             for(i=0; i<y; i++)
 476             {
 477                 ir_try[YY] /= div[0];
 478             }
 479         }
 480         for(i=0; i<x; i++)
 481         {
 482             ir_try[XX] /= div[0];
 483         }
 484     }
 485 }
 486
 487 static real optimize_ncells(FILE *fplog,
 488                             int nnodes_tot,int npme_only,
 489                             gmx_bool bDynLoadBal,real dlb_scale,
 490                             gmx_mtop_t *mtop,matrix box,gmx_ddbox_t *ddbox,
 491                             t_inputrec *ir,
 492                             gmx_domdec_t *dd,
 493                             real cellsize_limit,real cutoff,
 494                             gmx_bool bInterCGBondeds,gmx_bool bInterCGMultiBody,
 495                             ivec nc)
 496 {
 497     int npp,npme,ndiv,*div,*mdiv,d,nmax;
 498     gmx_bool bExcl_pbcdx;
 499     float pbcdxr;
 500     real limit;
 501     ivec itry;
 502
 503     limit  = cellsize_limit;
 504
 505     dd->nc[XX] = 1;
 506     dd->nc[YY] = 1;
 507     dd->nc[ZZ] = 1;
 508
 509     npp = nnodes_tot - npme_only;
 510     if (EEL_PME(ir->coulombtype))
 511     {
 512         npme = (npme_only > 0 ? npme_only : npp);
 513     }
 514     else
 515     {
 516         npme = 0;
 517     }
 518
 519     if (bInterCGBondeds)
 520     {
 521         /* For Ewald exclusions pbc_dx is not called */
 522         bExcl_pbcdx =
 523             (IR_EXCL_FORCES(*ir) && !EEL_FULL(ir->coulombtype));
 524         pbcdxr = (double)n_bonded_dx(mtop,bExcl_pbcdx)/(double)mtop->natoms;
 525     }
 526     else
 527     {
 528         /* Every molecule is a single charge group: no pbc required */
 529         pbcdxr = 0;
 530     }
 531     /* Add a margin for DLB and/or pressure scaling */
 532     if (bDynLoadBal)
 533     {
 534         if (dlb_scale >= 1.0)
 535         {
 536             gmx_fatal(FARGS,"The value for option -dds should be smaller than 1");
 537         }
 538         if (fplog)
 539         {
 540             fprintf(fplog,"Scaling the initial minimum size with 1/%g (option -dds) = %g\n",dlb_scale,1/dlb_scale);
 541         }
 542         limit /= dlb_scale;
 543     }
 544     else if (ir->epc != epcNO)
 545     {
 546         if (fplog)
 547         {
 548             fprintf(fplog,"To account for pressure scaling, scaling the initial minimum size with %g\n",DD_GRID_MARGIN_PRES_SCALE);
 549             limit *= DD_GRID_MARGIN_PRES_SCALE;
 550         }
 551     }
 552
 553     if (fplog)
 554     {
 555         fprintf(fplog,"Optimizing the DD grid for %d cells with a minimum initial size of %.3f nm\n",npp,limit);
 556
 557         if (inhomogeneous_z(ir))
 558         {
 559             fprintf(fplog,"Ewald_geometry=%s: assuming inhomogeneous particle distribution in z, will not decompose in z.\n",eewg_names[ir->ewald_geometry]);
 560         }
 561
 562         if (limit > 0)
 563         {
 564             fprintf(fplog,"The maximum allowed number of cells is:");
 565             for(d=0; d<DIM; d++)
 566             {
 567                 nmax = (int)(ddbox->box_size[d]*ddbox->skew_fac[d]/limit);
 568                 if (d >= ddbox->npbcdim && nmax < 2)
 569                 {
 570                     nmax = 2;
 571                 }
 572                 if (d == ZZ && inhomogeneous_z(ir))
 573                 {
 574                     nmax = 1;
 575                 }
 576                 fprintf(fplog," %c %d",'X' + d,nmax);
 577             }
 578             fprintf(fplog,"\n");
 579         }
 580     }
 581
 582     if (debug)
 583     {
 584         fprintf(debug,"Average nr of pbc_dx calls per atom %.2f\n",pbcdxr);
 585     }
 586
 587     /* Decompose npp in factors */
 588     ndiv = factorize(npp,&div,&mdiv);
 589
 590     itry[XX] = 1;
 591     itry[YY] = 1;
 592     itry[ZZ] = 1;
 593     clear_ivec(nc);
 594     assign_factors(dd,limit,cutoff,box,ddbox,mtop->natoms,ir,pbcdxr,
 595                    npme,ndiv,div,mdiv,itry,nc);
 596
 597     sfree(div);
 598     sfree(mdiv);
 599
 600     return limit;
 601 }
 602
 603 real dd_choose_grid(FILE *fplog,
 604                     t_commrec *cr,gmx_domdec_t *dd,t_inputrec *ir,
 605                     gmx_mtop_t *mtop,matrix box,gmx_ddbox_t *ddbox,
 606                     gmx_bool bDynLoadBal,real dlb_scale,
 607                     real cellsize_limit,real cutoff_dd,
 608                     gmx_bool bInterCGBondeds,gmx_bool bInterCGMultiBody)
 609 {
 610     int  npme,nkx,nky;
 611     real limit;
 612
 613     if (MASTER(cr))
 614     {
 615         if (EEL_PME(ir->coulombtype))
 616         {
 617             if (cr->npmenodes >= 0)
 618             {
 619                 if (cr->nnodes <= 2 && cr->npmenodes > 0)
 620                 {
 621                     gmx_fatal(FARGS,
 622                               "Can not have separate PME nodes with 2 or less nodes");
 623                 }
 624             }
 625             else
 626             {
 627                 if (cr->nnodes <= 10)
 628                 {
 629                     cr->npmenodes = 0;
 630                 }
 631                 else
 632                 {
 633                     cr->npmenodes = guess_npme(fplog,mtop,ir,box,cr->nnodes);
 634                 }
 635             }
 636             if (fplog)
 637             {
 638                 fprintf(fplog,"Using %d separate PME nodes\n",cr->npmenodes);
 639             }
 640         }
 641         else
 642         {
 643             if (cr->npmenodes < 0)
 644             {
 645                 cr->npmenodes = 0;
 646             }
 647         }
 648
 649         limit = optimize_ncells(fplog,cr->nnodes,cr->npmenodes,
 650                                 bDynLoadBal,dlb_scale,
 651                                 mtop,box,ddbox,ir,dd,
 652                                 cellsize_limit,cutoff_dd,
 653                                 bInterCGBondeds,bInterCGMultiBody,
 654                                 dd->nc);
 655     }
 656     else
 657     {
 658         limit = 0;
 659     }
 660     /* Communicate the information set by the master to all nodes */
 661     gmx_bcast(sizeof(dd->nc),dd->nc,cr);
 662     if (EEL_PME(ir->coulombtype))
 663     {
 664         gmx_bcast(sizeof(ir->nkx),&ir->nkx,cr);
 665         gmx_bcast(sizeof(ir->nky),&ir->nky,cr);
 666         gmx_bcast(sizeof(cr->npmenodes),&cr->npmenodes,cr);
 667     }
 668     else
 669     {
 670         cr->npmenodes = 0;
 671     }
 672
 673     return limit;
 674 }