pmcimg.c

蒙特卡罗模拟光子成像C语言版,代码简洁专业

  /* MAKE SURE THE SOURCES ARE AT INTERFACES */

  if (syscfg.flags.source_move)
    move_sources(&syscfg, (const TISSUE ***)tissueType);

  /* MAKE SURE THE DETECTORS ARE AT INTERFACES */

  if (syscfg.flags.detector_move)
    move_detectors(&syscfg, (const TISSUE ***)tissueType);

  /* RUN THE MONTE CARLO SIMULATION */

  run_simulation(&syscfg, (const TISSUE ***)tissueType, 
				II, JJ, II0, JJ0, argv[1]);

  /* DONE */

  endParallel();

  nvox = syscfg.nIxstep * syscfg.nIystep * syscfg.nIzstep * syscfg.nTstep;

  free_1d_array(II,  nvox);
  free_1d_array(JJ,  nvox);
  free_1d_array(II0, nvox);
  free_1d_array(JJ0, nvox);

  return 0;
}

/* ------------------------ RUN THE SIMULATION ---------------------- */

static void run_simulation(const struct Config *syscfg, 
                           const TISSUE ***tissueType,
                           TWOPT_T *II, TWOPT_T *JJ, 
                           TWOPT_T *II0, TWOPT_T *JJ0, const char *basefile)
{
  FILE *fp = NULL;		  /* FILE POINTERS FOR SAVING THE DATA */
  int i, ns, nvox, seed;
  unsigned long N;
#if (REALLY_RANDOM)
  FILE *rp;
#endif
  unsigned long NPhotons;
  int n;

#if (REALLY_RANDOM)
  /* /dev/random is a much better seed for parallel code than seed+n */

  if ((rp = fopen("/dev/random", "rb")) == NULL)
    {
      if ((rp = fopen("/dev/urandom", "rb")) == NULL)
	{
	  fread(&seed, sizeof(seed), 1, rp);
	  fclose(rp);
	}
      else
	{
	  printf("Warning: can't open /dev/random, using less random seed\n");
	  seed += myNode();
	}
    }
  else
    {
      fread(&seed, sizeof(seed), 1, rp);
      fclose(rp);
    }
#else
  seed = syscfg->seed;
#endif

  nvox = syscfg->nIystep * syscfg->nIxstep * syscfg->nIzstep * syscfg->nTstep;

  if (syscfg->nSrcs > 1)
    {
      if (syscfg->nDets == 0)
	printf("No detectors: writing all exiting photons to %s.his.aN\n",
	       basefile);
      else if (syscfg->flags.exact_exit_time)
	printf("Writing exact exit times to %s.his.tN\n", basefile);
      else
	printf("Writing exiting tgate number to %s.his.N\n", basefile);
    }
  else
    {
      if (syscfg->nDets == 0)
	printf("No detectors: writing all exiting photons to %s.his.a\n",
	       basefile);
      else if (syscfg->flags.exact_exit_time)
	printf("Writing exact exit times to %s.his.t\n", basefile);
      else
	printf("Writing exiting tgate number to %s.his\n", basefile);
    }

  /* The fraction of the total that I'm responsible for.  The number of
   * nodes is so vastly smaller than the number of photons, we won't
   * worry about rounding problems.  Subtract 1 since the master
   * thread doesn't itself run any photons */

  if (getNodeCount() > 1)
    NPhotons = syscfg->NT / (getNodeCount() - 1);
  else
    NPhotons = syscfg->NT;

#if (! REALLY_RANDOM)
  /* Diddle my random number seed, or there's no point to running
   * in parallel! */

  seed += myNode();
#endif

  /* FOR EVERY SOURCE, RUN THE SIMULATION */

  for (ns = 0; ns < syscfg->nSrcs; ns++)
    {
      /*******************************************************************
       * START MIGRATING THE PHOTON GENERATING PHOTONS UNTIL NUMBER OF
       * PHOTONS EXECUTED (N) IS EQUAL TO THE NUMBER TO BE GENERATED (NT)
       *******************************************************************/

      MPI_Barrier(MPI_COMM_WORLD);       /* Start all the nodes together */

      if (isMaster())
	{
          int nleft = 0;

	  if ((fp = openhis(ns, basefile, syscfg)) == NULL)
	    {
	      fprintf(stderr, "Can't open history file %s\n", basefile);
	      MPI_Abort(MPI_COMM_WORLD, 1);
	    }

	  /* As photons leave the system at the different nodes, write
	   * them to the history file.  Only the master copy should
	   * write to the history file to preserve data integrity. */

	  /* nleft does not include myself */
	  nleft = getNodeCount() - 1;

	  while (nleft > 0)
	    {
	      MPI_Status status;
	      float buffer[MAX_HISTENT];

	      /* Gather in photons as they leave the system
	       * and write history to disk (blocking call) */

	      MPI_Recv(buffer, MAX_HISTENT, MPI_FLOAT, MPI_ANY_SOURCE,
		       TAG_HISFILE, MPI_COMM_WORLD, &status);

	      MPI_Get_count(&status, MPI_FLOAT, &n);

	      /* We use empty records to signify that a node
	       *  has finished it's allocation of photons */

	      if (n == 0)
		{
#if (DEBUG)
		  printf("Node %d finished\n", status.MPI_SOURCE);
#endif
		  nleft -= 1;
		}
	      else
		{
		  int nw;

#if (PARANOID)
                  if (syscfg->nDets == 0)
                    assert(n == 4+syscfg->Ntissue);
                  else
                    assert(n == 2+syscfg->Ntissue);
#endif

		  if ((nw = fwrite(buffer, sizeof(float), n, fp)) != n)
		    {
		      perror("Error writing to history file");
		      MPI_Abort(MPI_COMM_WORLD, 1);
		    }
		}
	    }

          /* CLOSE HISTORY FILE */
          fclose(fp);
	}
      else			  /* not master process */
	{
	  float foo;		  /* Dummy memory for MPI_Send */

          for (N = 0; N < NPhotons; N++)
            run_photon(syscfg, ns, tissueType, II, JJ, fp, &seed);

	  /* Signal that I have no more photons to deliver */
	  MPI_Send(&foo, 0, MPI_FLOAT, 0, TAG_HISFILE, MPI_COMM_WORLD);

          /* Send any buffered messages and deallocate the buffer */

	  if (mpibuff != NULL)
	    {
              extern int mpiblen;
              void *p;
              int   s;

	      /* Detach buffer, blocks until all messages have been sent */
	      MPI_Buffer_detach(&p, &s);

	      free_1d_array(mpibuff, 100 * mpiblen);

	      mpibuff = NULL;
              mpiblen = 0;
	    }
	}

      /* SAVE THE PHOTON DENSITY (2-POINT GREENS FUNCTION) */

      if (syscfg->flags.gen_twopt)
	{
	  /* Sum all the nodes into the master's copy */

	  if (II != NULL)
	    {
#if (PARANOID)
	      assert(II  != 0);
	      assert(II0 != 0);
#endif
	      MPI_Reduce(II, II0, 
                         nvox, MPI_TWOPT, MPI_SUM, 0, MPI_COMM_WORLD);
	    }

	  if (JJ != NULL)
	    {
#if (PARANOID)
	      assert(JJ  != 0);
	      assert(JJ0 != 0);
#endif
	      MPI_Reduce(JJ, JJ0, 
                         nvox, MPI_TWOPT, MPI_SUM, 0, MPI_COMM_WORLD);
	    }

	  if (II != NULL && II0 != NULL)
	    memcpy(II, II0, nvox * sizeof(TWOPT_T));

	  if (JJ != NULL && JJ0 != NULL)
	    memcpy(JJ, JJ0, nvox * sizeof(TWOPT_T));

	  /* Save the final (summed) version of the 2pt file */

	  if (isMaster())
	    {
	      if (syscfg->nSrcs == 1)
		save2pt(basefile, -1, II, JJ, nvox, NULL);
	      else
		save2pt(basefile, ns, II, JJ, nvox, NULL);
	    }

	  /* Reset memory for next pass through */

          if (II != NULL)
            {
              for (i = 0; i < nvox; i++)
                II[i] = 0.0;

              if (II0 != NULL)
                memcpy(II0, II, nvox * sizeof(TWOPT_T));
            }

          if (JJ != NULL)
            {
              for (i = 0; i < nvox; i++)
                JJ[i] = 0.0;

              if (JJ0 != NULL)
                memcpy(JJ0, JJ, nvox * sizeof(TWOPT_T));
            }
	}
    }

  return;
}

/* Call the parallel savehis() routine */

void write_photon(const struct Config *syscfg, FILE *fp,
		  double x, double y, double z,
		  double kx, double ky, double kz, 
		  double Ttot, float *lenTiss)
{
#if (PARANOID)
  assert(!isMaster());
#endif

  parallel_savehis(syscfg, x, y, z, kx, ky, kz, Ttot, lenTiss);

  return;
}