xmdsmomentgroup.cc

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

  xmdsVector* nextVector = createxmdsVector();
  nextVector->setName("raw");
  myPostProcessingVectorNamesList.push_back("raw");
  if(complexRawVector) {
    nextVector->setVectorType(COMPLEX);
  }
  else {
    nextVector->setVectorType(DOUBLE);
  }
  nextVector->setInitialSpace(myRawSpace);
  nextVector->setComponents(mySamplingMomentsList);

  nextVector = createxmdsVector();
  nextVector->setName("fullstep");
  myPostProcessingVectorNamesList.push_back("fullstep");
  nextVector->setVectorType(DOUBLE);
  nextVector->setInitialSpace(myPostSpace);
  nextVector->setComponents(myPostMomentsList);

  if(simulation()->parameters()->stochastic) {
    nextVector = createxmdsVector();
    nextVector->setName("fullstep_sd");
    myPostProcessingVectorNamesList.push_back("fullstep_sd");
    nextVector->setVectorType(DOUBLE);
    nextVector->setInitialSpace(myPostSpace);
    nextVector->setComponents(myPostMomentsList);
  }

  if(simulation()->parameters()->errorCheck) {
    nextVector = createxmdsVector();
    nextVector->setName("halfstep");
    myPostProcessingVectorNamesList.push_back("halfstep");
    nextVector->setVectorType(DOUBLE);
    nextVector->setInitialSpace(myPostSpace);
    nextVector->setComponents(myPostMomentsList);
    if(simulation()->parameters()->stochastic) {
      nextVector = createxmdsVector();
      nextVector->setName("halfstep_sd");
      myPostProcessingVectorNamesList.push_back("halfstep_sd");
      nextVector->setVectorType(DOUBLE);
      nextVector->setInitialSpace(myPostSpace);
      nextVector->setComponents(myPostMomentsList);
    }
  }

  processVectors(myPostProcessingVectorNamesList,myPostSpace);

};

// ******************************************************************************
void xmdsMomentGroup::addSamples(
				 const unsigned long& n2add) const {
  if(debugFlag) {
    printf("xmdsMomentGroup::addSamples\n");
  }

  if(verbose())
    printf("Adding %li samples to moment group %li\n",n2add,myGroupNumber+1);

  nSamples += n2add;
};

// ******************************************************************************
void xmdsMomentGroup::finaliseGeometry() {
  if(debugFlag) {
    printf("xmdsMomentGroup::finaliseGeometry\n");
  }

  dimensionStruct firstDimension;

  firstDimension.name = simulation()->parameters()->propDimName;
  firstDimension.lattice = nSamples;
  firstDimension.domain.begin = 0;
  firstDimension.domain.end = 1;

  tempGeometry.setDimension(0,firstDimension);

  setGeometry(tempGeometry);
};

// ******************************************************************************
bool xmdsMomentGroup::requiresIntegrations() const {
  if(debugFlag) {
    printf("xmdsMomentGroup::requiresIntegrations\n");
  }

  return myRequiresIntegrations;
};

// ******************************************************************************
bool xmdsMomentGroup::needscomplexRawVector() const {
  if(debugFlag) {
    printf("xmdsMomentGroup::needscomplexRawVector\n");
  }

  return complexRawVector;
};

// ******************************************************************************
// ******************************************************************************
//				xmdsMomentGroup private
// ******************************************************************************
// ******************************************************************************

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

  if(verbose()) {
    printf("Writing moment group %li defines ...\n",myGroupNumber+1);
  }

  const char* mgFieldName = name()->c_str();

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// moment group %li defines\n",myGroupNumber);
  fprintf(outfile,"\n");
  fprintf(outfile,"#define _%s_ndims %li\n",mgFieldName,geometry()->nDims());
  fprintf(outfile,"#define _%s_lattice0 %li\n",mgFieldName,geometry()->dimension(0)->lattice);
  fprintf(outfile,"#define _%s_dx0 ((_%s_x0[_%s_lattice0-1] - _%s_x0[0])/(double)%li)\n",
	  mgFieldName,mgFieldName,mgFieldName,mgFieldName,geometry()->dimension(0)->lattice);
  fprintf(outfile,"#define _%s_dk0 (2*M_PI/(_%s_x0[_%s_lattice0-1] - _%s_x0[0]))\n",mgFieldName,mgFieldName,mgFieldName,mgFieldName);

  for(unsigned long i=1; i<geometry()->nDims(); i++) {

    const dimensionStruct* dimI = geometry()->dimension(i);

    fprintf(outfile,"#define _%s_lattice%li %li\n",mgFieldName,i,dimI->lattice);
    fprintf(outfile,"#define _%s_xmin%li %e\n",mgFieldName,i,dimI->domain.begin);
    fprintf(outfile,"#define _%s_dx%li ((%e - %e)/(double)%li)\n",mgFieldName,i,dimI->domain.end,dimI->domain.begin,dimI->lattice);
    fprintf(outfile,"#define _%s_dk%li (2*M_PI/(%e - %e))\n",mgFieldName,i,dimI->domain.end,dimI->domain.begin);
  }
  fprintf(outfile,"\n");

  fprintf(outfile,"#define _%s_raw_ncomponents %li\n",mgFieldName,(long)mySamplingMomentsList.size());
  fprintf(outfile,"#define _%s_fullstep_ncomponents %li\n",mgFieldName,(long)myPostMomentsList.size());

  if(simulation()->parameters()->errorCheck) {
    fprintf(outfile,"#define _%s_halfstep_ncomponents %li\n",mgFieldName,(long)myPostMomentsList.size());
  }

  if(simulation()->parameters()->stochastic) {

    fprintf(outfile,"#define _%s_fullstep_sd_ncomponents %li\n",mgFieldName,(long)myPostMomentsList.size());

    if(simulation()->parameters()->errorCheck) {
      fprintf(outfile,"#define _%s_halfstep_sd_ncomponents %li\n",mgFieldName,(long)myPostMomentsList.size());
    }
  }
  fprintf(outfile,"\n");

};

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

  if(verbose()) {
    printf("Writing moment group %li globals ...\n",myGroupNumber+1);
  }

  const char* mgFieldName = name()->c_str();

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// moment group %li globals\n",myGroupNumber);
  fprintf(outfile,"\n");
  xmdsField::writeGlobals(outfile);
  fprintf(outfile,"unsigned long _%s_sample_pointer;\n",mgFieldName);
  fprintf(outfile,"\n");
  fprintf(outfile,"double *_%s_x0 = new double[_%s_lattice0];\n",mgFieldName,mgFieldName);
  fprintf(outfile,"\n");
};

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

  if(verbose()) {
    printf("Writing moment group %li prototypes ...\n",myGroupNumber+1);
  }

  const char* mgFieldName = name()->c_str();

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// moment group %li prototypes\n",myGroupNumber);
  fprintf(outfile,"\n");
  xmdsField::writePrototypes(outfile);
  fprintf(outfile,"void _%s_sample();\n",mgFieldName);
  fprintf(outfile,"\n");
  fprintf(outfile,"void _%s_process();\n",mgFieldName);
  fprintf(outfile,"\n");
  fprintf(outfile,"void _%s_write_out(\n",mgFieldName);
  fprintf(outfile,"	FILE*);\n");
  fprintf(outfile,"\n");
};

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

  const unsigned long nDims = simulation()->field()->geometry()->nDims();
  const unsigned long myNDims = geometry()->nDims();
  const char* fieldName = simulation()->field()->name()->c_str();
  const char* mgFieldName = name()->c_str();

  unsigned long i;
  unsigned long j;

  if(verbose()) {
    printf("Writing moment group %s routines ...\n",mgFieldName);
  }

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"//		moment group %li routines\n",myGroupNumber);
  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"\n");

  xmdsField::writeRoutines(outfile);

  // ***************************************
  // ********  sample routine  *************
  // ***************************************

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _%s_sample() {\n",mgFieldName);
  fprintf(outfile,"\n");

  for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin(); pXMLString != mySamplingMomentsList.end(); pXMLString++) {
    if(complexRawVector) {
      fprintf(outfile,"complex %s;\n",pXMLString->c_str());
    }
    else {
      fprintf(outfile,"double %s;\n",pXMLString->c_str());
    }
  }
  fprintf(outfile,"\n");

  simulation()->field()->vectors2space(outfile,mySamplingSpace,myVectorNamesList,"");

  list<unsigned long>::const_iterator psamplingLatticeI = mySamplingLatticeList.begin();

  bool swapped=false;
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic){
    if(samplingSpace(0)&samplingSpace(1)) {
      swapped=true;
    }
  }

  if(!myRequiresIntegrations&(simulation()->parameters()->stochastic|(!simulation()->parameters()->usempi))) {
    {
      fprintf(outfile,"unsigned long _%s_raw_index_pointer=_%s_sample_pointer",mgFieldName,mgFieldName);
      for(i=1;i<myNDims;i++) {
	fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
      }
      fprintf(outfile,"*_%s_raw_ncomponents;\n",mgFieldName);
      fprintf(outfile,"\n");
    }
  }

  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
      if(swapped) {
	  int firstlattice = *psamplingLatticeI;
	  psamplingLatticeI++;
	  int secondlattice = *psamplingLatticeI;
	  psamplingLatticeI++;
        
	  if(secondlattice==0) {
	    // integration over this dimension
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"unsigned long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"\n");
            
	    fprintf(outfile,"double k%s = local_y_start_after_transpose*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName);
	    fprintf(outfile,"if(k%s>((_%s_lattice1-1)/2 + 0.1)*_%s_dk1)\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);	
	    fprintf(outfile,"	k%s -= _%s_lattice1*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);
                                    
	    fprintf(outfile,"\n");

	    fprintf(outfile,"for(long _i1=0; _i1<local_ny_after_transpose; _i1++) {\n");
	    fprintf(outfile,"\n");
	  }
	  else if(secondlattice==1) {
	    // cross-section in this dimension
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
	    }

	    fprintf(outfile,"\n");

	    fprintf(outfile,"double k%s = 0;\n",simulation()->field()->geometry()->dimension(1)->name.c_str());
	    fprintf(outfile,"unsigned long _i1 = 0;\n");                    
	    fprintf(outfile,"if(local_y_start_after_transpose==0) {\n");                    
	  }
	  else {
	    // normal sampling
	    fprintf(outfile,"long first_x_pointer, last_x_pointer;\n");
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"long _%s_%s_index_pointer;\n",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"double k%s;\n\n",simulation()->field()->geometry()->dimension(1)->name.c_str());
                    
	    fprintf(outfile,"if(local_y_start_after_transpose<(_%s_lattice%li-1)/2+1) {\n",mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   first_x_pointer=local_y_start_after_transpose;\n");
	    fprintf(outfile,"   k%s=local_y_start_after_transpose*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName);
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"   _%s_%s_index_pointer=0;",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"   }\n");
	    fprintf(outfile,"else if(local_y_start_after_transpose>(_%s_lattice%li-1)/2+_%s_lattice1-_%s_lattice%li) {\n"
		    ,mgFieldName,main2PostDim(1),fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   first_x_pointer=local_y_start_after_transpose-_%s_lattice1+_%s_lattice%li;\n",fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   k%s=(local_y_start_after_transpose-_%s_lattice1)*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"   _%s_%s_index_pointer=0;",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"   }\n");
	    fprintf(outfile,"else {\n");
	    fprintf(outfile,"   first_x_pointer=(_%s_lattice%li-1)/2+1;\n",mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   k%s=((_%s_lattice%li-1)/2-_%s_lattice%li+1)*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),mgFieldName,main2PostDim(1),mgFieldName,main2PostDim(1),fieldName);
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		fprintf(outfile,"   _%s_%s_index_pointer=(first_x_pointer+_%s_lattice1-_%s_lattice%li-local_y_start_after_transpose)",fieldName,pXMLString->c_str(),fieldName,mgFieldName,main2PostDim(1));
		fprintf(outfile,"*_%s_lattice0",fieldName);
		for(i=2;i<nDims;i++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,i);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	    fprintf(outfile,"   }\n\n");
                    
	    fprintf(outfile,"if(local_y_start_after_transpose+local_ny_after_transpose-1<(_%s_lattice%li-1)/2+1)\n",mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   last_x_pointer=local_y_start_after_transpose+local_ny_after_transpose-1;\n");
	    fprintf(outfile,"else if(local_y_start_after_transpose+local_ny_after_transpose-1>(_%s_lattice%li-1)/2+_%s_lattice1-_%s_lattice%li)\n"
		    ,mgFieldName,main2PostDim(1),fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   last_x_pointer=local_y_start_after_transpose+local_ny_after_transpose-1-_%s_lattice1+_%s_lattice%li;\n",fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"else \n");
	    fprintf(outfile,"   last_x_pointer=(_%s_lattice%li-1)/2;\n\n",mgFieldName,main2PostDim(1));
    
	    fprintf(outfile,"for(long _i1=first_x_pointer;_i1<last_x_pointer+1;_i1++) {\n");
	    fprintf(outfile,"\n");
	  }
                

	  const char* dimName = simulation()->field()->geometry()->dimension(0)->name.c_str();

	  fprintf(outfile,"	");
	  fprintf(outfile,"double k%s = 0;\n",dimName);
	  fprintf(outfile,"\n");

	  if(firstlattice==0) {
	    // integration over this dimension
	    fprintf(outfile,"	");
	    fprintf(outfile,"for(unsigned long _i0=0; _i0<_%s_lattice0; _i0++) {\n",fieldName);