xmdsintegrate.cc

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

    fprintf(outfile,",\n	const double *const _noise_vector");
  }
  fprintf(outfile,");\n");
  fprintf(outfile,"\n");
};

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

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

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_calculate_delta_a(\n",segmentNumber);
  fprintf(outfile,"	const double& _step");
  fprintf(outfile,",\n	const unsigned long cycle");
  if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
    fprintf(outfile,",\n	const double *const _noise_vector");
  }
  fprintf(outfile,") {\n");
  fprintf(outfile,"\n");

  const xmdsVector* mainVector;

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrate::writexSpaceRoutine: cannot find 'main' vector");
  }

  const char* typeName="";
  if(mainVector->vectorType()==COMPLEX) {
    typeName="complex";
  }
  else if(mainVector->vectorType()==DOUBLE) {
    typeName="double";
  }

  for(unsigned long i=0;i<mainVector->nComponents();i++) {
    fprintf(outfile,"%s d%s_d%s;\n",typeName,mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
  }
  fprintf(outfile,"\n");

  // add cross vectors to total vectors to use, but no duplicating names!

  list<XMLString> myTotalVectorsList = *vectorNamesList();

  for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
    list<XMLString>::const_iterator pXMLString2 = myTotalVectorsList.begin();
    while((pXMLString2 != myTotalVectorsList.end()) && (*pXMLString2 != *pXMLString)) {
      pXMLString2++;
    }
    if(*pXMLString2 != *pXMLString) {
      myTotalVectorsList.push_back(*pXMLString);
    }
  }

  simulation()->field()->vectors2space(outfile,0,myTotalVectorsList,"");

  if(crossVectorNamesList()->size() > 0) {
    fprintf(outfile,"_segment%li_calculate_cross_field();\n",segmentNumber);
    fprintf(outfile,"\n");
  }

  if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
    fprintf(outfile,"const double* _noises = _noise_vector;\n");
    fprintf(outfile,"\n");
  }


  list<XMLString>::const_iterator nextNLPElement = myNonLoopPropagationCodeList.begin();
  list<integrateMomentGroup>::const_iterator nextMGElement = myIntegrateMomentGroupList.begin();

  list<XMLString> theCodeList = *codeElementList();
  if(!(theCodeList.size()==1+numIntegrateMomentGroups()+numNonLoopPropagation())) {
    throw xmdsException("The list of integrate code elements is the wrong length!");
  }
    
  long unsigned int whichMG = 0;

  list<long> deleteMGArrayList;

  for(list<XMLString>::const_iterator codeElement = theCodeList.begin(); codeElement != theCodeList.end(); codeElement++) {
    if(!strcmp(codeElement->c_str(),"vectors")) {
      simulation()->field()->openLoops(outfile,0,myTotalVectorsList);
        
      fprintf(outfile,"\n");
        
      char indent[64];
      for(unsigned long i=0;i<nDims;i++) {
	indent[i]=0x09;
      }
      indent[nDims]=0;

      // integrate moment group pointer definition

      long unsigned int tempWhichMG = 0;
      for(list<integrateMomentGroup>::const_iterator tempNextMGElement = myIntegrateMomentGroupList.begin(); 
	  tempNextMGElement != myIntegrateMomentGroupList.end(); tempNextMGElement++) {

	long unsigned int nDimsRemaining = simulation()->field()->geometry()->nDims();
	tempWhichMG++;
			  
	for(list<bool>::const_iterator nextIntegrateDimension = tempNextMGElement->integrateDimensionList.begin(); 
	    nextIntegrateDimension != tempNextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	  if(*nextIntegrateDimension) {
	    nDimsRemaining--;
	  }
	}
				
	if(nDimsRemaining>0) {
	  fprintf(outfile,"        long _img%li_pointer = ",tempWhichMG);
	  for(long unsigned int j=0; j<nDimsRemaining-1; j++) {
	    fprintf(outfile,"(");
	  }
	  long unsigned int whichMoment = 0;
	  list<bool>::const_iterator nextIntegrateDimension = tempNextMGElement->integrateDimensionList.begin();

	  while(*nextIntegrateDimension) {
	    nextIntegrateDimension++;
	    whichMoment++;
	  }
	  fprintf(outfile,"_i%li",whichMoment);  
                  
	  nextIntegrateDimension++;
	  whichMoment++;
	  for(long unsigned int k=0; k<nDimsRemaining-1; k++) {
	    while(*nextIntegrateDimension) {
	      nextIntegrateDimension++;
	      whichMoment++;
	    }
	    fprintf(outfile,"*_%s_lattice%li+_i%li)",fieldName,whichMoment,whichMoment);
	    nextIntegrateDimension++;
	    whichMoment++;
	  }

	  fprintf(outfile,"*%i;\n",tempNextMGElement->momentNameList.size());
	}
		
      }
        
      fprintf(outfile,"// ************** propagation code **************\n");
      fprintf(outfile,"%s\n",propagationCode()->c_str());
      fprintf(outfile,"// **********************************************\n");
      fprintf(outfile,"\n");
      for(long unsigned int i=0;i<mainVector->nComponents();i++) {
	fprintf(outfile,"%s_active_%s_main[_%s_main_index_pointer + %li] = d%s_d%s*_step;\n",
		indent,fieldName,fieldName,i,mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
      }
      fprintf(outfile,"\n");
        
      if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
	fprintf(outfile,"%s_noises += _n_noises;\n",indent);
	fprintf(outfile,"\n");
      }
        
      simulation()->field()->closeLoops(outfile,0,myTotalVectorsList);
      fprintf(outfile,"\n");
    }
    else if(!strcmp(codeElement->c_str(),"functions")) {       
      fprintf(outfile,"// ************** propagation code **************\n");
      fprintf(outfile,"%s\n",nextNLPElement->c_str());
      fprintf(outfile,"// **********************************************\n");
      fprintf(outfile,"\n");
      nextNLPElement++;
    }
    else if(!strcmp(codeElement->c_str(),"moment_group")) {  
      whichMG++;  
              
      long unsigned int nDimsRemaining = simulation()->field()->geometry()->nDims();

      for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	  nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	if(*nextIntegrateDimension) {
	  nDimsRemaining--;
	}
      }

      if(nDimsRemaining>0) {
	deleteMGArrayList.push_back(whichMG);
      }

      fprintf(outfile,"// ************** integrate moment group code **************\n");
                
      //setup defines for the integrate_moment_groups if they are arrays
      long unsigned int whichMoment=0;
      if(nDimsRemaining>0) {
	for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin(); 
	    nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++){
	  fprintf(outfile,"#define %s _img%li_array[_img%li_pointer+%li]\n",nextMomentName->c_str(),whichMG,whichMG,whichMoment);
	  whichMoment++;
	}
      }
      fprintf(outfile,"\n");

      //setup actual arrays for the integrate_moment_groups (MPI and non-MPI)
      if(nDimsRemaining>0) {
	fprintf(outfile,"complex *_img%li_array = new complex[%i",whichMG,nextMGElement->momentNameList.size());
	if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	  whichMoment=0;
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	    if(!*nextIntegrateDimension) {
	      if(whichMoment==0){
		fprintf(outfile,"*local_nx");
	      }
	      else {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	    }
	    whichMoment++;
	  }
	}
	else {
	  whichMoment=0;
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	    if(!*nextIntegrateDimension) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	    }
	    whichMoment++;
	  }
	}
	fprintf(outfile,"];\n");
      }
      else {
	for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin(); 
	    nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++) {
	  fprintf(outfile,"complex %s = rcomplex(0,0);\n",nextMomentName->c_str());
	}
      }
      fprintf(outfile,"\n");
              
      //zero the arrays and setup a local environment for the sum
      if(nDimsRemaining>0) {
	fprintf(outfile,"for(long _i0=0; _i0<%i",nextMGElement->momentNameList.size());
	if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	  whichMoment=0;
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	    if(!*nextIntegrateDimension) {
	      if(whichMoment==0){
		fprintf(outfile,"*local_nx");
	      }
	      else {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	    }
	    whichMoment++;
	  }
	}
	else {
	  whichMoment=0;
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	    if(!*nextIntegrateDimension) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	    }
	    whichMoment++;
	  }
	}
	fprintf(outfile,"; _i0++) {\n");
	fprintf(outfile,"  _img%li_array[_i0] = rcomplex(0,0);\n",whichMG);
	fprintf(outfile,"}\n");
      }

      fprintf(outfile,"{\n");
      // open loops
      simulation()->field()->openLoops(outfile,0,myTotalVectorsList);
      // pointer definition
      if(nDimsRemaining>0) {
	fprintf(outfile,"        long _img%li_pointer = ",whichMG);
	for(long unsigned int j=0; j<nDimsRemaining-1; j++) {
	  fprintf(outfile,"(");
	}
	whichMoment = 0;
	list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();

	while(*nextIntegrateDimension) {
	  nextIntegrateDimension++;
	  whichMoment++;
	}
	fprintf(outfile,"_i%li",whichMoment);  
                  
	nextIntegrateDimension++;
	whichMoment++;
	for(long unsigned int k=0; k<nDimsRemaining-1; k++) {
	  while(*nextIntegrateDimension) {
	    nextIntegrateDimension++;
	    whichMoment++;
	  }
	  fprintf(outfile,"*_%s_lattice%li+_i%li)",fieldName,whichMoment,whichMoment);
	  nextIntegrateDimension++;
	  whichMoment++;
	}

	fprintf(outfile,"*%i;\n",nextMGElement->momentNameList.size());
      }
      // code
      fprintf(outfile,"%s\n",nextMGElement->integrateMomentGroupCode.c_str());

      // close loops
      simulation()->field()->closeLoops(outfile,0,myTotalVectorsList);
			  
      fprintf(outfile,"}\n");

      // if there was any integration, multiply by the relevant volume element
      if(nDimsRemaining<simulation()->field()->geometry()->nDims()) {
	if(nDimsRemaining>0) {
	  fprintf(outfile,"for(long _i0=0; _i0<%i",nextMGElement->momentNameList.size());
	  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	    whichMoment=0;
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	      if(!*nextIntegrateDimension) {
		if(whichMoment==0){
		  fprintf(outfile,"*local_nx");
		}
		else {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
		}
	      }
	      whichMoment++;
	    }
	  }
	  else {
	    whichMoment=0;
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	      if(!*nextIntegrateDimension) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	      whichMoment++;
	    }
	  }
	  whichMoment=0;
	  fprintf(outfile,"; _i0++) {\n");
	  fprintf(outfile,"   _img%li_array[_i0] = _img%li_array[_i0]",whichMG,whichMG);
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	    if(*nextIntegrateDimension) {
	      fprintf(outfile,"*_%s_dx%li",fieldName,whichMoment);
	    }
	    whichMoment++;
	  }
	  fprintf(outfile,";\n");
	  fprintf(outfile," }\n");
	}
	else {
	  for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin(); 
	      nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++) {
	    whichMoment = 0;
	    fprintf(outfile,"  %s = %s",nextMomentName->c_str(),nextMomentName->c_str());
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	      if(*nextIntegrateDimension) {
		fprintf(outfile,"*_%s_dx%li",fieldName,whichMoment);
	      }
	      whichMoment++;
	    }
	    fprintf(outfile,";\n");
	  }
	}
      }

      // if mpi and first variable integrated then mpireduce
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	if(* nextMGElement->integrateDimensionList.begin()) {
	  if(nDimsRemaining>0) {
	    fprintf(outfile,"\n  complex *_temp_img%li_array = new complex[%i",whichMG,nextMGElement->momentNameList.size());
	    whichMoment=0;
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	      if(!*nextIntegrateDimension) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	      whichMoment++;
	    }
	    fprintf(outfile,"];\n");
	    fprintf(outfile,"  MPI_Reduce(_img%li_array,_temp_img%li_array,2*%i",whichMG,whichMG,nextMGElement->momentNameList.size());
	    whichMoment=0;
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	      if(!*nextIntegrateDimension) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	      whichMoment++;
	    }
	    fprintf(outfile,",MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n");
	    fprintf(outfile,"\n  delete[] _temp_img%li_array;\n",whichMG);					
	  }
	  else {
	    fprintf(outfile,"\n  complex *_temp%li_mpireduce = new complex;\n",whichMG);
	    for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin(); 
		nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++) {
	      fprintf(outfile,"  MPI_Reduce(&%s,_temp%li_mpireduce,2,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n",nextMomentName->c_str(),whichMG);
	      fprintf(outfile,"\n  delete _temp%li_mpireduce;\n",whichMG);
	    }
	  }
	}
      }

      fprintf(outfile,"\n");
      fprintf(outfile,"// *********************************************************\n");
      fprintf(outfile,"\n");
      nextMGElement++;           
    }
    else {
      throw xmdsException("Unknown code element in the integrate block!");
    }
  }

  // delete any arrays that were generated
  for(list<long>::const_iterator deleteMGArray = deleteMGArrayList.begin(); 
      deleteMGArray != deleteMGArrayList.end(); deleteMGArray++) {
    fprintf(outfile,"delete[] _img%li_array;\n",*deleteMGArray);					
  }
  
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrate private
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
void xmdsIntegrate::outputSampleCount() const {
  if(debugFlag) {
    printf("xmdsIntegrate::outputSampleCount\n");
  }

  simulation()->output()->addSamples(mySamplesList);
};

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

  if(usesKOperators()) {
    if(constantK()) {
      fprintf(outfile,"_segment%li_calculate_k_operator_field();\n",segmentNumber);
    }
  }
};