xmdssimulation.cc

XMDS is a code generator that integrates equations. You write them down in human readable form in a

  if(verbose()) {
    printf("Writing simulation includes ...\n");
  }

  fprintf(outfile,
          "// ********************************************************\n"
  "// simulation includes\n"
  "\n"
  "#include <stdio.h>\n"
  "#include <math.h>\n"
  "#include <string>\n"
  "#include <fstream>\n"
  "#include <iostream>\n"
  "#include <sstream>\n"
  "#include <cstdlib>\n"
          "#include <ctime>\n");
  if(simulation()->argStruct()->nameList.size() != 0) {
    fprintf(outfile,"#include <getopt_xmds.h>\n");
  }
  fprintf(outfile,"#include <vector>\n");
  if(myParameters.usempi) {
    fprintf(outfile,"#include <mpi.h>\n");
  }
  if(myParameters.nThreads>1) {
    fprintf(outfile,"#include <fftw_threads.h>\n");
  }
  else if((!myParameters.stochastic)&(myParameters.usempi)) {
      fprintf(outfile,
              "#include <fftw_mpi.h>\n"
              "#include <fftw.h>\n");
    }
  else
    fprintf(outfile,"#include <fftw.h>\n");
  fprintf(outfile,"#include <xmdscomplex.h>\n");
  if (myParameters.binaryOutput) {
    fprintf(outfile,"#include <xmdsconfig.h>\n");
  }
  fprintf(outfile,
          "using namespace std;\n"
          "\n");
};

// ******************************************************************************
void xmdsSimulation::writeDefines(
                                  FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsSimulation::writeDefines\n");
  }

  if(verbose()) {
    printf("Writing simulation defines ...\n");
  }

  fprintf(outfile,
          "// ********************************************************\n"
          "// simulation defines\n"
          "\n");
  if(myParameters.nThreads>1) {
    fprintf(outfile,"#define _num_threads %li\n\n",myParameters.nThreads);
  }

  if(myParameters.stochastic) {
    fprintf(outfile,"#define _n_noises %li\n",myParameters.nNoises);
    for(long i=0;i<myParameters.nNoises;i++) {
      fprintf(outfile,"#define n_%li _noises[%li]\n",i+1,i);
    }
    if(myParameters.nPaths>1) {
      fprintf(outfile,"#define _n_paths %li\n",myParameters.nPaths);
    }
  }
  fprintf(outfile,"#define d%s _step\n",myParameters.propDimName.c_str());
  fprintf(outfile,"\n");

  xmdsElement::writeDefines(outfile);
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsSimulation::writeGlobals(
                                  FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsSimulation::writeGlobals\n");
  }

  if(verbose()) {
    printf("Writing simulation globals ...\n");
  }

  fprintf(outfile,
          "// ********************************************************\n"
          "// simulation globals\n"
          "\n");
  fprintf(outfile,"double %s;\n",myParameters.propDimName.c_str());
  fprintf(outfile,"\n");
  if(myParameters.errorCheck) {
    fprintf(outfile,"bool _half_step;\n\n");
  }
  if(myParameters.usempi){
    fprintf(outfile,"int rank, size;\n\n");
  }

  if(myParameters.stochastic) {
    if (myParameters.noiseKind == "poissonian") {
      fprintf(outfile, 
	      "double _noiseMean = %g;\n",myParameters.noiseMean);
    }
    if(myParameters.errorCheck) {
      fprintf(outfile,
              "unsigned short _gen1[3];\n"
              "unsigned short _gen2[3];\n");
    }
    else
      fprintf(outfile,"unsigned short _gen[3];\n");
    fprintf(outfile,"\n");
  }
  else if(myParameters.usempi){
    fprintf(outfile,
            "int local_nx, local_x_start, local_ny_after_transpose, local_y_start_after_transpose, total_local_size;\n\n");
  }

  xmdsElement::writeGlobals(outfile);
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsSimulation::writePrototypes(
                                     FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsSimulation::writePrototypes\n");
  }

  if(verbose()) {
    printf("Writing simulation prototypes ...\n");
  }

  if(myParameters.stochastic) {
    if (debugFlag) {
	printf("xmdsSimulation::writePrototypes noiseKind = %s\n",myParameters.noiseKind.c_str()); 
    }
    if (myParameters.noiseKind == "gaussian") {
      fprintf(outfile,
	      "// ********************************************************\n"
	      "// simulation prototypes\n"
	      "\n"
	      "void _make_noises(\n"
	      "   unsigned short *const _generator,\n"
	      "   const double& _var,\n"
	      "   double *const _noise_vector,\n"
	      "   const unsigned long& _n);\n"
	      "\n");
    }
    else if (myParameters.noiseKind == "gaussFast") {
      fprintf(outfile,
	      "// ********************************************************\n"
	      "// simulation prototypes\n"
	      "\n"
	      "void _make_noises(\n"
	      "   unsigned short *const _generator,\n"
	      "   const double& _var,\n"
	      "   double *const _noise_vector,\n"
	      "   const unsigned long& _n);\n"
	      "\n");
    }
    else if (myParameters.noiseKind == "uniform") {
      fprintf(outfile,
	      "// ********************************************************\n"
	      "// simulation prototypes\n"
	      "\n"
	      "void _make_noises(\n"
	      "   unsigned short *const _generator,\n"
	      "   const double& _var,\n"
	      "   double *const _noise_vector,\n"
	      "   const unsigned long& _n);\n"
	      "\n");
    }
    else if (myParameters.noiseKind == "poissonian") {
      fprintf(outfile,
	      "// ********************************************************\n"
	      "// simulation prototypes\n"
	      "\n"
	      "void _make_noises(\n"
	      "   unsigned short *const _generator,\n"
	      "   const double& _var,\n"
	      "   double *const _noise_vector,\n"
	      "   const unsigned long& _n);\n"
	      "\n"
	      "double gammln(double _xx);\n"
	      "\n"
	      "double poidev(double xm, unsigned short *const _generator);\n"
	      "\n");
    }
    else {
      sprintf(errorMessage(),"Incorrectly specified noise kind (%s) in xmdsSimulation::writePrototypes",myParameters.noiseKind.c_str());
      throw xmdsException(errorMessage());
    }
  }

  xmdsElement::writePrototypes(outfile);
};

// ******************************************************************************
void xmdsSimulation::writeRoutines(
                                   FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsSimulation::writeRoutines\n");
  }

  if(verbose()) {
    printf("Writing simulation routines ...\n");
  }

  // NOTE:  it turns out the erand48() is more likely to be thread safe, and generate
  // decent random numbers.  rand() isn't necessarily thread safe, and the thread safe
  // version (rand_r()) isn't as good a random number generator.
  if(myParameters.stochastic) {
    if (debugFlag) {
	printf("xmdsSimulation::writeRoutines noiseKind = %s\n",myParameters.noiseKind.c_str());
    }
    if (myParameters.noiseKind == "gaussian") {
      fprintf(outfile,
              "// ********************************************************\n"
              "//               simulation routines\n"
              "// ********************************************************\n"
              "\n"
              "// *************************\n"
              "void _make_noises(\n"
              " unsigned short *const _generator,\n"
              " const double& _var,\n"
              " double *const _noise_vector,\n"
              " const unsigned long& _n) {\n"
              "\n"
              "for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
              " const double _mag = sqrt(-2*log(erand48(_generator))*_var);\n"
              " const double _arg = 2*M_PI*erand48(_generator);\n"
              " _noise_vector[_i0 + 0] = _mag*cos(_arg);\n"
              " _noise_vector[_i0 + 1] = _mag*sin(_arg);\n"
              " }\n"
              "}\n"
              "\n");
    }
    else if (myParameters.noiseKind == "gaussFast") {
      fprintf(outfile,
              "// ********************************************************\n"
              "//               simulation routines\n"
              "// ********************************************************\n"
              "\n"
              "// *************************\n"
              "void _make_noises(\n"
              " unsigned short *const _generator,\n"
              " const double& _var,\n"
              " double *const _noise_vector,\n"
              " const unsigned long& _n) {\n"
              "\n"
              "for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
              "  double _v1, _v2, _rsq;\n"
              "  do {\n"
              "    _v1 = 2.0*erand48(_generator)-1.0;\n"
              "    _v2 = 2.0*erand48(_generator)-1.0;\n"
              "    _rsq = _v1*_v1 + _v2*_v2;\n"
              "  } while(_rsq >= 1.0 || _rsq == 0.0);\n"
              "  const double _fac = sqrt(-2.0*_var*log(_rsq)/_rsq);\n"
              "  _noise_vector[_i0 + 0] = _v1*_fac;\n"
              "  _noise_vector[_i0 + 1] = _v2*_fac;\n"
              " }\n"
              "}\n"
              "\n");
    }
    else if (myParameters.noiseKind == "uniform") {
      fprintf(outfile,
              "// ********************************************************\n"
              "//               simulation routines\n"
              "// ********************************************************\n"
              "\n"
              "// *************************\n"
              "void _make_noises(\n"
              " unsigned short *const _generator,\n"
              " const double& _var,\n"
              " double *const _noise_vector,\n"
              " const unsigned long& _n) {\n"
              "\n"
              "for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
              "  const double _fac = sqrt(_var);\n"
              "  _noise_vector[_i0 + 0] = _fac*erand48(_generator);\n"
              "  _noise_vector[_i0 + 1] = _fac*erand48(_generator);\n"
              " }\n"
              "}\n"
              "\n");
    }
    else if (myParameters.noiseKind == "poissonian") {
      fprintf(outfile,
              "// ********************************************************\n"
              "//               simulation routines\n"
              "// ********************************************************\n"
              "\n"
              "// *************************\n"
              "double gammln(double _xx) {\n"
              "  double _x, _y, _tmp, _ser;\n"
              "  const double _cof[6] = {76.18009172947146,-86.50532032941677,\n"
              "                          24.01409824083091,-1.231739572450155,\n"
              "                          0.1208650973866179e-2,-0.5395239384953e-5};\n"
              "  int _j;\n"
              "  _y = _x = _xx;\n"
              "  _tmp = _x + 5.5;\n"
              "  _tmp -= (_x+0.5)*log(_tmp);\n"
              "  _ser = 1.000000000190015;\n"
              "  for (_j=0; _j<=5; _j++) _ser += _cof[_j]/++_y;\n"
              "  return -_tmp+log(2.5066282746310005*_ser/_x);\n"
              "}\n"
              "\n"
              "double poidev(double xm, unsigned short *const _generator) {\n"
              "  static double sq, alxm, g;\n"
              "  double em,t,y;\n"
              "  if (xm < 12.0) {        // Use direct method \n"
              "    g=exp(-xm);\n"
              "    em = -1;\n"
              "    t=1.0;\n"
              "    // Instead of adding exponential deviates it is equivalent\n"
              "    // to multiply uniform deviates.  We never actually have to\n"
              "    // take the log, merely compare to the pre-computed exponential\n"
              "    do {\n"
              "      ++em;\n"
              "      t *= erand48(_generator);\n"
              "    } while (t > g);\n"
              "  } else {                  // Use rejection method\n"
              "    sq=sqrt(2.0*xm);\n"
              "    alxm=log(xm);\n"
              "    g=xm*alxm-gammln(xm+1.0);\n"
              "    do {\n"
              "      do { // y is a deviate from a Lorenzian comparison function\n"
              "        y=tan(M_PI*erand48(_generator));\n"
              "        em=sq*y+xm;        // em is y, shifted and scaled\n"
              "      } while (em < 0.0);  // Reject if in regime of zero probability\n"
              "      em=floor(em);        // The trick for integer-valued distributions\n"
              "      t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);\n"
              "      // The ratio of the desired distribution to the comparison\n"
              "      // function; we reject by comparing it to another uniform\n"
              "      // deviate. The factor 0.9 so that t never exceeds 1.\n"
              "    } while (erand48(_generator) > t);\n"
              "  }\n"
              "  return em;\n"
              "  }\n"
              "\n"
              "void _make_noises(\n"
              " unsigned short *const _generator,\n"
              " const double& _var,\n"
              " double *const _noise_vector,\n"
              " const unsigned long& _n) {\n"
              "\n"
              "for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
              "  const double _fac = sqrt(_var);\n"
              "  _noise_vector[_i0 + 0] = _fac*poidev(_noiseMean, _generator);\n"
              "  _noise_vector[_i0 + 1] = _fac*poidev(_noiseMean, _generator);\n"
              " }\n"
              "}\n"
              "\n");
    }
    else {
      sprintf(errorMessage(),"Incorrectly specified noise kind (%s) in xmdsSimulation::writeRoutines",myParameters.noiseKind.c_str());
      throw xmdsException(errorMessage());
    }
  }

  xmdsElement::writeRoutines(outfile);
};

// ******************************************************************************
xmdsGlobals* xmdsSimulation::createxmdsGlobals() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsGlobals\n");
  }

  xmdsGlobals* newxmdsGlobals =
    new xmdsGlobals(this,verbose());
  addChild((xmdsElement*) newxmdsGlobals);
  return newxmdsGlobals;
};

// ******************************************************************************
xmdsField* xmdsSimulation::createxmdsField() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsField\n");
  }

  xmdsField* newxmdsField =
    new xmdsField(this,verbose());
  addChild((xmdsElement*) newxmdsField);
  return newxmdsField;
};

// ******************************************************************************
xmdsArgv* xmdsSimulation::createxmdsArgv() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsArgv\n");
  }

  xmdsArgv* newxmdsArgv = new xmdsArgv(this,verbose());
  addChild((xmdsElement*) newxmdsArgv);
  return newxmdsArgv;
};

// ******************************************************************************
xmdsOutput* xmdsSimulation::createxmdsOutput() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsOutput\n");
  }

  xmdsOutput* newxmdsOutput =
    new xmdsOutput(this,verbose());
  addChild((xmdsElement*) newxmdsOutput);
  return newxmdsOutput;
};

// ******************************************************************************
xmdsSequence* xmdsSimulation::createxmdsSequence() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsSequence\n");
  }

  xmdsSequence* newxmdsSequence =
    new xmdsSequence(this,verbose());
  addChild((xmdsElement*) newxmdsSequence);
  return newxmdsSequence;
};